Unique Transcriptome and Gene Expression Analysis of Rice Seedling Reveals Different Cadmium Response Regulatory Mechanisms Between Indica and Japonica Rice

Background: In general, the Cd content in indica rice is usually higher than that in japonica. However, the mechanism for this discrepancy is unclear. Thus, understanding the genetic and molecular basis of Cd stress between indica and japonica is extremely important for rice improvement programs. Results: In this study, two varieties of rice, japonica 02428 and indica CH891, were continuously exposed to Cd and seedlings of the two varieties at two critical stages (3 rd and 5 th day) were selected for the dynamic genes analysis by transcriptome method. The results showed that CH891 was more sensitive to Cd than 02428, and a total of 7,204 and 6,670 differently expressed genes (DEGs) associated with Cd stress were detected at 3 rd day and 5 th day. Furthermore, we divided these DEGs into three categories: SCR (sensitive variety with Cd-responsive), RCR (resistant variety with Cd-responsive) and CCR (common Cd-responsive). The enriched metabolic pathways analysis of DEGs preferentially expressed in a stage-specic and cultivars-specic manner, and secondary metabolic processes were enriched in SCR while protein metabolism and plant hormone were enriched in RCR. The diverted metabolic pathways might be the major cause for different response mechanism of Cd in indica and japonica rice. Conclusion: These results provide a novel insight into the Cd response mechanism in rice seedlings between different varieties, and these important Cd-responsive DEGs were frequently involved in specic biological processes and metabolic pathways that might provide a novel insight over indica and japonica rice Cd response mechanism difference. Phenylpropanoid biosynthesis, glutathione metabolism, diterpenoid biosynthesis, brassinosteroid biosynthesis, alpha-linolenic acid metabolism and ABC transporters pathways were enriched in Cd5_CH891 vs. CK5_CH891. Protein processing in endoplasmic reticulum, porphyrin and chlorophyll metabolism, photosynthesis, phenylpropanoid biosynthesis and circadian rhythm-plant pathways were enriched in Cd5_02428 vs. CK5_02428.

the genetic and molecular basis of Cd stress between indica and japonica is extremely important for rice improvement programs.
Recently, many studies have used different rice cultivars to study the Cd stress [11], and several genes were identi ed to be associated with Cd response, such as P 1B -ATPase genes, natural resistanceassociated macrophage protein genes, cation diffusion facilitator genes, ATP-binding cassette transporter and other low-a nity cation transporter genes [12][13][14][15]. However, it was con rmed that genes encoding GST, heat shock protein and cytochrome P450 were also strongly induced under Cd stress [16]. Although functional analyses of these individual genes are helpful for understanding the regulatory mechanisms of Cd response in rice, the genetic basis of mechanisms for different subspecies is not clearly identi ed. Therefore, a comprehensive understanding of molecular mechanisms regulating seed development of different subspecies in Cd stress is required to facilitate the development of new insight into Cd response mechanism. Access to transcriptome and next generation RNA sequencing (RNA-seq) technology provides an opportunity to reveal genetic diversity among various genotypes/cultivars under Cd stress.
To identify genetic and molecular basis of Cd stress between indica and japonica, an elite japonica variety, 02428, and an elite indica variety, CH891, were continuously treated with Cd stress, and the seedling of the two varieties at two critical stages were selected for the analysis of the dynamic gene by transcriptome method. In total, we identi ed 7204 and 6584 DEGs at 3rd and 5th day. To investigate the different Cd response mechanisms in different cultivars and stages, we study enriched metabolic pathways of the DEGs and found that the response mechanisms preferentially were in a stage-speci c and cultivar-speci c manner. The results highlighted a new insight for difference in Cd tolerance between indica and japonica.

Plant materials
Two rice (oryza sativa L.) varieties were used in the study. Changhui 891 (CH891), excellent indica restore line in south China, generated from Key laboratory of crop genetics and breeding, Jiangxi Agricultural University, Jiangxi province, China, whose variety rights number in China Rice Data Center is CNA20161213.9 [17]. 02428 is a wide a nity japonica variety, selected and bred by Institute of genetic physiology, Jiangsu academy of agricultural sciences, Jiangsu province, China. It is a hybrid of radiation offspring of crab rice in Yunnan province and radiation offspring of Jibang rice in Shanghai [18], also being collected in the China Rice Data Center.

Plant growth conditions and treatments
Rice seedling were grown in culture dishes. Brie y, 30 seeds were selected by getting rid of shriveled and empty seed and surface-disinfected by immersing in 4% solution of sodium hypochlorite for twenty minutes and drenched for three times with distilled water. Thirty germinative seeds were placed in each dish, and the treatment groups (Cd) were treated with 100 mg/L CdCl 2 while the control groups (CK) with distilled water. Seedling were maintained in the incubator at 28℃ ± 2 for 16/8 h light/dark condition. Appropriate CdCl 2 solution and distilled water were added every day according to the requirement to avoid desiccation. Three biological replicates were arranged for each line.

Samples harvest and phenotype characterization
Seedling of two varieties (CH891 and 02428) were harvested and measured the seedling length of the population at the third day (Cd3 and CK3) and fth day (Cd5 and CK5). Seedlings were wrapped with tin foil paper, immediately frozenin liquid nitrogen, and stored at -80℃. Three biological replications and three technical replications were performed for the measurement of seedling length. mRNA library construction and sequencing Total RNA was extracted for three biological and technical replicates in each of 3rd and 5th day after Cd stress (DAS) and 3rd and 5th day after control (DAC) using Trizol reagent (Invitrogen, CA, USA) following the manufacturer's procedure. The total RNA quantity and purity were analysis of Bioanalyzer 2100 and RNA 6000 Nano LabChip Kit (Agilent, CA, USA) with RIN number > 7.0. Approximately 10 ug of total RNA representing a speci c adipose type was subjected to isolate Poly (A) mRNA with poly-T oligoattached magnetic beads (Invitrogen). Following puri cation, the mRNA is fragmented into smallpieces using divalent cations under elevated temperature. Then the cleaved RNA fragments werereverse-transcribed to create the nal cDNA library in accordance with the protocol for the mRNASeqsample preparation kit (Illumina, San Diego, USA), the average insert size for the paired-endlibraries was 300 bp (± 50 bp). And then we performed the paired-end sequencing on an IlluminaHiseq4000 at the (LC Sceiences,USA) following the vendor's recommended protocol.
Normalization of gene expression levels and identi cation of differentially expressed genes Sequencing reads were mapped to the reference sequences. The mapped reads of each sample were assembled using StringTie. Then, all transcriptomes from samples were merged to reconstruct a comprehensive transcriptome using perl scripts. After the nal transcriptome was generated, StringTie and edgeR was used to estimate the expression levels of all transcripts. StringTie was used to perform expression level for mRNAs by calculating fragments per kilobase of exon model per Million mapped reads (FPKM). The differentially expressed mRNAs and differentially expressed genes (DEGs) were selected with log2 (fold change) > 1 or log2 (fold change) <-1 and with statistical signi cance (p value < 0.05) by R package.

Quantitative real-time PCR (qRT-PCR) validation
In order to validate the RNA-seq results, different expression patterns of several genes was con rmed by quantitative real-time RT-PCR (qRT-PCR). For qRT-PCR, 1 µg of total RNA was used to synthesized cDNA using PrimeScript™ RT reagent Kit (Perfect Real Time) (TaKaRa). The qRT-PCR was carried out using SYBR® Premix Ex Taq II (Tli RNaseH Plus; TAKARA BIO Inc., Shiga, Japan) and determined in LightCycler 480 (Roche, Basel, Switzerland) according to the manufacturer's instructions. The qRT-PCR reactions were ampli ed for 95 •C for 30 s, followed by 40 cycles of 95•C for 5 s, 55•C for 30 s and 72•C for 30 s. All reactions were performed with three independent biological replicates for each sample and three technical replicates for each biological replicate were analyzed. The relative gene expression was calculated by the software of ABI7500 Real-Time PCR System using the 2 −∆∆Ct method.

Functional annotation and GO and KEGG classi cation
All expressed genes and signi cant genes were functional annotated against databases NR, KEGG, KOG, Pfam and Swiss prot, respectively. For the gene matched to multiple protein sequences, the protein with the highest similarity score was considered as the optimal annotation.

Phenotypic variation between CH891 and 02428 under Cd stress
Following exposure of seedlings of rice strains CH891 and 02428 to Cd, large phenotypic variations. The seedling lengths of CH891 and 02428 were shorter under Cd stress compared with the control condition ( Fig. 1A-D). Furthermore, the seedling length of CH891 was signi cantly longer than 02428 in control condition, while there was no signi cant differences between the two genotypes with Cd stress (Fig. 1E, F). It indicated that CH891 was more sensitive to Cd than 02428. The experiment carried out continuously for 7 days and found that the top two seedling length inhibition rate were on 4th and 6th day (Fig. 1G). It indicated that genes had higher expression in the 3rd and the 5th day than other days. Therefore, seedlings were harvested at the 3rd and 5th day and then sequenced RNAs.
RNA sequencing of seedling transcriptome of the two genotypes A total of 1,070 million reads with average of 46 million reads per sample were generated. While the total of valid reads is 1,057 million and the average valid read per sample is 46 million. The ratio of Q20 for each was above 99%, and the Q30 base percentage was above 95% (Supplementary Table S1). Therefore, the quality of the data was very high and meet the requirements for further analysis. To con rm the accuracy and reproducibility of the RNA-Seq results, 16 genes were selected for qRT-PCR using the speci c primers (Supplementary Table S2). The validation results for the 16 genes are shown in Fig. 5A, B. The qRT-PCR results were all consistent with the RNA-Seq data (Supplementary Table S2). To conclude our transcriptome sequencing results were credible.

Identi cation of differentially expressed genes (DEGs)
By comparing samples of the same rice cultivar in different conditions (control and stress) and different rice cultivar (CH891 and 02428) in the same condition at different stages, we assigned four comparison groups at the two stages and eight comparison groups in total were obtained. And then selected DEGs in different comparision respectively by restricting pval ≤ 0.05. It indicated that not only the cultivar but also the treatments and time points affect the gene expression level. The number of up-regulated and downregulated genes in each comparison group is shown in Fig. 3A.
Functional enrichment analysis was performed for all these DEGs ( Supplementary Fig. S1). KEGG enrichment analysis showed that plant-pathogen interaction, phenylpropanoid biosynthesis, tryptophan metabolism, avonoid biosynthesis and diterpenoid biosynthesis pathways were enriched in Cd3_CH891 vs. Cd3_02428. Plant-pathogen interaction, tryptophan metabolism, phenylpropanoid biosynthesis, avonoid biosynthesis, diterpenoid biosynthesis and amino sugar and nucleotide sugar metabolism pathways were enriched in CK3_CH891 vs. CK3_02428. Starch and sucrose metabolism, phenylpropanoid biosynthesis, glycolysis, carbon metabolism, alpha-linolenic acid metabolism, fatty acid degradation and brassinosteroid biosynthesis pathways were enriched in Cd3_CH891 vs. CK3_CH891. Plant-pathogen interaction and plant hormone signal transduction pathways were enriched in Cd3_02428 vs. CK3_02428. On the 5th day, phenylpropanoid biosynthesis, glutathione metabolism, avonoid biosynthesis, brassinosteroid biosynthesis, amino sugar and nucleotide sugar metabolism and ABC transporters pathways were enriched in Cd5_CH891 vs. Cd5_02428. Tryptophan metabolism, sesquiterpenoid and triterpenoid biosynthesis and phenylpropanoid biosynthesis pathways were enriched in CK5_CH891 vs.
All of the results revealed that there were more up-regulated DEGs under Cd stress than those under control condition in CH891, while there were more down-regulated DEGs under Cd stress than those under control condition in 02428 on either the 3rd or 5th day. Moreover, there was a higher number of DEGs obtained for Cd stress in CH891 (Cd-sensitive cultivar) than those in 02428 (Cd-resistant cultivar).

Classi cation of differentially expressed genes (DEGs)
A total of 7204 unique DEGs were detected in four groups on the 3rd day and 6670 unique DEGs obtained on the 5th day, these DEGs divided into 15 subgroups (Fig. 3B, C). Excluding these DEGs of the groups (CK3_CH891 vs. CK3_02428 at the third day and the groups contained CK5_CH891 vs. CK5_02428 at the fth day) irrelevant to Cd stress, the DEGs of the rest seven subgroups divided into three categories: genes from sensitive variety with Cd-responsive (SCR), genes from resistant variety with Cd-responsive (RCR), and common Cd-responsive (CCR) DEGs. There were 849, 676 and 770 DEGs in SCR, RCR and CCR at the 3rd day, respectively (Supplementary Table S3). Whereas 959, 592 and 1516 DEGs were in SCR, RCR and CCR at the 5th day, respectively (Supplementary Table S3).
In addition, we selected important Cd-responsive (ICR) genes in the three categories with |log2 (fold change)|≥2 and FPKM ≥ 2 in at least one group to select representative DEGs. At the 3rd day, there were 554 ICR genes were screened, among which 194 for SCR, 154 for RCR and 206 for CCR (Supplementary  Table S3). Whereas a total of 941 ICR genes were screened at the 5th day, of these 389 for SCR, 216 for RCR and 336 for CCR (Supplementary Table S3). Subsequently, there were 17, 9 and 21 ICR genes were detected at the two stages in SCR, RCR and CCR, respectively (Supplementary Table S3).
Gene functional annotation analysis of SCR Gene functional annotation analysis was conducted for ICR of SCR. Cytoplasm and DNA binding were especially enriched in SCR3, oxidation-reduction process and protein binding were particularly enriched in SCR5. Transcription factor activity, sequence-speci c DNA binding regulation of transcription, DNAtemplated, integral component of membrane and plasma membrane were enriched in both SCR3 and SCR5 ( Supplementary Fig. S2). In order to gain more biological information and regulatory network, the Kyoto Encyclopedia of Genes and Genomes (KEEG) enrichment pathways was performed for understanding the molecular mechanism of Cd response in rice seedlings. KEGG analysis results showed that, iso avonoid biosynthesis, inositol phosphate metabolism, avonoid biosynthesis, and anthocyanin biosynthesis pathways were enriched in SCR3 ( Fig. 2A). While iso avonoid biosynthesis and anthocyanin biosynthesis pathways were signi cantly enriched. One DEG involved in anthocyanin biosynthesis was down-regulated while one was up-regulated in Cd3_CH891. And two DEGs involved in iso avonoid biosynthesis were up-regulated in Cd3_CH891. Phenylpropanoid biosynthesis, alpha-linolenic acid metabolism, linoleic acid metabolism, glutathione metabolism, anthocyanin biosynthesis, ABC transporters, brassinosteroid biosynthesis and diterpenoid biosynthesis pathways were enriched in SCR5 (Fig. 2B) Transcription factor activity, sequence-speci c DNA binding protein binding, regulation of transcription and DNA-templated were enriched in both RCR3 and RCR5 (Supplementary Fig. S2).
In terms of KEGG pathway enrichment analysis, ribosome, regulation of autophagy and insulin resistance pathways were enriched in ICR of RCR3 (Fig. 2C). Among which ribosome and regulation of autophagy pathways were most enriched. One DEG involved in regulation of autophagy was down-regulated and the other two were up-regulated in Cd3_02428. One DEG involved in ribosome was up-regulated and six were down-regulated in Cd3_02428. Protein processing in endoplasmic reticulum, plant hormone signal transduction, thiamine metabolism, glucosinolate biosynthesis and base excision repair pathways were enriched in ICR of RCR5 (Fig. 2D). And protein processing in endoplasmic reticulum and plant hormone signal transduction pathways were most enriched. Four DEGs involved in protein processing in endoplasmic reticulum were up-regulated and other four were down-regulated in Cd5_02428. One DEG involved in plant hormone signal transduction was up-regulated and the other ten were down-regulated in Cd5_02428.
Gene functional annotation analysis of CCR As far as ICR of CCR, chloroplast and zinc ion binding were enriched on the third day while ATP binding was enriched in the fth day. And protein binding regulation of transcription, DNA-templated, cytoplasm, plasma membrane, integral component of membrane, transcription factor activity, sequence-speci c DNA binding and DNA binding were enriched in both CCR3 and CCR5 (Supplementary Fig. S2).
In terms of KEGG pathway analysis, alpha-linolenic acid metabolism, phagosome, limonene and pinene degradation and fatty acid degradation pathways were enriched in CCR3 (Fig. 2E). Moreover, alphalinolenic acid metabolism and fatty acid degradation pathways were markedly enriched in CCR3. One DEG involved in alpha-linolenic acid metabolism was down-regulated and two DEGs were up-regulated. And one DEG involved in fatty acid degradation was up-regulated while the other one was downregulated. Plant-pathogen interaction, iso avonoid biosynthesis, vitamin B6 metabolism, diterpenoid biosynthesis, and glutathione (GSH) metabolism pathways were enriched in the fth day (Fig. 2F).
Among which plant-pathogen interaction, iso avonoid biosynthesis, and glutathione (GSH) metabolism pathways were observably enriched. One DEG involved in iso avonoid biosynthesis was down-regulated, two were up-regulated. Two DEGs involved in glutathione metabolism were down-regulated, three were up-regulated. Seven DEGs involved in plant − pathogen interaction were down-regulated and twenty were up-regulated.
To understand the interaction of all the enriched GO terms at all the stages of SCR, RCR and CCR in the general, we constructed a network of signi cantly enriched GO terms (biological process) (Fig. 4). It showed that the biological process was focus on transport, response to different stimulus, metabolic of hormone, gene express, cellular activity, and secondary metabolic process. It suggested that different stress reactions were activated to respond to Cd stress.

DEGs both in 3rd and 5th day after stress (DAS)
In total, we obtained 46 common ICR in three categories at two stages and they were involved in several different pathways (Fig. 5C, D, E). Except for the pathways listed above, some other pathways were detected, such as inositol phosphate metabolism, cutin, suberine, and wax biosynthesis, biosynthesis of amino acids, starch and sucrose metabolism, ascorbate and aldarate metabolism, base excision repair, aminoacyl-tRNA biosynthesis, galactose metabolism, and circadian rhythm-plant. It indicated that not only the enriched pathways above but also insigni cantly enriched pathways had interaction to regulate Cd stress in rice.

Genetic basis of Cd stress in rice
The Cd stress extremely affects rice growth and development [19][20], such as photosynthesis, transpiration and other physiological processes of plants and produces excessive oxygen free radicals. In recent years, Cd has fascinated great attention due to its harmful effects on plant productivity. Several studies have shown that different cultivars showed different responses under Cd stress [11,21]. However, the regulation mechanism between indica and japonica under Cd stress have not yet been reported. In the present study, we evaluated an elite japonica variety, 02428 (Cd sensitive cultivar), and one elite indica variety, CH891 (Cd resistant cultivar), by exposing seedling were continuously treated with Cd solution 100 mg/L for one week, seedling of the two varieties at two critical stages, when it determines the length of the seedling of the two varieties, were selected for the dynamic genes analysis under Cd stress by transcriptome method. The results showed that the metabolic pathways between 3rd DAS and 5th DAS were different not only for 02428 but also for CH891, implying that distinct genetic systems might be responsible for Cd stress at different stages in rice, a thorough and dynamic gene analysis for rice Cd stress was necessary. Otherwise, the metabolic pathways between CH891 and 02428 were also different at the same stage (3rd or 5th DAS), these metabolic pathways might be the cause that indica and japonica showed different response to Cd stress, therefore, the exploitation genes in these metabolic pathways might fully explain the differences in Cd stress between indica and japonica. All, the DEGs analysis showed that the difference in metabolic pathway not only exists in different genotypes but also in different stages.
The stage speci c expression of DEGs between indica and japonica KEGG pathway analysis was performed for all DEGs with signi cant differences. It showed that the enriched pathways of all DEGs are not signi cantly different from those of the following three categories (SCR, RCR, CCR), except for some amino sugar metabolism, amino acid metabolism, sugar metabolism and nucleotide sugar metabolism. In general, metabolic pathways of all DEGs mainly focus on basic life activities, secondary metabolic processes, plant-pathogen interaction and plant hormone signal transduction pathways. This suggests that the following analysis would not be too disturbed by the background.
To give more detailed information for the genes response to Cd stress in rice, we divided DEGs into three categories (SCR, RCR, and CCR) at the two critical stages. The enriched metabolic pathways analysis of DEGs indicated that these DEGs preferentially expressed in a stage-speci c and cultivars-speci c manner, some signi cant metabolic pathways were identi ed for SCR, RCR and CCR.
Expressed Gene special in SCR The iso avonoid biosynthesis and anthocyanin biosynthesis pathways were enriched in SCR3. Iso avones were found to play an important role in plant defense reactions [22]. Under Cd stress, the cell wall of rice was damaged, which made it susceptible to pathogen infection and then induced the expression of iso avones to improve the disease resistance in plants [23]. Previous studies showed that iso avones are usually abundant in legumes but low in other plants [24]. Combined with our result, it may provide new insight into the response of rice to Cd stress at the iso avone level.
Moreover, Anthocyanin is a kind of antioxidant, which can scavenge oxygen free radicals in plant cells and relieve the toxicity of reactive oxygen [25]. Roychoudhury reported that anthocyanin biosynthesis was induced under CdCl 2 [11], which was consistent with our result. It suggested that anthocyanin, as an antioxidant, relates to Cd stress response in rice and may relieve the toxicity of Cd.
Phenylpropanoids biosynthesis was one of the enriched KEGG pathways in SCR5. Phenylpropanoids participated in the antioxidant activity of cell walls and the biosynthesis of lignin, which plays an important role in the response of plants to abiotic stress [26]. Several reserchers have found that phenylpropanoids metabolite was signi cantly altered by Cd stress in rice [27] and Arabidopsis helleri [28]. Our results con rmed that phenylpropanoids contribute to Cd response in rice. Besides, another enriched pathway in SCR5 was ABC transporters pathway, which was metabolic pathway authenticated to be important under Cd in rice [29]. Using the energy released by hydrolyzed ATP, ABC transporters can actively transport a variety of heavy metals through the cytomembrane, which is associated with plant stress tolerance and detoxi cation of heavy metals [30]. However, the result showed that ABC transporters was found to be enriched in SCR5 may put forward that ABC transporters participate in Cd response in rice and have possibility to inhibit Cd stress. From our results and previous study, ABC transporters actually play important role in response to Cd stress. Diterpenoids, as secondary metabolites, was reported to inhibit intracellular ROS production and lipid peroxidation, and then enhanced the antioxidant defense system in Sideritis [31]. Beside, our study found that diterpenoid biosynthesis pathway was enriched in SCR5, it suggested that diterpenoid metabolism also be induced by Cd in rice and plays an important role in response to Cd stress in rice. Brassinosteroid (BR) has a wide range of biological functions and plays an important role in the adaptability of plants under stress. Bajguz and Hayat found that external application of BR can improve the salt tolerance of plants [32]. Our result indicated that BR was also related to Cd response in rice, and DEGs involved in this pathway were all related to cytochrome P450 and were up-regulated, which is consistent with the fact that cytochrome P450 participates in detoxi cation metabolism [16].

Expressed Gene special in RCR
In terms of RCR, the KEGG pathways were also different at 3rd and 5th DAS. Ribosome-related DEGs play key roles in cold stress signal transduction [33]. Therefore, ribosome pathway was identi ed to be enriched in this study. In addition to ribosome, the regulation of autophagy was another enriched pathway in RCR3. Autophagy can transport oxidative proteins and injured intracellular substances caused by stress into vacuoles for degradation, thereby reducing the accumulation of toxic substances in the cytoplasm [34]. It plays an important role in plant development and response to stress [35]. Thus, our study found that regulation of autophagy pathway was induced in RCR3 further veri ed that it played a role in Cd stress in rice to a certain extent, and these two pathways were newly identi ed in rice under Cd stress here, and which identi ed only at 3rd DAS.
As far as RCR5, plant hormone signal transduction was enriched. When in the unfavorable environment, plants must respond quickly and accurately to activate the necessary physiological responses which are usually mediated by plant hormones [36]. Phytohormones have been demonstrated to play an important role in mediating plant growth plasticity in response to metal stress include Cd [37]. Consistently, we detected plant hormone signaling transduction pathway in RCR5. Our result was consistent with the previous study in rice, however, it was only identi ed at RCR5, and protein processing in endoplasmic reticulum pathway was also enriched in RCR5. Cellular exposure to Cd is known as strongly induce the unfolded protein response, which suggests that the endoplasmic reticulum (ER) is preferentially damaged by Cd in yeast [38]. In our study, it was newly identi ed in rice under Cd stress. It suggested that Cd can lead to protein alteration among ER in rice.

Expressed Gene special in CCR
Differences also existed in CCR in the 3rd and 5th DAS. Alpha-linolenic acid metabolism and fatty acid degradation can be summarized as fatty acid metabolism, which was related to Cd response in CCR3. Previous studies revealed that heavy metal can speci cally alter fatty acid metabolism in plants, Cd can enhance the level of lipid peroxides and lead to strong changes in fatty acid content in rice [39][40]. Fatty acid metabolism was found to be induced by heavy metal Cr in rice [41]. In this study, our ndings provided further evidence that fatty acid metabolism is a stress response to Cd stress in rice. As far as the KEGG pathway analysis in CCR5, plant-pathogen interaction, GSH metabolism and iso avonoid biosynthesis pathway were enriched. There is a very complex interaction between the endogenous bacteria or fungi in the plant. These microorganisms can produce a variety of biological effects in plants, they can form a special growth state to better adapt to the environment [42]. Several previous studies have found that plant pathogen interaction is induced by plant pathogen invasion [43]. The nding of this study puts forward the possibility that Cd damages the function and structure of cells and provides the conditions for pathogen attack. Many studies have reported that GSH alleviates heavy metal toxicity in rice [28,41]. However, the binding of GSH to heavy metal was catalyzed by glutathione S-transferase (GST), which is an important class of antioxidant and detoxi cation proteins in plant cells. In this study among the DEGs involved in GSH metabolism, all of the six were GST genes. It suggested that GST make a difference in GSH metabolism, which participates in Cd response in rice. Also, similar to pathway analysis in SCR3, iso avonoid biosynthesis pathway reacted on Cd stress in CCR5. It implied that iso avonoid biosynthesis pathway generally exist in Cd response in rice.
It can be observed that not only the response mechanisms between 3rd and 5th DAS in the same category but also mechanisms between different categories were different. Enriched pathways in SCR were associated with secondary metabolic processes while protein metabolism and plant hormone signal transduction were enriched in RCR. It indicated that CH891 might respond to Cd by inducing secondary metabolic processes and then can develop its immunity and resistance to adversity. Differently, the process of protein production and degradation was affected in 02428, so as to protect itself from Cd. It may prove that the response mechanism of Cd in indica and japonica rice was different in one way.

The common DEGs between 3rd and 5th DAS in indica and japonica
With the exception of metabolic pathways discussed above, the common DEGs in three categories at two stages were detected, which were enriched in various pathways (Fig. 4C, D, E). Some pathways were discussed above and others were not signi cantly enriched in this study, which may involve the minor genes and conspire to Cd response in rice. Of these, DEGs, LOC_Os03g57240 (DST), a common DEG in RCR at 3rd and 5th DAS, was reported to be a zinc nger transcription factor and relate to DNA replication, which negatively regulates drought and salt tolerance of rice [44]. Therefore, it was detected in this study, which implied DST regulated Cd stress in rice. We selected four genes to conduct qRT-PCR, and they were veri ed compared with RNA-Seq (Fig. 2B).

Conclusions
In this study, high-throughput sequencing was used to sequence the rice seedling transcriptome treated with Cd at different stages, which highlighted the transcriptional variations among two different rice varieties under Cd conditions. Statistical analysis of 7204 and 6670 DEGs revealed three categories for a total of 554 and 941 ICR DEGs in rice at two stages, respectively. Furthermore, these important Cdresponsive DEGs were frequently involved in speci c biological processes and metabolic pathways that    The comparison of RNA-seq results and qRT-PCR analysis of gene expression levels and the KEGG pathways of common ICR in SCR, RCR, CCR at two stages. (A) LOC_Os03g08470, LOC_Os04g39880 and LOC_Os09g36700 were in Cd3_CH891 vs. CK3_CH891, LOC_Os10g40720 was in Cd3_02428 vs.