Comparative Transcriptomic Analysis of Root Cadmium Responses in Two Chinese Rice Cultivars Yuzhenxiang and Xiangwanxian 12

Long Ping Branch of Graduate School, Central South University, Changsha, Hunan 410125, China Guizhou Rice Research Institute, Guizhou Academy of Agricultural Sciences, Guiyang, Guizhou 550006, China Hunan Rice Research Institute, Hunan Academy of Agricultural Sciences, Changsha, Hunan 410125, China Hunan Agricultural Biotechnology Research Institute, Hunan Academy of Agricultural Sciences, Changsha, Hunan 410125, China


Introduction
Rice (Oryza sativa L.) is the primary source of calorie intake in plenty of countries, including China, and the staple food of over half of the world's population [1,2]. To meet the food demand of increasing world populations, it is imperative to improve the rice yield [3]. However, the sustainability and productivity of the rice production system are constantly threatened by multiple factors, such as environmental pollution, climate change, and excessive use of pesticides and fertilizers [2].
Heavy metal pollution is the predominant type in all soil contaminants, and cadmium (Cd) contamination ranks the first among heavy metal contaminants in China [4,5]. Cd contamination, mainly caused by mining, e-waste, and overuse of nitrogen and phosphate fertilizers, has been found in large-scale agricultural soil in China [6,7]. Cd is readily absorbed by rice plants and easily transferred to food chains [8]. Moreover, it is often accumulated in rice grains and human bodies due to its long half-life of up to 25-30 years [8,9]. Excessive Cd exposure not only can inhibit the growth of rice plants but also can decrease the quality/ nutrients/yields of rice grains [7]. Additionally, the intake of Cd-contaminated rice can markedly increase the risks of multiple diseases, such as cancer, osteoporosis, and liver/ kidney injury, as well as nervous system diseases [9][10][11]. To reduce the potential harm of Cd on human health, it is imperative to screen out the genes that play vital roles in Cd response in rice plants [7].
It has been reported that there is a notable difference in Cd content in the grains of different rice cultivars, which might be caused by the genotypic and environmental diversities among rice cultivars [12,13]. More specifically, rice Cd concentration is closely associated with heading time, soil pH, mutation/dysregulation/diversity of multiple genes, or quantitative trait locus related to Cd response (absorption, translocation, and accumulation) among different rice subspecies and cultivars [12,13]. Over the past decades, high-throughput RNA sequencing (RNA-seq) in combination with function annotation/ enrichment and bioinformatics prediction analysis has been widely used to decipher or speculate the essential genes/biological processes/pathways under different conditions at the transcriptomic level and gain a deep and comprehensive understanding of molecular basis underlying the phenotypic/biological differences in plants including rice [14][15][16].
In this study, RNA-seq-based transcriptomic analysis of root samples of YZX and XWX with or without Cd treatment at the booting stage was performed to identify key genes in root Cd response in Yuzhenxiang (YZX, a high Cd accumulation rice cultivar) and Xiangwanxian 12 (XWX, a low Cd accumulation rice cultivar) and genes associated with root Cd accumulation difference between YZX and XWX. RNA-seq results showed that 341 and 161 transcripts were differentially expressed in the roots of YZX or XWX rice after Cd treatment. Moreover, the protein-protein interaction (PPI) networks of dysregulated transcripts in response to Cd exposure in YZX or XWX were established by the STRING database and 7 genes that might play vital roles in the response to Cd pollution were screened out based on the node degrees of proteins in the PPI networks. Additionally, 41 transcripts that function as crucial players in Cd response were filtered out based on GO annotation analysis. Furthermore, 257 transcripts that might be associated with the difference in root Cd response between YZX or XWX were screened out by comparative transcriptomic analysis. Also, we further examined the expression patterns of 10 genes (ferredoxin-1 (ADI1), fructose-1, 6-bisphosphatase (CFBP1), glycine cleavage system H protein (GDCSH), laccase-6 (LAC6), peroxiredoxin Q (Os06g0196300), ribosome-recycling factor (Os07g0570700), peroxisomal membrane protein 11-4 (PEX11-4), probable serine/threonineprotein kinase WNK3, chaperone protein CLPB1, and heat stress transcription factor B-2c (HSFB2C)) by RT-qPCR assay and filtered some key genes related to Cd stimulation based on RNA-seq/RT-qPCR/interaction/GO enrichment/ GO annotation analysis together with previous study outcomes.

Plant Materials and Treatment.
Two Indica rice cultivars (i.e., YZX and XWX) were cultivated in the LT-36VL climatic growth chamber (Percival, Perry, IA, USA) at the Hunan Agricultural Biotechnology Research Institute (Changsha, China). In the preliminary tests, YZX and XWX rice at 4 different growth phases (i.e., seedling, tillering, booting, and grain-filling stages) were cultivated in a hydroponic system containing 2 mg/L of Cd under low-temperature (15-20°C), moderate-temperature (22-27°C), and high-temperature (30-35°C) condition for 48 h. Cd concentration was determined according to the instructions in the Chinese National Standard GB 5009.15-2014.

RNA Sequencing.
In the RNA-seq and real-time quantitative reverse transcription PCR (RT-qPCR) validation experiments, YZX and XWX rice were grown at 22-27°C (moderate temperature) in the hydroponic system with or without 2 mg/L of Cd exposure for 48 h. Next, rice root samples were collected. e samples were divided into 4 groups: AMD group (YZX group with Cd treatment), AMK group (YZX group without Cd treatment), BMD group (XWX group with Cd addition), and BMK group (XWX group without Cd exposure). Each group contains 3 root samples.
RNA was extracted from the root samples using the TRIzol reagent ( ermo Scientific, Waltham, MA, USA) according to the protocols of the manufacturer. e concentration, purity, and quality of RNA were examined by NanoDrop 2000 spectrophotometer ( ermo Scientific). RNA sequencing was carried out as previously described [17]. e cDNA library was constructed using the TruSeq RNA sample preparation kit (Illumina, San Diego, CA, USA) following the protocols of the manufacturer. e mRNA with poly-A structure was enriched by Poly-T oligo-attached magnetic beads (Illumina). e cDNA library was sequenced using the PE150 strategy by Biomarker Technologies Co., Ltd. (Beijing, China). Raw sequencing data were processed and filtered into high-quality clean data. Next, the clean data were aligned to the reference genome of Oryza sativa L. subsp. Indica using HISAT2 and then assembled and quantified using StringTie. e genome version used in the sequence assembly was ASM465v1. (http://plants.ensembl. org/Oryza_indica/Info/Index).

Differential Expression
Analysis. Differential expression analysis was performed using the DESeq R package. Genes were defined to be differentially expressed at the fold change ≥2 and FDR <0.01.

Interaction Network Construction.
e interaction networks of gene-coding proteins were constructed by the STRING database (https://www.string-db.org/).

RT-qPCR Assay.
RNA was reversely transcribed into first-strand cDNA without genomic DNA contamination using TransScript One-Step gDNA Removal and cDNA Synthesis SuperMix (TransGen Biotech Co., Ltd, Beijing, China). Briefly, RNA was coincubated with Oligo (dT) 18 primer for 5 min at 65°C and then immediately placed on an ice bath for 2 min. en, the reaction mixture was coincubated with TS Reaction Mix, TransScript RT/RI Enzyme Mix, and gDNA Remover at 42°C for 30 min and then terminated at 85°C for 5 min. Next, cDNA was amplified and quantified using the ChamQ Universal SYBR qPCR Master Mix (Vazyme Biotech Co., Ltd., Nanjing, China) and specific quantitative primers under the following thermocycling conditions: 95°C for 3 min, 40 cycles of 95°C for 10 s, and 60°C for 30 s. 18sRNA served as the housekeeping gene to normalize the expression of other genes. RT-qPCR was performed on Roche LightCycler 480 II Real-Time PCR Detection System (Basel, Sweden). e quantitative primer sequences are presented in Table 1. e RNA samples in the RT-qPCR assay were identical to those used in RNA-seq experiments.

Statistical
Analysis. Data analysis was performed using GraphPad Prism 7 software (GraphPad Software, Inc., San Diego, CA, USA). Results were displayed as mean ± standard deviation. e difference among groups was examined using one-way ANOVA and Dunnett post hoc test. Statistically significant differences were set at P < 0.05.

Cd Accumulation Difference in the Roots of YZX and XWX at Four Different Growth Stages.
Previous studies have reported that there is a higher Cd accumulation in the first node under panicle, panicle node, and grains of YZX cultivar compared to XWX cultivar [16,18]. We supposed that the difference in grain Cd accumulation in these two rice cultivars might be caused by the discrepancy of root Cd absorptive ability. Moreover, root Cd uptake capacity has been found to be varied at different growth temperatures and growth stages in rice [19][20][21]. To identify the optimal conditions under which the difference of Cd concentration in the roots of YZX and XWX was the most prominent, the concentrations of Cd in the roots of these two rice cultivars at the seedling/tillering/booting/grain-filling stage under low/moderate/high temperature and 2 mg/L of Cd stress were measured. Results showed that there was a noticeable difference in root Cd content at the booting stage when these two rice cultivars were cultured under the moderate-temperature condition ( Figure 1). Hence, YZX and XWX rice were cultivated at moderate temperature and 2 mg/L of Cd stress for 48 h to explore the gene expression alterations in response to Cd treatment in the following experiments.

Identification of Differentially Expressed Transcripts after
Cd Treatment in YZX and XWX. In this study, RNA-seq technology was used to identify Cd response-related transcripts in YZX and XWX and investigate genes that might play vital roles in governing the difference in root Cd accumulation between YZX and XWX. e statistics information of clean data after filtration of RNA-seq raw data is shown in Table 2. e GC content of clean data was approximately 50% (Table 2). Moreover, more than 94% of bases have a recognition accuracy of over 99.9% in the clean data ( Table 2)   Journal of Chemistry 6, respectively. Genes with greater node degrees might play central roles in the response to Cd. Among genes with higher node degrees, 7 genes (Os06g0196300 (OsJ_019618), Os07g0570700 (OsJ_24808), ADI1, GDCSH, HSFB2C, PEX11-4, and CLPB1) were screened out for further RT-qPCR validation. Moreover, interaction analysis of dysregulated genes in the AMD vs. AMK group revealed that there was a complex interaction among OsJ_019618, OsJ_24808, ADI1, GDCSH, and PEX11-4. For instance, OsJ_019618 could interact with OsJ_24808, ADI1, GDCSH,   Tables 3 and 5). Moreover, there was a potential interaction between OsJ_24808 and OsJ_019618 or GDCSH (Supplementary Tables 3 and 5).

Comparative Analysis of Cd-Responsive Genes in XWX
and YZX. GO annotation analysis of the 341 differentially expressed transcripts in the AMD vs. AMK group showed that 37 transcripts were involved in the response to Cd (Supplementary Table 7). Moreover, 23 of the 161 differentially expressed transcripts in the BMD vs. BMK group were identified to be implicated in Cd response based on GO annotation analysis of dysregulated transcripts in the BMD vs. BMK group (Supplementary Table 8). Among these Cdresponsive transcripts, 19 common elements were found in AMD vs. AMK and BMD vs. BMK groups (Figure 2(a) and Supplementary Table 9). ese common transcripts had the same alteration trends in the AMD vs. AMK and BMD vs. BMK groups (Supplementary Table 9). And 18 or 4 transcripts were found to be differentially expressed only in the AMD vs. AMK group or BMD vs. BMK group, respectively (Figure 2(a) and Supplementary Table 9). ese 22 transcripts might have crucial roles in regulating the difference in root Cd accumulation between YZX and XWX at the booting stage. Additionally, we noticed that 43 transcripts were differentially expressed only in the BMD vs. BMK group, but not notably altered in the AMD vs. AMK group (Supplementary Table 10). Also, 214 transcripts were found to be markedly upregulated or downregulated in the AMD vs. AMK group, whereas the difference in the expression of these 214 transcripts was not significant in the BMD vs. BMK group (Supplementary Table 11). e 257 transcripts in Supplementary Tables 10 and 11 are integrated into Supplementary Table 12. We believed that these transcripts might also be related to the difference in root Cd accumulation between XWX and YZX. Also, transcripts with notably different upregulated or downregulated ratios in the BMD vs. BMK and AMD vs. AMK groups might be implicated in root Cd accumulation difference between two rice cultivars, which were not analyzed in our project. KEGG enrichment analysis showed that the transcripts in Supplementary Table 12 mainly participated in the regulation of metabolic pathways, biosynthesis of secondary metabolites, carbon fixation in photosynthetic organisms, carbon metabolism, glutathione metabolism, fructose and mannose metabolism, glycolysis/gluconeogenesis, starch and sucrose metabolism, and glycerophospholipid metabolism (Supplementary Table 13 and Figure 2(b)). For instance, 23, 16, 5, 7, or 5 differentially expressed genes were significantly enriched in KEGG pathways related to metabolism, biosynthesis of secondary metabolites, carbon fixation in photosynthetic organisms, carbon metabolism, or glutathione metabolism, respectively (Supplementary Table 13). e top 20 KEGG pathways are shown in Figure 2(b).

Expression Analysis of 10 Interested Genes by RT-qPCR
Assay. Next, 10 genes including the 7 abovementioned genes with greater node degrees in the PPI networks, 1 interested gene (LAC6), and 2 genes related to heavy metal stress tolerance (WNK3 and CFBP1) (detailed information is shown in the Discussion section) were selected for further RT-qPCR validation. RT-qPCR results showed that the expression levels of ADI1, CFBP1, GDCSH, LAC6, Os06g0196300, Os07g0570700, PEX11-4, and WNK3 were notably downregulated in both AMD vs. AMK and BMD vs. BMK groups (Figure 3(a), 3(b), 3(d), and 3(f )-3(j)). CLPB1 and HSFB2C were highly expressed in the AMD vs. AMK group, but not notably changed in the BMD vs. BMK group (Figures 3(c) and 3(e)). Comparative analysis of RT-qPCR and RNA-seq outcomes disclosed that the alteration trends of these 10 genes were consistent in response to Cd treatment in YZX and XWX rice cultivars, although some differences were not statistically significant (Figures 3(a)-3(j)).

Interaction and Enrichment Analysis of ese 10 Genes.
e relationships of these 10 genes were further deciphered by PPI analysis via the STRING database. Results suggested that 7 gene-encoded proteins (ADI1, CFBP1, PEX11-4, OsJ_019618, OsJ_24808, GDCSH, and CLPB1) could form a potential regulatory network (Figure 4). e interaction Notes. Clean reads: the total number of paired-end reads in clean data; clean bases: the total number of bases in clean data; GC content: the percentages of G and C bases in clean data; ≥Q30%: the percentage of bases with the quality value ≥ 30 (base identification accuracy ≥ 99.9%).
analysis of these 10 proteins is shown in sheet 1 of Supplementary Table 14. Moreover, GO enrichment analysis revealed that 7 genes (LAC6, CFBP1, OsJ_019618, WNK3, OsJ_24808, ADI1, and GDCSH) were enriched in the metabolic process (Supplementary Table 14, sheet 2). Additionally, 6 genes (LAC6, CFBP1, CLPB1, WNK3, OsJ_24808, and ADI1) had a potential ion binding activity (Supplementary Table 14, sheet 3). Furthermore, annotation analysis revealed that 3 genes (LAC6, OsJ_019618, and PEX11-4) were involved in the regulation of the response to oxidative stress (Supplementary Table 14, sheet 4).  Table 6). To explore the association of these genes and Cd response, Venn analysis of genes that could interact with the 10 abovementioned genes and GO-annotated Cd-related genes in Supplementary Tables 7 and 8 was performed. Results showed that OsJ_019618 could interact with 4 genes related to the Cd response (i.e., APX2, HCF136, GLN2, and ALDP). ADI1 could interact with 2 Cd-related genes (i.e., ALDP and GLN2). OsJ_24808, GDCSH, or CLPB1 could interact with Cd-responsive gene TPI, GLN2, or HSP17.4, respectively. ese data further suggested the close association of these genes and Cd response.

Discussion
Rice is a big reservoir of potentially toxic heavy metals in paddy soil-rice systems [22,23]. e contamination of heavy metals in soil not only influences rice growth and grain yields but also seriously threatens the health of animals and humans that consume rice [22][23][24]. In this study, genes related to Cd response and root Cd content difference between YZX and XWX were further examined in the roots of these two rice cultivars with different grain Cd accumulation.
Oxidative stress has been found to be closely linked with Cd tolerance and Cd-induced injury and toxicity [38][39][40]. For instance, molybdenum could relieve Cd-induced toxicity and potentiate Cd tolerance by restricting Cd uptake, reducing oxidative stress, and improving antioxidant defense responses [41]. ese data suggested that these genes might function as crucial players in the response to Cd stress. Cao et al. showed that 706 transcripts were specially expressed in the roots of the mutant Indica rice (Icd1 group), which had an extremely low Cd accumulation in root and shoot under Cd exposure [42]. Moreover, 987 transcripts were specially expressed in the roots of wild-type rice (WT group) under Cd exposure [42]. ese transcripts might be implicated in the difference in Cd accumulation in wild-type and mutanttype Icd1 rice. Combined with our data and Cao's outcomes, we found that 4 transcripts (BGIOSGA004865 (CFBP1),

Conclusion
Taken together, a host of Cd response-related transcripts were identified in the roots of YZX and XWX at the booting stage. Some transcripts that might be associated with the difference in the root Cd-responsive ability of these two rice cultivars were screened out. Moreover, our data suggested that ADI1, CFBP1, PEX11-4, OsJ_019618, OsJ_24808, GDCSH, CLPB1, and WNK3 might play vital roles in regulating Cd response and Cd tolerance. An in-depth insight into the genetic and molecular mechanisms underlying root Cd stress response difference in YZX and XWX at the booting stage might contribute to the development of new strategies that can improve rice yield and reduce Cd harm to rice and humans.
Data Availability e data supporting these results are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest. Figure 4: Interaction network of 10 interested genes. e yellow line represents that the potential interaction between proteins was obtained by text mining. e black line denotes that the potential interaction between proteins was identified by coexpression analysis. e green line means that the potential interaction between proteins was predicted based on gene neighborhood relationships. PPI analysis suggested that PEX11-4 could interact with CFBP1 (FBP1), ADI1, and Os06g0196300 (OsJ_019618). ere was a potential interaction between CFBP1 (FBP1) and PEX11-4/ADI1/Os06g0196300 (OsJ_019618)/GDCSH. GDCSH coud interact with CFBP1 (FBP1)/ADI1/Os06g0196300 (OsJ_019618)/Os07g0570700 (OsJ_24808). Os07g0570700 (OsJ_24808) could interact with GDCSH/Os06g0196300 (OsJ_019618). Os06g0196300 (OsJ_019618) could interact with ADI1/PEX11-4/CFBP1 (FBP1)/GDCSH/Os07g0570700 (OsJ_24808)/CLPB1. ADI1 could interact with PEX11-4/CFBP1 (FBP1)/GDCSH/Os06g0196300 (OsJ_019618). CLPB1 could interact with Os06g0196300 (OsJ_019618).  Table 7: differentially expressed transcripts related to Cd responses in the AMD vs. AMK group, which was annotated by the GO database. Supplementary Table 8: GOannotated Cd-responsive transcripts, which were also differentially expressed in the BMD vs. BMK group. Supplementary Table 9: Venn analysis outcomes for differentially expressed transcripts related to Cd responses in the AMD vs. AMK and BMD vs. BMK groups. Supplementary Table 10: transcripts that were markedly upregulated or downregulated in the BMD vs. BMK group, but not in the AMD vs. AMK group. Supplementary Table 11: transcripts that were notably upregulated or downregulated in the AMD vs. AMK group, but not in the BMD vs. BMK group. Supplementary Table 12: the information of transcripts in Supplementary Tables 10 and 11. Supplementary Table 13: KEGG enrichment analysis for genes in Supplementary  Table 12. Supplementary Table 14