Identification and Comparative Analysis of Cadmium Tolerance-Associated miRNAs and Their Targets in Two Soybean Genotypes

MicroRNAs (miRNAs) play crucial roles in regulating the expression of various stress responses genes in plants. To investigate soybean (Glycine max) miRNAs involved in the response to cadmium (Cd), microarrays containing 953 unique miRNA probes were employed to identify differences in the expression patterns of the miRNAs between different genotypes, Huaxia3 (HX3, Cd-tolerant) and Zhonghuang24 (ZH24, Cd-sensitive). Twenty six Cd-responsive miRNAs were identified in total. Among them, nine were detected in both cultivars, while five were expressed only in HX3 and 12 were only in ZH24. The expression of 16 miRNAs was tested by qRT-PCR and most of the identified miRNAs were found to have similar expression patterns with microarray. Three hundred and seventy six target genes were identified for 204 miRNAs from a mixture degradome library, which was constructed from the root of HX3 and ZH24 with or without Cd treatment. Fifty five genes were identified to be cleaved by 14 Cd-responsive miRNAs. Gene ontology (GO) annotations showed that these target transcripts are implicated in a broad range of biological processes. In addition, the expression patterns of ten target genes were validated by qRT-PCR. The characterization of the miRNAs and the associated target genes in response to Cd exposure provides a framework for understanding the molecular mechanism of heavy metal tolerance in plants.


Introduction
Heavy metal contamination in soil is a serious environment issue over the world, originating mainly from a set of anthropogenic activities, such as mining, industrial activities and utilization of phosphate fertilizers [1,2]. In soil cadmium (Cd) can be easily taken up by plants, which results in various toxicity symptoms, such as reduced biomass, leaf chlorosis, inhibition of root growth, and morphological alterations, even plant death at excessive exposure [3]. In order to survive from stress environment, higher plants possess six possible ways to resist heavy metal exposure at the cellular level: (1) reduce metal bioavailability; (2) control metal influx; (3) chelate metals; (4) promote metal efflux; (5) compartmentalize and sequester metals; and (6) detoxify metal-induced reactive oxygen species (ROS) [4,5]. Efforts have been made to identify genetic elements that are involved in Cd detoxification in plants. For example, the natural resistance-associated macrophage protein (NRAMP) transporter is reported to function in Cd uptake in Arabidopsis thaliana and rice (Oryza sativa) [6,7]. On the contrary, heavy metal transporter 3 (HMA3) plays a role in sequestration of Cd into vacuoles in Arabidopsis, rice and soybean (Glycine max) [8][9][10]. In addition, two major quantitative trait loci for seed Cd accumulation, Cda1 and cd1, have been well verified by Jegadeesan et al. [11] and Benitez et al. [10,12], respectively. However, it is still unknown how the Cd tolerance is controlled by soybean genes or loci identified from these studies.
MiRNAs are one of the major types of endogenous, small, nonprotein coding single-stranded RNA, around 22-24 nucleotides in higher plants, which regulate gene expression at the posttranscription and translation levels [13][14][15]. Numerous studies have demonstrated that miRNAs implicate in most of the essential physiological processes in plants, including development regulation, signal transduction, and stress responses [16,17]. MiRNAs are important for plant responses to heavy metal stress [18][19][20]. For instance, miR393 and miR171 respond to heavy metals in Brassica (B. napus) [21,22] and Medicago (M. truncatula) [23]. Additionally, 19 potential novel miRNAs responsive to Cd were isolated by using conventional sequencing approaches [24]. Recently, the genome-wide discovery of new miRNAs and analysis of the expression patterns becomes easier by using new small RNA sequencing and microarray technologies. Ding et al. [25] verified that miR166, miR171, miR390, miR156 and miR168 responded to Cd stress in rice, and Zhou et al. [26] reported that miR396, miR397, miR398 and miR408 were related to Cd exposure in Brassica.
To thoroughly characterize the function of each miRNAs, it is necessary not only to accurately identify the miRNAs, but also to validate the targets and to explore their interactions. Generally, according to the perfect sequence complementarity between a miRNA and its target(s) or the conservation of miRNA targets among different plant species, computational target prediction has been globally adopted to evaluate targets of miRNAs in plant [27,28]. Nevertheless, due to the existence of the higher mismatch in miRNA-target pairing, computational target prediction is often questionable on the authenticity of predicted target transcripts [29]. Therefore, all the prediction targets need to be verified by experimental approaches. The modified 59 RACE (rapid amplification of cDNA ends) has been widely applied in target confirmation and cleavage site mapping [30]. However, this method is used only in small-scale due to time and labor consuming and higher cost. Fortunately, degradome sequencing, a high-throughput approach known as PARE (parallel analysis of RNA ends), has been adopted to validate the miRNA spliced target transcripts [31][32][33]. Recently, this technology has been successfully used to identify miRNA targets in Arabidopsis [31,32], rice [33], and soybean [34].
Soybean is one of the most important agriculture crops in the world. However, Cd pollution in soil affects the yield and quality of soybean seeds. Several soybean miRNAs have been found to play roles in response to biotic and abiotic stresses [35][36][37]. However, little is known about the Cd-responsive miRNAs in soybean, even less understanding on the targets of Cd-responsive miRNAs for the regulation. To understand the regulatory mechanism of miRNAs in response to Cd treatment, a custombuilt miRNAs microarray chip was designed and used to detect the expression of miRNAs in HX3 and ZH24 roots under Cd stress or Cd-free. A total of 26 miRNAs were identified in response to Cd exposure in soybean roots. To evaluate the target transcripts of the miRNAs, a high-throughput degradome sequencing was adopted using a small RNA library consisting of four samples described above. Fifty five targets cleaved by 14 Cd-responsive miRNAs were identified. In addition, a number of Cd-responsive miRNAs and target mRNAs in soybean have been validated by quantitative RT-PCR (qRT-PCR).

Plant materials and stress treatment
Two contrasting soybean cultivars, HX3 (Cd-tolerant) and ZH24 (Cd-sensitive), were used in this research. Seeds of each cultivar were sown in sterile quartz sand and left for 4 days (d) to germinate. Then, seedlings were transferred into 1/2 full-strength nutrient solution containing 2.5 mM KNO 3

Microarray assay analysis
The miRNA microarray assays were performed by a service provider, LC Sciences in Houston, TX. The array design was based on 694 conserved miRNA mature sequences of legumes downloaded in miRBase Release 18.0 (http://www.mirbase.org/) and 259 well characterized miRNAs published in papers [38][39][40][41][42][43] ( Table S1). Each probe was repeated three times on the chip to ensure reproducibility of microarray. The 5S rRNA was designed as the inner positive control; blank and non-homologous nucleic acids served as negative controls.
For each hybridization experiment, 2 to 5 mg total RNA sample, which was size fractionated using a YM-100 Microcon centrifugal filter (from Millipore) and the small RNAs (,300 nt) isolated were 39-extended with a poly(A) tail using poly(A) polymerase. An oligonucleotide tag was then ligated to the poly(A) tail for later fluorescent dye staining; Hybridization was performed overnight on a mParaflo microfluidic chip using a micro-circulation pump (Atactic Technologies) [44]. On the microfluidic chip, each detection probe consisted of a chemically modified nucleotide coding segment complementary to target microRNA and a spacer segment of polyethylene glycol to extend the coding segment away from the substrate. The detection probes were made by in situ synthesis using PGR (photogenerated reagent) chemistry. The hybridization melting temperatures were balanced by chemical modifications of the detection probes. Hybridization used 100 mL 6xSSPE buffer (0.90 M NaCl, 60 mM Na 2 HPO 4 , 6 mM EDTA, pH 6.8) containing 25% formamide at 34uC. Hybridization images were collected using a laser scanner (GenePix 4000B, Molecular Device) and digitized using Array-Pro image analysis software (Media Cybernetics).
Data were analyzed by first subtracting the background and then normalizing the signals using a LOWESS filter (Locallyweighted Regression) [45]. To identify Cd-induced miRNAs, a criterion of signal .500, fold change .2 and P,0.01 was used.

Microarray Data Deposition
All microarray data discussed in this publication have been deposited in NCBI's Gene Expression Omnibus and are accessible through GEO Series accession number GSE50388 (http://www. ncbi.nlm.nih.gov/geo/query/acc. cgi?acc = GSE50388).

Validation of miRNAs expression profile via stem-loop qRT-PCR
The miRNA expression was examined by qRT-PCR experiments. Total RNA was isolated from roots of contrast and treated plants at 6,12,24,48, 96 and 144h after Cd exposure and with TRIZOL reagent following the manufacturer's instructions (Invitrogen). RNA samples were treated with RNase-free DNase I (TaKaRa) to avoid amplification from genomic DNA. These samples were collected at the same time as those for miRNA microarray and degradome sequencing. The first cDNA strand was synthesized from total RNA using the MMLV-reverse transcriptase (Invitrogen) and miRNA specific stem loop primers, which were designed mostly according to the method described by Chen et al. [46] and Varkonyi-Gasic et al. [47]. Briefly, six to seven nucleotide tips pairing with the mature miRNA 39 end were linked to a self-looped sequence (GTCGTATCCAGTGCGTGTCGT-GGAGTCGGCAATTGCACTGGATACGAC) to make up the stem-loop reverse transcription primer. QRT-PCR analysis was carried out using the SsoFast EvaGreen Supermix kit (BIO-RAD) in a CXF96 (BIO-RAD) qRT-PCR System. Reaction conditions for thermal cycling were: 95uC for 3 min, 40 cycles of 95uC for 10 s, 57.0-63.3uC for 10 s and 72uC for 30 s. The annealing temperature (57.0-63.3uC) was adjusted to suit the amplification of the individual transcript. The sequences of stem-loop reverse transcriptase primers and miRNA-specific PCR primers are listed in Table S2. In each real qPCR experiment, each gene was run in triplicate with different cDNAs synthesized from three biological replicates. Relative fold changes of gene expression were calculated using the comparative DDC t method [48], and F-box was used as the reference gene [49]. Standard errors and standard deviations were calculated from replicates and significance was measured through Student's t-test at the level of 0.01,P#0.05 and P#0.01.

Sequencing of degradome library and bioinformatics analysis
Construction of degradome library of soybean roots differed considerably from past efforts [31,50] and followed as Ma et al. [51] with some modification. (1) Approximately 150 ng of poly(A)+ RNA was used as input RNA and annealing with Biotinylated Random Primers; (2) Strapavidin capture of RNA fragments through Biotinylated Random Primers; (3) 59 adaptor ligation to only those RNAs containing 59-monophosphates; (4) Reverse transcription and PCR; (5) Library was sequenced using the 59 adapter only, resulting in the sequencing of the first 36 nucleotides of the inserts that represented the 59 ends of the original RNAs. A Public software package, CleaveLand3.0 was used for analyzing sequencing data generated. All targets were classified into four categories based on the abundance of the resulting mRNA tag relative to the overall profile of degradome reads that matched the target [31,50]. In category I, the abundance of cleavage sequences was equal to the most abundant degradome sequences on the transcript, and there was only one maximum on the transcript; in category II, the abundance of the degradome sequences at the cleavage site was equal to the maximum abundance on the transcript, and there was more than one maximum on the transcript; in category III, the abundance of cleavage sequences was less than the maximum but higher than the median for the transcript; in category IV, the abundance at cleavage site was equal to or less than the median for the transcript. The optimized score thresholds were set to 4.5 for category I, 4 for category II, 3.5 for category III, and 3 for category IV. These thresholds were used to select the resulting output. The gene ontology (GO) analysis of the candidate target transcripts of the identified miRNAs in this research was performed using the AgriGO toolkit [52].

Degradome sequencing Data Deposition
All degradome sequencing data discussed in this publication have been deposited in NCBI's Gene Expression Omnibus and are accessible through GEO Series accession number GSE50063 (http://www.ncbi.nlm.nih.gov/geo/query/acc. cgi?acc = GSE50063).

Validation of the miRNA target gene expression profile by qRT-PCR
The expression patterns of the target genes identified by degradome sequencing were evaluated by qRT-PCR. Total RNA was isolated from soybean root as the method described above. One mg of total RNA was used for initiating the reverse transcription reaction, and incubated with oligo(dT) primer. The reverse transcription reaction was performed by MMLV-reverse transcriptase (Invitrogen) following the supplier's manual. Then, target gene primers were added to carry out the PCR array. A soybean housekeeping gene F-box [49] was used for the reference gene for qRT-PCR. Supplementary Table S3 shows the primer sequencing of the target gene and the F-box gene. In each qPCR experiment, each gene was run in triplicate with different cDNAs synthesized from three biological replicates. The DDC t method was used to determine the expression level differences among samples of root with Cd treatment and without Cd treatment [48]. Standard errors and standard deviations were calculated from replicates and significance was measured through Student's t-test at the level of 0.01,P#0.05 and P#0.01.

Identification and expression patterns of Cd-responsive miRNAs in soybean
In this study, miRNA microarrays containing 953 miRNA mature sequences were used to investigate the expression patterns of miRNAs in soybean roots under Cd stress (microarray data were deposited into the NCBI-GEO with accession no. GSE50388). There were 34 miRNAs in HX3 and 48 miRNAs in ZH24 which differentially expressed between Cd treatment and control with at least one signal .500 and P,0.01. To identify Cdresponsive miRNAs, a further strict criterion with fold change .2 was used. A total of 26 miRNAs were identified as Cd-responsive miRNAs ( Table 1). Among those miRNAs, 16 miRNAs from G. max, two miRNAs from wild soybean (Glycine soja), two miRNAs from Arachis hypogaea, and one miRNA from Vigna unguiculata were released in miRBase, while one miRNA (Gma-m040-5p) from G. max [42], three miRNAs (PN-miR397a_L-1, PN-miR1509b_R+1 and PC-15-5p) from G. soja [43] and one miRNA (Vun78330_1521_100) from V. unguiculata [41] were newly identified as novel miRNA in each species, respectively. According to the similarity of the sequences, Gma-m040-5p was characterized as one member of miRNA390 family, Vun78330_1521_100 was identified as one member of miRNA319 family, while, PC-15-5p, having 25 nucleotides, shared 14 common sequences with gma-miR6300 from 3 to 18.

Confirmation of differentially expressed miRNAs by qRT-PCR
To validate the expression profile of the soybean miRNAs obtained through the microarray assays, qRT-PCR was performed for 16 Cd-responsive miRNAs (gma-miR3522, gma-miR397a, gma-miR408, gma-miR408b-5p, gma-miR4996, gma-miR396a-3p, gma-miR398c, gma-miR1509b, gma-miR5037b, PC-15-5p, gma-miR319c, gma-miR1535b, Gma-m040-5p, gma-miR4403, gma-miR396b-5p, and Vun78330_1521_100) in the same samples used in the microarray ( Figure S1 and S2). The results from qRT-PCR were compared with those from the microarray (Figure 1). The expression patterns of all 16 Cd-responsive miRNAs obtained by qRT-PCR were mostly similar to the results from the microarray, except for those of Gma-m040-5p and gma-miR4403. The expression patterns of Gma-m040-5p and gma-miR4403 obtained by qRT-PCR in HX3 were not consistent with that fromthe microarray but were similar with that in ZH24. Moreover, gma-miR1535b was confirmed to be up-regulated in HX3 and downregulated in ZH24, and miR396b-5p was verified to be downregulated in both cultivars and showed different trends compared with miR396a-3p. Although the fold change of expression valued by qRT-PCR was not exactly identical to that calculated by the microarray, the differential expression trends were the same by both methods. The qRT-PCR results ultimately reflected consistency with the microarray data.

Identification of targets of Cd-responsive miRNAs in soybean
MiRNA target confirmation is a prerequisite to get deep understanding of the functional role of miRNAs. Only several targets of miRNAs from soybean have been predicted [35][36][37][38]53]. To identify more targets in soybean, particularly specific targets of Cd-stress-responsive miRNAs, high-throughput degradome sequencing was used (the degradome sequencing data were deposited into the NCBI-GEO with accession no. GSE50063). In total, we obtained 8913111 raw reads from the library which was constructed from a mixture of four samples (HX3-CK, HX3-Cd-treatment, ZH24-CK and ZH24-Cd-treatment). After removing the reads without the CAGAG adaptor, 5430126 unique raw-reads were obtained. The unique sequences were aligned to the G. max genome database, and 6516276 reads were mapped to the genome. The mapped reads from the library represented 51481 annotated G. max genes. The sliced target transcripts were categorized into four groups according to the relative abundance of the tags at the target mRNA sites.
In total, 376 targets that could potentially be cleaved by 204 miRNAs were identified (Table S5). There were 193, 22, 210 and 4 targets in categories I, II, III and IV respectively. Most of miRNAs cleaved two or more different transcription targets, while 35 miRNAs (17.16%) were detected to cleave only one transcription target. Twenty five transcriptions were identified to be cleaved by mtr-miR169m, which was the highest amount of transcriptions cleaved by the same miRNA in this study. Furthermore, gma-miR396e was identified to cleave 19 members of growth-regulating factor families and a gene (Gly-ma03g36370.1) encoding a protein of unknown function. Most targets for conserved miRNAs were conserved, but still some conserved miRNAs had non-conserved or novel transcripts. For instance, a transcript encoding GATA type zinc finger transcription factor family protein was identified for gma-miR398 (T-plots of the four of the targets are presented in Figure 2).
Fourteen out of the 26 Cd-stress-responsive miRNAs (9 miRNA families) were identified to cleave 55 transcripts by the degradome sequencing (Table 2). Gene ontology (GO) analysis was carried out by using the AgriGO toolkit [52]. The results revealed that 54 genes were annotated as being involved in 12 biological processes ( Figure 3). More than 72% of targets took part in cellular and metabolic processes. It is notable that 20 targets (37.04%) were involved in regulation of biological process which was more enriched in the targets of Cd-responsive miRNA than in the soybean genes as a whole. In contrast, the proportion of the Cdresponsive miRNA targets in response to stimulus (16.67%) was less than that of the whole soybean genes.
Seven out of 14 miRNAs which were all up-regulated in HX3 and ZH24 were identified to cleave 25 targets belonging to nine The dark grey shading bar represents the relative expression level of miRNAs in HX3 measured by microarray in responsive to Cd stress. The dark grey diagonal stripe shading bar represents the relative expression level of miRNAs in ZH24 measured by microarray in responsive to Cd stress. The french grey shading bar represents the relative expression level of miRNAs in HX3 measured by qRT-PCR in responsive to Cd stress. The french grey diagonal stripe shading bar represents the relative expression level of miRNAs in ZH24 measured by qRT-PCR in responsive to Cd stress. The results of qRT-PCR are mean 6 SD of duplicates of three biological replicates. Significance of the changes between HX3 and ZH24 at the levels of Cd/CK was checked with Student's t-test at the level of P#0.01 (shown as ''**''). doi:10.1371/journal.pone.0081471.g001 families as blue-copper-binding protein family, CDPK-related kinase family, copper/zinc superoxide dismutase family, cupredoxin superoxide dismutase family, GATA type zinc finger transcription factor family, laccase (LAC) family, plantacyanin family, prolyl oligopetidase family and uclacyanin family. Six down-regulated miRNAs in HX3 and ZH24 were found to slice 27 transcriptions belonging to ARM repeat superfamily, ATPase E1-E2 type family, CYS, MET, PRO, and GLY protein family, growth-regulating factor family, myb domain protein family, HCO3-transporter family and cycloidea and PCF transcription factor family and to cleave one target (Glyma18g03980.2) with unknown function. Two members of isopentenyltransferase gene family were found to be cleaved by gma-miR1535b which was up-regulated in HX3 but down-regulated in ZH24.

Discussion
Heavy metal, especially Cd, is toxic to plants and plants initiate a variety of tolerant strategies to overcome the Cd toxicity. Recently, increasing evidences have revealed that miRNAs played the crucial role on regulation of plant genes at the posttranscriptional in responding to metals stresses. Under Cd stress, miR397, miR408 and miR398 were reported to be up-regulated and miRNA396 and miRNA319 to be down-regulated in response to Cd exposure in the root of rice or Brassica [25,26]. In this study, a total of 26 Cd-responsive miRNAs were identified; almost half of which, such as gma-miR397a, PN-miR397a_L-1, ahy-miR408-3p, gma-miR408b-5p, gma-miR408b-5p, gma-miR396b-5p, ahy-miR398, gma-miR398a, gma-miR398c, vun-miR319b, gma-miR319c and Vun78330_1521_100, have the similar alteration patterns in HX3 and ZH24 as that of rice and Brassica mentioned above (Table 1 and Figure 1). In addition, we also found that two other miR396 family members, gma-miR396a-3p and gma-miR396i-3p, showed up-regulated in the both cultivars, which was contrary with that of gma-miR396b-5p and its homologous members. The expression patterns of gma-miR396b-5p in two soybean genotypes were very similar with that of miR396 in rice [25] and Brassica [26]. For gma-miR390a-5p, an apparently down-regulated miRNA in ZH24, it has the same alteration with that in the root of rice [25] and shoot of Brassica [26], but was contrary with that in root tissue of Brassica [26]. It has been found so far that some miRNAs, including gma-miR3522, gso-miR3522a, gso-miR3522b, gma-miR4996, gma-miR1535b, Gma-m040-5p, gma-miR4403, gma-miR1509b, PN-miR1509b_R+1, gma-miR5037b and PC-15-5p, were responsive to Cd only in soybean, in which, gso-miR3522a, gso-miR3522b (released in miRBase 18.0), PN-miR1509b_R+1 and PC-15-5p [43], were found before only in wild soybean.
Of the 26 miRNAs identified in our study, almost all of them showed similar expression patterns in HX3 and ZH24, except for gma-miR1535b, which was detected as being up-regulated in HX3 and down-regulated in ZH24. The similarly of miRNA regulation may represent the fundamental mechanism of adapting to Cd exposure, and the difference of regulated miRNAs, including those which had distinct expression trends or had the same alteration but showed apparent expression level, may explain the distinct Cd sensitivities between the two soybean cultivars.
Usually, miRNAs regulate the expression of target genes at the post-transcriptional level by mediating gene silencing in plants. Large amount of target genes were also predicted by computer arithmetic, while quite a few of which have been validated by experiment, like 59-RACE assay [53]. In this study, degredome library was prepared from root samples of two genotypes with or without Cd stress. Sequencing and analysis resulted in the identification of 376 targets for 204 miRNAs. Among these, 55 genes cleaved by 14 Cd-stress-responsive miRNAs were verified by degradome sequencing (Table 2). qRT-PCR results revealed that the expression patterns of five miRNAs were significantly negatively correlated with their corresponding target genes. ( Figures 4 and 5).These results suggest that miRNAs play crucial roles in response to Cd by mediating target genes silencing.
Many abiotic conditions including heavy metal result in oxidative stress in plants [54]. Several miRNAs are involved in the regulation of genes responsible for antioxidation. MiR398 is the first miRNA identified to regulate plant responses to oxidative stress [55]. The MIR398 promoter contains the GTAC motif that has an important feature in Cu responsiveness in Arabidopsis [56]. This motif is recognized by the SPL7 (SQUAMOSA promoter binding protein-like 7) transcription factor that directly binds to the promoter and activate transcription of MIR398 gene. SPL7 is homologous to copper response regulator1, the transcription factor required for Cu deficiency response in Chlamydomonas reinhardtii [57]. In addition, SPL7 is required for the expression regulation of other Cu-related miRNAs such as miR397, miR408, and miR857 [56]. In this study, miR397a, miR408 and miRNA398c showed almost the similar up-regulated alteration to response to Cd exposure, which might imply that SPL7 involved in regulation of Cu deficiency and Cd response in soybean (Figures 1, S1 and S2).
Three most abundant copper proteins, copper/zinc superoxide dismutase (Cu, Zn-SOD; CSD), plastocyanin (PC) and cytochrome c oxidase (COX) are essential for oxidative stress alleviation, photosynthesis and respiratory electron transport, respectively [58]. CSD which is one target gene of miR398, as ROS scavenger and a major copper protein, is important for stress tolerance and survival in plant [59]. The mode of Curegulating ROS signal transduction was well studied in A. thaliana [60,61]. Under high Cu stress, oxidative stress suppresses miR398 expression, which is essential for the accumulation of CSD1 and CSD2 mRNA levels for catalyzing the dismutation of superoxide radicals into H 2 O 2 and sink Cu. On the contrary, decreased CSD1 and CSD2 transcripts abundance allows efficient delivery of limited Cu to PC, which is essential for photosynthesis when Cu is limited [60,61]. It is likely with the concentration of Cu limited, excess Cd induces miR398 expression, which inhibits the function of Cu/ Zn/SOD (CSD) and further induces ROS accumulation [4]. In this study, miRNA398a showed up-regulation in response to Cd treatment in HX3 and ZH24 (Table 1). Glyma03g40280.2 (coding CSD1), one target of miRNA398a (Table S5), was identified to be down-regulated in two cultivars and the decreasing trend in HX3 was more obvious than that in ZH24 (data not shown). This result may imply HX3 could save more Cu for photosynthesis and a higher tolerance than that of ZH24 when exposed to Cd. In addition, miRNA398c, one homologues member of miRNA398a, was confirmed to be up-regulated in two cultivars by qRT-PCR. Glyma06g19680.1 and Glyma14g39910.1, two targets of miR-NA398c, was confirmed to be down-regulated in HX3 and ZH24 (Figures 4 and 5). Glyma06g19680.1 encodes zinc finger transcription factor protein, which plays an important role in stress tolerance [62]. As shown in Table S7, the expression of Glyma06g19680.1, however, was down-regulated in the two cultivars (Figures 4 and 5), and in HX3 its expression was higher than that in ZH24 under Cd treatment. Glyma14g39910.1 encodes Prolyl oligopeptidase (prolyl endopeptidase), an atypical serine protease, which hydrolyses peptides and peptide homones after proline in peptides up to around 30 residues long [63]. When suffer from Cd exposure, the expression of Glyma14g39910.1 decreased sharply in HX3 in comparing with that in ZH24 (Figures 4 and 5). This could help alleviate the harmful effect caused by the hydrolysis of protease.
MiRNA397a and miRNA408 are two other important Cdresponsive miRNAs, which are also regulated by SPL7. Glyma18-g42520.1encodes laccase, which have been provided involved in the lignification of A. thaliana stem with the lac17 mutant [64]. Several studies reported on an increased lignin synthesis upon metal treatment [65,66], and lignification is one of defense mechanism under Cd exposure in soybean root [67,68]. It is well known that lignin provides structural support, and contributes to plant defense mechanism against both biotic and abiotic stresses [69]. As shown in Table S7, Glyma18g42520.1, down-regulating by miRNA397a in two cultivars (Figures 4 and 5), showed higher expression level in ZH24 than that in HX3 under Cd stress. This result seems inconsistent with that HX3 which possess a higher Cd tolerance than ZH24. As one kind of copper protein, LAC should be considered. And Cd treatment may lead to Cu deficiency, which have been elucidated by Chmielnicka and Sowa [70] in rats. So, even higher expression level of laccase in ZH24 with Cd treatment, compared with that in HX3 with Cd treatment, might not form functional protein when it encounters Cu deficiency caused by Cd exposure. This may explain why HX3 possess a higher tolerance. In addition, miRNA408 was identified to cleaved cupredoxin, laccase and plantacynin by degradome sequencing. Glyma03260.2 and Glyma08g13510.1, encoding cupredoxin and planacynin (another two copper protein), separately, showed almost the same expression patterns in four mixture samples with that of laccase, though the expression of Glyma03260.2 showed lower in ZH24 than that in HX3 with Cd treatment (Figures 4, 5 and Table S7).
One of the primary symptoms of Cd toxicity is root growth inhibition, with the root apex being the most sensitive part of the root [71][72][73]. Several conserved miRNAs, such as miR160, miR164, miR167, miR390 and miR396, have been well characterized playing functional roles in root development [74][75][76][77][78]. Glyma15g19460 and Glyma17g35090.1, cleaved by miR396b-5p, encode for growth-regulating factor (GRF), which is known to act in a functionally redundant fashion to positively control cell proliferation and size in leaves in A. thaliana [79][80][81]. Consistent with the fact that miR396 acts as a negative regulator of GRF gene expression, overexpression of miR396 negatively impacted cell proliferation in leaves and meristem size [78,82]. The miR396b-5p, Glyma15g19460.1 and Glyma17g35090.1 were down-regulated in two cultivars (Figures 1, 4, 5, S1 and S2). This deviation may be the result of a less dominant regulation by miRNAs and the involvement of other regulations, or of dependent of the balance between the transcription rate and the post-transcriptional regulation of the target when a certain stress stimulates the transcription of a miRNA and its target. However, the expression levels of Glyma15g19460.1 and Glyma17g35090.1 in ZH24 were lower than that in HX3 with or without Cd treatment (Tables S6 and S7). This interesting result may explain why the primary root of HX3 was longer than that of ZH24 (data not shown).
Other conserved gma-miR390a-5p, which identified silencing TAS3 (trans-acting siRNA) in rice or Medicago with Cd [25], Hg [40] and Al [39], performed over 2 fold down-regulation in ZH24 and less than 1.1 fold up-regulation in HX3 (Table 1). The downregulation of miR390 would lead the accumulation of intact TAS3 transcript and the decrease of tasiARFs and ARFs resulting in the inhibition of lateral root growth [40,83]. ARFs (auxin response factors) play critical roles in lateral root development [84]. This result might explain the growth inhibition in ZH24 was serious and sensitive to Cd exposure.
As one novel soybean Cd-responsive miRNA, miR1535b was identified to cleave Glyma07g38620.1. And Glyma07g38620.1 encoding isopentyl transferase (IPT), which catalyses the ratelimiting first step in de novo cytokinin (CK) biosynthesis and promotes the formation of isopentenyladenosine-59-monophosphate (iPa) [85,86]. CK has been known to inhibit primary root elongation in A. thaliana [87], and overexpression of ipt in leaves and roots can promote stress tolerance in Agrostis stolonifera [88]. In our study, Glyma07g38620.1 displayed an apparent up-regulation in HX3 and slight down-regulation in ZH24 under Cd exposure (Table S8), so does the CK, which may explain why HX3 possessed higher tolerance and distinctly primary root elongation inhibition than that of ZH24 under Cd stress. In addition, ARFs were up-regulated in HX3 and down-regulated in ZH24 due to the regulation of gma-miR390a-5p. CK has been known to inhibit primary root elongation and suggested to act as an auxin antagonist in the regulation of lateral formation. And the upregulation of ARFs might lead to the promotion of lateral root growth. But in what rats of ARFs and CK it regulated the growth of root and stress tolerance in soybean under Cd exposure need further study. Figure S1 The relative expression levels of 16 Cdresponsive miRNAs between HX3-CK and HX3-Cd in the root of soybean by qRT-PCR. The expression levels of miRNAs were normalized to the level of F-box. The french grey shading bar represents the relative expression level of miRNAs in HX3-CK. The dark grey shading bar represents the relative expression level of miRNAs in HX3-Cd. The results are averages 6 SD of the duplicates of three biological replicates. Significance of the changes between HX3-CK and HX3-Cd was checked with Student's t-test at the level of 0.01,P#0.05 (shown as ''*'') and P#0.01 (shown as ''**''). (TIF) Figure S2 The relative expression levels of 16 Cdresponsive miRNAs between ZH24-CK and ZH24-Cd in the root of soybean by qRT-PCR. The expression levels of miRNAs were normalized to the level of F-box. The french grey shading bar represents the relative expression level of miRNAs in ZH24-CK. The dark grey shading bar represents the relative expression level of miRNAs in ZH24-Cd. The results are averages 6 SD of the duplicates of three biological replicates. Significance of the changes between ZH24-CK and ZH24-Cd was checked with Student's t-test at the level of 0.01,P#0.05 (shown as ''*'') and P#0.01 (shown as ''**'').