The genome-wide impact of cadmium on microRNA and mRNA expression in contrasting Cd responsive wheat genotypes

Background Heavy metal ATPases (HMAs) are responsible for Cd translocation and play a primary role in Cd detoxification in various plant species. However, the characteristics of HMAs and the regulatory mechanisms between HMAs and microRNAs in wheat (Triticum aestivum L) remain unknown. Results By comparative microRNA and transcriptome analysis, a total three known and 19 novel differentially expressed microRNAs (DEMs) and 1561 differentially expressed genes (DEGs) were found in L17 after Cd treatment. In H17, by contrast, 12 known and 57 novel DEMs, and only 297 Cd-induced DEGs were found. Functional enrichments of DEMs and DEGs indicate how genotype-specific biological processes responded to Cd stress. Processes found to be involved in microRNAs-associated Cd response include: ubiquitin mediated proteolysis, tyrosine metabolism, and carbon fixation pathways and thiamine metabolism. For the mRNA response, categories including terpenoid backbone biosynthesis and phenylalanine metabolism, and photosynthesis - antenna proteins and ABC transporters were enriched. Moreover, we identified 32 TaHMA genes in wheat. Phylogenetic trees, chromosomal locations, conserved motifs and expression levels in different tissues and roots under Cd stress are presented. Finally, we infer a microRNA-TaHMAs expression network, indicating that miRNAs can regulate TaHMAs. Conclusion Our findings suggest that microRNAs play important role in wheat under Cd stress through regulation of targets such as TaHMA2;1. Identification of these targets will be useful for screening and breeding low-Cd accumulation wheat lines. Electronic supplementary material The online version of this article (10.1186/s12864-019-5939-z) contains supplementary material, which is available to authorized users.


Background
Cadmium (Cd) contamination of soils is rapidly increasing from both human activities and other environmental causes [1]. Deleterious effects of Cd are dramatic and varied in both humans and plants. In humans, Cd ingestion leads to osteoporosis and emphysema, among other diseases [2,3]. In plants, high Cd results in a wide range of biochemical and physiological disorders, inhibiting yield and quality [4]. Unfortunately, even moderate levels of Cd in soils can result in dramatic Cd accumulation in leaves or grains [5,6]. Contamination in the food chain with Cd thus poses a major widespread hazard to human health [7][8][9]. Wheat, one of the most important worldwide crops, is grown on nearly 20% of cultivated land and serves as a staple food source for 30% of humanity [10,11]. Previous studies indicated that compared with other cereals, wheat accumulates Cd primarily via roots, eventually accumulating Cd in edible grains [12,13]. Therefore, it is very important to understand wheat responses to Cd in order to improve wheat growth, grain quality, and human welfare.
MicroRNAs (miRNAs) are typically 21 nucleotide endogenous non-coding RNAs that play key roles in regulating gene expression post-transcriptionally, usually via cleavage or translational repression of target mRNAs [14,15]. In plants, miRNAs are involved in various physiological processes and play important roles in growth and development [16,17], and numerous miR-NAs participate in responses to biotic and abiotic stresses [18][19][20][21]. Deep-sequencing of developing wheat grains revealed 605 miRNAs [22] and 57 conserved miRNA families have been identified in wheat hybrid necrosis [23]. Hundreds of miRNAs were differentially expressed in leaf and roots tissues in wheat under drought stress, respectively [24]. Individual analysis of miRNAs or mRNAs have been widely reported under various stressors. However, the interplay of Cd, regulatory miR-NAs, and the transcriptome is not at all understood, so there is still an obvious gap in our understanding of their interplay in a major crop. P 1B -ATPases, also known as heavy metal ATPases (HMAs), are proteins that use ATP to pump a wide range of cations across membranes against electrochemical gradients [25]. AtHMA2 and AtHMA4 act as pumps to efflux Zn/Cd out of cells, playing critical role in rootshoot Zn/Cd translocation through xylem loading [26]. The HMA3 proteins, such as AtHMA3, TcHMA3 and OsHMA3, are typically located in tonoplasts and are not regulated through Zn or Cd; these proteins play a role in pumping Zn and/or Cd into vacuoles [27][28][29]. Much available evidence indicates that HMAs play an important role in heavy metal transmembrane transport. However, there is little understanding of HMAs in the wheat genome and their roles in wheat under Cd stress.
Therefore, to gain an understanding of the global transcriptional response and its interplay with miRNAs and HMAs in wheat during Cd stress, we used next generation transcriptome and miRNA sequencing to profile the responses of low-and high-cadmium-accumulating wheat cultivars in response to Cd. Additionally, we paid special attention to regulatory dynamics of the HMA gene family. We further verified by qRT-PCR significantly differentially expressed of miRNAs, mRNAs and TaHMAs. Finally, we performed for gene ontology and pathway analysis, as well as depicted general functional landscapes of miRNAs-TaHMAs expression. We expect that these findings will provide a meaningful resource for the understanding the global transcriptional response to Cd and how these are modulated by miRNAs.

Results and discussion
Differences in cd accumulation in grains among wheat germplasms We first screened 185 wheat germplasms for grain Cd levels using ICP-MS (Additional file 1: Table S1). By group, the average concentrations of Cd in grain decreased in the order of CMC > Chuanyu > CIMMYT, with significant differences between group averages (Fig. 1a). The tolerance limit of Cd in food is 0.1 mg/kg set by Ministry of Health of PRC (People's Republic of China). Thus, there were 15 (28.85%) cultivars that exceeded the tolerance limit in Chuanyu group, 1 (1.79%) in the CIMMYT group, and 56 (76.71%) in the CMC group (Fig. 1b). Our results revealed that the average grain Cd concentration of four durum wheat genotypes was 0.224 mg/kg, and ranged from 0.191 to 0.371 mg/kg. These four durum cultivars contained higher grain Cd concentrations compared with bread wheat (Additional file 2: Figure S1), in agreement with previous studies [30][31][32]. Previous work showed that durum wheat cultivars grown in a variety test in South Dakota contained similar average grain Cd concentrations (0.22 mg/kg, ranged from 0.13 to 0.25 mg/kg) [33], indicating Fig. 1 a Grain Cd concentration among three major wheat groups. The Chuanyu group included 52 wheat cultivars, the CIMMYT group included 56 wheat cultivars, and the CMC group included 73 wheat cultivars. b Percentage of low Cd concentration (≤ 0.0579 mg/kg), moderate Cd concentrations (< 0.1 mg/kg) and exceeding Cd concentration standard (≥ 0.1 mg/kg) wheat in each of the three groups that our results are good representatives of world-wide sampling.
Differentially expressed miRNAs and mRNAs between L17Cd and L17CK, H17Cd and H17CK To illuminate the underlying molecular mechanism, we chose cultivars strongly contrasting for Cd content: a low-Cd-accumulating cultivar (Chuanyu 17, hereafter L17) and a high-Cd-accumulating cultivar (a Chuanyu 17 sib-line, hereafter H17) to determine differentially expressed miRNAs and mRNAs in roots. We generated 12 small RNA sequencing libraries from three biological replicates each for L17 and H17 roots under control and CdCl 2 treatment and sequenced them on an Illumina HiSeq TM 2500. The major characteristics of these libraries are summarized in Table 1. Between 18.8 and 32.4 million total raw reads per library were obtained. After quality filtration, between 12.6 and 20.4 million clean reads were obtained per replicate.
After comparing the distributions of miRNA expression values in the samples after normalization and correlation analysis of miRNA expression levels between samples and samples (Additional file 2: Figure S2), we next performed miRNA expression analysis and found many clear contrasts in miRNA expression between genotypes and groups. Two miRNAs were upregulated in L17Cd (treated with 100 μM CdCl 2 for 24 h Cd) and one was downregulated compared with L17CK (CdCl 2 -free). In the high-Cd-accumulating genotype (H17), 12 miR-NAs in H17Cd (treated with 100 μM CdCl 2 for 24 h Cd) were downregulated compared with H17CK (CdCl 2 -free) ( Table 2). Regarding novel miRNAs in L17, 10 miRNAs were upregulated and nine were downregulated in L17Cd compared with L17CK, while in H17, 49 miRNAs were upregulated and eight were downregulated in H17Cd compared with H17CK (Table 3). In addition, one uniquely expressed miRNA between L17Cd and L17CK, one uniquely expressed miRNA between H17Cd and H17CK, and two expressed miRNAs were commonly expressed. These results indicated that tae-miR9664-3p and tae-miR159a were upregulated in L17Cd group compared with that in L17CK group, while they were downregulated in H17Cd group compared with that in H17CK group. Moreover, according to our RNA sequencing data, the putative target mRNAs (such as Traes_7BL_66C9AAA60, Traes_1BS_1E568BC28, Traes_ 5DL_BDE25A669, Traes_2DS_47DDC83DC and Traes_ 5BS_C3B48D8CC) for upregulated tae-miR9664-3p and the putative target mRNAs (such as Traes_3AL_ F09C7CA19, Traes_7AS_2831C6473, Traes_7DS_0FFF2 FB02 and Traes_4AL_03723C140) for upregulated tae-miR159a were downregulated in L17Cd group compared with that in L17CK group. Previous work demonstrated that tae-miR-9664 was expressed in the flag leaf and developing seed of wheat (Triticum aestivum L.) [34]. Tae-miR159a played positive roles in wheat response to Puccinia striiformis f.sp.tritici through the regulation of TaMyb3 expression [35]. Our results suggested that contrasting Cd-accumulating ability of L17 and H17 might associated with the up-regulation or down-regulation between tae-miR9663-3p, tae-miR159a and their target genes. We also found that tae-miR9774 was uniquely down-regulated between L17Cd and L17CK, while tae-miR398 was specifically down-regulated between H17Cd and H17CK. A previous study demonstrated that tae-miR-9774 was expressed in wheat leaves, stems, roots and young spikes [36]. miR398 was induced by UVB light and heat stress in Arabidopsis but was inhibited by ABA, salinity, cold, and oxidative stress [37,38]. miR398 can regulate copper superoxide dismutase (CSD) 1 and 2 that were both induced by salinity treatment [37,39]. These results suggested that miRNAs could participate in the growth and development of wheat and also function in wheat response to both biotic and abiotic stress. To our knowledge, this is the first study to identify differentially expressed miRNAs in response to Cd accumulation in multiple wheat genotypes contrasting for Cd accumulation.
After checking the distributions of expression values after normalization and correlation analysis of expression levels between samples (Additional file 2: Figure  S3), we next performed differential expression analysis of mRNA transcripts. In the low-Cd-accumulating genotype (L17), a total of 1561 mRNAs were detected to be In each wheat samples, data is shown as the mean of the three libraries generated ± standard deviation differentially expressed with fold change > 2.0, P < 0.05 (volcano plot filtering) and FDR < 0.05 (Additional file 2: Figure S4, Fig. 2a). Among them, 1208 and 353 mRNAs were upregulated and downregulated (fold change > 2.0, P < 0.05 and FDR < 0.05) in L17 roots compared with controls, respectively. In the high-Cd-accumulating genotype (H17), a total of 297 mRNAs were found to be differentially expressed with a fold change > 2.0, P < 0.05 (volcano plot filtering) and FDR < 0.05 (Additional file 2: Figure S4, Fig. 2b). Among them, 204 and 93 mRNAs were upregulated and downregulated (fold change > 2.0, P < 0.05 and FDR < 0.05) in H17 roots compared with controls, respectively.

GO and KEGG enrichment analysis
It is well known that miRNAs play essential roles in plant growth and development and also control plant responses to biotic and abiotic stresses by regulating targeted gene expression [23,40]. Therefore, Gene Ontology (GO) analysis of targeted mRNAs of miRNAs can real the role of differentially expressed miRNAs. Enriched GO categories between L17Cd and L17CK groups, H17Cd and H17CK groups, according to the gene number enriched in GOs, the top 10 biological process (BP), molecular function (MF) and cellular component (CC) for downregulated miRNAs and upregulated miRNAs were shown in (Additional file 2: Figure  S5). Overlaying gene sets (Additional file 2: Figure S6) revealed that down-regulated miRNAs in L17 consisted of genes involved in ubiquitin-dependent protein catabolic processes and cell wall organization, whereas upregulated miRNAs in L17 consist of those involved in protein glycosylation and oxidoreductase activity. Downregulated miRNAs in H17 consist of those involved in vesicle-mediated transport, and response to heat, whereas up-regulated miRNAs in H17 consist of those involved in photosynthetic electron transport chain and iron ion binding. For mRNAs, GOs between L17Cd and L17CK groups, H17Cd and H17CK groups, according to the gene number enriched in GOs, the top 10 BP, MF and CC for upregulated mRNAs and downregulated mRNAs were shown in (Additional file 2: Figure S7). Based on the  results of Venn diagrams, we found that up-regulated mRNAs in L17 prefer to cell wall macromolecule catabolic process and response to stress, whereas downregulated mRNAs in L17 prefer to cell wall organization and DNA replication; up-regulated mRNAs in H17 prefer to photosynthesis and nicotianamine biosynthetic process, whereas down-regulated mRNAs in H17 prefer to defense response and metal ion transport (Additional file 2: Figure S6). Analysis of pathway enrichment for miRNAs between L17Cd and L17CK groups revealed a diverse set of processes responsive to Cd (Fig. 4a, b). Among these pathways, starch and sucrose metabolism was the most prominent pathway targeted by downregulated miNAs (Fig. 4a), whereas glycolysis / gluconeogenesis was the top enriched pathway for up-regulated miRNAs (Fig. 4b).
In the H17Cd and H17CK contrasts we found that carbon metabolism was the most prominent pathway targeted by downregulated miNAs (Fig. 4c), whereas photosynthesis was the most prominent pathway targeted by upregulated miNAs (Fig. 4d). For mRNAs differentially regulated between L17Cd and L17CK groups, the top 10 pathways associated with up-or down-regulated mRNAs are given in Fig. 5. Among  these pathways, phenylpropanoid biosynthesis was the top pathway in upregulated mRNAs (Fig. 5a), whereas starch and sucrose metabolism was the top one in downregulated mRNAs (Fig. 5b). Pathways enriched in the H17Cd and H17CK groups are given in Fig.  5c, where glutathione metabolism was most prominent. On the other hand, phenylpropanoid biosynthesis was the most prominent in down-regulated mRNAs (Fig. 5d).
Our results are complement previous work in other systems. Feng et al. [41] found enrichment of genes differentially-expressed in response to Cd in two sweet sorghum genotypes. Enriched in the high Cdaccumulating sweet sorghum genotype were genes that Sb04g033320, Sb04g033330, and Sb01g012440. Enriched in the low Cd-accumulating sweet sorghum genotypes were genes that Sb02g006960, Sb01g043400, and Sb10g022390. Other studies have indicated that transporters for essential elements such as Fe 2+ , Zn 2+ and Ca 2+ are involved in Cd uptake and transport [42]. Numerous heavy metal transporters, such as ATPbinding cassette transporters (ABC), natural resistanceassociated macrophage proteins (Nramp) and heavy metal ATPases (HMA) participate in the acquisition, distribution, and homeostasis of Cd in plants [27,43,44].
The MAPK signaling pathway is also part of the miRNA-target gene network. MAPK signaling cascades are one of the most conserved pathways in plants and are known to modulate plant response to abiotic stress [45]. Taken together, these studies and ours indicate that the interaction effects of multiple miRNAs are involved in the regulation of diverse physiological systems to relieve the toxicity induced by Cd.

Genes related to transport
Between L17Cd and L17CK groups 223 up-regulated genes and 183 down-regulated genes are classified as transporters based on functional annotation. Between H17Cd and H17CK, 162 up-regulated genes and 134 down-regulated genes were transporter genes according to the functional annotation of the differentially expressed genes (Additional file 3). Of these, Traes_ 2AL_9B175F3DA, Traes_2AS_95611CAD2, Traes_2DL_ D5F2272A9, Traes_2BS_3212EB7DF, Traes_7BL_39745 BF3F, Traes_7DL_44DA42073, Traes_2DL_2A1D45D02, Traes_3DL_4EE343C68, Traes_5BL_4BCE2DBEF, Traes_ 2BS_3212EB7DF, Traes_3AS_CBB5F7EB6, Traes_4DL_ 560FD0832, Traes_7AL_3E508AA9C, Trae_7DL_A5269 C73F, Trae_7AL_8304348B7 and Trae_7BL_A4BE6BD10 are related to metal ion. Through the availability of the wheat genome sequence, it is now possible to identify all the HMA family members in wheat. Here we identify a total of 32 HMA members in the wheat genome (Table 4) based on the genome sequence. We further performed a BLASTN search against the wheat expressed sequence tag (EST) using these HMAs queries to verify the existence of these wheat HMAs. Results showed that most of the TaHMAs' we found in the genome sequence were supported by EST hits. Among the supported TaHMA genes, TaHMA1;4 has the largest hits of ESTs, with 61, followed by TaHMA6;3 with 53. The lengths of these predicted TaHMA proteins varied from 419 aa to 1702 aa with the putative molecular weights (Mw) ranging from 44.07 to 185.21 kDa and theoretical isoelectric points (PI) ranging from 5.23 to 8.21. Eight HMA proteins are known in A.thaliana and nine in O.sativa, which range from 542 aa to 1171aa and from 792 aa to 1060 aa, respectively [46]. Using full-length HMA protein sequences from Triticum aestivum, A.thaliana and O.sativa, we constructed a phylogentic tree to determine the phylogenetic relationship between HMA proteins from Triticum aestivum and those from other species. The HMA gene families could be divided into three groups: the Cu/Ca/Zn/Cd/Co-ATPases group, the Zn/ Cd/Pb/Co-ATPase group, and the Cu-ATPase group [47]. TaHMA1;1-TaHMA1;5 belong to the Cu/Ca/Zn/ Cd/Co-ATPases group, TaHMA2;1-TaHMA2;6, TaHM A3;1-TaHMA3;2 belong to the Zn/Cd/Pb/Co-ATPase group, TaHMA4;1-TaHMA4;3, TaHMA5;1-TaHMA5;3, TaHMA6;1-TaHMA6;3, TaHMA7;1-TaHMA7;3, TaHMA 8;1-TaHMA8;3, TaHMA9;1-TaHMA9;3 belong to the Cu-ATPase group (Fig. 6). To our knowledge, it is first genome-wide analysis of the HMA family in wheat.

Tissue-specific expression patterns of TaHMA genes
Using available RNA-seq data (the wheat expression database, http://wheat.pw.usda.gov/WheatExp/) for five different tissues, the tissue specificity of the TaHMA genes was investigated to focus on the temporal and spatial expression patterns and putative functions of HMA genes in the growth and development. According to FPKM values, we found that the expression levels of the TaHMAs varied significantly in different tissues (Fig. 7). Most HMAs were found to be expressed in five   The red or blue colors indicate higher or lower relative abundance, respectively. Z represent Zadoks scale, a decimal code for the growth stages of cereals. P-value< 0.05 were regarded as statistically significant 1, TaHMA9;3 were up-regulated (Fig. 8a). In H17, the expression levels of TaHMA1 (Fig. 8b).

Response of selected TaHMA genes to exogenous cd stress
HMAs have been implicated in the transmembrane transport of metals, such as Zn, Cd, Pb and Cu in different plant species [46,[49][50][51]. Therefore, we were interested to determine the response of select HMA genes under Cd stress. We thus selected eight TaHMA genes (TaHMA1;4-TaHMA3;1) belonging to the Cu/Ca/Zn/Cd/Co-ATPase group and the Zn/Cd/ Pb/Co-ATPase group to determine their relative expression levels in roots under Cd stress using qRT-PCR. Results showed that relative expression levels of TaHMA1;4, TaHMA1;5, TaHMA2;1 and TaHMA3;1 were significantly down-regulated under Cd treatment compared with those under control conditions in high-accumulating wheat genotype (H17), while relative expression levels of TaHMA1;4, TaHMA1;5, TaHMA2;1 andTaHMA2;2 were significantly upregulated under Cd treatment compared with those under control conditions in low-accumulating wheat genotype (L17) (Fig. 10). Relative expression of TaHMA1;5 in H17, TaHMA2;1 in L17 was consistent with our transcriptome sequencing data. Interestingly, expression data from other TaHMA genes in H17 and L17 treated with Cd were not consistent with sequencing data. A possible explanation for this inconsistence is the number of samples used on this experiment. Previous work indicated that AtHMA1 is responsible for Cu to the stroma, exporting Zn2+ form the chloroplast, or as a Ca 2+ /heavy metal transporter to the intracellular organelle [52][53][54]. AtHMA2 and AtHMA4 in A.thaliana are involved in xylem loading of Zn as well as Cd [55][56][57][58][59]. Compared with wild-type plants, plants overexpressing AtHMA3 exhibited a 2-3-fold increase in Cd accumulation [28]. AtHMA5 is responsible for the Cu translocation from roots to shoots or Cu detoxification of roots [60,61]. AtHMA6 (also known as PAA1) is involved in transporting Cu over the chloroplast envelope, whereas AtHMA8 (PAA2) most likely transports Cu into the thylakoid lumen to supply plastocyanin [62,63]. AtHMA7 has been proposed to delivery Cu to ethylene receptors and Cu homeostasis in the seedlings [64,65]. Among rice HMA characterized, OsHMA1 transports Zn and Cd [66], OsHMA2 is responsible for Cd accumulation [67], OsHMA5 is involved in loading Cu to the xylem of the roots and other organs [68], and OsHMA9 is responsible for transporting Zn, Cu, Pb and Cd [69]. In wheat, TaHMA2 is a plasma membranelocalized Zn/Cd transporter that pumps Zn/Cd into apoplast [26]. Our results indicated that Cd stress induced the up-regulation of TaHMA2;1 in L17.

Conclusions
We present here a broad range of genome-wide data surrounding Cd response in contrasting wheat lines. A total three known and 19 novel differentially expressed microRNAs (DEMs), and 1561 differentially expressed genes (DEGs) were found in L17 after Cd treatment. In H17 we found 12 known and 57 novel DEMs, 297 Cdinduced DEGs. Moreover, we present the identification of 32 TaHMA genes in wheat. Finally, we provide microRNA-TaHMA expression networks that suggest that miRNAs can regulate TaHMAs which need to be validated experimentally in future. Taken together, our results indicate that microRNAs may play important roles in wheat under Cd stress through regulating target genes such as TaHMA2;1. This information will be useful for screening and breeding low-Cd accumulation wheat genotypes.

Study site and wheat germplasms
The Shifang heavy industry and Cd-polluted area was selected as our study site in the northeast of Sichuan province. Low grain Cd accumulative wheat cultivars were selected from the Cd-polluted area as suitable for planting in Cd pollution farmland. A total of 185 wheat germplasms (provided by Chengdu Institute of Biology, Chinese Academy of Sciences) were chosen to include wide genetic diversity. Detailed information of wheat cultivars were shown in Additional file 1:

Plant growth conditions and cd treatment
A low-Cd-accumulating cultivar (Chuanyu 17, namely L17) and its sib-line H17 (a high-Cd-accumulating cultivar) from Chuanyu group were used in the study. Seeds were sterilized by soaking in 2% H 2 O 2 for 10 min and fully rinsed with deionized water. After sterilizing, the seeds were soaked in deionized water at room temperature for 6 h, and then planted hydroponically in Hoagland solution in growth chambers at 24 ± 2°C and 50% relative humidity with a photoperiod of 16 h light/8 h dark. One-week-old seedlings were treated with 0 (CK) and 100 μM CdCl 2 for 24 h (Cd). Roots from the plants of similar size were harvested separately and washed three times with deionized water. Roots from three biological replicates were frozen in liquid nitrogen immediately and stored at − 80°C for small RNA and transcriptome sequencing.

Determination of grain and root cd concentrations
For wheat grain, whole grain was ground into flour (< 1 mm). Flour was dried at 80°C for 12 h and then digested by HNO 3 and H 2 O 2 in a microwave digester. Roots from three plants as replicates for each of CK and Cd treatments were dried at 70°C to a constant weight and then also digested by HNO 3 and H 2 O 2 in a microwave digester. All grain samples were processed together with quality controls. Cd concentration was measured using inductively couple plasma-mass spectrometry (ICP-MS, NexIONTM 300, USA) following the manufacturer's instructions [70]. A Certified Reference material (CRM; GBW10020, provided by the National Research Center for CRM, China) was applied to assess the precision of the analytical procedures for plant material.

RNA extraction
One week-old L17 and H17 seedlings were treated with 0 (CK) and 100 μM CdCl 2 (Cd) for 24 h. Then, roots of 12 seedlings (each treatment for three biological replicates) were collected for RNA extraction. Total RNA were extracted using the TRIzol Reagent (Invitrogen, USA) following the manufacturer's protocol. RNA integrity was evaluated using the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA), and each sample had an RNA integrity number (RIN) > 8.0.
Small RNA sequencing and analysis cDNA libraries were built using SuperScript II Reverse Transcriptase basing on the manufacturer's instructions.
Resultant libraries were then sequenced on the Illumina HiSeqTM 2500 sequencing platform (OEbiotech Company, Shanghai, China). Raw reads were converted into sequence data by base calling. After reads with lowquality, 5′ primer contaminants and ploy (A), without 3'adapter and insert tag, shorter than 15 nt and longer than 41 nt from the raw reads were filtered, the clean reads were obtained and mapped against the wheat genome. Clean reads were annotated with miRBase (version 21.0, http://www.mirbase.org/) to identify know miRNAs using Bowtie software. To identify novel miRNAs, the remaining unannotated clean reads that could be aligned to the genome were analyzed by miRDeep2 and RNAfold. The criteria for designation as a novel miRNA were as follows: length of 18-24 nt, precursors with a signature hairpin structure and the formation of maturation achieved by Dicer, minimum free energy of precursors of less than 18 kcal/mol, a minimum support number for the maturity sequence of precursors of at least 5, and potential miRNA sequence with less than 3 nt mispairing in the sequence of the mature and perfectly matched middle sequence. Differentially expressed miRNAs (DEMs) were identified with the threshold of p value < 0.05. While the p value was calculated with the DEG algorithm in the R package. The targets of DEMs were predicted by using targetfinder [48]. per kilobase of transcript sequence per million base pairs sequenced (FPKM) value of each gene was calculated using cufflinks, and the read counts of each gene were obtained by htseq-count. Differentially expressed genes (DEGs) were identified using the DESeq. R package functions estimateSizeFactors and nbinomTest. P value < 0.05, fold Change > 2 and false discovery rate (FDR) < 0.05 were set as the threshold for significantly differential expression. Hierarchical cluster analysis of DEGs was performed to explore genes expression pattern.

Gene ontology (GO) and KEGG pathway analysis
For miRNAs, GO enrichment and KEGG pathway enrichment analysis of different expressed miRNA-target-Gene were performed using R based on the hypergeometric distribution. For mRNAs, GO enrichment and KEGG pathway enrichment analysis of DEGs were performed using R based on the hypergeometric distribution [72].

Identification of HMA gene family in wheat
The HMA gene family was identified following the method as described by Li  Gene structure construction, protein domain and motif analysis Gene structure information was obtained from the Ensemble plants database (http://plants.ensembl.org/index. html). All full-length amino acid sequences of the TaH-MAs were also used to identify conserved domain motifs by the Multiple Em for Motif Elicitation (MEME) tool [75]. The parameters were set as follows: maximum numbers of different motifs, 10; minimum motif width, 6; maximum motif width, 50.

Phylogenetic analysis
The HMA protein sequences from A.thaliana, O.sativa and T.aestivum L. were performed for multiple alignments by CLUSTALW and the results of alignment were used to construct phylogenetic tree using the NJ method in MEGA (version 6.0) [76]. Bootstrap test method was emplyed and the replicate was set to 1000.
The TaHMA gene expression analysis by RNA-seq data To study the expression of TaHMA genes in different organs, the wheat expression database (http://wheat.pw. usda.gov/WheatExp/) was used to investigate the differential expression of TaHMAs. The FPKM (fragments per kilobase of transcript per million fragments mapped) value was calculated for each HMA gene, the log2 transformed FPKM value of the TaHMA genes were used for heat map generation. And p-value < 0.05 were taken as statistically significant threshold [77].

Quantitative real time polymerase chain reaction (qRT-PCR) analysis
Eight miRNAs and six genes with different expression patterns and eight TaHMA genes were selected for validation by quantitative real-time RT-PCR (qRT-PCR). Total RNA was extracted from roots of L17 and H17 seedlings treated with 0 (CK) or 100 μM CdCl2 (Cd) for 24 h. For mRNA, first strand cDNA was synthesized using HiScript II Q RT SuperMix (Vazyme, R223-1). For miRNA, the first strand cDNA was synthesized using miScript II Reverse Transcription Kit (Qiagen, 218,161). The qRT-PCR was carried out using Quanti-Fast® SYBR® Green PCR kit (Qiagen, 204,054) according to the manufacturer's instructions. The primers used in the qRT-PCR analyses were listed in (Additional file 4). GAPDH and 5S were used as internal genes for mRNAs and miRNAs, respectively. Three technical replicates were performed for each gene. The expression levels were calculated from the 2 -ΔΔCt value.

Statistical analysis
All data were represented as mean ± SD and analyzed with SPSS20.0 software. For multiple group comparison, ANOVA with least significant difference was applied. P < 0.05 was considered statistically significant.

Additional files
Additional file 1: Grain Cd concentration among wheat germplasms.
(XLSX 16 kb) Additional file 2: Figure S1. Grain Cd concentration between bread wheat and durum wheat. Figure S2. A Box plot of miRNA. B Correlation analysis of expression levels of miRNAs between samples and samples. Figure S3. A Box plot of mRNA. B Correlation analysis of expression levels of mRNAs between samples and samples. Figure S4. A The volcano plot illustrates the distribution of mRNA expression fold changes vs P values between L17Cd and L17CK. B The volcano plot illustrates the distribution of the data in mRNA profiles between H17Cd and HL17CK. Red, green and black points in the volcano plot represents significantly upregulated mRNAs, significantly downregulated mRNAs and not differential expressed mRNAs, respectively. Figure S5.