Integration of small RNAs, degradome and transcriptome sequencing in hyperaccumulator Sedum alfredii uncovers a complex regulatory network and provides insights into cadmium phytoremediation

Summary The hyperaccumulating ecotype of Sedum alfredii Hance is a cadmium (Cd)/zinc/lead co‐hyperaccumulating species of Crassulaceae. It is a promising phytoremediation candidate accumulating substantial heavy metal ions without obvious signs of poisoning. However, few studies have focused on the regulatory roles of miRNAs and their targets in the hyperaccumulating ecotype of S. alfredii. Here, we combined analyses of the transcriptomics, sRNAs and the degradome to generate a comprehensive resource focused on identifying key regulatory miRNA‐target circuits under Cd stress. A total of 87 721 unigenes and 356 miRNAs were identified by deep sequencing, and 79 miRNAs were differentially expressed under Cd stress. Furthermore, 754 target genes of 194 miRNAs were validated by degradome sequencing. A gene ontology (GO) enrichment analysis of differential miRNA targets revealed that auxin, redox‐related secondary metabolism and metal transport pathways responded to Cd stress. An integrated analysis uncovered 39 pairs of miRNA targets that displayed negatively correlated expression profiles. Ten miRNA‐target pairs also exhibited negative correlations according to a real‐time quantitative PCR analysis. Moreover, a coexpression regulatory network was constructed based on profiles of differentially expressed genes. Two hub genes, ARF4 (auxin response factor 4) and AAP3 (amino acid permease 3), which might play central roles in the regulation of Cd‐responsive genes, were uncovered. These results suggest that comprehensive analyses of the transcriptomics, sRNAs and the degradome provided a useful platform for investigating Cd hyperaccumulation in S. alfredii, and may provide new insights into the genetic engineering of phytoremediation.


Introduction
With the food safety issues caused by heavy metals being highly emphasized worldwide, more and more intensive studies on the decontamination of polluted environments have been carried out in recent years (J€ arup, 2003;Schmitt, 2009). Among the several types of ubiquitous heavy metal pollutants, cadmium (Cd) is known as a harmful agent, which is toxic to all living cells by interfering with DNA repair machines, generating reactive oxygen species (ROS) and disturbing regular gene expression (Jomova and Valko, 2011;Valko et al., 2005). In most edible plants, Cd, as a nonessential element, can be competitively accumulated against the uptake of some essential metals, impeding vegetative growth and decreasing crop quality (Seth et al., 2008;Varotto et al., 2013). However, a special kind of plant, called hyperaccumulating, evolved strategies not only for extreme tolerance to Cd stress but also for hyperaccumulating Cd ions from polluted soils into its above-ground component (Verbruggen et al., 2009). Taking advantage of these plants for the remediation of heavy metals in contaminated soils, namely phytoremediation, has been regarded as an environmentally friendly and economical method (Lombi et al., 2001). Among the more than 400 identified hyperaccumulating plants, the hyperaccumulating ecotype of Sedum alfredii Hance is a Cd/zinc/lead co-hyperaccumulating species of Crassulaceae in China (Yang et al., 2002). Possessing a high-performance cellular detoxification system, this plant is a promising phytoremediation candidate that accumulates substantial heavy metal ions in shoots without obvious signs of poisoning (Yang et al., 2006). However, more studies at the molecular level are needed to illustrate the mechanism underlying its hyperaccumulating ability.
Cd toxicity in plants can indirectly lead to the overproduction of ROS, such as hydrogen peroxide (H 2 O 2 ), superoxide (O 2˙À ) and hydroxyl radical (OH˙) (Gill and Tuteja, 2010;Jaspers and Kangasj€ arvi, 2010), and subsequently causes lipid peroxidation, enzyme inactivation, cell membrane integrity damage, ion leakage and eventually cell death (Boominathan and Doran, 2003). The biochemical efficient enzymatic and nonenzymatic antioxidant defence systems protect plant cells from oxidative stress damages by preventing dangerous elevations of ROS levels (Gill and Tuteja, 2010;Petrov et al., 2015). Additionally, the increased ROS concentrations may induce plant programmed cell death (PCD) when the plants were exposed to Cd (Arasimowicz-Jelonek et al., 2012). Plants exposed to Cd can also induce signalling perception through an external signal, and the signal is rapidly transmitted to the responsible transcription factors (Lin and Aarts, 2012). Plant hormones together with some signalling molecules can regulate plant PCD (Gadjev et al., 2008). Thus, plants regulate the response to Cd exposure through a complex mechanism.
In recent years, the ever-accelerated updating of sequencing platforms has made it possible to conduct a large-scale association analysis of both mRNAs and small RNAs (sRNAs) in nonmodel plants (Morozova and Marra, 2008;Mutz et al., 2013). Although the whole transcriptome sequencing of the hyperaccumulating ecotype of S. alfredii was described in previous research (Gao et al., 2013), a more detailed comparative transcriptome analysis among different Cd stress durations was required to further elucidate the mechanism underlying its hyperaccumulation and hypertolerance of Cd. In the past few years, research has suggested that plant miRNAs and their targets play important regulatory roles in their adaptation to different heavy metal stresses (Gupta et al., 2014;Yang, 2013). Highthroughput analysis of these endogenous conserved miRNAs under heavy metal stress might provide a new viewpoint on stress-response mechanisms in plants (Hsieh et al., 2009). At present, many potential heavy metal-responsive miRNA candidates were identified in rice, medicago, Brassica juncea, radish and many other plants (Liu and Zhang, 2012;Wang et al., 2015;Zhou et al., 2012a,b). A number of miRNAs regulated by aluminium stress were detected in wild soya bean (Zeng et al., 2012). Brassica napus miRNAs and their targets were profiled on a genomewide scale in response to Cd (Zhou et al., 2012a). Furthermore, a genomewide identification of radish miRNAs and their targets was conducted under the stress of Pb exposure . However, little is known regarding miRNA expression profiles in S. alfredii in response to heavy metal stress. Therefore, there is also a need for a comprehensive identification of miRNAs involved in heavy metal absorption and resistance in the hyperaccumulating ecotype of S. alfredii.
The objective of this study is to systematically identify potential heavy metal-responsive miRNAs and their targets in the hyperaccumulating ecotype of S. alfredii. Thus, the combination of the transcriptomics, sRNAs and degradome generated a comprehensive resource focused on identifying key regulatory miRNAtargeted circuits under Cd stress. Our results will provide valuable information regarding Cd-responsive mechanism in the hyperaccumulating ecotype of S. alfredii, and help to improve the genetic engineering of phytoremediation.

Transcriptome sequencing in hyperaccumulating S. alfredii under Cd treatments
The transcriptome library was constructed from mixed RNA pools consisting of root, stem and leaf samples from plants treated with 400 lM CdCl 2 for eight durations. Then, 44 495 628 raw reads were generated, accounting for approximately 8.8 Gb of sequencing data. After quality control, the remaining high-quality reads were assembled into 87 721 unigenes, having an average length of 586 bp. The summary of the Illumina transcriptome sequencing of S. alfredii is shown in

Sequencing and identification of known and novel miRNAs
The small sRNA library was assembled from eight batches of sequencing data sets to increase coverage (Table S1). After removing the low-quality sequences, 17-to 25-nt-long sequences were obtained. The size distribution of the unique sRNAs is summarized in Figure S1, and no significant differences were found among the eight libraries. The displayed lengths of these S. alfredii sRNAs ranged from 18 to 25 nt, and the 24-nt sRNAs displayed the highest redundancy size class among all eight libraries.
All clean and unique sRNAs read sequences were queried using the BLAST algorithm to search against the relevant noncoding RNAs deposited in the Rfam and NCBI GenBank databases to remove rRNA, tRNA, snRNA and snoRNA sequences. To predict novel miRNAs, the remaining unannotated sRNA unique reads that did not match the transcriptome of any library were retrieved and subjected to a secondary structure prediction. All of the locigenerating sRNAs that could be folded into a secondary structure were considered as potential novel miRNA candidates. In total 356 miRNAs were found in S. alfredii and categorized into four different groups based on their abundance in the miRNA database and sequencing reads ( Figure 1, Table S2). Of these miRNAs, 260 pre-miRNAs corresponding to 299 known unique mature miRNAs were identified as having high similarity compared with known plant miRNAs. These miRNAs belonged to 82 already documented miRNA families ( Figure S2). The 21 nt length occupied 41.47% of the known unique miRNAs (Table S3). In total, 40 pre-miRNAs, corresponding to 57 mature miRNAs, were identified as novel miRNA candidates that were not registered in miRBase. The lengths of these miRNAs varied from 18 to 25 nt, with 47.37% being 24 nt long (Table S3). Based on Mfold calculations, the Minimum Free Energy (MFE) and MFE index of these predicted pre-miRNAs ranged from À23.9 to À147.1 kcal/ mol and 0.9 to 1.12 kcal/mol, respectively (Table S4). These characteristics meet the requirements to maintain the stability of the hairpin structures of miRNAs. To further identify authentic miRNAs without a reference in the Sedum genome, the secondary structures of pre-miRNAs and the reads coverage in mature and star areas were analysed. We found 128 highly confident miRNAs (Table S5), of which 30 miRNAs were not found in other plant species. Cd-responsive miRNAs in S. alfredii To identify miRNAs in S. alfredii that respond to Cd, the differential expression of miRNAs in the eight libraries was analysed and compared using the read counts generated from the high-throughput sequencing. In total, 153 miRNAs (P-value <0.05) showed differential expression patterns. When a stricter criterion of total expression abundance >10 was used, 79 differentially expressed miRNAs, including 67 known miRNAs and 12 novel miRNAs, were found among the eight libraries. The heatmap of potential differentially expressed miRNAs is illustrated in Figure 2a. Most miRNA members of the same family had similar expression profiles. For example, three miR408 and three miR171 family members were significantly up-regulated by Cd exposure. Five miR2916, three miR6173 and three miR8005 family members were found to be significantly down-regulated by Cd treatments. Among the 12 novel miRNAs, 10 miRNAs were up-regulated by Cd exposure, while the other two miRNAs (PC-5p-98092_120 and PC-5p-207477_71) were down-regulated ( Figure 2a). In addition to the validation of differentially expressed miRNAs, we also examined the distribution of 153 differential miRNAs between control and treatments ( Figure 2b). Among these miRNAs, the largest number (45 miRNAs) of up-regulated miRNAs was found with a CdCl 2 treatment of 48 h compared with the control, while 32 miRNAs were significantly down-regulated with a CdCl 2 treatment of 0.5 h compared with the control (Figure 2b). These data indicated that the duration of Cd treatments might cause short-term as well as long-term effects on the expression levels of miRNAs.

Target prediction of the known and novel miRNAs by degradome sequencing
Through degradome sequencing, 9 649 268 raw reads representing 4 721 650 unique reads, were generated from the mixed degradome pools. After removing the reads lacking the adaptor, 6 812 869 (70.6% of all reads) sequences were successfully mapped to the 48 446 unigenes (55.22% of the 87 721 input cDNA sequences) of the hyperaccumulating S. alfredii. The CleaveLand 3.0 pipeline was used to identify degraded targets for each of the miRNA families (Addo-Quaye et al., 2009). All of the potential cleaved transcripts were classified into five categories, based on the signature abundance at each occupied transcript position (Xu et al., 2013;Yang et al., 2013): category 0 -over one raw read at the position, abundance at the position is equal to the maximum on the transcript, and there is only one maximum on the transcript; category 1-over one raw read at the position, abundance at the position is equal to the maximum on the transcript, and there is more than one maximum position on the transcript; category 2-over one raw read at the position, abundance at position is less than the maximum but higher than the median for the transcript; category 3-over one raw read at the position, abundance at the position is equal to or less than the median for the transcript; and category 4-only one raw read at the position. Among these identified targets, there were 171, 16, 216, 8 and 346 in categories 0, 1, 2, 3 and 4 respectively (Table S6, Figure S3).
A total of 722 targets for 180 known miRNAs were identified from the mixed degradome (Table S6). We also identified the targets of 57 novel newly discovered miRNAs by degradome sequencing and analysis. Finally only 14 novel miRNAs targeted 32 transcripts ( Table 2). Most of the miRNAs (187) cleaved a single transcript target, while 29 miRNAs were detected to cleave two or more different transcript targets (Table S6). Forty-seven targets were identified to be cleaved by ptc-miR171i-p3_1ss21TG, which was the highest amount of transcripts cleaved by the same miRNA; while 11 miRNAs were detected to cleave only one transcript target, which was annotated to an unknown protein with hydrolase activity.

Annotation and enrichment analysis of targets for miRNAs under cadmium stress
An analysis using the BLASTX algorithm showed that the identified target transcripts detected in our degradome library shared homology with other plant proteins. A gene ontology (GO) functional classification analysis was conducted to further understand the functions of these 754 identified targeted genes ( Figure 3a). These target genes predominantly participated in biological processes, cellular component and molecular function. Fifteen transcripts were identified in biological processes, with the two most frequent categories being 'regulation of transcription' and 'transcription'. Of the 15 molecular component categories, the two most highly represented were 'nucleus' and 'integral to membrane'. Finally, there were 15 molecular function categories, with the most abundant being 'DNA binding' (Figure 3a). The functional annotations indicated that there were five transcription factors, NAC82, NAC90, HD-Zip III,  Wangeline et al., 2004). This analysis suggested that the miRNA targets were concentrated in the transcription factors and respond to stimuli.
An analysis of the GO term enrichment indicated that key biological pathways were involved in the Cd responses. Significantly enriched GO terms (Fisher's test P-value <0.05) were visualized using the reviGO tool (Supek et al., 2011) and plotted ( Figure 3b). Some interesting GO terms related to metal transport (GO: 0000287 magnesium ion binding and GO: 0005234 extracellular-glutamate-gated ion channel activity), redox-related secondary metabolism (GO: 0016706 oxidoreductase activity, GO: 0000822 inositol hexakisphosphate binding, GO: 0008295 spermidine biosynthetic process), and phytohormone auxin (GO:

Correlation analysis of miRNAs expression profiles and their target genes
To investigate the variable trend of the global transcriptome in the hyperaccumulating S. alfredii under different Cd stress treatment durations, the previously assembled unigenes were applied to a DEG analysis. In total, 342 target genes presented differential expression levels under Cd stress. There were 45 differentially expressed targets from 40 miRNAs (Table 3) when the targets corresponded to differentially expressed miRNAs, and 39 miRNA-target pairs showed a reverse expression pattern (negatively correlated; Figure 4). For example, cme-miR858_L-1R+1 was significantly down-regulated by Cd exposure, while its target MYB12 (BMK. 57978) was significantly up-regulated.

Gene coexpression network analysis
To determine the genes related to the coexpressed target genes, the DEG data set was used to construct a network. After excluding nonexpressed and lowly expressed genes, we identified 60 gene coexpression modules containing 31 985 unigenes ( Figure S4). Ultimately, the resulting network was composed of 31 distinct gene modules. The constructed gene modules ranged in size from 4238 (turquoise module) to 136 genes (plum1 module). To investigate the network connections for the target genes, we focused on the two target genes ARF4 (BMK. 50544, auxin response factor 4, PC-5p-3760_2402 target gene) and AAP3 (BMK. 21215, amino acid permease 3, han-miR3630-3p_L-2R-1_1ss21TA target gene), which play a prominent roles in heavy metal stress (Gielen et al., 2012;Xu et al., 2012). A complete list of module assignments and network metrics for the two genes is included in Table S7. In this network, the two hub genes ARF4 and AAP3 were directly connected with 1754 edges and 1212 edges. Among them, 1115 edges were coregulated by the two hub genes, suggesting that the coexpressed genes were most likely coregulated.
We also created a subnetwork containing 91 genes by identifying five categories that were related to signal transduction, transcription factor, antioxidant-related genes, PCD and metal transport (Figure 6). Of the 91 genes, 32 and 36 were in the signal transduction and transcription factor categories respectively when testing the functional assignments (Table S8). The two categories may trigger transport processes for detoxification, activate biochemical defence reactions, and PCD under Cd stress. Seventeen metal transport-related genes might be involved in Cd uptake, transport and sequestration processes. Five members of the peroxidase family found in the subnetwork may remove H 2 O 2 formed because of Cd tress, and lead to the activation of other defence mechanisms. When the plants experience heavy metal tress, plant PCD, as a positive and negative aspect of environment adaptation, provides survival benefits for the whole plant. More importantly, this network can be exploited in future studies to identify novel genes that interact with ARF4 and AAP3.

Discussion
In comparison with miRNAs identified from the other plants, little research has been conducted on the role of miRNAs in Crassulaceae plants. Only 146 EST sequences of S. alfredii have been deposited in GenBank, which is insufficient to analyse miRNA behaviour under Cd stress. Here, a transcriptome data set was used as reference sequence for miRNA and degradome sequencing analyses in S. alfredii. The approach highlighted the benefits of providing more qualitative and quantitative descriptions of the miRNA functions during stress tolerance. Furthermore, 356 miRNAs, including 128 highly confident miRNAs, were identified based on the transcriptome data in the hyperaccumulating S. alfredii. Of these miRNAs, more than half (68.35%) had relatively low expression abundances, indicating that highthroughput sequencing is a most powerful strategy to identify miRNAs with low expression levels in plants. Additionally, the sRNA length distribution patterns peaked at 24 nt ( Figure S1), which is consistent with previous results for most angiosperms, such as peanut (Zhao et al., 2010), soya bean (Song et al., 2011), Medicago truncatula (Szittya et al., 2008), Paulownia australis (Niu et al., 2014) and cotton . Therefore, the sRNA length distribution pattern in S. alfredii is similar to that in other plants.
In this study, we identified 79 miRNAs differentially regulated by Cd exposure (Figure 1). Most miRNA members of the same family showed similar expression profiles. For instance, three miR396 and three miR171 family members were significantly upregulated by Cd exposure. In M. truncatula, it was reported that miR171 was also up-regulated by Cd stress (Zhou et al., 2008). However, miR171 in rice was reported to be down-regulated in response to Cd treatment (Ding et al., 2011;Tang et al., 2014); miR171 and miR396 were also down-regulated by Cd exposure in B. napus (Xie et al., 2007;Zhou et al., 2012a). Thus, the differential expression patterns of miR171 and miR396 might be due to the different Cd treatment conditions of plant species. We also discovered that some previously known miRNAs were regulated by Cd stress. Five miR2916, three miR6173 and three miR8005 family members were found to be significantly downregulated by Cd treatments (Figure 2a). In plants, miR2916 was found to be responsive to abiotic stresses, such as salinity (Qin et al., 2015) and boron stress (Ozhuner et al., 2013). In Saccharina japonica, miR6173 was reported to regulate the response to heat stress . These findings implied that the known miRNAs were most likely involved in cross-adaptation to regulate plant tolerance to abiotic stresses, including heavy metal stress. However, their exact functions remain to be verified in  future investigations. Different miRNA members within the same family were present at different expression levels after Cd exposure. Similar events also occured with miRNA family members in other plants. In soya bean, miR396a-3p and miR396i-3p were significantly up-regulated, while miR396b-5p was downregulated by Cd exposure (Fang et al., 2013). Thus, a complex mechanism may regulate the expression profiles of members of the same family. A transcriptome-wide analysis of the degradome was performed, and numerous target transcripts for the known and novel miRNAs were determined. For the miRNAs, 754 targets were identified by degradome sequencing (Table S6). Most targets of known miRNAs are conserved, including transcription factors, and those involved in signal transduction and plant response to heavy metals, all of which are found in other plant species (Cao et al., 2013;Zhang et al., 2012Zhang et al., , 2014. Recent studies showed that a number of miRNAs regulate the expression of their targets through transcript cleavage or translation repression (Brodersen et al., 2008). An integrated analysis of miRNA expression profiles and their targets can help to identify the functional miRNA-target modules involved in regulating specific biological processes (He et al., 2014;Pei et al., 2013). In this study, we identified 45 differentially expressed targets from 40 differentially expressed miRNAs by DEG analysis and RT-qPCR (Table 3), and most of these miRNA-target pairs showed reverse expression patterns (Figure 4). For example, cme-miR858_L-1R+1 was significantly down-regulated by Cd exposure, while its target MYB12 (BMK. 57978) was significantly up-regulated. However, miR858 was also reported to be up-regulated in roots exposed to different Pb concentrations in cotton (He et al., 2014). However, many miRNA-target pairs involved in heavy metal stress in previous reports were not detected in our study, such as miR604/LTP (lipid transfer protein) (Huang et al., 2009) and miR167/NRAMP (Burklew et al., 2012;Srivastava et al., 2013). S. alfredii may have evolved some unique miRNA-target regulatory mechanisms in response to Cd stress.
It was revealed that two hub genes, ARF4 and AAP3, might play central roles in the regulation of Cd-responsive genes by a coexpression regulatory network construction. Amino acid transport is highly regulated by environmental signals, such as  light, cold, high salt and drought, and mechanical signals (Liu and Bush, 2006). In Populus, the expression levels of most AAP genes were down-regulated by salt stress (Wu et al., 2015). Additionally, three amino acid transporters showed downregulated expression patterns in S. nigrum roots with a Cd treatment (Xu et al., 2012). Similarly, our investigation indi-cated that the expression of AAP3 in S. alfredii exhibited a down-regulated pattern under Cd stress (Figure 5i). Furthermore, many ARF transcription factors, as critical components of auxin signalling, play important roles in regulating various abiotic stressors (Jain and Khurana, 2009;Wang et al., 2010). The expression levels of numerous ARF genes change when  plants respond to abiotic stresses in soya bean (Ha et al., 2013), Sorghum bicolor (Sekhwal et al., 2015) and banana (Hu et al., 2015). In this study, we found that the expression of ARF4 was up-regulated by Cd stress (Figure 5j). Thus, AAP3 and ARF4 are potentially relevant to Cd responses in S. alfredii.
In conclusion, our study reported an integrated analysis of the transcriptomics, sRNAs and degradome to generate a comprehensive resource focused on identifying key regulatory miRNAtarget circuits under Cd stress. A total of 87 721 unigenes and 356 miRNAs were identified, and 79 miRNAs were significantly differentially expressed and identified as Cd-responsive miRNAs. The target genes for Cd-responsive miRNAs functioned in signal transduction, response to heavy metal stress and secondary metabolite pathways. Furthermore, 45 differentially expressed targets were screened from 40 differentially expressed miRNAs using an integrated analysis of miRNA-target expression profiles. Ninety-one coexpressed genes formed a coregulation subnetwork in which two hub targets, ARF4 and AAP3, might play a central role in controlling the transcriptomic regulation in response to Cd. These results will facilitate a comprehensive understanding of Cd hyperaccumulation in S. alfredii, and help elucidate miRNA-mediated molecular mechanisms underlying plant responses to Cd stress and provide insights into Cd phytoremediation.
All of the samples were frozen in liquid nitrogen immediately, and stored at -80°C. Total RNA was extracted using the Total RNA Purification Kit (NORGEN, Thorold, Canada). All RNA samples were treated with DNase I (TaKaRa, Dalian, China) to avoid genomic DNA contamination. RNA quality and purity were checked using denaturing 1.0% (p/v) agarose gel electrophoresis and a NanoDrop 2000 spectrophotometer (Thermo, Wilmington, DE) at 260/280 nm (ratio >2.0). The total mixed RNA samples, consisting of roots, stems and leaves, treated with 400 lM CdCl 2 for eight different durations, were prepared for de novo transcriptome and degradome sequencing. Then, equal parts of RNA roots, stem and leaf samples for the same time point were mixed together for sRNA sequencing and DEG analysis.

Transcriptome sequencing and de novo assembly analysis
The sequencing library for transcriptome analysis was prepared according to Illumina's kit. mRNA with Poly (A) was purified from the total RNA, using oligo(dT) magnetic beads, and then fragmented with an RNA fragmentation kit. First-strand cDNA was synthesized using six random hexamer primers and a template of short fragments. Second-strand cDNA was synthesized using the buffer, dNTPs, RNase H and DNA polymerase I and then ligated to sequencing adapters. Following the isolation of suitable cDNA fragments by gel electrophoresis analysis, and PCR enrichment, the final products were loaded onto an Illumina HiSeq TM 2000 platform for transcriptome sequencing. Total reads were trimmed of low-quality and adapter sequences using the Illumina Pipeline, and then the remaining high-quality reads were assembled into nonredundant unigenes using Trinity (http://trinityrnaseq.github.io/). Unigenes were tentatively identified based on the best hits against known sequences in the database.

sRNA sequencing and miRNAs identification
Equal amounts of RNA samples from roots, stems and leaves at each of the eight different Cd treatment durations were pooled to generate eight sRNA libraries. The raw RNA reads generated by NGS from the eight sRNA libraries were processed to remove 5' and 3' adapters, and contaminated and low-quality sequences, as well as those smaller than 18 nt, through Illumina's Genome Analyzer Pipeline V1.5. Then, the filtered reads were subjected to a further filtration step, to remove the common RNA families (rRNA, tRNA and snRNA) with a proprietary pipeline script ACGT101-miR v4.2 (LC Sciences, Houston, TX). Then, the remaining clean and unique reads were aligned against the latest miRBase database, version 20.0 (http://www.mirbase.org/), using the BLAST algorithm to identify known miRNAs. The stem-loop hairpin structures were aligned to sequencing reads and mature miRNAs from miRBase using Bowtie (Langmead et al., 2009). Four groups of miRNAs were revealed through bioinformatics analysis of sRNA sequencing based on the classification method (Ambady et al., 2012). The reads distribution was checked to meet the principles for miRNA prediction and authentic miRNAs were regarded as described previously (Tarver et al., 2012;Taylor et al., 2014). All of the sequencing data set, including transcriptome sequencing and sRNA sequencing results, were deposited into the NCBI SRA database under accession number SRP058333.

Degradome sequencing, target identification and analysis
Equal amounts of all 24 RNA samples were mixed together to generate one degradome library. Then, the eight sRNA mixed samples and one degradome mixed sample were sent to Hangzhou LC-Bio Co., Ltd (Hangzhou, China) for cDNA library construction and Illumina Genome Analyzer GA-I Sequencing (Illumina, San Diego, CA). The extracted sequencing reads generated by degradome sequencing having 20 and 21 nucleotides were used to identify potentially cleaved targets by the CleaveL and 3.0 pipeline. Then, the degradome reads were mapped to the hyperaccumulating S. alfredii transcriptome data. The targets selected were categorized as 0, 1, 2, 3 and 4 as in a previous study (Xu et al., 2013;Yang et al., 2013). Based on the signatures (and abundances) along the hyperaccumulating S. alfredii transcriptome data, t-plots were built for the highefficiency analysis of the potential miRNA targets. Finally, all of the identified potential target genes were subjected to an NCBI search using the BLASTX algorithm and a GO analysis.

Differentially expressed target gene analysis
To discover the expression profiles of the target genes, eight independent libraries were constructed from each of the RNA samples from eight different Cd treatment durations. For each library, all of the sequences were processed to filter out adaptor sequences and low-quality sequences. Then, all of the clean tags were mapped to the assembled unigenes of the hyperaccumulating S. alfredii for annotation. The reads per kb per million reads method was used to calculate the gene expression level. Then, a rigorous algorithm method was performed to identify the differentially expressed genes between two samples. The false discovery rate method was used to determine the P-value threshold in multiple tests and analyses. The significantly differentially expressed genes among all of the different samples were judged by the threshold as follows: P-value <0.005, false discovery rate ≤0.001 and the absolute value of log2 ratio ≥1. The heatmap of the differentially expressed miRNAs and their corresponding targets was constructed using the ggplot2 package in R (version 3.1.3).

RT-qPCR validation
To validate the high-throughput sequencing results, RT-qPCR on miRNAs and their targets was performed on a 7300 Real-Time PCR System (Bio-Rad, Hercules, CA). We chose 12 miRNA-target pairs, including 12 differentially expressed miRNAs and 20 targets, for the RT-qPCR analysis. The RNA samples used for RT-qPCR were the same as for the experiments mentioned above. The reverse transcription reactions for miRNAs and cDNA were carried out using the SYBR â PrimeScript TM miRNA RT-PCR Kit (TaKaRa, Dalian, China) and the Superscript III First-Strand Synthesis system followed by RNase H treatment (Invitrogen, Carlsbad, CA), respectively. The specific miRNA forward primers were designed based on the sequences of the respective miRNAs, and the primers for the target genes were designed using the online Primer3 program (http://frodo.wi.mit.edu/ primer3/). Beta-tubulin (TUB) was selected as a reference gene (Sang et al., 2013). The gene names, sequences and the primers used for RT-qPCR analysis are available in File S1. RT-qPCR reactions for miRNAs and their targets were performed in 96well plates using a SYBR â Premix Ex Taq TM Kit (TaKaRa, Dalian, China) and SYBR â PrimeScript TM miRNA RT-PCR Kit (TaKaRa, Dalian, China). The amplification procedure and further data analysis were performed according to a previous study (Han et al., 2012). Each PCR reaction was repeated three times independently.

Coexpression network construction
Construction of the network was performed using WGCNA to identify modules of highly correlated genes based on the DEG data (Langfelder and Horvath, 2008;Zhang and Horvath, 2005). Expression values with more than four zeros for each gene were selected and log2 transformed before being processed through the WGCNA package (version 1.46) in R (version 3.1.3). An extensive overview of WGCNA explaining the analysis steps can be found on the webpage: http://labs.genetics.ucla.edu/horvath/ CoexpressionNetwork/Rpackages/WGCNA/Tutorials/index.html (Langfelder and Horvath, 2008). For quality control, we performed a clustering and principal components analysis based on the gene expression levels of the DGE data. The correlation networks were produced using the WGCNA package with the default power of six (Venables et al., 2012). All of the other parameters were used with the default values. Eigengenes were calculated for each gene coexpression module to visualize the gene expression patterns. The coexpressed genes with strong interconnection were defined as hub genes. Subsequently, we selected the differentially expressed target genes from the hub genes to analyse their edges. Results of coexpression analysis were visualized in Cytoscape (version 3.2.1) as described by Saito et al. (2012).

Supporting information
Additional Supporting information may be found in the online version of this article: Figure S1 Length distributions of the unique sRNAs. Figure S2 Numbers of identified miRNAs in known miRNA families in Sedum alfredii. Graphical representation of the different members of conserved miRNA families by sequencing and bioinformatics prediction. 'Other' represents 34 of the known miRNA families containing only one member. Figure S3 Five categories based on target cleavage positions. Figure S4 WGCNA coexpression modules based on DEG data. Table S1 Overview of sRNA sequencing. Table S2 The list of miRNAs identified through the bioinformatics pipeline. Table S3 Length distribution of unique miRNA. Table S4 Profiles of novel miRNAs originating from predicted RNA hairpins. Table S5 miRNA genes supported by miRBase and reads coverage. Table S6 The predicted targets of known and novel miRNAs. Table S7 A complete list of ARF4 and AAP3 module assignments. Table S8 The edge genes annotated in the coexpression subnetwork of ARF4 and AAP3. File S1 miRNA and their targets names, sequences and primers for RT-qPCR.