Cadmium (Cd) is a toxic heavy metal that is harmful to the environment and human health. Cd pollution threatens the cultivation of rice (Oryza sativa L.) in many countries. Improving rice performance under Cd stress could potentially improve rice productivity.
In this study, 9 growth traits of 188 different cultivated rice accessions under normal and Cd stress conditions were found to be highly variable during the seedling stage. Based on ~3.3 million single nucleotide polymorphisms (SNPs), 119 Cd-mediated growth response (CGR) quantitative trait loci (QTL) were identified by a genome-wide association study (GWAS), 55 of which have been validated by previously reported QTL and 64 were new CGR loci. Combined with the data from the GWAS, transcriptome analysis, gene annotations from the gene ontology (GO) Slim database, and annotations and functions of homologous genes, 148 CGR candidate genes were obtained. Additionally, several reported genes have been found to play certain roles in CGRs. Seven Cd-related cloned genes were found among the CGR genes. Natural elite haplotypes/alleles in these genes that increased Cd tolerance were identified by a haplotype analysis of a diverse mini core collection. More importantly, this study was the first to uncover the natural variations of 5 GST genes that play important roles in CGRs.
The exploration of Cd-resistant rice germplasm resources and the identification of elite natural variations related to Cd-resistance will help improve the tolerance of current major rice varieties to Cd, as well as provide raw materials and new genes for breeding Cd-resistant varieties.
Cadmium (Cd) is a well-known heavy metal element that has a long decomposition cycle, easy migration, high toxicity, and is difficult to degrade, thereby posing a major threat to environmental safety and harming human health through the food chain in the form of biological enrichment . Cd results in harmful pathogenic, carcinogenic, and mutagenesis effects on the human body, including bone pain caused by chronic Cd poisoning [2, 3]. Cd is not an essential element for plant growth and has adverse effects on growth. When Cd accumulates to a certain extent, plants exhibit toxic symptoms, stunted root growth, and inhibited absorption of water and nutrients, leading to a series of physiological metabolic disorders, including blocked chlorophyll, sugar and protein synthesis, decreased photosynthesis, and altered enzymatic activities, which ultimately lead to reduced yield . Cd destroys plant roots by altering RNA synthesis and proton pump activity, and can replace essential elements, including sulfur, calcium, and magnesium, thereby leading to a shortage of these essential elements . Cd mainly binds to proteins in plants and affects their growth and development by interfering with enzymatic activities, resulting in growth disorders . However, the genetic basis of Cd effects on rice growth remains poorly studied.
Recently, several rice genes have been reported to be involved in the uptake and transport of Cd, including OsIRT1 , OsIRT2, OsNRAMP1 , OsNramp5 , OsHMA3 , OsHMA2 , OsLCT1 , CAL1 , OsCd1 , OsZIP1 , LCD , and OsCCX2 . Among them, OsNramp5 are responsible for transporting Cd from the apoplast to root cells. OsHMA3 in the tonoplast selectively sequestrates Cd into vacuoles and OsHMA2 and CAL1 transport Zn and Cd from the roots to shoots. OsLCT1, a transporter gene for Cd transport in the phloem, regulates grain Cd accumulation. However, to date, only 3 quantitative trait loci (QTL) (i.e., OsHMA3, CAL1, and OsCd1) have been cloned [10, 13, 14]. Three genes (i.e., OsPCS1, OsPCS2, and OsCLT1) have been found to be involved in the chelation of Cd [18,19,20] and OsPCS2 was more strongly activated by Cd than by As(III) [21,22,23]. Additionally, the auxin transporter, OsAUX1, is involved in primary root and root hair elongation and Cd stress responses in rice . Although our knowledge of the genetic control of Cd uptake and accumulation in rice has gradually increased, the molecular characteristics of these genes and other genes remain to be identified.
Phytoremediation techniques for Cd contamination in soil has received increasing attention [25, 26]. It has been reported that different rice varieties have different Cd tolerance [27, 28]. In order to cultivate Cd-tolerant rice varieties, it is important to explore the genetic factors that control Cd-mediated growth responses (CGRs). Genome-wide association study (GWAS) is an effective tool for exploiting elite alleles that control important agronomic traits in germplasm resources [29,30,31,32]. Wu et al.  performed a GWAS using 100 barley accessions and identified 9 QTL for root Cd, 21 for shoot Cd, 15 for grain Cd, and 14 for root-to-shoot Cd translocation. Zhao et al.  detected 14 QTL for Cd accumulation in rice grain from a collection of 312 rice accessions by a GWAS. Recently, the grain Cd accumulation QTL, OsCd1, which belongs to the major facilitator superfamily, was identified by a GWAS based on 127 rice cultivars . Thus, GWAS for CGRs would aid us to identify more excellent natural variations related to Cd-tolerant.
This study investigated the responses of rice growth to high Cd stress. Based on a GWAS and transcriptome analysis, Cd-tolerant germplasm resources were explored and elite natural variations associated with Cd-resistance were identified. Collectively, these results serve as a resource for potentially improving the tolerance of current major cultivars to Cd and provide useful information for cloning candidate genes in future studies.
Phenotypic variation of 9 seedling growth traits under normal and Cd stress conditions
A collection of 188 cultivated rice accessions, consisting of 80 indica, 83 temperate japonica, and 25 tropical japonica, obtained from a mini-core collection of rice in China and other regions of the world, were used in this study [35, 36] (Table S1). This collection represents a large variation of geographical origins and genetic diversity. Various root and shoot traits at the seedling stage were investigated under normal and Cd stress conditions (Table S1), including the sum of root length (SRL), root area (RA), superficial area of roots (SA), root volume (RV), root diameter (RD), maximum root length (MRL), shoot length (SL), shoot weight (SW), and root weight (RW). Compared to normal conditions, 8 traits (i.e., SRL, RA, SA, RV, MRL, SL, SW, and RW) decreased in varying degrees after Cd treatment (Table 1), while 5 traits (i.e., SRL, RA, SA, RV, and RW) were more sensitive, especially RV. These results were consistent with the findings of previous studies that found that Cd stress may inhibit the growth of seedlings [37, 38]. RD did not change under Cd stress, indicating that the growth of root thickness was not affected by Cd. MRL was also less affected by Cd. In order to study the effects of Cd on rice growth, the ratios of the values under Cd stress versus the values under normal conditions (G/N) were used to measure CGRs. Considering the large genetic differences between subspecies, CGR indices among indica, temperate japonica, and tropical japonica rice were compared. The SW and RW indices for temperate japonica were higher than indica (P = 0.038 and 0.040, respectively), while the SL index for temperate and tropical japonica was less than indica (P = 0.0006 and 0.007, respectively). Indica and tropical japonica had similar RA (P = 0.59), SA (P = 0.66), and RV (P = 0.98) values, and the RV value of indica were lower than temperate japonica (P = 0.035). There were no differences in SRL values between subspecies (P > 0.05).
The CGRs of germplasm resources exhibited a great deal of variation (Fig. 1). The variation ranges of SA, SRL, RA, RV, SW, and RW were large, while the variation ranges of RD, MRL, and SL were relatively small. High correlation coefficients were detected for SRL, RA, SA and RV indices (Fig. S1). High correlation coefficients (> 0.65) were also detected between RW and SW, RA, SA, and RV indices. Low correlation coefficients (< 0.2) were detected between RD, MRL and SL indices, suggesting that there were distinct genetic architectural differences among these indices and that a high MRL index did not indicate enhanced SL index.
Identification of QTL for CGRs by GWAS
A GWAS was performed to identify QTL for CGR traits under an expedited mixed linear model approach (EMMAX program) based on ~3.3 million single nucleotide polymorphisms (SNPs) with missing rates of ≤ 30% and a minor allele frequency (MAF) > 0.05 covering the whole rice genome (MAF) [31, 39]. Thirteen, 9, 12, 14, 22, 8, 5, 13 and 23 CGR QTL were obtained for SA, SRL, RA, RV, SW, RD, MRL, SL and RW indices, respectively (Fig. 2, S2 and S3; Table S2). There were several common QTL intervals among the different CGR indices (Table S2). Eight common QTL were detected among SA, RA and SRL, indicating that they were pleiotropic (Table S2). There were twenty-two common QTL between two CGR indices, and 1 common QTL were detected between six CGR indices. However, there were no common QTL detected between MRL and RD indices, which was consistent with a low correlation (Fig. S1).
Comparison of GWAS results with reported QTL/genes
The genomic positions of known Cd-related functional genes with the QTL intervals obtained in this study were compared. Seven genes were co-localized with the associated sites (Figs. 2 and 3; Table S2). OsAUX1, an auxin transporter involved in primary root and root hair elongation, and in rice Cd stress responses , was located in the linkage disequilibrium (LD) region where the qSRL1.1, qRA1.1, qSA1.1, qRV1.3, and qSL1.3 loci were detected (Table S3). CAL1 is a gene that encodes a defensin-like protein acted on by chelating Cd in the cytosol and facilitates Cd secretion to extracellular spaces. It was previously found to lower cytosolic Cd concentration, while driving long-distance Cd transport via xylem vessels . In this study, CAL1 was detected in the LD region of the qSRL2.1, qRA2.1, qSA2.1, qRV2.1, and qSW2.4 loci. OsHMA2, a major transporter of Zn and Cd from roots to shoots, was located in the region of the qSRL6.1/ qRA6.1/ qSA6.1 loci. OsCLT1 encodes a CRT-like transporter 1 and functions as an important component of glutathione homeostasis and Cd tolerance in rice roots . In this study, it was located in the region of the qSW1.1 and qRW1.2 loci. OsCd1, a major facilitator superfamily gene involved in root Cd uptake that contributes to grain accumulation in rice, was located in the region of the qSW3.1 and qRW3.1 loci. OsLCT1, a low-affinity cation transporter for phloem Cd transport in plants, was located in the region of the qRD6.1 and qSL6.6 loci. OsHMA3, a P1B-type heavy metal ATPase that transports Cd into the vacuoles for sequestration , was located in the region of the qRW7.1 locus.
Moreover, the localization of associated sites detected in this study were compared with previously detected Cd-related QTL from the linkage mapping and GWAS of previous studies [14, 27, 42,43,44,45,46]. A total of 55 associated sites were co-localized with 28 previously reported QTL (Fig. 3; Table S2). Specifically, qSRL3.1, qRA3.2, qSA3.2, qRV3.3, and qSW3.2 were all located in the regions of qSRR3 (SRR, shoot/root ratio of Cd concentration) ; qSRL8.1, qRA8.1, qSA8.1, and qRV8.1, were located in the region of qRDW8.2 (RDW, root dry weight); qSRL9.1, qRA9.1, and qSA9.1, were located in the region of qLR-9; qRA3.4, qSA3.4, and qRV3.5, were located in the region of qCd3 (Cd, grain Cd content). These results indicate the important roles of repeatable QTL detected by linkage mapping and GWAS, and these QTL may be related to each other.
Elite alleles in 7 cloned genes for CGR
To better understand the natural variation of 7 known genes, haplotype analyses were performed using all of the non-synonymous SNPs within their ORFs and promoters. Twelve SNPs in the promoter and 2 non-synonymous SNPs (i.e., Chr1_36999122 and Chr1_36999123) were detected at the OsAUX1 locus, which play important roles in response to Cd stress and root development . Four distinct haplotypes, including 2 major haplotypes, Hap.1 and Hap.2, were identified based on the 14 SNPs in cultivated rice and exhibited a large genetic difference between indica and japonica (Fig. 4a). Hap.1 was mostly present in temperate and tropical japonica, while Hap.2 was mostly present in indica. Significant differences in CGR indices (i.e., RV, SW and RW) were detected between Hap.1 and Hap.2 (Fig. 4e). Accessions with the Hap.1 genotype had higher RV, SW and RW indices than accessions of other haplotypes, indicating Cd tolerance in seedlings.
A total of 8 SNPs (MAF ≥ 0.05) in OsCd1 identified 4 haplotypes (Fig. 4b). The vast majority of indica accessions were Hap.1, while japonica accessions were Hap.2 or Hap.3. Accessions carrying Hap.1 showed lower RV, SW and RW indices than japonica carrying Hap.2 or Hap.3 (Fig. 4f). Furthermore, recent studies revealed that accessions with OsCd1D449 have higher grain Cd concentrations compared to those with OsCd1V449 , which also exhibit a lower Cd transport ability, suggesting that the missense mutation, V449D, is responsible for the divergence of rice CGRs between indica and japonica. Interestingly, 2 indica accessions with higher CGR indices were Hap.3, indicating that elite japonica alleles had been introgressed into indica accessions through breeding (Fig. 4b).
Recently, Liu et al.  reported that variations in OsHMA3 contributed to differential grain Cd accumulation between indica and japonica. Here, 6 non-synonymous SNPs (i.e., Chr7_7408467, Chr7_7408183, Chr7_7407112, Chr7_7406920, Chr7_7406728, and Chr7_7406505) and 8 SNPs in the promoter of OsHMA3 revealed 2 major haplotypes, Hap.1 and Hap.2 (Fig. 4c). Hap.1 contained most indica rice and was associated with lower RV, SW, and RW indices. In contrast, Hap.2 contained most temperate and tropical japonica rice and was associated with higher RV, SW and RW indices (Fig. 4g). Similarly, 3 major haplotypes were detected in OsCLT1 with Hap.1 was associated with higher RV, SW, and RW indices (Fig. 4d and h); its frequency was high in temperate and tropical japonica, but low in indica.
Three of 16 SNPs (i.e., Chr6_29481363, Chr6_29481366 and Chr6_29481398) in the promoter region of OsHMA2 distinguished 2 distinct major haplotypes, Hap.1 and Hap.2 (Fig. S4a). The RV, SW and RW indices of Hap.1 were significantly higher than Hap.2 (Fig. S4d). Additionally, a premature termination codon mutation, K153* of Hap.3 was detected in OsHMA2. Hap.3 also had higher RV, SW and RW indices, but lower frequency.
Nine SNPs in the promoter and 1 non-synonymous SNP (i.e., Chr2_25190881) of CAL1 revealed 3 major haplotypes (Fig. S4b). Most japonica carried Hap.2 and the difference in the number of indica rice of three major haplotypes was small. Hap.1 had higher CGR indices than Hap.3 in indica (Fig. S4e). However, there was no significant phenotypic difference between Hap.3 and Hap.4, indicating that the non-synonymous SNP could not be the cause of the variation affecting RV, SW and RW indices. These results suggest that the codon sequences in CAL1 for maintaining protein function were highly conserved and that the phenotypic differences among haplotypes could be caused by differences in the expression levels of germplasm resources, which was consistent with a previous study on CAL1 that found that it positively regulated the Cd contents in the leaves and xylem sap .
Only 2 haplotypes were detected by 13 non-synonymous SNPs and 10 SNPs in the promoter of OsLCT1. All of the indica accessions were belonged to Hap.2 (Fig. S4c). Both Hap.1 and Hap.2 had temperate japonica and tropical japonica rice accessions. RV and SW indices of Hap.1 were higher than Hap.2 (Fig. S4f).
Differentially expressed rice genes in response to Cd stress from RNA-Seq data
To investigate the transcriptomic response of rice to Cd stress, an RNA-Seq analysis was performed on 3 Cd-tolerant varieties (CTVs) and 3 Cd-sensitive varieties (CSVs) under normal and Cd stress conditions (Fig. S5; Table S1). RNA-Seq data of 3 biological replicates were combined to screen the common DEGs in each variety. RNA-Seq data of 3 Cd-tolerant varieties and 3 Cd-sensitive varieties were combined to screen the common DEGs in all CTVs and CSVs, respectively. A total of 2,528 differentially expressed genes (DEGs) that respond to Cd stress were identified. More DEGs were identified in CTVs (1,068 up-regulated and 773 down-regulated genes) than CSVs (712 up-regulated and 729 down-regulated genes) under Cd stress (Fig. S5). According to the gene ontology (GO) enrichment analysis, the up-regulated DEGs in CSVs were highly enriched in oxidation reduction, zinc ion transmembrane transport, and tetrapyrrole, heme and iron ion binding, while down-regulated DEGs were enriched in oxidoreductase, carbohydrate metabolic processes, hydrolyzing O-glycosyl compounds and apoplast (Fig. S6; Table S4). Genes enriched in the GO terms oxidation reduction, aminoglycan catabolic process, heme binding and iron ion binding were up-regulated, while genes enriched in the GO terms photosynthesis, oxidoreductase activity, thylakoid and photosynthetic membrane were down-regulated in CTVs (Fig. S7; Table S5). Interestingly, the down-regulated genes in CSVs enriched in hydrolase activity, hydrolyzing O-glycosyl compounds, endopeptidase inhibitor activity, and peptidase inhibitor activity were up-regulated in CTVs (Figs. S6 and S7).
There were 476 co-up-regulated and 278 co-down-regulated genes among CTVs and CSVs (Fig. 5a and b; Table S6). The GO analysis revealed that these common genes were enriched in oxidation reduction, tetrapyrrole binding, heme binding, metal ion transmembrane transport, and phenylpropanoid metabolic processes (Fig. 5c and d). For CTV-specific DEGs, up-regulated genes were enriched in oxidation reduction, iron ion binding, and heme binding, while down-regulated genes were enriched in photosynthesis, thylakoid, and photosynthetic membrane (Fig. S8; Table S7). Chitin has selective permeability in material transport and plays a role in attracting cations [49,50,51,52]. Interestingly, genes enriched in the GO term, chitin metabolic process, were up-regulated in CTVs (Fig. S8), suggesting that chitin may be associated with CGRs in CTVs. The term oxidation reduction in biological process and the term oxidoreductase activity in molecular function were the most significant GO terms in CTV-specific up-regulated terms and they contained 29 cytochrome P450, 8 peroxidase precursor, 7 dehydrogenase, 4 flavin-containing monooxygenase family protein, and 3 NADP-dependent oxidoreductase genes, which could help to reduce reactive oxygen species (ROS)-relevant oxidative damage (Fig. S8; Table S7). The terms photosynthesis in biological process and the term membrane in cellular component in CTV-specific down-regulated terms contained 25 genes related to photosynthesis and 13 transporters genes, which may function in plant growth and Cd transport.
Determination of candidate genes within CGR QTL by integrated genomic and transcriptomic analyses
The CGR loci identified by GWAS provided important clues for understanding the genetic architecture of the observed variations in rice growth under Cd stress. To identify the candidate genes responsible for each CGR locus, based on the LD decay values in indica and japonica (123 kb and 167 kb, respectively), all of the genes within 200 kb of the most significant SNPs were extracted and the data derived from the RNA-Seq analysis, the gene ontology (GO) Slim analysis , and their annotations and functions of homologous genes were considered. Firstly, there were 1,594 genes in intervals of 119 CGR QTL, including 921 clearly annotated genes (Table S8). Eighty-eight candidate genes from 74 CGR QTL were obtained by the integrated genomic and transcriptomic analyses (Table S9). Enriched GO terms of common genes between the GWAS and RNA-Seq analysis included response to stimulus (GO:0050896), transporter activity (GO:0005215), electron carrier activity (GO:0009055), and iron ion binding (GO:0005506) (Fig. S9). Secondly, in order to identify Cd-related membrane transporters in rice, a GO Slim analysis was conducted and the genes related to membrane and transport were selected as important candidate genes. Altogether, 53 candidate genes from 48 CGR QTL located on 12 chromosomes were found to be associated transport and membrane (Fig. S10; Table S10). Thirdly, according to the function of gene annotations and homologous genes, 12 candidate genes in 13 QTL intervals were identified (Table S11). By applying these approaches, a total of 148 genes from 85 CGR QTL were identified, being selected as potential candidate genes for each of the loci controlling CGR traits in rice (Table S11). Among them, 7 cloned Cd-related genes (i.e., OsCd1, OsAUX1, OsHMA2, OsHMA3, OsCLT1, CAL1 and OsLCT1) were identified for the CGR genes. More importantly, some reported genes also identified in this study (i.e., OsATX1, OsBOR3, OsTIP2, OsZIP8, OsVAMP714, OsHMA5, OsCAX1b, QHB, OsABI5, UbL402, OsGA20ox1/qEPD2/GNP1 and OsUGE1) were found to play certain roles in CGRs.
Haplotype analyses of 6 QTL genes for CGRs
Among the closely associated candidate genes, several reported genes had no evidence regarding their functions in the control of CGR traits. Six of these genes (i.e., OsGSTU31 (LOC_Os10g38189), OsGSTU6.1 (LOC_Os10g38340), OsGSTU6.2 (LOC_Os10g38360), OsGSTU21 (LOC_Os10g38150), OsGSTU32 (LOC_Os10g38314), and OsATX1) were selected for functional analyses as case studies.
One QTL (i.e., qSW10.2/qRW10.2) controlling both SW and RW on chromosome 10 was identified by GWAS and 38 genes were in the interval. Among them, 5 candidate genes, all encoding glutathione S-transferase (GST), were up-regulated under Cd stress according to the transcriptomic data (Table S11). Using the MSU Rice Genome Annotation (Osa1) Release 7 for annotated genes, 20 GST genes in or around this QTL interval were identified, 12 of which were co-up-regulated under Cd stress (Tables S4 and S5). The tolerance of plant cells to toxic elements is highly dependent on glutathione metabolism. First, GST proteins indirectly act on Cd accumulated reactive oxygen species (ROS) by maintaining the antioxidant flavonoid pool . GSTs detoxify cytotoxic substrates and ameliorate their toxicity by catalyzing the conjugation of glutathione to substrates . These findings suggest that this gene cluster could play an important role in the CGRs of rice. Further investigations were conducted to analyze the haplotypes of 5 GST genes within this QTL. Four SNPs in the promoter and 6 non-synonymous SNPs (i.e., Chr10_20449102, Chr10_20449187, Chr10_20449235, Chr10_20449241, Chr10_20449359, and Chr10_36999971) were detected in OsGSTU31. Three distinct haplotypes, including 2 major haplotypes, Hap.2 and Hap.3, were identified based on the 10 aforementioned SNPs and exhibited large genetic differences between indica and japonica (Fig. 6a). Hap.2 contained most temperate and tropical japonica, while Hap.3 contained the most indica. Significant differences in CGR indices were detected between Hap.2 and Hap.3 (Fig. 6e). Accessions with the Hap.2 genotype had higher RV, SW, and RW indices than accessions of other haplotypes and exhibited Cd tolerance in seedlings. Four non-synonymous SNPs and 11 SNPs in the promoter of OsGSTU6.1 revealed 2 major haplotypes, Hap.1 and Hap.3 (Fig. 6b). Hap.1 contained most japonica and was associated with higher RV, SW, and RW indices. In contrast, Hap.3 contained most Indica and was associated with lower RV, SW, and RW indices (Fig. 6f). Similarly, 2 major haplotypes were detected in OsGSTU6.2, OsGSTU21, and OsGSTU32 with Hap.1, Hap.3, and Hap.2 associated with higher RV, SW, and RW values, respectively (Fig. 6c, g and S11). Their frequency was high in temperate and tropical japonica, but low in indica.
There were 14 annotated genes in the qSRL8.2 QTL interval (Table S5). A heavy metal-associated domain containing protein, OsATX1 (LOC_Os08g10480), was detected, which was reported to exhibit the heterologous expression of OsATX1 in a Cd-sensitive mutant of yeast (Saccharomyces cerevisiae), Δycf1, which increased the tolerance to Cd by decreasing their respective concentrations in transformed yeast cells . Interestingly, no polymorphism was detected in the code region of OsATX1 among rice accessions, suggesting that the sequences in OsATX1 for maintaining protein function were highly conserved. Twelve SNPs in the promoter region of OsATX1 revealed 2 major haplotypes and Hap.1 was predominant and associated with higher RV, SW, and RW indices than Hap.4 (Fig. 6d, h). These results suggest that phenotypic differences among different haplotypes could be caused by differences in the expression levels of OsATX1 among rice accessions.
The growth and development of rice are inhibited under Cd stress. Whether low grain Cd rice or high shoot Cd rice varieties for phytoremediation, both of these varieties were required to grow well under Cd stress. Using rice mini-core germplasms to systematically study the differences in rice growth responses to Cd, Cd-tolerant varieties and Cd-sensitive varieties were screened in this study in order to identify Cd-tolerant genes. To the best of our knowledge, this was the first attempt to conduct a GWAS for CGR traits. In total, 119 QTL for CGRs were identified, 55 of which overlapped with previously reported Cd-related QTL. Based on an integrated analysis strategy, a total of 148 candidate genes for CGRs, including 7 cloned Cd-related genes, were identified. Elite alleles of 13 genes were investigated and will serve as potential candidates for the genetic improvement of Cd-tolerant rice.
The complexity of genetic control of CGR traits and exploration of candidate genes
Cd is a toxic heavy metal that inhibits the growth of roots and shoots, reduces leaf and tiller number, physiologically impairs photosynthesis and mitochondrial respiration, and results in DNA degradation and cell death [44, 56]. Various traits exhibited different responses to Cd. The SRL, RA, SA, RV, and RW were considerably different between normal and Cd stress conditions, while RD and MRL exhibited minor change (Table 1). Different responses to Cd were detected among accessions from different subgroups. Both temperate and tropical japonica had higher RW and SW indices than indica, while RV index of indica and tropical japonica were lower than temperate japonica (Table 1). In this study, > 100 QTL were identified for CGR traits (Table S2), suggesting that the genetic regulation of CGRs is very complex and many genes play an important role in the growth of rice under Cd stress. However, only 3 Cd-related QTL (i.e., OsHMA3, CAL1, and OsCd1) were cloned. The findings of this study provide important information for cloning novel CGR genes in future studies.
Considering the relatively low LD decay of rice, 1 associated locus in this study was defined as a 200 kb region containing > 10 genes ; therefore, it was rather difficult to pinpoint the causal genes for these loci. However, the combined use of QTL information, expression profiles, GO Slim analysis, and prediction of gene functions could help uncover candidate genes, just as candidate genes were uncovered in this study. Based on the RNA-Seq data, GO Slim analysis, and gene annotations, 88 candidate genes for 74 CGR QTL, 53 for 48 CGR QTL, and 12 for 13 CGR QTL were identified (Tables S5, S6, and S7). Furthermore, many known genes were located in these CGR QTL. Collectively, this study provided a relatively comprehensive analysis of the genetic architecture of CGR in rice.
Favorable natural haplotypes/alleles for the improvement of Cd tolerance in rice
The mining of more favorable alleles of CGR genes is required in order to achieve ideal Cd tolerance in rice. At the single-gene level, favorable haplotypes for 13 CGR genes were identified (Figs. 4 and 6, S4, and S11), including 7 known Cd-related genes and 6 novel genes. Most japonica accessions carried superior OsHMA3Hap2, OsHMA2Hap1, OsCLT1Hap1, OsCd1Hap3, OsAUX1Hap1, OsATX1Hap1, LOC_Os10g38150Hap1, LOC_Os10g38189Hap2, LOC_Os10g38314Hap2, LOC_Os10g38340Hap1, and LOC_Os10g38360Hap1 haplotypes. Additionally, elite japonica alleles of these genes can be considered as primary alternatives for improving Cd-tolerance in indica. Because the expression of OsHMA3, CAL1, and OsAUX1 were induced by Cd [13, 24, 41, 48], natural variations in the promoter region of rice accessions were likely important functional SNPs associated with CGR traits. The results of the haplotype analysis indicated that major haplotypes, which consist of non-synonymous SNPs in the coding sequence (CDS) and/or promoter regions within single loci, represent important allelic diversity of QTL underlying the variation of CGRs in rice populations.
OsHMA3 is a cadmium transporter located in the vacuolar membrane of rice roots, belonging to the heavy metal ATPase (HMA) family. OsHMA3 could transport Cd2+ into vacuoles to isolate it and reduce the transport of Cd2+ to the aboveground, thus reducing cadmium toxicity [41, 58]. Four haplotypes were identified in rice diversity populations. There were five non synonymous SNPs and 8 SNPs on the promoter between haplotype 1 and 2 (Fig S4). Among them, F229L, V323G were located in A-Domain, V550I, S614G, C678R in ATP- binding domain, and V752A in Metal- binding domain of OsHMA3 (Yan et al. 2016). The SNPs in the promoter had been proved to affect the transcriptional activity . Therefore, the differential functions of OsHMA3 between Hap.1 and Hap.2 could be attributed to the eight nucleotide changes occurring in the promoter region. OsHMA2 is a major Zn and Cd transporter in rice roots and shoots. It is homologous with OsHMA3, a heavy metal ATPase. Its C-terminal region is essential for Cd transport in shoot . In rice accessions, Hap.1 had significantly higher RV, SW and RW indices than Hap.2, indicating that the three SNPs (i.e., Chr6_29481363, Chr6_29481366 and Chr6_29481398) could play important role in response to Cd for rice. Hap.3 of OsHMA2 had a premature termination codon mutation and resulted in a truncated protein product, which could affect its cadmium transport function. Ten indica rice and 4 temperate japonica rice in Hap.3 were primary improved varieties, only one was recently improved varieties (Fig S4g and Table S12). Compared with Hap.1 and Hap.2, Hap.3 had relatively fewer varieties. These results suggested that Hap.3 was just found recently and had already been used in recent rice breeding, but it was rarely used at present. Future breeding with Hap. 3 is a potential alternative for improving Cd tolerance in current elite varieties.
Natural variations in 6 QTL genes served important roles in CGR
ATX1 of Saccharomyces cerevisiae encodes an 8.2-kDa peptide, which has significant similarity with many bacterial metal transporters. ATX1 is involved in the transport and distribution of copper and protects cells from the toxicity of superoxide anion and hydrogen peroxide . The resistance of Saccharomyces cerevisiae atx1 deletion strain to Cd2+ was higher than that of wild type and ATX1 can specifically bind Cd2+ . The two copper chaperones of Arabidopsis thaliana, namely Antioxidant Protein1 (ATX1) and ATX1-Like Copper Chaperone (CCH) (CCH), share high sequence homology . Arabidopsis Antioxidant Protein1 (ATX1) plays an essential role in copper (Cu) homeostasis, conferring tolerance to both excess and subclinically deficient Cu . The high affinity of Cd for thiols might be the reason that Cd2+ also can bind the Cu-binding motif MXCXXC, which was required for the physiological function of ATX1 . Knockout of OsATX1 resulted in an increase of Cu concentration in roots, while overexpression of OsATX1 decreased root Cu concentration but increased shoot Cu accumulation. The concentrations of Cu in developing tissues, including upper nodes and internodes, younger leaf blades, and leaf sheaths, were increased significantly in OsATX1-overexpressing plants and decreased in osatx1 mutants compared with the wild type , indicating that rice varieties with high OsATX1 expression might show more sensitive to Cd. In fact, overexpression of OsATX1 increased Cd accumulation in the shoots . Haplotype analysis showed that OsATX1 showed significant indica-japonica differentiation. The high RV, SW, and RW indices of japonica rice might be resulted from their low OsATX1 expression (Fig S12). Thus, gene editing of OsATX1 promoter could improve cadmium tolerance of cultivated rice.
The induction of oxidative stress by Cd was one of the major alterations in plant cells . When redox imbalance occurs, membrane lipids, proteins, and nucleic acids were oxidized, which in turn affects plant metabolism. Glutathione functioned as an antioxidant and moderated the redox imbalance induced by toxic metal accumulation in Arabidopsis . Moreover, OsGST4 played an important role during oxidative stress by ROS-scavenging in rice . In this study, the RNA-Seq data revealed that many genes enriched in the GO term, oxidation reduction, responded to Cd stress (Fig. 4c and d). Cd may be formed as a complex with phytochelatins or glutathione, and is subsequently transported to the vacuoles through an unidentified ABC transporter. Glutathione is used to detoxify xenobiotics through GSTs. GST family genes were previously found to play roles in Cd resistance and accumulation of pak choi . In this study, a GST gene cluster on chromosome 10 was identified by the integrated genomic and transcriptomic analyses (Table S11). Twelve of these genes were up-regulated under Cd stress. The haplotype analysis revealed that all 5 OsGSTUs in the QTL interval showed indica-japonica differentiation and their japonica haplotypes had higher CGR indices (Figs. 6a–c and S11). These results suggest that the GST gene cluster played an active role in detoxification and the japonica alleles of the 5 GSTs enabled rice to grow better under Cd stress.
The CGRs of germplasm resources exhibited a great deal of variation and the influence of Cd on the growth of indica rice was greater than that of japonica rice. A total of 148 genes from 85 CGR QTL were obtained by comprehensive analyses. Natural elite haplotypes/alleles of 13 genes, including 7 cloned Cd-related genes and 6 novel genes, are identified and will serve as potential candidates for the genetic improvement of Cd-tolerant rice. The cultivation of novel Cd-tolerant varieties also helps to ensure a stable rice yield.
Plant materials and phenotyping
A total of 188 rice accessions from around the world were used for evaluating CGRs in this study. Hydroponic experiments were performed at the greenhouse of Agricultural Genomics Institute in Shenzhen, China during the summer of 2018. All of the seeds were soaked in deionized water at 37°C in the dark for 2 d, then transferred to a net floating on deionized water for 5 d. Seedlings were cultured in a half-strength Kimura B nutrient solution (pH, 5.4) with the following composition (μM): 90 KH2PO4, 270 MgSO4, 180 (NH4)2SO4, 90 KNO3, 180 Ca(NO3)2, 3 H3BO3, 0.5 MnCl2, 1 (NH4)6Mo7O24, 0.4 ZnSO4, and 20 Fe(III)-EDTA. Solutions were changed 3 times per week and the pH was adjusted to 5.4 every day. Plants were grown in a greenhouse with natural sunlight at 30°C during the day and 25°C at night . In order to compare the growth of normal and Cd-treated seedlings, 15-d-old plants were exposed to Cd stress in a 1/2 Kimura B nutrient solution containing 50 μM CdCl2 for 7 d, solutions were renewed every 2 d. The experiment was conducted three times. Five plants per variety were sampled and the sum of root length (SRL), root area (RA), superficial area of root (SA), root volume (RV), root diameter (RD), maximum root length (MRL), shoot length (SL), shoot weight (SW), and root weight (RW) were measured.
The sequence data of all of the rice accessions for GWAS were obtained from the 3,000 Rice Genomes Project . The SNP data were filtered out with minor allele frequencies (MAF) < 0.05 and missing rates > 30%. The efficient mixed model analysis (EMMA) feature of the EMMA eXpedited (EMMAX) software was employed for GWAS . The significance threshold was calculated using the formula “-log10(1/the effective number of independent SNPs)” as described previously , and effective numbers of independent SNPs were determined by PLINK to be 398107 in this population . The suggestive P values was 2.5 × 10−6. Finally, the threshold was set at −log(P) = 5.6 to identify significant association signals. Based on the LD decay values in indica and japonica rice (123 kb and 167 kb, respectively), several SNPs passing the threshold on the same chromosome were clustered as one associated locus with a region of < 200 kb. All genes located within the candidate region were extracted .
RNA-Seq and GO enrichment analyses
Three Cd-tolerance cultivars (CTVs) (CX47, Yungeng 23 and IRIS_313_9050) and 3 Cd-sensitive cultivars (CSVs) (GUI630, ALBANIA_SPECIES and BAXIANG) based on phenotyping results were selected for the RNA-Seq analysis (Fig. S5; Table S1). 15-d-old plants of the 6 varieties were planted under normal and Cd stress conditions for 12 h, and then the roots were sampled for RNA extraction. RNA was extracted by preparing samples using a Micro RNA Extraction kit (Axygen, NY, USA) and reverse transcribed into cDNA using a ReverTra® Ace qPCR-RT kit (TOYOBA, Osaka, Japan). RNA-Seq libraries were prepared using 3 biological replicates for each variety and sequenced separately using a Hiseq Xten sequencer. TOPhat2 software  was used to align the cleanup data to the reference genome MSU V7.0  and gene expression was quantified by fragment per kilobase million (FPKM) using the Cufflinks  default parameters. A false discovery rate (FDR) < 0.05 and absolute value of log2 ratio ≥ 1 were used to identify differentially expressed genes (DEGs) as previously described . GO enrichment analyses were conducted using agriGO v2, an online GO analysis toolkit and database for agricultural communities .
Correlation coefficients between the measured traits were calculated using the R package PerformanceAnalytics  as described in Note S1. The violin map for haplotype analysis was also constructed in R. Data were statistically analyzed and multiple comparisons were made using Duncan’s multiple range test as described . P values of less than 0.05 were considered to indicate statistical significance. Different letters denote significant differences. Statistical calculations were performed using Microsoft Excel 2010.
Nakanishi H, Ogawa I, Ishimaru Y, Mori S, Nishizawa NK. Iron deficiency enhances cadmium uptake and translocation mediated by the Fe2+ transporters OsIRT1 and OsIRT2 in rice. Soil Sci Plant Nutr. 2006;52(4):464–9.
Satoh-Nagasawa N, Mori M, Nakazawa N, Kawamoto T, Nagato Y, Sakurai K, et al. Mutations in rice (Oryza sativa) heavy metal ATPase 2 (OsHMA2) restrict the translocation of zinc and cadmium. Plant Cell Physiol. 2012;53(1):213–24.
Uraguchi S, Kamiya T, Sakamoto T, Kasai K, Sato Y, Nagamura Y, et al. Low-affinity cation transporter (OsLCT1) regulates cadmium transport into rice grains. Proc Natl Acad of Sci U S A. 2011;108(52):20959–64.
Liu XS, Feng SJ, Zhang BQ, Wang MQ, Cao HW, Rono JK, et al. OsZIP1 functions as a metal efflux transporter limiting excess zinc, copper and cadmium accumulation in rice. BMC Plant Biol. 2019;19(1):283.
Das N, Bhattacharya S, Bhattacharyya S, Maiti MK. Identification of alternatively spliced transcripts of rice phytochelatin synthase 2 gene OsPCS2 involved in mitigation of cadmium and arsenic stresses. Plant Mol Biol. 2017;94(1-2):167–83.
Shimpei H. Masato, Kuramata, Tadashi, Abe, Hiroki, Takagi, Kenjirou, Ozawa: Phytochelatin synthase OsPCS1 plays a crucial role in reducing arsenic levels in rice grains. Plant J. 2017;91(5):840-8.
Shimpei U, Nobuhiro T, Christian H, Kaho A, Naoko O-O. Phytochelatin Synthase has Contrasting Effects on Cadmium and Arsenic Accumulation in Rice Grains. Plant Cell Physiol. 2017;58(10):1730-42.
Shinichi Y, Yosuke U, Aya M, Kumiko O, Toru M. Rice phytochelatin synthases OsPCS1 and OsPCS2 make different contributions to cadmium and arsenic tolerance. Plant Direct. 2018;2(1):e00034.
Yu C, Sun C, Shen C, Wang S, Liu F, Liu Y, et al. The auxin transporter, OsAUX1, is involved in primary root and root hair elongation and in Cd stress responses in rice (Oryza sativa L.). Plant J. 2015;83(5):818–30.
Santos MS, Pedro CA, Goncalves SC, Ferreira SM. Phytoremediation of cadmium by the facultative halophyte plant Bolboschoenus maritimus (L.) Palla, at different salinities. Environ Sci Pollut R. 2015;22(20):15598–609.
Xiuyan L, Sunlu C, Mingxue G, Zheng Y, Peng X. Association Study Reveals Genetic Loci Responsible for Arsenic, Cadmium and Lead Accumulation in Rice Grain in Contaminated Farmlands. Front Plant Sci. 2019;10:61.
Wang Q, Zeng X, Song Q, Sun Y, Lai Y. Identification of key genes and modules in response to Cadmium stress in different rice varieties and stem nodes by weighted gene co-expression network analysis. Sci Rep. 2020;10(1):9525.
Huang X, Wei X, Sang T, Zhao Q, Feng Q, Zhao Y, et al. Genome-wide association studies of 14 agronomic traits in rice landraces. Nat Genet. 2010;42(11):961.
Yano K, Yamamoto E, Aya K, Takeuchi H, Lo PC, Hu L, et al. Genome-wide association study using whole-genome sequencing rapidly identifies new genes influencing agronomic traits in rice. Nat Genet. 2016;48(8):927-34.
Yu J, Xiong H, Zhu X, Zhang H, Li H, Miao J, et al. OsLG3 contributing to rice grain length and yield was mined by Ho-LAMap. BMC Biol. 2017;15(1):28.
Zhao J, Yang W, Zhang S, Yang T, Liu Q, Dong J, et al. Genome-wide association study and candidate gene analysis of rice cadmium accumulation in grain in a diverse rice collection. Rice. 2018;11(1):61.
Jianping Y, Jinli M, Zhanying Z, Haiyan X, Xiaoyang Z. Alternative splicing of OsLG3b controls grain length and yield in japonica rice. Plant Biotechnol J. 2018;16(9):1667-78.
Chen Z, Feng Y, Wang Y, Li Y, Liu Q, Xu S, et al. Study on the growth and photosynthetic characteristics of wheat seedlings under [C4mim][OAc] (1-butyl-3-methyl-imidazolium acetate) with Cd2+ stress. B Environ Contamin Tox. 2015;94(5):627–32.
Kollarova K, Kamenicka V, Vatehova Z, Liskova D. Impact of galactoglucomannan oligosaccharides and Cd stress on maize root growth parameters, morphology, and structure. J Plant Physiol. 2018;222:59–66.
Miyadate H, Adachi S, Hiraizumi A, Tezuka K, Nakazawa N, Kawamoto T, et al. OsHMA3, a P1B-type of ATPase affects root-to-shoot cadmium translocation in rice by mediating efflux into vacuoles. New Phytol. 2011;189(1):190–9.
Huang Y, Sun C, Min J, Chen Y, Tong C, Bao J. Association mapping of quantitative trait loci for mineral element contents in whole grain rice (Oryza sativa L.). J Agric Food Chem. 2015;63(50):10885–92.
Liu CL, Gao ZY, Shang LG, Yang CH, Ruan BP, Zeng DL, et al. Natural variation in the promoter of OsHMA3 contributes to differential grain cadmium accumulation between Indica and Japonica rice. J Int Plant Biol. 2020;62(3):314-29.
Bhatnagar A, Sillanp M. Applications of chitin- and chitosan-derivatives for the detoxification of water and wastewater--a short review. Adv Colloid Interfac. 2009;152(1-2):26–38.
Hernández LE, Sobrinoplata J, Monteropalmero MB, Carrascogil S, Florescáceres ML, Ortegavillasante C, et al. Contribution of glutathione to the control of cellular redox homeostasis under toxic metal and metalloid stress. J Exp Bot. 2015;66(10):2901.
Lin SJ, Culotta VC. The ATX1 gene of Saccharomyces cerevisiae encodes a small metal homeostasis factor that protects cells against reactive oxygen toxicity. Proc Natl Acad of Sci U S A. 1995;92(9):3784–8.
Xu N, Chu Y, Chen H, Li X, Wu Q, Jin L, et al. Rice transcription factor OsMADS25 modulates root growth and confers salinity tolerance via the ABA–mediated regulatory pathway and ROS scavenging. PLoS Genet. 2018;14(10):e1007662.
Wu X, Chen J, Yue X, Wei X, Zou J, Chen Y, et al. The zinc-regulated protein (ZIP) family genes and glutathione s-transferase (GST) family genes play roles in Cd resistance and accumulation of pak choi (Brassica campestris ssp. chinensis). Ecotox Environ Safe. 2019;183:109571.
Li MX, Yeung JMY, Cherny SS, Sham PC, Li MX, Yeung JM, et al. Evaluating the effective numbers of independent tests and significant p-value thresholds in commercial genotyping arrays and public imputation reference datasets. Hum Genet. 2011;131(5):747–56.
Wang X, Ning Y, Zhang P, Li C, Zhou R, Guo X. Hair multi-bioelement profile of Kashin-Beck disease in the endemic regions of China. J Trace Elem Med Biol. 2019;54:79-97.
Zhao Y, Zhang H, Xu J, Jiang C, Yin Z, Xiong H, et al. Loci and natural alleles underlying robust roots and adaptive domestication of upland ecotype rice in aerobic conditions. PLoS Genet. 2018;14(8):e1007521.
We thank LetPub (www.letpub.com) for its linguistic assistance during the preparation of this manuscript. We also thank Haiyan Xiong (University of Cambridge) and Muhammad Abdul Rehman Rashid (University of Agriculture Faisalabad, Pakistan) for critical reading and suggested revision of the manuscript.
This work was supported by grants from the Shenzhen Science and Technology Program (Nos. KQTD2016113010482651, 2017050414212249, JCYJ20170303154319837 and JCYJ20170303154506881), National Postdoctoral Program for Innovative Talents (BX201600151), China Postdoctoral Science Foundation (No. 2018M631616 and 2019T120150), National Natural Science Foundation of China (Grant Nos. 31901514). The funding numbers provided the financial support to the research programs, but didn’t involve in work design, data collection, analysis and preparation of the manuscript.
Jianping Yu, Chaolei Liu and Hai Lin contributed equally to this work.
Authors and Affiliations
Shenzhen Branch, Guangdong Laboratory of Lingnan Modern Agriculture, Genome Analysis Laboratory of the Ministry of Agriculture and Rural Affairs, Agricultural Genomics Institute at Shenzhen, Chinese Academy of Agricultural Sciences, Shenzhen, China
Jianping Yu, Hai Lin, Bin Zhang, Xiaoxia Li, Qiaoling Yuan, Tianjiao Liu, Huiying He, Zhaoran Wei, Chao Zhang, Hongsheng Gao, Quan Wang, Qian Qian & Lianguang Shang
Key Laboratory of Crop Heterosis and Utilization, Ministry of Education/ Beijing Key Laboratory of Crop Genetic Improvement, China Agricultural University, Beijing, 100193, China
Jianping Yu & Longbiao Guo
State Key Laboratory of Rice Biology, China National Rice Research Institute, Chinese Academy of Agricultural Sciences, Hangzhou, 310006, China
L.S. and Q.Q. designed the research, and together with J.Y. and C.L. performed most of experiments and analyzed the data. J.Y. analyzed data; H.L., B.Z., X.L., Q.Y., H.H., Z.W., C.Z., S.D., H.G. and Q.W. helped with characterizing the phenotypes; H.L. and T.L. helped with the transcriptomic analysis; L.G. helped to revise the manuscript. J.Y. and L.S. and C.L. conceived the experiments and wrote the manuscript. All of the authors read and approved the final manuscript.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
Yu, J., Liu, C., Lin, H. et al. Loci and natural alleles for cadmium-mediated growth responses revealed by a genome wide association study and transcriptome analysis in rice.
BMC Plant Biol21, 374 (2021). https://doi.org/10.1186/s12870-021-03145-9