Identification of drought-responsive and novel Populus trichocarpamicroRNAs by high-throughput sequencing and their targets using degradome analysis

Illumina sequencing of P. trichocarpaleaves under control and drought conditions

The size distribution of all unique sRNAs is summarized in Figure 1A. The displayed length of P. trichocarpa sRNA ranged from 16 to 27 nt, and the two major size classes were 24 nt (41.06% in CL and 43.94% in DL) and 21 nt (16.64% in CL and 16.05% in DL). This is in agreement with previous studies on sRNAs of P. trichocarpa [26] and M. truncatula [17] using high-throughput sequencing. To analyze the average abundance of each length between sRNAs of CL and DL further, we measured the ratio of raw and unique reads (Figure 1B). The redundancies of sRNAs varied widely in length, and the 20 and 21 nt sRNAs displayed the highest redundancies. The average ratio of redundant and unique sequences of sRNAs of the two libraries showed obvious changes in 21 nt sRNAs; the redundancy of DL was 37.16% greater than that of CL. This may be why drought stress strongly induced the expression of these 21 nt sRNAs; most conserved miRNAs belong to this group.

After genomic annotation of the P. trichocarpa sRNAs, small interfering RNA (siRNA) and miRNA with various important post-transcription regulating functions were the largest of our acquired sequences. The siRNA is a 22 to 24 nt double-strand RNA, each strand of which is 2 nt longer than the other on the 3’ end [31]. These aligned sequences might represent siRNA candidates. In total, deep sequencing obtained 577,393 and 956,979 siRNA candidates after the control and drought stress treatments, respectively (Additional file 1: file S1). Interestingly, the ratio of siRNA reads to all sRNAs reads increased sharply from 2.20% (CL) to 3.17% (DL). To obtain the annotation of known miRNAs, sRNAs were aligned to the miRBase 18.0 of P. trichocarpa. In total, 10,784,410 and 15,674,365 sequencing reads were identified as known poplar miRNAs in the two libraries. Thirty-four families from 207 known miRNAs were found, which accounted for about 87.3% of the total members. The remaining 30 miRNAs were not detected (Additional file 2: S2), possibly because of the tissue specificity of expression in poplar.

Novel non-conserved miRNAs from P. trichocarpa

These RNA structures were predicted by MFOLD software and manually checked according to the criteria of Meyers [32]. The lowest minimum free energy (MFE) of all hairpin structures of the novel miRNAs precursors was −31.9 kcal/mol (Table 2), which is slightly lower than the threshold of −30 kcal/mol reported in a previous study [35]. All precursors of novel miRNAs had regular hairpin structures (Additional file 4: S4), and four of these (Ptc-miRn5, Ptc-miRn11, Ptc-miRn38, and Ptc-miRn50) are presented in Figure 2.

Differential expression of miRNAs in P. trichocarpa

To identify drought-responsive miRNAs from P. trichocarpa, the number of normalized miRNA reads of CL and DL were compared. Based on the sequencing results, the differential expression of miRNAs greater than two-fold were chosen for experimental validation by quantitative real time polymerase chain reaction (RT-qPCR) (Additional file 5: S5). As shown in Figure 3, the expression patterns of the sequencing and RT-qPCR results of drought-responsive miRNAs were consistent, both indicating that four miRNAs (Ptc-miR159a-c, Ptc-miR472a, Ptc-miR472b, and Ptc-miR473a) were upregulated after drought treatment, and that five miRNAs (Ptc-miR160a-d, Ptc-miR164a-e, Ptc-miR394a/b-5p, Ptc-miR408, and Ptc-miR1444b-c) were downregulated by drought stress [25].

We further analyzed the expressions of the 65 new miRNAs under the two treatments. The drought-responsive miRNAs are listed in Figure 3; all were confirmed by the sequencing and RT-qPCR results. Among the 65 miRNAs, two novel miRNAs (Ptc-miRn6a-d and Ptc-miRn16) were downregulated by drought stress, and only miRn5 was upregulated in response to drought stress (Additional file 5: S5).

Target analysis of novel and conserved miRNAs by degradome sequencing

The previously known miRNA targets also identified in this study are available on the PopGenIE site (http://bioinformatics.cau.edu.cn/PMRD/adjunct/ptc_miR_target.txt). For new miRNAs whose targets were not known, we predicted their targets using the plant target prediction pipeline by the P. trichocarpa genome V2.0. The rules used for target prediction were based on those suggested by Allen et al. (2005) [36] and Schwab et al. (2005), as follows: (i) no more than four mismatches between the sRNA and the target (G-U bases count as 0.5 mismatches); (ii) no more than two adjacent mismatches in the miRNA/target duplex; (iii) no adjacent mismatches in positions 2–12 of the miRNA/target duplex (5’ of miRNA); (iv) no mismatches in positions 10–11 of the miRNA/target duplex; (v) no more than 2.5 mismatches in positions 1–12 of the miRNA/target duplex (5’ of miRNA); and (vi) the minimum free energy (MFE) of the miRNA/target duplex should be equal or greater than 74% of the MFE of the miRNA bound to its perfect complement [37]. We predicted 281 targets for 53 miRNA families; the other six were not found (Additional file 6: S6).

High-throughput sequencing of Populus

In a comparison of six Populus miRNA studies (Table 1) [11, 15, 25, 26, 29, 30], two used traditional Sanger sequencing [25, 29], two others used 454-pyrosequencing [26, 30], and the remaining two used the latest Illumina sequencing technology (as in the present study) [11, 15]. Along with the rapid development of sequencing technology, CL and DL can result in more sequences and greater sequencing depths than those reported in previous publications, due to the high throughput of the Illumina sequencer. In our study, because of the in-depth search, a large number of novel non-conserved miRNAs were found. The P. trichocarpa genome of Version 2.0 was used in this study; the transcript assemblies of the P. trichocarpa genome Version 2.0 are more meticulous than those of Version 1.1. This can increase the likelihood of finding more new miRNAs in general and drought-induced novel miRNAs in particular.

Novel miRNAs

Compared to six previous studies of Populus plants [10, 11, 15, 18, 19, 26], we identified 28 novel miRNAs have been identified (Table 2). Eleven of these were found at least once. On comparing the miRNA counts, 24 had counts greater than 100. Interestingly, two of the members of the Ptc-miRn54 family are the most frequently and robustly miRNAs identified in poplar high-throughput sequencing studies. Furthermore, the counterparts of Ptc-miRn40, Ptc-miRn52, Ptc-miRn54a, and Ptc-miRn54b in P. beijingensis were verified by RT-qPCR [11]. This provides more, strong evidence for the novel miRNAs identified from P. trichocarpa.

Drought-responsive miRNAs in P. trichocarpa

We further studied the target genes of these drought-responsive miRNAs by sequencing of the degradome library and comparing our work to previous studies [25, 29]. We found two upregulated miRNAs (Ptc-miR472 and Ptc-miRn5) that were both predicted to target putative disease resistance proteins in P. trichocarpa (Additional file 5: S5) [25]. The cross adaptation between disease resistance and drought stress tolerance in plants exists through unknown mechanisms. Ptc-miR159 is another upregulated miRNA; its Arabidopsis homolog targets an MYB transcription factor. The ABA-induced accumulation of the miR159 homolog makes the MYB transcript degradation desensitize hormone signaling during seedling stress responses in Arabidopsis [40]. According to our degradome sequencing results, the Ptc-miR159 was confirmed to target a methionine sulfoxide reductase (MSR). The homologs of MSR were induced by biotic and abiotic stresses in plants [47–50]. They catalyze the reduction of methionine sulfoxide to methionine [47] and play a major role in regulating the accumulation of reactive oxygen species (ROS), which can damage proteins in plant cells [50]. Regulation of the MSR gene by Ptc-miR159 may occur through a homeostatic mechanism in response to drought stress in P. trichocarpa.

Ptc-miR473 was also upregulated in drought stress. It targets a member of a plant-specific GRAS transcription factor gene family [29]. Another member of this family (PeSCL7) from P. euphratica was confirmed to play key roles in salt and drought stress tolerance [51]. In the present study, Ptc-miR473 was confirmed to be targeted to Vein Patterning 1 (VEP1), which belongs to a short-chain dehydrogenase/reductase (SDR) superfamily [52]. The homolog of VEP1 in Arabidopsis was confirmed to be required for vascular strand development and to be upregulated by osmotic stress [52, 53]. Ptc-miR473 regulates the expression the GRAS protein and VEP1, both of which were responsive to drought stress, this may be the drought tolerance mechanism in P. trichocarpa.

The number of downregulated miRNAs was larger than the number of upregulated miRNAs. The two downregulated miRNAs (miR160 and miR164) were both identified to be cold-responsive miRNAs in P. trichocarpa [25]. TMV-Cg virus infection in Arabidopsis causes the accumulation of miR160 and miR164 [54]. Three auxin responsive factor (ARF) genes (ARF10, ARF16, and ARF17) are the targets of miR160 [55]. Repression of ARF10 by miR160 is critical for the seed germination and post-germination stages [56]. MiR164 has been predicted targete six NAC-domain proteins (PNAC041, PNAC042, PNAC151, PNAC152, PNAC154, and PNAC155) from subfamily NAC-a [57], and NAC-domain proteins have been confirmed to be important in drought stress tolerance [58, 59]. These mechanisms may also be at work in drought-stress tolerance in P. trichocarpa for these two miRNAs.

Two downregulated miRNAs (Ptc-miR408 and Ptc-miR1444) have been reported to be Cu-responsive miRNAs in P. trichocarpa. Their targets include miR408-targeted plastocyanin-like proteins and miR1444-targeted all plastid polyphenol oxidases [60, 61]. Drought treatment may increase the relative concentration of Cu ion in the cytoplasm. When the Cu supply is sufficient, it is envisaged that the conjunction between mature miRNAs and their precursors will be suppressed, leading to the upregulation of miRNA-targeted Cu proteins [60]. Accordingly, the balance of Cu ion contributes to the healthy growth and development of poplars during stress. In P. trichocarpa, Ptc-miR1444a is reportedly downregulated by dehydration [25], and Ptc-1444b/c was also found to be downregulated by drought in this study. MiR408 is reportedly downregulated by drought stress in rice [13] and has been experimentally identified to target an early responsive dehydration-related (ERD) protein in P. trichocarpa. Drought stress might induce the expression of ERD protein by downregulating the expression of miR408 in P. trichocarpa. This may be one of the mechanisms of regulation of drought-stress tolerance [25].

Other downregulated miRNA is Ptc-miR394, whose predicted targets are annotated as dehydration-responsive protein (POPTR_0002s07760.1) and F-box proteins (POPTR_0001s13770.1 and POPTR_0003s16980.1), which were recently reported to be differentially regulated by stress conditions and to play significant roles in the abiotic stress-response pathway. In Arabidopsis, salt-induced miR394 targets the mRNA of F-box proteins [12, 56].

From the analysis of predicted targets to downregulated Ptc-miRn6, a CCCH-type zinc finger protein and two trichome birefringence-like proteins (TBLPs) were functionally predicted. Although a cotton CCCH-type zinc finger protein has been identified to enhance abiotic stress tolerance in tobacco [62], we did not find any additional possible regulatory mechanisms between CCCH-type zinc finger protein and drought tolerance in P. trichocarpa. The homolog of TBLP in Arabidopsis is important to the formation of crystalline cellulose in trichomes [63]. As previous studies have reported, trichome density increases with water shortage [64], and the thick trichome layer could prevent water loss [65]. This may be the mechanism by which miRn6 regulates the expression of TBLP to adapt to drought stress.

Degradome analysis of non-drought-responsive miRNAs

In Arabidopsis, miR390 was reported to target TAS genes [66], while in P. trichocarpa, no TAS homologs have been found [26]. From our study, the degradome sequencing data proved the adjustment mechanism of Ptc-miR390 and lipoxygenases (LOXs). The activity of LOX protein can partially reduce the production of radicals and ROS [67]. This may explain the regulatory mechanism of miR390 in poplars. Four UDPGs were found to be targeted by Ptc-miR482, and all were classified as category I. The UDP-glucosyltransferases (UDPGs) are enzymes that attach a sugar molecule to a specific acceptor in plants [68]. As in Arabidopsis, the UDPG is a key regulator of stress adaption through auxin IBA [69] and plays a role in fine-tuning nitrogen assimilation in cassava [70]. This is a novel mechanism by which miR482 regulates the UDPG gene family in P. trichocarpa.

The degradome sequencing results imply that the miRNAs with no detected targets may silence genes by repressing translation. However, we could not obtain information about translation repression by miRNA through degradome sequencing. Only 19 targets of new miRNAs were identified. The targets of these non-conserved miRNAs are difficult to detect, possibly because of low abundance or a spatial expression pattern. More studies are needed to shed light on the regulation network of these miRNAs in P. trichocarpa. Over-expressing or repressing expression of these miRNAs in P. trichocarpa may help to elucidate the regulation mechanism.

Plant materials and total RNA extraction

P. trichocarpa seedlings of the same size (~5 cm) from tissue culture were planted in individual pots (15 L) containing loam soil and placed in a greenhouse at Beijing Forestry University. They were well irrigated and grown under control conditions (25°C day/20°C night, 16-h photoperiod) for three months, the heights of them were about 45 cm. During the period of drought-stress treatment, P. trichocarpa seedlings were sustained at two RSMC levels (70–75% and 15–20%) for 1 month according to a previous publication [28]. The mature leaves were used as drought materials. Mature leaves from soil with sufficient irrigation (RSMC at 70–75%) were used as a control, and a relatively modest dehydration level (RSMC at 15–20%) was chosen for the drought treatment. Each treatment contained three repeat individuals. Leaf water potential (WP) was measured by PsyPro WP data logger (Wescor) (Additional file 8: S9). Photosynthetic rate, water conductance, intercellular CO2 concentration, and transpiration rate were measured by Li-6400 Photosynthesis System (Li-Cor) (Additional file 9: S10). For material harvest, mature leaves from the same position of different individual plants were collected and frozen immediately in liquid nitrogen for RNA extraction. The total RNA was extracted by the standard CTAB method for plants [71]. Then they were used for sequencing and RT-qPCR.

High-throughput sequencing and bioinformatics analysis

Illumina sequencing on sRNAs (ranged from 18 nt to 30 nt) was conducted using an Illumina Genome Analyzer, following the Illumina protocol [72]. After removing contaminants, low-quality sequences, and <18 nt sequences, clean reads were obtained and aligned against the P. trichocarpa genome (version 2.0) using SOAP software [73]. tRNA, rRNA, snRNA, snoRNA, and some other repeat sequences were removed from the sequences with a perfect match to the genome through a search of the NCBI Genbank and Rfam databases [74]. The remaining unique sequences were divided into known miRNAs and candidate miRNAs by alignment with the miRbase 18.0 [75]. The candidate miRNAs were further analyzed by MFOLD software on the RNA secondary structure of the miRNA::miRNA* and pre-miRNA hairpin energy [76]. Parameters were set to meet the criteria of plants [32].

Differentiatial expression analysis of miRNAs between the two treatments

The sequence reads of the two libraries were normalized to 1 million by the total number of sRNA reads in each sample. The calculation of the p-value for comparison of the miRNA expression between the two libraries was based on previously established methods [77, 78]. Specifically, the log₂ ratio formula was: log₂ ratio = log₂ (miRNA reads in drought treatment/miRNA reads in control).

where N₁ is the total number of reads in the sequencing library of the control, N₂ is the total number of reads in the sequencing library of the drought treatment, x is the number of reads for an miRNA in the control library, and y is the number of reads for an miRNA in the drought treatment library.

All calculations were performed on a BGI Bio-Cloud Computing platform (http://www.genomics.cn/en/navigation/show_navigation?nid=4143). Normalized miRNAs of <1 were filtered in both libraries.

RT-qPCR of mature miRNAs

To validate the results of miRNAs from high-throughput sequencing, RT-qPCR was performed. The RNAs were extracted from leaves using the CTAB method [71]. A poly (A) was added to the 3’ end, and reverse transcription was begun. In particular, a known sequence at the 5’ end of the oligo-dT primer was designed to be a communal reverse primer of the RT-qPCR. The One Step Prime-Script miRNA cDNA Synthesis Kit and SYBR Premix ExTag II (TaKaRa) were used. All primers used in this study are listed in Additional file 10: S8. The 5.8S ribosomal RNA was used as the internal control [25]. RT-qPCR was performed using an ABI StepOnePlus instrument. Calculation of RT-qPCR results were revised as follow: Sample cycle threshold (Ct) values were determined and then standardized based on the 5.8S gene control primer reaction, and the 2^-ΔΔCT method was applied to calculate the relative changes in gene expression from RT-qPCR experiments [79].

Target prediction and confirmation by degradome sequencing

New P. trichocarpa miRNA targets were predicted as described before [36, 80–82]. During the prediction, a penalty score (alignment score) criterion was induced according to the alignment between the miRNA and its potential target. Our cut-off values in both prediction and degradome sequencing data analysis were also set to <2.5 as used in previous studies on poplar miRNA target prediction. The biological function of the predicted targets was retrieved from the Universal Protein Resource (http://www.uniprot.org).

Degradome sequencing following the PARE protocol was used [38]. Only miRNA-cleaved mRNA and other degraded mRNA could be ligated by a 5’ RNA adapter because the 5’-phosphate and intact mRNAs were protected by the 5’ cap. First, adapters and low-quality nucleotide reads were removed from raw reads using the Fastx-Toolkit. Then the clean reads were further analyzed by Cleaveland 2.0 software [83]. Briefly, the reads were first mapped to the P. trichocarpa transcripts database from JGI Phytozome 2.0. At this step, a target plot was also created to distinguish the true miRNA cleavage site from background noise. We ran Cleaveland 2.0 with default parameters using 100 randomized sequencing shuffles. The NCBI database was used to predict functions of targets that were not annotated in JGI Phytozome 2.0. The cleaved target transcripts were categorized into three categories according to the following criteria: I, the abundance of reads in its cleavage site is the maximum on the transcript; II, the abundance of reads in its cleavage site is not the maximum, but is equal to or higher than the median for the transcript; and III, the abundance of reads in its cleavage site is less than the median for the transcript.