Comparative Small RNA Profiling and Functional Exploration on Wheat With High- and Low-Cadmium Accumulation

Cadmium is a toxic metal widely found in workplaces and plant soil because of extensive industrialization. Wheat is an important source of food generated from plant soil. The different responses of wheat against different omic levels of cadmium have been observed and widely studied worldwide. With the development of high-throughput sequencing, micro-level biological research has extended to the microRNA level. In this study, high-cadmium-accumulating wheat cultivars (Annong9267) and low-cadmium-accumulating wheat cultivars (Qian 102032) were used as experimental models. The two cultivars were treated by Cd for 2 h to explore the microRNA profiles in root and leaf tissues through small RNA sequencing. Important small RNAs, such as tae-miR9663-5p and tae-miR6201, and potential small RNA-mediated mechanisms associated with cadmium accumulation were identified by summarizing specific microRNA profiling patterns and their respective target genes. At the wheat roots and leaves, differentially expressed small RNAs related to cadmium accumulation in different plant tissues (roots or leaves) were identified, and functional enrichment analyses on target genes of differentially expressed miRNAs in low- and high-cadmium-accumulating wheat cultivars in different plant tissues (roots or leaves) obtained some known mature miRNAs and new miRNAs. The identified miRNA will be regarded as a potential screening biomarker for low-cadmium-accumulating wheat.

Cadmium is a toxic metal widely found in workplaces and plant soil because of extensive industrialization. Wheat is an important source of food generated from plant soil. The different responses of wheat against different omic levels of cadmium have been observed and widely studied worldwide. With the development of high-throughput sequencing, micro-level biological research has extended to the microRNA level. In this study, high-cadmium-accumulating wheat cultivars (Annong9267) and low-cadmiumaccumulating wheat cultivars (Qian 102032) were used as experimental models. The two cultivars were treated by Cd for 2 h to explore the microRNA profiles in root and leaf tissues through small RNA sequencing. Important small RNAs, such as tae-miR9663-5p and tae-miR6201, and potential small RNA-mediated mechanisms associated with cadmium accumulation were identified by summarizing specific microRNA profiling patterns and their respective target genes. At the wheat roots and leaves, differentially expressed small RNAs related to cadmium accumulation in different plant tissues (roots or leaves) were identified, and functional enrichment analyses on target genes of differentially expressed miRNAs in low-and high-cadmium-accumulating wheat cultivars in different plant tissues (roots or leaves) obtained some known mature miRNAs and new miRNAs. The identified miRNA will be regarded as a potential screening biomarker for low-cadmium-accumulating wheat.

INTRODUCTION
Heavy metal pollution caused by extensive industrialization, urbanization, and improper wastewater treatment can affect the physical processes of plants (Ashraf et al., 2019(Ashraf et al., , 2020Fischer et al., 2020). Approximately 2.35 × 10 12 m 2 of arable land worldwide is contaminated by heavy metals (Li et al., 2017). Among the nonessential heavy metals, cadmium (Cd) is perhaps the metal that has attracted the most attention. Cd concentrations in durum wheat grain harvested in some areas of the northern Great Plains in the United States and adjoining regions of Canada have been reported to exceed 100 µg kg −1 dry weight (Guo et al., 2010). In China, approximately 2.786 × 10 9 m 2 of agricultural soil was polluted with Cd (Li et al., 2017). Several cereals uptake cadmium from soil and water and accumulate this metal in their shoots and edible parts. In general, cadmium enters from the root and then transfers to the stem, leaves, and grains by metal transporters. Wheat, a conventional and worldwide staple food source that accounts for 30% of calorie consumption, can also uptake cadmium and often exceeds the maximum standard exposure level (Alaux et al., 2018). Cadmium concentration in durum wheat grain harvested on Canadian prairies has reached 300 µg kg −1 , thereby exceeding the limit set by the Codex Alimentarius Commission (200 µg kg −1 ) . Decreasing the cadmium uptake in wheat to ensure grain quality and improve wheat growth is necessary for human health. Thus, the interaction between cadmium and wheat must be thoroughly understood.
Investigating microRNAs (miRNAs) can provide insights into plant biological processes. MiRNAs are a group of endogenous, single-stranded, short non-coding RNAs with partially selfcomplementary stem-loop structures that play key roles in gene expression regulation, predominantly in post-transcription via cleavage or translational repression of target mRNAs. In plants, these molecules regulate cellular differentiation, development, and metabolism, which are involved in responses to environmental changes, such as light, nutrition, temperature, and abiotic and biotic stresses (Axtell and Meyers, 2018). Several miRNAs actively play roles in wheat biological development. Debernardi et al. (2017) proved that miRNA172 plays a crucial part in wheat spike morphogenesis and grain threshold. Zhao et al. (2016) found that tae-miR408 can regulate the heading time in wheat via controlling the TaTOC1 gene. Guo et al. (2018) showed that miR9678 affects seed germination in wheat. Han et al. (2014) found that miR164 and miR169 show different expression patterns during seeding, suggesting that they perform different kinds of regulation in various developmental stages of wheat seeding. With high-throughput sequencing, different groups of miRNAs are found to participate in different biological processes, such as cell division, carbohydrate metabolism, and stress response. Meng et al. (2013) identified 104 miRNAs that are potentially involved in grain filling. Ravichandran et al. (2019) identified 104 miRNAs that regulate heat stress response, 36 of which are differentially expressed after heat stress. For metal toxicity response, Wang et al. (2018) found that miR319 can regulate cadmium, mercury, aluminum, and manganese uptake in Medicago truncatula. Zhou et al. (2012) reported that miR393 and miR396 regulate cadmium, mercury, and aluminum uptake in Brassica napus. Zhou et al. (2019) identified several miRNAs targeting heavy metal ATPases, which are responsible for Cd root-to-shoot translocation. Zhou et al. (2017) identified and characterized miRNAs in two cultivars: low (SJ19) and high-Cdaccumulating cultivars (CX4) of Brassica parachinensis. However, miRNAs related to wheat cadmium uptake remain to be explored.
Understanding wheat responses to cadmium will help to decrease the toxic metal uptake, thereby increasing food safety and quality. However, the yield penalty for harboring low cadmium in wheat and acreage either in China or in other countries, which grows low-cadmium accumulators, has not been reported. Based on field trials and cluster analysis in previous studies, 60 wheat cultivars were preliminary screened and divided into three categories: low-Cd accumulators, middle-Cd accumulators, and high-Cd accumulators in grain. In the following studies, the high-Cd accumulator (Qian 102032) and low-Cd accumulator (Annong9267) were verified again. The high-and low-Cd-accumulating wheat types have good reproducibility and stability in the annual accumulation, which are suitable for later cultivation and variety identification. In this work, two different cadmium-accumulating wheat (Triticum aestivum L.) cultivars, namely, Annong9267 and Qian 102032, were used as plant materials to identify miRNAs that respond to cadmium uptake. Then, next-generation miRNA sequencing was applied to profile the different responses of miRNA expression in low-and high-cadmium-accumulating wheat cultivars. Known and novel miRNAs associated with cadmium accumulation were identified, and their functions were analyzed to understand the genes related to cadmium uptake. The findings provided insights into the regulation of cadmium uptake and an opportunity to improve grain quality and reduce food cadmium contamination.

Wheat Growth Conditions and Cadmium Treatment
Two different cadmium-accumulating wheat (T. aestivum L.) cultivars, namely, Annong9267 (low-cadmium accumulator) and Qian 102032 (high-cadmium accumulator), were used as plant materials. Uniform seeds of wheat were surface sterilized in 10% NaClO and germinated at 25 • C for 3 days in the dark. The germinated seeds were grown in Petri dishes (diameter 18 cm, each containing 50 seeds) floating on half-strength 800 mmol m −2 s −1 Hoagland's solution in a growth chamber under a 12 h photoperiod, 70% relative humidity, and 25 • C/18 • C (day/night). The plants were grown hydroponically for 10 days and then transferred to the same nutrient solution containing 200 mmol L −1 CdCl 2 .
Based on previous studies, after 2 h of cadmium treatment, the content of cadmium accumulated in the two cultivars reached the peak. The leaves and roots of the seedlings treated with 200 mmol L −1 of CdCl 2 hydroponic culture for 2 h on the 10th day of growth were separately harvested, immediately frozen in liquid nitrogen, and stored at 80 • C. Six root samples and six leaf samples were collected and classified into four groups for future small RNA sequencing: Group 1 (samples 1-3), Annong9267 leaf samples; Group 2 (samples 4-6), Annong9267 root samples; Group 3 (samples 7-9), Qian 102032 leaf samples; and Group 4 (samples 10-12), Qian 102032 root samples. Each group contained three replications.

RNA Extraction
The frozen root and leaf samples were subjected to RNA extraction. Total RNA samples were extracted using the mirVana miRNA Isolation Kit (Cat#. AM1561, Austin, TX, United States) following the manufacturer's protocol. RNA integrities were evaluated using the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, United States), and all qualified RNA samples should meet the criteria of RNA integrity number (RIN) ≥ 7.0.

Small RNA Library Construction
Ligate 3 was added following the general workflow for small RNA library construction. In brief, 1 µg of total RNA was subjected to connect 3 adapter, heated using a polymerase chain reaction (PCR; Applied biosystems, 9700) instrument at 70 • C for 2 min, added to a ligation buffer, RNase inhibitor, and T4 RNA Ligation2, mixed, and reacted at 25 • C for 1 h. Stop solution was added and reacted at 28 • C for 15 min. The solution was immediately placed on ice, and 5 adapters were added. The 5 adapter mix system was placed in a PCR instrument at 70 • C for 2 min and then added into the solution to react at 25 • C for 1 h. After the reaction, the solution was immediately placed on ice.
For reverse transcription-polymerase chain reaction (RT-PCR), First-strand Master Mix, Super Script II (Invitrogen), dNTP mix, and 5 and 3 adapter RNA and RNase inhibitors were mixed and heated at 50 • C for 1 h. PCR Primer Cocktail and Master Mix were then added, and the solution was subjected to 98 • C for 30 s, followed by 11 cycles of 98 • C for 10 s, 60 • C for 30 s, and 72 • C for 15 s, maintained at 72 • C for 10 min, and held at 4 • C. The PCR products were purified by electrophoresis, and the gel stripe between 147 and 157 nt was cut off and recollected. The DNA was precipitated and dissolved in 10 µL of 10 mM Tris-HCl (pH 8.5).
The library was quantified by a Qubit 2.0 Fluorometer, and its quality was validated using an Agilent Technologies 2100 Bioanalyzer (Agilent Technologies). Cluster generation was then conducted, and 10 µL of 2 nM template DNA was subjected to cluster generation, hybridized to a paired end flow cell, and amplified by bridge amplification on the Illumina cBot. The flow cells were then loaded onto the HiSeqTM 2500 platform and sequenced.

Small RNA Sequencing and Analysis
Base calling was conducted to convert raw reads into sequence data. Impurities in raw data were then filtered. Reads with low quality and reads shorter than 10 nt, 5 primer contaminants and ploy (A), and reads without 3 adapter and insert tag were eliminated. After filtering, the remaining clean reads were stored in FASTQ format.
Clean reads with lengths between 18 and 40 nt were then mapped with wheat genome [international wheat genome sequencing consortium (IWGSC), version number miRBaseV22] (Wang et al., 2014) and annotated with miRBase 1 to identify known miRNAs using Bowtie software (version: v0.12.7) (Griffiths-Jones et al., 2007;Langmead and Salzberg, 2012). For novel miRNAs, the unannotated reads that could be aligned to the genome were analyzed using miRCat (version: V4.4.Alpha.D) (Moxon et al., 2008). When analyzing the folding model, if the sequence is located in the stem loop structure, then it is preliminarily determined as a new miRNA candidate.
The expression level of known annotated miRNAs and novel miRNAs was analyzed. The number of reads for each miRNA was normalized by the trimmed mean of M values (TMM) and then converted into transcripts per million (TPM, the number of reads on a miRNA × 10 6 /total read number) to standardize the expression of miRNA and allow the comparison for miRNA expression levels among different miRNAs and samples.
Differentially expressed miRNAs (DEMs) were identified by edgeR (version: v3.2.5) analysis and calculated p-values. After p-value calculation, multiple-hypothesis testing, and correction, the p-value threshold was determined by controlling the false discovery rate (FDR) to obtain the corrected p-value, namely, q-value. Multiples of differential expression were also calculated on the basis of the TPM value to obtain fold change. In the analysis, the p-value threshold was set to ≤0.05, and the fold change was set to ≥2. MiRanda (version: v3.3a) was used to predict whether the miRNA could bind to the 3 UTR of mRNAs to identify the targets of these DEMs (Turner, 1985).

Gene Ontology and KEGG Pathway Analysis
The differentially expressed miRNA target genes were analyzed by Gene Ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment using R studio based on the hypergeometric distribution (Kanehisa and Goto, 2000;Consortium, 2004). Target genes corresponding to the three levels of the biological process, namely, cellular component and molecular function, were first counted for GO enrichment analysis. The specific principle is to map the target genes of different miRNAs to the terms of the GO database, calculate the number of genes in each entry, and apply the hypergeometric test to screen for significantly enriched target genes compared with the entire genome background. For GO entries, the calculation formula is presented as follows: where N is the number of genes annotated with GO in all genes, n is the number of target genes in N, M is the number of genes annotated as a specific GO term in all genes, and m is the number of target genes annotated as a specific GO term. After the calculated p-value was corrected by multiple hypothesis testing, q ≤ 0.05 was used as the threshold. GO terms that met this condition were defined as significantly enriched in the differentially expressed genes. For GO enrichment analysis, emphasis was given for the rich factor, which was calculated as follows: (the number of target genes in a GO term/the number of all target genes correspond to the GO database)/(the number of genes contained in a GO term/the total number of genes correspond to the GO database). A high rich factor indicated a high degree of enrichment. Q-value was the corrected p-value after multiple hypothesis testing correction: when the Q-value is small, the enrichment is significant.
KEGG was used to identify significantly enriched metabolic or signal transduction pathways in DES target genes. The KEGG enrichment formula was similar to that of GO analysis. p-value was corrected by the Bonferroni method, and the threshold was set to q-value ≤ 0.05. KEGG terms that met this condition were defined as significantly enriched.

Deep Sequencing of Small RNA in Different Cadmium-Accumulating Wheat Cultivars
The roots and leaves from two different cadmium-accumulating wheat cultivars, Annong9267 (low-cadmium accumulator) and Qian 102032 (high-cadmium accumulator), were collected, and 12 sRNA libraries were created and sequenced to investigate different sRNA profiles in different cadmium-accumulating wheat cultivars.
The sequencing results showed that 19,755,465 to 48,544,084 raw reads were generated per library ( Table 1). After preprocessing, low-quality reads were removed. A total of 18,630,273 to 48,056,605 clean reads were filtered out from each library. The clean ratio ranged from 94.3 to 99.6%, thereby showing excellent sequencing quality.

Annotation of Small RNAs
After filtering, the clean reads were mapped to IWGSC and annotated by miRBase. Figure 1 and Table 2 show the distribution of catalogs of small RNA of mapped reads. Several known mature miRNAs and a number of unannotated reads were successfully obtained. This result would lead to the further identification of novel miRNAs. A total of 126 known mature miRNAs (Supplementary Table 1) and 63 pri-miRNAs (Supplementary Table 2) were identified.

Prediction of Novel Small RNAs
After annotation, the unannotated tags and predicted novel miRNAs were further analyzed by miRCat. A total of 751 novel miRNAs were predicted (Supplementary Table 5). Most of the novel miRNAs had a relatively low expression, but some had great abundance. In comparison with root small RNAs (Group 2 vs. Group 4), miRNAs such as novel.504 (only detected in group 2, p-value of 1.87E-19 and FDR of 3.30E-17) and novel.411 (log2 transformed fold change of 2.077441, p-value of 9.97E-06, and FDR of 0.00706) exhibited differential expression levels in the two aforementioned groups. In comparison with the leaves (Group 1 vs. Group 3), novel.526 was differentially expressed with log2 transformed fold change of 3.137454561, p-value of 5.30E-09, and FDR of 2.19E-07. These results indicated that these novel miRNAs were conservatively expressed in certain wheat cultivars, which might play important roles in cadmium uptake regulation.
In Group 2/4, nine DEMs were identified (Supplementary Table 4). Among them, six miRNAs were upregulated, of which tae-miR9664-3p was the most abundant. Novel.574 showed the largest fold change of 4.79. Three miRNAs were downregulated, of which tae-miR9663-5p was the most abundant, which showed the largest fold change of -3.58. Tae-miR9663-5p was expressed in the roots (with log2 transformed fold change of -3.58114 in the leaves, p-value of 2.41E-15, and FDR of 2.85E-13) and leaves (with log2 transformed fold change of -3.85717031, p-value of 8.84E-25, and FDR of 1.74E-22), and it showed similar expression change pattern in the two comparison groups. This finding indicated the potential downregulation role in cadmium uptake regulation. Figure 2B shows the heatmap of the expression level of DEMs in Group2/4. The two top candidates (novel.504 and novel.409) were expressed in both the leaves and roots of group 1.

Target Gene Functions of DEMs in Different Cadmium-Accumulating Wheat Cultivars
Some target genes of DEMs were also identified and subjected to gene function analysis. Thousands of genes targeted by FIGURE 1 | sRNA count distribution. Count distribution of different RNA categories in all 12 samples from four groups. sRNAs were categorized into multiple groups according to their biological classification. The first category "miRbase" indicates annotated some RNAs using the miRbase database. differentially expressed miRNAs were previously determined; however, analyzing their individual target gene function and biological effects would be impossible and laborious. Here, GO and KEGG terms were applied to describe the functional distribution of these targeted genes and understand the detailed correlations between the targeted genes of DEMs and their responsibility in cadmium uptake regulation. Gene Ontology assignments identified 3318 possible target genes for 89 DEMs. The functions were analyzed according to three categories: "biological process, " "cellular component, " and "molecular function." GO functional classifications are listed to illustrate the distribution of gene functions from the macro and biological effect levels.
In group 1/3, all target genes were largely distributed in nine GO terms (Figure 3A), of which the most frequent were "cell, " "cell part, " and "organelle" all clustered in the "cellular component" category. In "biological process, " the most frequent terms were "cellular process" and "metabolic process." In "molecular functions, " the most frequent terms were "binding" and "catalytic activity." In group 2/4, target genes were largely distributed in nine GO terms ( Figure 3B). Similar to root samples, the most frequent terms were all clustered in the "cellular component" category. In "biological process, " "single-organism process" also showed high frequency. In "molecular functions, " the most frequent terms were "binding" and "catalytic activity, " similar to that in group 1/3 results.
Gene Ontology enrichment analysis was conducted to understand the details of target gene functions. The top 30 enrichment terms of each comparison are shown in Figure 4. The rich factors were "AMP-activated protein kinase activity" and "histone H3-K9 deacetylation" in the root samples and "response to 1-aminocyclopropane-1-carboxylic acid" and "peptidyl-serine dephosphorylation" in the leaf samples. The term "spermidine synthase activity" was shared in both tissues, indicating its possible role in Cd uptake. Two major GO terms were differently enriched as "TORC1 complex" (significant in root samples) and "response to 1-aminocyclopropane-1 carboxylic acid" (significant in leaf samples), and one GO term was significantly shared and enriched as "spermidine synthase activity" in both tissues.
Pathway enrichment analysis based on KEGG database and generated reports for target genes in each pairwise was conducted to understand the interactions among target genes. The functions were analyzed according to six categories: "cellular processes, " "environmental information processing, " "genetic information processing, " "human diseases, " "metabolism, " and "organismal systems." KEGG classifications are listed to illustrate the distribution of pathway functions at the biological effect level. A bar plot for the statistics of KEGG term classification is shown in Figure 5, and a scatter plot to display the top 20 KEGG enrichment results is illustrated in Figure 6.
For KEGG term classification, the two comparisons showed similar results. The terms with most abundant gene numbers were "environmental adaptation, " "global and overview maps, " "carbohydrate metabolism, " "translation, " "folding, sorting, and degradation, " and "signal transduction." For KEGG enrichment analysis, the richest factors were "biosynthesis of the siderophore group non-ribosomal peptides, " "atrazine degradation" in root samples, and "atrazine degradation" in leaf samples. "Atrazine degradation" was shared in both tissues, indicating that the regulation of Cd uptake and atrazine degradation might share some common steps in the roots and leaves. The "AMPK signaling pathway" is a top shared pathway in both tissues. In particular, it was found that "porphyrin and chlorophyll II metabolism" were enriched in the roots, and "necroptosis" was enriched in the leaves, indicating potential tissue-specific cadmium intake and accumulation-associated pathways.

DISCUSSION
A comprehensive analysis was conducted on the miRNA expression patterns in low-and high-cadmium-accumulating  wheat cultivars with tissue-specific exploration on the roots and leaves by using small RNA sequencing techniques. New small RNA properties associated with cadmium accumulation and tissue specificity in wheat were identified at two major levels: (1) differentially expressed small RNAs associated with cadmium accumulation in different plant tissues (roots or leaves) and (2) functional enrichment analyses on target genes of DEMs in lowand high-cadmium-accumulating wheat cultivars in different plant tissues (roots or leaves). These results were validated on the basis of previously reported publications in wheat.

Differentially Expressed Small RNAs in Different Plant Tissues
A group of functional small RNAs were differentially expressed in low-and high-cadmium-accumulating wheat cultivars either in the roots or in the leaves. In the roots, emphasis was given on the two levels of functional small RNAs with either high abundance or differentially expressed patterns comparing low-and highcadmium-accumulating plants. Three typical small RNAs in the roots with relatively high expression level in the high-cadmium accumulation group and significant fold change comparing lowand high-cadmium-accumulating wheat cultivars were selected to filter and identify the significant small RNAs that contributed to cadmium accumulation.
In the roots, tae-miR9663-5p, a typical miRNA identified in wheat in 2014, was found to be a significant miRNA with TPM of 4158.93719 in high-cadmium-accumulating roots and 286.9844726 in low-cadmium-accumulating roots (log transformed fold change of -3.86, Supplementary Table 3). According to recent publications in 2014, this miRNA was upregulated during the initiation of wheat seeds, particularly during the pre-seedling stage, thereby confirming its functional role in wheat (Han et al., 2014). Another independent research  identified this miRNA as a typical miRNA biomarker for the identification of high-cadmium-accumulating genotypes, thereby validating the present result (Zhou et al., 2020). The abovementioned research selected the root as the target plant organ for studying cadmium accumulation (Zhou et al., 2020). Therefore, the first identified biomarker, tae-miR9663-5p, is a root-specific cadmium accumulation indicator in wheat. Another upregulated miRNA in high-cadmium accumulation roots is tae-miR6201 with TPM of 661.3927583 in low-cadmium accumulation roots and 107.1818898 in highcadmium accumulation roots (log transformed fold change of 2.625446091). According to recent publications, tae-miR6201 participates in the initiation of roots in wheat Li et al., 2019), thereby confirming its tissue specificity. For its correlations with cadmium accumulation, only a few publications have discussed the detailed functional correlations between tae-miR6201 and cadmium accumulation. Early in 2014, a systematic analysis on new wheat miRNAs with high-throughput sequencing techniques confirmed that such a miRNA was associated with the intake and metabolism of cadmium (Kurtoglu et al., 2014), thereby corresponding to our prediction. Therefore, the top miRNA, namely, tae-miR6201, may also be a cadmium accumulation-associated miRNA with root specificity.
In addition to root-specific miRNAs, a group of effective miRNAs were identified in the leaves of high-and low-cadmium accumulation groups. The number of typical miRNAs with differentially expressed levels in the leaves of high-and lowcadmium accumulation groups was lower than that in the roots. One of the significant miRNAs in the leaves had a similar fold change of -3.581, p-value of 2.41E-15, and FDR of 2.85E-13. Such a miRNA is functional in the roots for cadmium accumulation. However, the same publication confirming the presence of tae-miR9663-5p in wheat roots also validated a relatively high expression level of this miRNA in the flag leaves of wheat (Han et al., 2014). In 2019, an independent study (Zhou et al., 2019) on the effects of cadmium on the miRNA of the entire plant confirmed that plants with different Cd metabolism capacities might have different expression levels of typical miRNAs, including our target miRNA, namely, tae-miR9663. Therefore, as a shared typical miRNA biomarker for distinguishing high-and low-cadmium-accumulating plants in the roots and leaves, tae-miR9663-5p is an important potential biomarker for monitoring the cadmium accumulation capacity of plants in either the leaves or roots. Another top upregulated miRNA comparing low-and high-cadmium accumulating leaves is tae-miR9664-3p, which is also highly expressed in the leaves (Han et al., 2014). For its detailed correlations with cadmium intake and accumulation, tae-miR9664 is upregulated in the flag leaves of wheat (Ramachandran et al., 2020), thereby corresponding to the tissue specificity of our findings. For its correlations with cadmium accumulation, an independent study (Zhou et al., 2019) on the effects of cadmium on miRNA also identified this miRNA as a candidate cadmium intake and metabolism-associated miRNA, thereby validating our results.
A group of shared or unshared small RNAs that contributed to the regulation of cadmium metabolism and accumulation were identified. These molecules might be potential biomarkers for detecting and monitoring the cadmium accumulation status of wheat to promote the development of relevant agricultural technologies against cadmium pollution.

Functional Enrichment Analyses on Target Genes of DEMs
A group of miRNAs that can distinguish low-and highcadmium-accumulating plants in either the roots or leaves were identified. Given the complex role of miRNAs in gene expression regulation, a single or several DEMs cannot describe the comprehensive functional differentiation. Here, GO and KEGG pathway enrichment analyses were conducted to present the functional distribution profiling of the target genes for our candidate miRNAs. The top GO and KEGG terms in either the roots or leaves were selected for detailed analyses.
In the roots, two effective GO terms, namely, AMP-activated protein kinase activity and histone H3-K9 deacetylation, were identified as the top GO terms derived from the target genes of differentially expressed miRNAs, indicating their potential relationships with cadmium accumulation. According to recent publications, early in 2012, researchers from Poland confirmed that a specific protein kinase SNF1-related protein kinase type 2 contributed to the accumulation of cadmium in the roots via the activation of AMP-activated protein kinase activity (Kulik et al., 2012). Therefore, the target genes of significant small RNAs could be enriched in such GO terms. Furthermore, histone H3-K9 deacetylation has been widely reported to be correlated with the pathological effects of cadmium in different kinds of plants, including wheat (Zheng et al., 2016). As for KEGG enrichment results, a specific term, namely, "Atrazine degradation, " has shown to be the richest factor. Few publications discussed the relationships between Atrazine and cadmium. Only one publication in 2005 confirmed that the uptake of Atrazine might affect the accumulation rate of cadmium (Su et al., 2005). Therefore, this KEGG pathway might be one of the top enriched functional pathways, although further studies in this field were needed to reveal the detailed relationships between the two important metabolites in wheat.
Similar to the functional enrichment results in target genes of DEMs in the roots, atrazine degradation as a KEGG pathway term has also been identified in the leaves, indicating the specific relationship between atrazine and cadmium (Su et al., 2005). For GO terms, response to 1-aminocyclopropane-1-carboxylic acid and peptidyl-serine de-phosphorylation are two top terms, in which the target gene of DEMs enrich. According to early research on tomatoes (Liu et al., 2008), 1-aminocyclopropane-1-carboxylic is correlated with the improvement of cadmium tolerance, presenting a correlation between such a GO term and cadmium metabolism. Similar correlations have also been further validated in different plants, including carrot (Di Toppi et al., 1998), pea (Safronova et al., 2006), soybean (Chmielowska-Bak et al., 2013), and wheat (Milone et al., 2003). Therefore, in wheat, the target genes of our candidate small RNAs are correlated with the metabolism of cadmium, thereby validating our prediction. As for peptidyl-serine dephosphorylation, although few publications report the functional role of such biological processes in plants, in 2019, a specific study on Caenorhabditis elegans confirmed that such biological processes are correlated with weak cadmium stress via complex microbiological regulation, thereby validating our prediction (Dölling et al., 2019).
The top GO and KEGG enrichment terms are correlated with high-and low-cadmium accumulation in their respective plant tissues (roots or leaves), as supported by recent literature. This study identified a group of candidate miRNAs that might play a regulatory role for cadmium accumulation in different plant tissues and predicted various candidate regulatory pathways as the potential candidate biological mechanisms of cadmium accumulation for further validation.

CONCLUSION
High-throughput sequencing was applied to compare the small RNA profiling of high-and low-cadmium-accumulating wheat types in two independent tissues of roots and leaves. The differentially expressed small RNAs were identified as candidate biomarkers to evaluate cadmium metabolism in roots and leaves. The target genes of differentially expressed small RNAs were identified, and such genes underwent GO and KEGG pathway enrichment to reveal the biological mechanisms for different cadmium accumulation capacities in wheat. This study selected a group of cadmium accumulation-associated small RNAs and took the initial step on revealing the biological mechanisms of cadmium accumulation in wheat. Further experimental validation and plant model-based experiments were necessary to verify the current results. In our further plans, the SNP variations of the precursor and mature miRNA sequences between highand low-cadmium-accumulating wheat types with significant differences in cadmium accumulation will be identified to develop a biomarker for cadmium uptake.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: NCBI (accession: PRJNA698286).

AUTHOR CONTRIBUTIONS
YL and YZ designed the study. YL and XW performed the experiments. LY, YL, and TS analyzed the results. YZ wrote the manuscript. All authors contributed to the research and reviewed the manuscript.