Natural variations in Cd resistance and morpho-physiological characterization of wheat genotypes
To identify genotypic differences regarding Cd resistance among natural wheat genotypes, a panel consisting of 442 cultivars were grown under 10 μM CdCl2 for 14 d using hydroponic culture system, and the chlorotic spots of leaves was used as a proxy for Cd resistance. It was observed that maximal lengths of leaves and roots ranged from 5.0 cm to 38.8 cm (with the variance coefficient of 27.2%) and from 6.0 cm to 33.6 cm (with the variance coefficient of 27.3%), respectively (Fig. 1A and 1B). Based on wide natural variation in Cd resistance occurring among different wheat genotypes, the lengths of chlorotic leaves caused by Cd toxicity were further investigated to identify the most sensitive wheat genotypes. The results showed that the chlorotic leaf lengths ranged from 1.83 cm to 15.8 cm, with the variance coefficient of 13.3% (Fig. 1C). Based on these investigations,
Next, to examine the biological basis of these natural variations, we selected ‘Zhengmai 1860’ (‘ZM1860’) and ‘Zhoumai 32’ (‘ZM32’) as the most Cd-resistant and Cd-sensitive wheat genotypes for the subsequent analysis. Under Cd-free (control) conditions, the two wheat genotypes showed no marked differences in growth performance (Fig. 1D). However, ‘ZM32’ showed obvious plant lodging compared with ‘ZM1860’ under Cd toxicity (Fig. 1D). After exposure to 10 μM Cd for 14 d, the leaves of ‘ZM32’ showed severe chlorosis than those of ‘ZM1860’, which showed few significant differences with the plants under control (Fig. 1E). Moreover, the roots of ‘ZM32’ also showed significant inhibition (Fig. 1F). To be accurate, the dry biomasses of shoots and roots of ‘ZM1860’ were 73.53% and 124.4% higher than those of ‘ZM32’ under Cd toxicity (Fig. 1G, H). Further, we also found that the ratio of root/shoot was significantly higher in ‘ZM1860’ than in ‘ZM32’ (Fig. 1I). In terms of the leaf morphology, both maximal leaf lengths and total leaf areas were tested to be 21.67% and 65.90% higher in ‘ZM1860’ than in ‘ZM32’ (Fig. 1J, K). However, the leaf angles of ‘ZM32’ were significantly larger than those of ‘ZM1860’ (Fig. 1L), the leaves of whichshowed significantly higher chlorophyll concentrations (Fig. 1M).In terms of the root system architecture, all of the maximal root lengths (Fig. 1N), root surface areas (Fig. 1O), root volumes (Fig. 1Q), root tip numbers and (Fig. 1R) were markedly larger in ‘ZM1860’ than in ‘ZM32’ (Fig. 1N-R). However, average root diameters were not identified to be significantly different between the two wheat genotypes (Fig. 1N).
To examine the cellular effects underlying the morphological differences, we characterized the leaf ultrastructure using electron microscopy. The stoma density were highly similarly between ‘ZM1860’ and ‘ZM32’ under both control and Cd conditions (Fig. 2A-E), whereas the trichome density of ‘ZM1860’ was significantly higher than that of ‘ZM32’ under Cd toxicity (Fig. 2J). Further, TEM analysis showed that the cell morphology of ‘ZM1860’ under Cd toxicity remained similar to that under Cd toxicity (Fig. 2F, H); in contrast, more leaf detached cells in leaves were examined in ‘ZM32’ than in ‘ZM1860’ (Fig. 2G, I). Moreover, more complete chloroplasts with better grana were identified in the leaf mesophyll cells of ‘ZM1860’ compared with ‘ZM32’ (Fig. 2K-O).In addition, we also found that ‘ZM1860’ had much thicker cell walls than did ‘ZM32’ under Cd toxicity although no obvious differences were tested between the genotypes under control condition (Fig. 2P-T). Further analysis of cell wall components showed that the shoot lignin concentration of ‘ZM1860’ was 17.91% higher than that of ‘ZM32’ (Figure S1E) although the concentrations of cellulose, semi-cellulose, and pectin were highly similar between the two wheat genotypes under Cd (Figure S1A, C, G). Meanwhile, we also investigated their cell wall components in the roots. A significantly higher pectin concentration was observed in ‘ZM32’ than in ‘ZM1860’ (Figure S1H) whereas the three other components did not differ between ‘ZM1860’ and ‘ZM32’ (Figure S1B, D, F).
Differential ROS and phytohormone homeostasis in wheat genotypes
Significantly higher concentrations of MDA, an indicator of membrane lipid peroxidation, were found in ‘ZM32’ than in ‘ZM1860’ under Cd toxicity (Fig. 3A, B). The pivotal role of sulfydryl in Cd detoxification led us to test the GSH differences between the Cd-resistant and Cd-sensitive genotypes. The result revealed a significantly higher of GSH in the shoots of ‘ZM32’ than in the shoots of ‘ZM1860’ (Fig. 3C) whereas obvious differences were not observed in the roots (Fig. 3D).
To identify the ROS catabolism differences between the wheat genotypes under Cd toxicity, we tested the ROS levels and antioxidant enzyme activity. The results showed that both H2O2 and O2- had obviously higher levels in ‘ZM32’ than in ‘ZM1860’ (Fig. 3E-H), particularly in the roots of these two wheat genotypes (Fig. 3F, H). Further DAB and NBT staining confirmed the ROS over-accumulation in ‘ZM32’ relative to ‘ZM1860’ (Fig. 3I, J). Subsequently, the activity assays of antioxidant enzymes showed that different antioxidant enzymes showed distinct roles in the ROS homeostasis between shoots and roots. To be specific, the CAT, APX, and GPX enzymes presented significantly higher activity in the shoots of ‘ZM32’ than in the shoots of ‘ZM1860’ (Fig. 3K-T). However, we found that ‘ZM1860’ had much higher activity ofbothSOD and CAT in the roots (Fig. 3K-T).
Although several phytohormones were found to be different between the Cd-resistant and Cd-sensitive genotypes under control (Fig. 3U), more attention was paid to their differences in these two genotypes under Cd toxicity. In the shoots, an obviously higher concentration of SA was observed in ‘ZM32’ than in ‘ZM1860’ whereas no significant differences were not identified in other six phytohormones (Fig. 3U). In the roots, both IAA and GA3 showed markedly higher levels in ‘ZM1860’ than in ‘ZM32’, whereas the concentrations of ABA and JA in ‘ZM32’ were almost twice of those in ‘ZM1860’ (Fig. 3U).
Ionomic basis of differential plant Cd resistance in wheat genotypes
Based on the significant morpho-physiological differences between the Cd-resistant and Cd-sensitive genotypes, ICP-MS was employed to reveal the ionomic basis of differential Cd resistance. Firstly, total Cdconcentrations of whole plants were found to be significantly higher in ‘ZM1860’ than in ‘ZM32’ (Fig. 4A), whereas ‘ZM1860’ accumulated 47.75% of Cd higher than ‘ZM32’ in whole plants (Fig. 4B). Furthermore, the Cd concentration of whole shoots was significantly higher in ‘ZM1860’ than in ‘ZM32’ (Fig. 4C). Further, we tested the Cd concentrations at different positions of leaves, and the Cd concentrations were gradually decreased from mature leaves to young leaves in both ‘ZM1860’ and ‘ZM32’ (Fig. 4D). Moreover, the Cd concentration of each leaf was significantly higher in ‘ZM1860’ than in ‘ZM32’ (Fig. 4D). Subsequently, we found that the root Cd concentration of ‘ZM32’ was 25% higher than that of ‘ZM1860’ (Fig. 4E). Based on the plant biomasses and Cd concentrations above, both shoot and root content were obtained. The results showed that more Cd was accumulated in the shoots of ‘ZM1860’ than in the shoots of ‘ZM32’ (Fig. 4F), which had a higher content of root Cd than ‘ZM1860’ did (Fig. 4G). This differential Cd allocation led to a higher translocation coefficient of Cd from roots to shoots in ‘ZM32’ than in ‘ZM1860’ (Fig. 4H). To further characterize the differential long-distance translocation of Cd, the vascular Cd concentrations were determined. The results showed that the Cd concentrations in the stems and xylems of ‘ZM32’ were more than twice and 1.5-fold of those in ‘ZM1860’ (Fig. 4I, J).
To further identify the differential long-distance translocation of Cd between the Cd-resistant and Cd-sensitive genotypes, the grains were harvested at the maturity stage of wheat plants that were grown under 10 μM Cd condition. Although obvious significant differences were not in the grain appearances between ‘ZM1860’ and ‘ZM32’ (Fig. 4K), a significant higher Cd concentration was identified in the grains of ‘ZM32’ than in the grains of ‘ZM1860’ (Fig. 4L). Moreover, we also found that the ‘ZM1860’ grains had much higher concentrations of K, Mg, Fe, Cu, Mn, and Zn than ‘ZM32’ did (Fig. 4M, O, P, Q, R, S) except for a higher level of grain Ca in ‘ZM32’ than in ‘ZM1860’ (Fig. 4N).
WGS-based identification of genomic variations between Cd-resistant and Cd-sensitive genotypes
To identify genomic variations between the Cd-resistant and Cd-sensitive genotypes, we performed deep WGS of both ‘ZM1860’ and ‘ZM32’, generating 170.5 Gb (about 10.21 ×) 192.5 Gb (about 11.51 ×) clean data for ‘ZM32’ and ‘ZM1860’, respectively (Supplementary Table S2). The high-quality sequencing parameters, including Q20 and Q30, were ≥ 92%, and the mapping and coverage rates of WGS data were ≥ 98% and ≥ 95% for these two genotypes, respectively (Supplementary Table S2). Such a broad coverage and read depth validated the high quality of genome assembly and sequencing data.
As shown in Fig. 5, thousands of genomic variants, including SNPs and InDels, between ‘ZM32’ and ‘ZM1860’, were randomly distributed across the 21 chromosomes of allohexaploid wheat. In total, 14,833,473 SNPs were identified across the 21 chromosomes of wheat (Fig. 5A), ranging from 30,532 (Chr. 6D) to 3,088,362 (Chr. 3B) with an average of 674,289 SNPs on each chromosome (Fig. 5A, B). Strong differences were observed in terms of SNP proportions among the A, B, and D subgenomes, with the B subgenome being the more polymorphic (8,109,845 SNPs; 54.67%), followed by the A-genome (5,779,020; 38.96%) then the D-genome (816,689; 5.507%) (Fig. 5B). The frequency of SNPs on each chromosome varied from 1 per 24,585 bp (Chr. 6D) to 1 per 241 bp (Chr. 3B), with an average value of 1 per 3,967 bp (Fig. 5B). The nucleotide diversity π (average SNP number per nucleotide) ranged from 4.00 × 10–5 (Chr. 6D) to 4.15 × 10−3 (Chr. 3B), with an average value of 1.27 × 10−3, 1.74 × 10−3, and 1.90 × 10−4 on the A, B, and D subgenomes, respectively (Fig. 5C). The annotations of reference (‘Chinese Spring’) genome were used to examine the distribution of SNPs within various genomic features. Among the genome-wide 14,833,473 SNPs between the Cd-resistant and Cd-sensitive genotypes, 71.04% of SNPs belonged to the transition type (A/G and C/T; Ts), which was more than that 28.96% of the trans-versions (A/C, A/T, C/G and G/T; Tv) (Supplementary Figure S2A). This led to an average Ts/Tv ratio of 2.45 (Supplementary Figure S2A), a value that is similar to what has been reported previously in wheat (Rimbert et al., 2018). There was 1.424% of the total SNPs detected in the genic regions (including exonic, intronic, and splicing regions), whereas a much larger proportion of SNPs (98.576%) were located in the 2.0-kb upstream (promoter), 1.0-kb downstream, and other intergenic regions (Supplementary Fig. S2B). Among the SNPs located in the genic areas, 40.38% of the variations occurred in the coding sequences (Supplementary Fig. S2B). Only a minor fraction (1.01%) of the SNPs located in the exonic regions were mapped onto the stop codons (Supplementary Fig. S2C). In general, the ratio of synonymous SNPs (46.62%) to non-synonymous SNPs (52.37%) was about 0.89 (Supplementary Fig. S2C).
Also, a total of 1,147,284 InDels were identified that were unevenly distributed over the wheat genome (Fig. 5D), ranging from 7,265 (Chr. 4D) to 183,650 (Chr. 3B) with an average of 54,633 InDels on each chromosome (Fig. 5D, E). The majority (97.6%) of total InDel variants occurred in the 2.0-kb upstream (promoter), 1.0-kb downstream, and other intergenic regions, whereas a little proportion (2.40%) of InDels were identified in the genic sequences (Supplementary Fig. 3A). Among the InDels occurring in the exons, only ~1.64% were identified in the stop codons. For the other InDels, 47% of SNPS caused frame-shift deletion or insertion (Supplementary Fig. 3B). In terms of the InDel length variation occurring in the genomic or genic regions, the mono-nucleotide type was the most (Supplementary Fig. 3C, D).
Based on the annotation results shown in Fig. 6A, there were 117,861 unigenes (94.70%)
annotated in NR, 117,090 unigenes (94.08%) annotated in Nt, 36,445 unigenes (29.28%) annotated in KO, 82,499 unigenes (66.28%) annotated in SwissProt; 84,666 unigenes (68.02%) annotated in PFAM; 84,661unigenes (68.02%) annotated in GO; 23,499 unigenes (18.88%) annotated in KOG. In summary, there were 15,305 unigenes (12.29%) annotated in all databases, and 122,930 unigenes (98.77%) annotated in at least one database.
KOG, GO, and KEGG are international standardized systems of gene functional classification that describe properties of genes and their products in any organism. The genome-wide unigenes with polymorphisms between the Cd-resistant and Cd-sensitive genotypes were categorized into 26 groups of KOG function clusters, among which the ‘general function prediction only’ cluster had the highest number (3,396; 14.45%) of unigenes, whereas the small clusters of ‘cell wall/membrane/envelope biogenesis’, ‘inorganic ion transport and metabolism’, and ‘defense mechanisms’ were paid more attention (Fig. 6B) due to their close relationships with Cd resistance. Regarding the GO analysis, the genome-wide DNA variants were found to be accumulated in the terms involving detoxification, response to stimulus, antioxidant activity, and transporter activity (Fig. 6C). In terms of the KEGG analysis, revealed that the pathways were identified to be implicated in processes of environmental adaptation, signaling transduction, membrane transport, transport and catabolism, and cell growth and death (Fig. 6D).
Transcriptomics-based variations between Cd-resistant and Cd-sensitive genotypes
Based on the morpho-physiologic, ionomic, and genomic differences, to further identify the DEGs regulating differential Cd resistance in wheat genotypes, a high-throughput genome-wide transcriptome sequencing was performed for both ‘ZM32’ and ‘ZM1860’. After removing adapter sequences and low-quality reads, approximately 4.7 × 107 clean reads were obtained from each sample. The total length of clean reads from the 24 samples reached about 1.67 × 1010 nt with Q20 > 97% and Q30 > 93%, and the mapping rate of all samples were more than 90% (Supplementary Table S3). Taken together, these parameters indicated that the transcriptome sequencing data were of high quality suitable for the following analysis.
Principal component analysis showed that compared with the control, the two genotypes and their tissues exhibited significantly different transcriptomic features under Cd toxicity (Fig. 7A), which indicated genotype- and organ-specific transcriptomic responses to Cd toxicity. Further, Pearson correlation coefficients between each RNA-seq group showed that the same organs (shoots/roots) had much higher correlation than the same treatments (control/Cd) irrespective of ‘ZM1860’ or ‘ZM32’ (Fig. 7B). In both shoots and roots, more Cd-responsive DEGs were identified in ‘ZM32’ than in ‘ZM1860’ (Fig. 7C). Moreover, a larger number of DEGs between ‘ZM1860 and ‘ZM32’ were found in both shoots and roots, and the DEGs had higher expression levels in ‘ZM32’ than in ‘ZM1860’ under both control and Cd conditions (Fig. 7D). Venn diagram analysis showed that 192 and 1,776 genes responsive to Cd toxicity were differentially expressed in the shoots and roots between ‘ZM1860’ and ‘ZM32’ (Fig. 7E, F), and some of these 1,776 DEGs were presumed to be the core regulators responsible the differential long-distance translocation of Cd from roots to shoots. GO enrichment analysis revealed that the Cd-responsive DEGs in the shoots were mainly involved in the categories of active transmembrane transporter activity, oxidoreductase activity and ATPase (Fig. 7G). KEGG analysis showed that the metabolism of fatty acid, ascorbate, sulfur (cysteine, methionine, and glutathione), and phytohormone signal transduction were over-represented under Cd toxicity (Fig. 7H).
Identification of core Cd transporter family members regulating differential Cd resistance
The differential ionomic profiling indicated that less root xylem Cd translocation causing a lower shoot Cd concentration, was the key factor for Cd resistance in ‘ZM1860’ relative to ‘ZM32’ (Fig. 4). Although many transporters are reported to be involved in the Cd resistance and detoxification, HMA2/4s have been the only two transporters essential for root-to-shoot Cd translocation to counter Cd toxicity in plants (Li et al. 2014). Based on these findings, we defined the candidate genes, responsible for differential wheat Cd resistance, as the DEGs only in the roots.
Among the nine TaHMA (HMA1-HMA9) subfamilies, only six of them were identified to be expressed in this study (Fig. 8A). Both HMA2 and HMA4 have a single member in both Arabidopsis ad rice, namely AtHMA2/OsHMA2 and AtHMA4/OsHMA4 (Berthomieu et al. 2003). However, all of the genome-wide six TaHMA2s and three TaHMA4s were identified in this study, however, we only examined the differential expression of TaHMA2b-7A (TraesCS7A02G420000) and TaHMA2b-7B (TraesCS7B02G320900) in the roots under Cd toxicity (Fig. 8A). Given the fact that the Cd-sensitive genotype ‘ZM32’ had a higher translocation coefficient of Cd from roots to shoots than ‘ZM1860’ did, therefore, TaHMA2b-7B was excluded from the two candidate genes due to its higher expression in ‘ZM1860’ than in ‘ZM32’. In summary, TaHMA2b-7A was proposed to be the key candidate responsible for differential Cd resistance between these two wheat genotypes. A genome-wide association analysis of a natural population consisting 206 wheat genotypes was performed based on the maximal root length under Cd toxicity, and we also found that the TaHMA2 homologs, including TaHMA2b-7A, were located in the candidate genomic regions for natural variations in Cd resistance (Fig. 8B).
Molecular characterization of TaHMA2b-7A
Subsequently, more in-depth analysis was conducted to characterize the molecular function of TaHMA2b-7A and the possible reasons for its differential expression between the Cd-resistant and Cd-sensitive genotypes. Gene sequence analysis showed that the full-length genomic region, including nine exons and eight introns, of TaHMA2b-7A was 6,535 bp, encoding 939 amino acids (Fig. 8C). WGS and PCR-based gene cloning showed that sequence differences were not identified in the exon regions, which led to identical amino acid sequences between ‘ZM1860’ and ‘ZM32’. Further sequence polymorphism analysis revealed a 25-bp deletion (-697GTGTGGCCCATGAATTCTAGGTACA-721) in the TaHMA2b-7A promoter region of ‘ZM32’ compared with ‘ZM1860’ (Fig. 8D). Cis-acting regulatory element analysis showed that a core Cu-responsive element (-716GTAC-720) existed within the InDel site (Supplementary Table S4), which might be the major reason for the differential expression of TaHMA2b-7A in the Cd-resistant and Cd-sensitive genotypes.
In addition to eight transmembrane domains with a very long C-terminal region (Fig. 8E), the conserved motifs of N-terminal metal-binding (GICCPSE), phosphatase (TGES), ion translocation (CPC), aspartyl phosphate formation (DKTGT), and ATP-binding (GDGMNDA) were also identified in the amino acid sequences of TaHMA2b-7A (Fig. 8F). Within the N-terminal metal-binding domain (GICCPSE), the two cysteine residues determining the HMA activity and affinity with metal substrates (Wong et al. 2009), were identified in TaHMA2b-7A, OsHMA2, and AtHMA2 (Fig. 8F). To further characterize the evolutionary pattern of TaHMA2b-7A, the six HMA2 homologs were physically mapped onto the wheat chromosomes (Fig. 8G). We found that the physical distances within both the TaHMA2a-7A and TaHMA2b-7A pair and the TaHMA2a-7D and TaHMA2b-7D pair were smaller than 100 kb (Fig. 8G), which indicated these four HMA homologs were formed as a result of tandem duplication, whereas the other two homologs (TaHMA2a-7B and TaHMA2b-7B) arose from segmental duplication. Evolutionary selection pressure analysis revealed that the non-synonymous (Ka)/synonymous (Ks) nucleotide substitution rate (0.85) of TaHMA2b-7A was < 1.0 (Fig. 8H), which suggesting a purifying selection to preserve its function. Furthermore, subcellular localization analysis showed that TaHMA2b-7A was localized in the plasma membrane (Fig. 8I).
Phylogeny analysis showed that the HMA2 homologs of diverse plant species were obviously classified into two clades: monocots and dicots (Fig. 9A), and this suggested the evolutionary divergence of HMA2s occurred after plant speciation. Functional studies and phylogeny analysis of HMAs have shown that these transporters can be divided into three subgroups based on their metal-substrate specificity: (i) a Cu/Ca/Zn/Cd/cobalt (Co) group, (ii) a Zn/ Co/Cd/lead (Pb) group, and (iii) a Cu group (Hermand et al. 2014; Fig. 9B), and all of AtHMA2, OsHMA2 and TaHMA2b-7A were clustered within the Zn/ Co/Cd/Pb-ATPase clade (Fig. 9B). Then, we investigated the spatio-temporal expression patterns of TaHMA2b-7A (Fig. 9C-H). After exposure to 12-h Cd toxicity, the TaHMA2b-7A expression was dramatically increased until 1 d, and its expression was kept steadily from 1-d to 5-d Cd toxicity, which showed periodic increase-steadiness changes (Fig. 9C). In general, the TaHMA2b-7A expression showed an increased trend with the solution Cd concentrations increasing (Fig. 9C). In addition to the transcriptional response of TaHMA2b-7A to Cd toxicity, we also examined the TaHMA2b-7A expression under low Fe, low Zn and salt stress. The results showed that TaHMA2b-7A was induced by the three nutrient stresses, particularly in the roots (Fig. 9E-G). Therefore, it led us to investigate the growth performance of the Cd-resistant ‘ZM1860’ and Cd-sensitive ‘ZM32’under low Fe and low Zn stresses. The results showed that ‘ZM1860’ presented higher resistance to low Zn but sensitivity to low Fe (Supplementary Fig. S4).
To further explore the roles of TaHMA2b-7A in the growth and development of wheat plants, we tested its expression at the stages of seedlings, tillering, flag leaves, booting, anthesis, grain filling and ripening. The expression profiling showed that at the stages of both tillering and flag leaf, TaHMA2b-7A showed relatively higher expression in the roots; however, when wheat plants enter the reproductive stages, TaHMA2b-7A wasmainly expressed in the internodes, spikelets, and grains (Fig. 9H).
Natural variations of TaHMA2b-7A in regulating wheat Cd resistance
To further understand the pivotal roles of TaHMA2b-7A in regulating wheat Cd resistance, an association analysis between their differential expression and wheat Cd resistance was performed in 14 wheat cultivars randomly selected from our collected wheat resources, including ‘ZM1860’ and ‘ZM32’. Based on the Cd resistance assessment, ‘ZM1860’ and other six cultivars (‘Lunxuan 715’, ‘Zheng 103’, ‘Shannongtedali 1’, ‘Zhongyuan 6’, ‘Fengdecun 5’, and ‘Xinmai 19’) were identified to be Cd-resistant with less leaf chlorosis under Cd toxicity (Fig. 10A, B), whereas ‘ZM32’ and other xin cultivars (‘Yangmai 20’, ‘Zengmai 518’, ‘Jingdong 1’, ‘Zhengmai 369’, ‘Han 6172’, and ‘Ruiquanmai 24’) were rated with high Cd sensitivity with more severe leaf chlorosis (Fig. 10A, B). In general, the chlorotic lengths of leaves in the seven Cd-sensitive genotypes were significantly larger than those in seven Cd-resistant genotypes (Fig. 10C). Further,we examined the Cd concentration in both shoots and roots of the 14 wheat genotypes, and calculated the translocation coefficients of Cd from roots to shoots. On the whole, the Cd-sensitive genotypes allocated a larger proportion of Cd to the shoots than the Cd-resistant genotypes did (Fig. 10D). TaHMA2b-7A showed significantly higher expression levels in the seven Cd-sensitive genotypes than in the seven Cd-resistant genotypes (Fig. 10E). Moreover, the leaf chlorosis of wheat plants was positively correlated with the TaHMA2b-7A expression under Cd toxicity (Fig. 10F).