High throughput deep degradome sequencing reveals microRNAs and their targets in response to drought stress in mulberry (Morus alba)

Abstract

MicroRNAs (miRNAs) play important regulatory roles by targeting mRNAs for cleavage or translational repression. Identification of miRNA targets is essential to better understanding the roles of miRNAs. miRNA targets have not been well characterized in mulberry (Morus alba). To anatomize miRNA guided gene regulation under drought stress, transcriptome-wide high throughput degradome sequencing was used in this study to directly detect drought stress responsive miRNA targets in mulberry. A drought library (DL) and a contrast library (CL) were constructed to capture the cleaved mRNAs for sequencing. In CL, 409 target genes of 30 conserved miRNA families and 990 target genes of 199 novel miRNAs were identified. In DL, 373 target genes of 30 conserved miRNA families and 950 target genes of 195 novel miRNAs were identified. Of the conserved miRNA families in DL, mno-miR156, mno-miR172, and mno-miR396 had the highest number of targets with 54, 52 and 41 transcripts, respectively, indicating that these three miRNA families and their target genes might play important functions in response to drought stress in mulberry. Additionally, we found that many of the target genes were transcription factors. By analyzing the miRNA-target molecular network, we found that the DL independent networks consisted of 838 miRNA-mRNA pairs (63.34%). The expression patterns of 11 target genes and 12 correspondent miRNAs were detected using qRT-PCR. Six miRNA targets were further verified by RNA ligase-mediated 5’ rapid amplification of cDNA ends (RLM-5’ RACE). Gene Ontology (GO) annotations and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis revealed that these target transcripts were implicated in a broad range of biological processes and various metabolic pathways. This is the first study to comprehensively characterize target genes and their associated miRNAs in response to drought stress by degradome sequencing in mulberry. This study provides a framework for understanding the molecular mechanisms of drought resistance in mulberry.

About the Authors:

Ruixue Li

Affiliations School of Biology and Technology, Jiangsu University of Science and Technology, Zhenjiang, Jiangsu, P. R. China, The Sericultural Research Institute, Anhui Academy of Agricultural Sciences, Hefei, Anhui, P. R. China

Dandan Chen

Affiliation: School of Biology and Technology, Jiangsu University of Science and Technology, Zhenjiang, Jiangsu, P. R. China

Taichu Wang

Affiliation: The Sericultural Research Institute, Anhui Academy of Agricultural Sciences, Hefei, Anhui, P. R. China

Yizhen Wan

Affiliation: School of Biology and Technology, Jiangsu University of Science and Technology, Zhenjiang, Jiangsu, P. R. China

Rongfang Li

Affiliation: School of Biology and Technology, Jiangsu University of Science and Technology, Zhenjiang, Jiangsu, P. R. China

Rongjun Fang

Affiliation: School of Biology and Technology, Jiangsu University of Science and Technology, Zhenjiang, Jiangsu, P. R. China

Yuting Wang

Affiliation: The Sericultural Research Institute, Anhui Academy of Agricultural Sciences, Hefei, Anhui, P. R. China

Fei Hu

Affiliation: The Plant Protection and Agro-products Safety Institute, Anhui Academy of Agricultural Sciences, Hefei, Anhui, P. R. China

Hong Zhou

Affiliation: School of Biology and Technology, Jiangsu University of Science and Technology, Zhenjiang, Jiangsu, P. R. China

Long Li

Affiliation: School of Biology and Technology, Jiangsu University of Science and Technology, Zhenjiang, Jiangsu, P. R. China

Weiguo Zhao

* E-mail: [email protected]

Affiliation: School of Biology and Technology, Jiangsu University of Science and Technology, Zhenjiang, Jiangsu, P. R. China

Introduction

MicroRNAs (miRNAs) are endogenous non-coding small RNAs containing 21–24 nucleotides. miRNAs regulate gene expression accurately and effectively at the post-transcriptional level by repressing translation or directly degrading target mRNAs [1,2]. In plant cells, primary miRNA (pri-miRNA) with a cap and a poly (A) tail is processed into precursors of miRNA (pre-miRNA) containing the distinctive stem-loop structure [3], and thereafter the pre-miRNA is cleaved into a miRNA/miRNA* duplex [4]. Then the miRNA/miRNA* duplex is methylated by Hua Enhancer 1 (HEN1) and transported into the cytoplasm [5]. Mature plant miRNAs are loaded into the ARGONAUTE 1 (AGO1) complex to suppress translation or to cleave the target transcripts [6]. Plant miRNAs have been largely implicated in degradation of their RNA targets by slicing precisely between the 10th and 11th nucleotides (nt) from the 5’ end of miRNAs. Recently, increasing evidence has indicated that miRNAs participate in various processes, such as plant growth, physiological and biochemical processes, signal transduction, cell apoptosis, and biotic and abiotic stress responses [7–11].

To date, hundreds of small RNAs have been isolated by direct cloning and by deep sequencing in higher plants, such as Arabidopsis thaliana, rice, soybean, chrysanthemum, and so on [12–14]. As negative regulators of gene expression, plant miRNAs can respond to biotic and abiotic stress [15–17]. By establishing small RNA libraries of A. thaliana under high salt, drought, ABA and cold treatment, 30 miRNAs associated with abiotic stress were identified [18]. More than 40 plant miRNA family genes are now associated with high salt and drought stress [19,20]. Upon drought treatment, expression of miR159, miR160, miR166, miR169, miR172, miR395, miR396, miR408, miR472, miR477, miR482, miR1858, miR2118, and miR5049 were found to be significantly differentiated in bread wheat [21]. Xie et al found that a series of miRNAs are associated with these top-ranked genes that combat drought and salinity stress in cotton, including miR164, miR172, miR396, miR1520, miR6158, -n24, -n56, and -n59 [22]. Salinity and drought stress induces elevated expression of miR319 in creeping bentgrass, resulting in downregulation of at least four putative target genes of miR319 (AsPCF5, AsPCF6, AsPCF8, and AsTCP14) as well as a homolog of the rice NAC domain gene AsNAC60, and therefore contributes to plant abiotic stress response [23]. Eldem et al found that the expression level of 262 (104 up-regulated, 158 down-regulated) of the 453 miRNAs changed significantly in leaf tissue, whereas 368 (221 up-regulated, 147 down-regulated) of the 465 miRNAs had expression levels that changed significantly in root tissue upon drought stress. The expression level of miR159, miR169, miR393, miR397, miR398 and miR395 had different changes between root and leaf in response to drought [24].

However, few studies have identified miRNAs involved in response to abiotic stress in mulberry [3]. Due to the high degree homology of sequence and function of miRNAs between the different species [25], miRNAs are likely to represent a more primitive physiological regulatory mechanism [26]. Illuminating the function of these tiny molecule nucleotides requires efficient approaches to identify their target genes. Originally, plant miRNA targets were studied via computational prediction, which is based on the perfect sequence complementarity between a miRNA and the target mRNA or sequence conservation among different species [27]. However, targets prediction often has a high level of mismatch in miRNA-mRNA pairs and every single predicted gene must be verified independently. This one-at-a-time isolation of target cleavage remnants is laborious, time-consuming and costly [28,29]. To overcome this limitation, a new method called degradome sequencing has been successfully established to identify small RNA target transcripts at a global scale [30–32]. The method includes deep sequencing, bioinformatic analysis and 5’-rapid amplification of cDNA ends (5’ RACE). It has been used for global identification of miRNA-target RNA pairs in numerous plant species, such as A. thaliana, Oryza sativa, Glycine max and Zea mays [29,30,33,34].

Drought is one of the major natural disasters affecting plants, and dry weather is becoming more frequent with the intensification of the greenhouse effect. This results in huge losses in industrial and agricultural production, and serious damage to the ecological environment.

Mulberry (Morus alba) is an important perennial economic tree, and has a broad ecological distribution in China. It is used in sericulture and has important economic and ecological value. It is highly adaptable to poor environments and can endure salt, drought, cold, waterlogging, and heavy metal ions [35,36]. Certain miRNAs (such as miR159, miR169 and miR398) are known to play an important role in response to drought stress in A. thaliana [12,16,30]. However, these miRNAs and their molecular roles have not been verified in mulberry. This study was conducted to: (1) comprehensively investigate miRNA targets in response to drought stress in mulberry; (2) analyze the regulatory relationship between miRNAs and their target genes in mulberry through the analysis of miRNA-target molecular network of DL and CL; (3) further verify six of the miRNA targets by RNA ligase-mediated 5’ rapid amplification of cDNA ends (RLM-5’ RACE); and (4) characterize differential expression patterns of 12 drought-responsive miRNAs and target mRNAs in mulberry by stem-loop qRT-PCR. The findings open a new pathway to greatly understand the molecular basis in response to drought stress in mulberry.

Materials and methods

Plant materials and drought stress treatment

The mulberry species (Morus alba), ‘Yu711’, was obtained from the National Mulberry Gene Bank in Zhenjiang, Jiangsu, China. The plants were grown in a greenhouse under a photo period of 14 h light/10 h dark at 25°C day/20°C night. Mulberry cuttings were grafted to rootstocks. The grafted nurseries were planted in 35 cm diameter pots containing loam soil with a nursery per pot. The grafted plants were randomly grouped when new shoots had grown to 20 cm in length. One group was used as the control, and the other one was treated with drought stress for 5 d, 10 d, and 15 d (Fig 1). Each grouping contained three replicates. Young leaves were collected from the same position of the plants and immediately frozen in liquid nitrogen and stored at -80°C for RNA extraction.

[Figure omitted. See PDF.]

Fig 1. Effects of drought stress on phenotypic alterations in mulberry sapling.

(A) After drought treatment for 15 days. (B) Control mulberry sapling.

https://doi.org/10.1371/journal.pone.0172883.g001

Construction of the degradome library and sequencing

All centrifuge tubes and pipette-tips were RNase-free or were treated with DEPC. All buffers were RNase free or prepared using DEPC-treated ddH₂O. Total RNA of each sample was extracted by Trizol (Life Technology, USA) according to the manufacturer’s instructions. For the same grouping of treatments, equal amounts of the RNA from each time point were mixed together. The RNA degradome library was constructed as previously described [5,37]. In brief, poly (A) RNA was isolated from approximately 200 μg of total RNA using the Oligotex mRNA mini kit (Qiagen, USA). A 5’-RNA adapter containing a Mme I recognition site in its 3’ terminus was ligated to the poly (A) RNA possessing a free 5’-monophosphate by T4 RNA ligase (Takara, China). The ligation products were purified using the Oligotex mRNA mini kit and reverse transcribed using the Oligo (dT)₁₈ primer and Superscript II reverse transcriptase (Invitrogen, USA). The first-strand cDNA was amplified for five cycles using Ex Taq DNA Polymerase (Takara, China) and the PCR products were digested with Mme I. Next, the digested products were ligated to a 3’-double-stranded DNA adaptor using T4 DNA ligase (Takara, China), and amplified by PCR for 20 cycles. The final PCR products were gel purified and subjected to SE51 sequencing using Illumina HiSeq2000 (Illumina Inc., San Diego, USA).

Bioinformatics analysis and identification of miRNA targets

Two mulberry degradome libraries were sequenced using Illumina HiSeq 2000. Low quality nucleotide reads and clip adapter sequences were removed from the raw data using the Fastx-Toolkit. High quality 20–21 nucleotide-long reads were subjected to the CleaveLand pipeline for small RNA targets identification as previously described [38]. Briefly, the raw tags were first normalized to ‘tags per 100 million’ (TP100M). After preprocessing, the non-coding RNAs were removed and clean tags were generated, saved, and classified by the alignment to the database and remove the ncRNAs (no coding RNAs). Clean tags were mapped to mulberry complete reference genome (http://morus.swu.edu.cn/morusdb/). The sequences of rRNAs, tRNAs, snoRNAs and snRNAs were retrieved from the Rfam database (http://www.sanger.ac.uk/Software/Rfam/). Degraded sequences in the clean reads were identified as polyN when over 70% of the sequence was a single base. Next, distinct tags that perfectly matched mulberry cDNA or mRNA sequences by SOAP2.20 based on Genbank of NCBI (http://www.ncbi.nlm.nih.gov/), the known miRNAs by PAREsnip from miRBase (miRBase21.0: http://www.mirbase.org/), or miRNAs identified in other studies [3,39] were used to analyze miRNA-mRNA pairs. Alignments were then scored according to a previously described scheme developed for plant miRNA/target pairings [40]. All alignments with scores up to seven and no mismatches at the cleavage site (between the 10th and 11th nucleotide) were considered as candidate targets. CL and DL data were analyzed separately. The identified targets were grouped into five categories based on the relative abundance of the degradome tags mapping at the miRNA target site through the height of the degradome peak at each occupied transcript position. Category 0: >1 raw tags at the position, abundance at the position was equal to the maximum on the transcript, and there was only one maximum on the transcript; category 1: >1 raw tags at the position, abundance at the position was equal to the maximum on the transcript, and there was more than one maximum position on the transcript; category 2: >1 raw tags at the position, abundance at the position was less than maximum but higher than the median for the transcript; category 3: >1 raw tags at the position, abundance at the position was equal to or less than the median for the transcript; category 4: only 1 raw tag at the position.

Construction of miRNA-targets network

Identified targets and related miRNAs based on the analysis of degradome sequencing in DL and CL the profiles were used for building and visualizing miRNA-mRNA interaction global network using Cytoscape 3.2.1 software as previously described [41,42,43]. Subsequently, in order to further understand the regulatory relationships between drought responsed miRNAs and their target genes, we further filtered the target genes with P ≤ 0.05 in DL. Using these filtered interaction pairs, the specific network of miRNA-mRNA which P ≤ 0.05 and homologous to drought-related miRNA targets in DL were also constructed. In the network, nodes represent miRNAs or target genes. The network structure is formed of basic elements (target genes and miRNAs; named nodes) and the connections representing miRNA-target interactions (named edges). If two genes were annotated to be related, an edge was added between them in the network.

Detection of the expression profiles of miRNAs and their targets by qRT-PCR

To determine the regulatory relationship between miRNAs and their targets, the expression profile of 12 drought responsive miRNAs and their targets were examined by stem-loop qRT-PCR as previously described [44,45]. Total RNA was isolated from the leaves of the contrast and treated plants at 5 d, 10 d, and 15 d under drought stress and with RNAiso plus (Takara, China) according to the manufacturer’s instructions. The first cDNA strand was synthesized from total RNA using the M-MLV-reverse transcriptase (RTase) (Takara, China) and miRNA specific stem-loop primers were designed according to the method described by Chen et al. [44]. Briefly, six nucleotides that paired with the 3’ end of the miRNA were linked to a stem-looped sequence (GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCACTGGATACGAC) to synthesize the stem-loop reverse transcription primer. The reactions contained 1 μg of total RNA, and each reaction was performed with 1 μ l20 μM gene-specific stem-loop primers. The RNA and primers were mixed with DEPC-treated water up to 6 μl, then After, incubated at 70°C for 10 min and immediately chilled on ice for more than 2 min. Then, 2.5 μl 5 × RT-Buffer, 0.5 μl 10 mM dNTP, 0.5 μl RNase Inhibitor and 0.5 μl M-MLV- RTase were added together to a total volume of 10 μl. Synthesis was performed at 42°C for 60 min and inactivation of the enzyme was performed at 72°C for 15 min. The product was then diluted 3-fold and 1 μl cDNA was used as the template to perform the qRT-PCR with each miRNA specific forward primer and universal primer. For target genes, 1 μg RNA and 1 μl 50 μM Oligo(dT)₁₈ primer were mixed with DEPC-treated water to 17 μl, then incubated at 70°C for 10 min and after, ice-cooled immediately for more than 2 min. Then, 5 μl 5× RT-Buffer, 1.25 μl 10 mM dNTP, 0.75 μl RNase Inhibitor and 1 μl M-MLV RTase were added. The next steps were identical to the reverse transcription of miRNA. The product was then diluted 4-fold and 1 μl cDNA was used as the template to perform the qRT-PCR with each target gene primer. The primers for miRNAs, targets and β-actin (selected as a reference gene for normalization) were mixed with the FastStart Universal SYBR Green Master Mix kit (Roche, USA) based on the handbook, and 20 μl of the reaction mix was added to each well. Reactions were performed in LightCycler^® 96 Real-Time PCR System (Roche, USA) with thermal cycling parameters at 95°C for 600 s followed by 45 cycles of 95°C for 10 s, 59°C for 10 s, and 70°C for 10 s. During amplification, melting curves were constructed. The sequences of stem-loop reverse transcriptase primers, miRNA-specific PCR primers and target-specific PCR primers were listed in S1 and S2 Tables, respectively. All reactions were assayed in three biological replicates with three technical replicates. The relative expression differences of miRNAs and their targets were calculated by 2^-ΔΔCt method. Standard errors and standard deviations were calculated from replicates and significance was measured through one-way ANOVA Duncan’s multiple range test at the level of 0.01 < P ≤ 0.05 and P ≤ 0.01.

Validation of targets by RLM-5’ RACE

In order to validate the cleavage sites of miRNA to target genes, a modified RLM-5’ RACE was performed using the FirstChoice RLM-RACE Kit (Ambion, USA) without calf intestine alkaline phosphatase (CIAP) and tobacco acid pyrophosphatase (TAP) treatments for six target genes in the drought stressed degradome library as previously described [46,47]. Total RNA was extracted from drought stressed mulberry leaves with RNAiso plus (Takara, China) as described by the manufacturer. Then, approximately 2 μg RNA was ligated with 5’ RACE adaptors (5’-GCUGAUGGCGAUGAAUGAACACUGCGUUUGCUGGCUUUGAUGAAA-3’) using T4 RNA ligase. The ligated mRNAs were then reverse-transcribed using random decamers via M-MLV RTase following the manufacturer’s instructions. Two rounds of 5’ RACE reactions were performed and the initial PCR was carried out using the RT reaction product, the 5’ RACE outer primer (5’-GCTGATGGCGATGAATGAACACTG-3’), and the gene-specific outer primer, Nested PCR was carried out using 1/10 of the initial PCR reaction product, the 5’ RACE inner primer (5’-CGCGGATCCGAACACTGCGTTTGCTGGCTTTGATG-3’) and gene-specific inner primer. The final PCR products were gel purified, cloned into the pMD18-T Vector (Takara, China), and sequenced. Each target was confirmed by at least four clones. Gene-specific PCR primers were designed from between 150 and 600 nucleotides from the 5’ end of the predicted target site and were listed in S3 Table.

Functional annotations of the miRNA targets

To investigate the putative biological functions of target genes and biological processes possibly regulated by miRNAs in mulberry, Gene Ontology (GO) annotations, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis (http://www.kegg.jp/kegg/pathway.html), and the Non-Redundant (NR) Protein Database [48–51] were employed to annotate and classify target genes using the DAVID gene annotation tool. For enrichment analysis, a hypergeometric distribution based statistical test (level of significance at 0.05) was used to reject the chances of randomness in the miRNA’s associations to target genes with their corresponding ontology term. The GO categorization results were listed as three independent hierarchies for biological process, cellular component, and molecular function. GO classification of the drought responsive miRNAs targets and the metabolic pathways of these genes were obtained based on KEGG analysis.

Results

Mulberry degradome library construction, sequencing, and sequence analysis

In order to identify the miRNA targets in mulberry plants at a global level, two degradome libraries (drought library and contrast library) which captured the cleaved mRNAs were constructed for sequencing by Illumina sequencer. Our data have been deposited in the NCBI’s Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo/) with the accession number GSE84889. A total of 17,193,823 raw tags for CL and 20,200,629 tags for DL were obtained. After trimming low quality tags, adapter tags and N tags, there were 17,188,651 (99.97%) clean tags in CL and 20,192,466 (99.96%) clean tags in DL. Of those, there were 26,818,165 total common clean tags in CL and DL. There were 3,772,569 and 6,790,383 specific clean tags in CL and DL, respectively (Table 1). Of the 17,188,651 clean tags, 948,369 were identified as fragments of rRNAs, 348 were identified as fragments of tRNAs, 398 were identified as fragments of small nuclear RNA (snRNAs), and 1026 were identified as fragments of rRNAs, tRNAs, small nuclear RNA (snRNAs) and small nucleolar RNAs (snoRNAs) in CL, respectively using a BLASTN search against the Rfam database. Of the 20,192,466 clean tags, 693,817 were identified as fragments of rRNAs, 546 were identified as fragments of tRNAs, 612 were identified as fragments of small nuclear RNAs (snRNAs), and 1516 were identified as fragments of small nucleolar RNAs (snoRNAs) in DL, respectively (Table 1). The priority rule, Rfam > Genbank > Poly N, was used to map every unique degradome to the corresponding annotation, and map unannotated tags to reference genes (cDNA), classification of clean tags in CL and DL were listed in Table 1. The chromosomes distribution of clean tags was obtained by mapping clean tags to the complete mulberry reference genome using SOAP2.20, as shown in Fig 2. Finally, 10,132,806 (58.95%) and 9,353,151 (46.32%) clean tags in CL and DL were mapped to the genome, respectively. These data indicated that our two degradome libraries were of high quality and recovered most of the degraded mRNA targets that contained the sequence profile of miRNA-mediated cleavage and allowed us to conduct further analysis.

[Figure omitted. See PDF.]

Table 1. Summary of mulberry degradome sequences by Illumina sequencing.

https://doi.org/10.1371/journal.pone.0172883.t001

[Figure omitted. See PDF.]

Fig 2. Chromosomes distribution of clean tags in CL and DL.

The x axis indicates the chromosomes. The y axis indicates the number of tags that locate on each chromosome: area above 0 shown on blue is the number of tags on the sense strand of chromosome, whereas area below 0 shown in red is the number of tags on the antisense strand of chromosome. (A) CL. (B) DL.

https://doi.org/10.1371/journal.pone.0172883.g002

Systematic identification of miRNA targets in mulberry

Confirmation of miRNA targets was a prerequisite to better understand the functional roles of miRNAs and may result in the discovery of new non-conserved miRNAs. In mulberry, conserved miRNA targets were previously investigated mainly by bioinformatics prediction and only a few miRNA targets have been experimentally validated [3,39]. High-throughput degradome sequencing was performed in this study to identify more miRNA targets in mulberry, particularly specific targets of drought-stress-responsive miRNAs. A characteristic feature of miRNA-guided cleaving was that the cleavage takes place precisely between the 10th and 11th nt from the 5’ of miRNA in the complementary region of the target transcript. Therefore, cleaved RNA targets should have distinct peaks of degradome tags at the predicted cleavage site relative to other regions of the transcript. Systematic identification of miRNA targets was accomplished by analyzing the 20 and 21 nt tags with the CleaveLand pipeline for miRNA target identification using the methods described above.

As listed in S4 and S5 Tables, 409 target genes of 30 conserved miRNA families and 990 target genes of 199 novel miRNAs were identified in CL. In DL, 373 target genes of 30 conserved miRNA families and 950 target genes of 195 novel miRNAs were identified. Based on the abundance of tags at the targets’ cleavage sites, the miRNA targets were categorized into five classes as described above. There were 10, 26, 522, 23 and 818 targets in categories 0, 1, 2, 3 and 4 respectively, in CL, and 26, 9, 708, 73 and 507 targets in categories 0, 1, 2, 3 and 4 respectively, in DL. Remarkably, a conserved miRNA family could target various numbers of genes ranging from 1 up to 54 and a target could be cleaved by 1–7 different conserved miRNA or miRNA family. Of those conserved miRNA families in DL, mno-miR156, mno-miR172 and mno-miR396 had the highest number of targets with 54, 52 and 41 unique transcripts, respectively (S5 Table), indicating that these three miRNA families may be in the center of the gene regulation networks. In contrast, fewer targets were identified for mno-miR2111, mno-miR4376, mno-miR4995, mno-miR827, mno-miR390 and mno-miR5523, indicating that these miRNAs might act in specialized pathways. The target gene XM_010104394.1, a squamosa promoter-binding-like protein, predicted to be involved in metabolic processes in response to stimulus and catalytic activity, was cleaved by six mno-miR156s and mno-miR535, indicating that this miRNA may play an important role in the process of growth, physiological metabolism and the corresponding stress of plants. Many of the target genes of the conserved miRNAs were classified as the transcription factors (TFs), such as NAC (a target of mno-miRn108-5p), ARF (a target of mno-miR160, mno-mi172, mno-miR156 et al.), NFY (a target of mno-miRn223-5p), MYB (a target of mno-miR159, mno-miR858, mno-miRn148-5p et al.), HD-ZIP (a target of mon-miR166), GRAS (a target of mno-miR171), and GRF (a target of mno-miR396), and so on. These TFs are known to regulate diverse aspects of plant growth, development, and biotic and abiotic stresses. Compared with A. thaliana and rice, many conserved miRNA targets were found to be conserved in mulberry, indicating that the miRNA-target relationship was evolutionarily conserved [27,30].

Of the 72 miRNA-mRNA pairs in CL and the 63 miRNA-mRNA pairs in DL with P ≤ 0.05, there were 29 and 27 target genes for 11 and 12 conserved miRNA families, respectively, and 36 and 27 target genes for 26 and 23 novel miRNAs, respectively (Tables 2 and 3). In CL, 72 miRNA-mRNA pairs were categorized into five classes as 6, 11, 14, 7 and 34 in category 0 to category 4, respectively. In miRNA-mRNA pairs of DL, 15 targets belonged to category 0, whereas the category 1–4 targets accounted for 6, 22, 8 and 12, respectively. The abundance of tags for each of them were presented in the form of target plots (T-plots) that plot the abundance of the tags relative to their position in the transcript [31]. For instance, six representative T-plots were shown, three from CL and another three from DL (Fig 3). In each of the six models, a clear peak for the absolute number of tags was found at the identified cleavage site for mno-miR166f, mno-miR408c, mno-miR396b, mno-miR319c, mno-miR156d and mno-miR828. Cleavage sites of all of the miRNAs were located in the coding DNA sequence (CDS) of the target genes identified by degradome sequencing in mulberry.

[Figure omitted. See PDF.]

Table 2. Target genes of miRNAs in the mulberry CL library.

https://doi.org/10.1371/journal.pone.0172883.t002

[Figure omitted. See PDF.]

Table 3. Target genes of miRNAs in mulberry DL library.

https://doi.org/10.1371/journal.pone.0172883.t003

[Figure omitted. See PDF.]

Fig 3. T-plots of the targets cleaved by miRNA in CL and DL.

The T-plots show the distribution of 3’ end of the degradome tags within the full-lengthof the target mRNA sequence (bottom). The red line represents the cleaved target tags and is shown in red arrow. The alignment along with the detected cleavage frequencies (absolute numbers) are shown beside the black arrow and it shows the miRNA with a portion of its target sequence (top). The two dots indicate matched RNA base pairs; one dot shows a GU mismatch whereas none dot represent other types of mismatch. The bold letter (nucleotide) on the target transcript indicates the cleavage site showed by the blank arrow. (A) XM_010100268.1, a category 0 target for mno-miR166f in CL. (B) XM_010105580.1, a category 2 target for mno-miR408c in CL. (C) XM_010099938.1, a category 4 target for mno-miR396b in CL. (D) XM_010109976.1, a category 0 target for mno-miR319c in DL. (E) XM_010099345.1, a category 2 target for mno-miR156d in DL. (F) XM_010089439.1, a category 4 target for mno-miR828 in DL.

https://doi.org/10.1371/journal.pone.0172883.g003

Identification of miRNA targets in response to drought stress

All reference genes containing degradation sites were listed in S6 Table. In total, 1,154 and 1,098 targets were investigated in CL and DL, respectively. And 651 (36.07%) targets were found only in DL (Fig 4). Of the target genes with P ≤ 0.05 in DL, the mno-miR156 family conferred the highest number of targets with 17 (45.95%) of the conserved miRNA families (Table 3), and only one miR156-mRNA pair belonged to CL (Table 2), indicating that the miR156 family and their target genes played a key important function in responding to drought stress in mulberry. mno-miR156 targeted a number of mulberry genes annotated as SPB, vesicle transport protein GOT1A, ARF, GRAM domain-containing protein, and chloroplast stem-loop binding protein. Other targete genes of miRNAs were annotated as TFs TCP4, tubulin gamma-1 chain, pentatricopeptide repeat-containing protein, NEDD8-activating enzyme E1 catalytic subunit, auxilin-related protein 2, arginine/serine-rich-splicing factor RSP31, polyamine transporter At3g13620, AT-rich interactive domain-containing protein, mitochondrial import inner membrane translocase subunit, and ATP binding protein.

[Figure omitted. See PDF.]

Fig 4. Venn chart of the target genes detected in the CL and DL libraries.

The number and percentage was the quantity and proportion of target genes in CL specific or DL specific or both CL and DL, respectively.

https://doi.org/10.1371/journal.pone.0172883.g004

GO, KEPP and NR analysis of miRNA target genes in mulberry

The identified targets for miRNAs in the two mulberry degradome libraries were annotated by searching NR databases. By blasting 1800 identified target genes against the NR database, 35,841 functional module structures with the domain sequences were obtained, as shown in S6 and S7 Tables. The highest percentage of the target genes (50.83%) had a global similarity of 60%-80% (Table 4), most of which were targeted by the conserved miRNAs. In the NR annotation by species distribution, the target genes in M. multicaulis had the highest percentage (45.94%) showing similar protein function with Prunus persica, followed by 13.35% in Vitis vinifera, suggesting that M. multicaulis had a similar molecular basis for drought stress with these two species (Table 5).

[Figure omitted. See PDF.]

Table 4. Global analysis of the homology in the target genes by NR annotation.

https://doi.org/10.1371/journal.pone.0172883.t004

[Figure omitted. See PDF.]

Table 5. Functional similarity in the targets with other species.

https://doi.org/10.1371/journal.pone.0172883.t005

To better understand the metabolic pathways of miRNA target genes in mulberry, KEGG was employed to predict the function of those targets, and 119 pathways were involved (S7 and S8 Tables). Importantly, most target genes were associated with metabolic pathways (318 genes), followed by biosynthesis of secondary metabolites (136 genes), plant hormone signal transduction (89 genes), RNA transport (78 genes), and plant-pathogen interaction (76 genes), as shown in S8 Table. The plant hormone signal transduction pathway includes many TFs and key enzymes, such as HD-ZIP and shaggy-related protein kinase (SK), which play important roles under stress response and abiotic stress.

A total of 1,561 miRNA targets identified in DL were classified into three categories according to GO analysis: biological process, cellular component and molecular function, which were classified into 24, 14 and 12 terms respectively, as shown in Fig 5. More than 1000 targets were found to be involved in the cell, cell part, cellular process, metabolic process and organelle. The target genes were mostly associated with stimulus (including 568 genes), regulation of biological process and biological regulation. The results suggested that the corresponding miRNAs regulated the expression of these target genes during drought stress in mulberry by affecting various TFs to induce or shut off specific metabolic networks during response to adversity.

[Figure omitted. See PDF.]

Fig 5. GO analysis of miRNA target genes identified in DL.

The digitals on the left y-axis and right y-axis shows the percentage and enrichment of miRNA targets, respectively. Only the identified target genes for miRNAs by degradome sequencing were considered.

https://doi.org/10.1371/journal.pone.0172883.g005

Overview of miRNA-target networks

The consecutive miRNA-target networks of each library were constructed by Cytoscape 3.2.1 using the miRNAs and target genes. The network structure consists of the basic elements (genes and miRNAs; named nodes) and the connections representing miRNA-target interactions (named edges). To reveal the details of the network, the global networks were categorized into overlapped network and independent networks in CL and DL (Fig 6). Analysis of the nodes of the network, revealed that the independent networks consisted of 910 miRNA-mRNA pairs (65.05%) and 838 miRNA-mRNA pairs (63.34%) in CL and DL, respectively. This finding indicated that the interaction of miRNA-mRNA pairs differed between CL and DL, suggesting that the expression of numerous miRNAs may play roles in regulating the expression of target genes in response to drought stress.

[Figure omitted. See PDF.]

Fig 6. The global profiles of miRNA-target networks in CL and DL trials.

The global profiles of CL and DL miRNA-target networks were constructed, respectively, and included the overlapped network and independent networks of CL and DL. The nodes represent miRNAs or target genes, and the edges represent the connection strength. (A) CL global network. (B) DL global network. (C) Overlapped network. (D) Independent network of CL. (E) Independent network of DL.

https://doi.org/10.1371/journal.pone.0172883.g006

To increase our understanding of the regulatory role of drought-related miRNAs and their target genes, we constructed specific profiles of miRNA mediated interaction network using the target genes with P ≤ 0.05 in DL (Fig 7). This network contained 43 miRNAs and 53 genes, and revealed that many nodes were connected through miRNA-mRNA interaction data and formed a complex network. In this network, most of these targets identified by degradome sequencing were TFs, including ARF, SPB, HD-ZIP, MYB, TCP, SCL, GRAS and SPL. Targets analyzed by Nr annotation included diverse and important enzymes, such as NEDD8-activating enzyme E1 catalytic subunit, copperion-binding protein, type 2A protein phosphatase, transaldolase and CSP41A, suggesting that these genes may play important roles in the stress response.

[Figure omitted. See PDF.]

Fig 7. The specific profiles of miRNA-target networks in DL trials.

The specific profiles of DL miRNA-target network were constructed using the target genes with P ≤ 0.05. The nodes represent miRNAs or target genes, and the edges represent the connection strength.

https://doi.org/10.1371/journal.pone.0172883.g007

Detection of the expression of miRNA and their targets by qRT-PCR

To investigate the regulatory function of miRNAs on their target genes, the expression profiles of 13 miRNA-mRNA pairs (mno-miR166f for XM_010099828.1, mno-miR166f for XM_010100268.1, mno-miR166c for XM_010099828.1, mno-miR171a for XM_010090594.1, mno-miR319c for XM_010109976.1, mno-miR535 for XM_010104394.1, mno-miR4376 for XM_010106326.1, mno-miR156d for XM_010092101.1, mno-miR172a for XM_010106115.1, mno-miRn120-3p for XM_010090594.1, mno-miRn144-5p for XM_010102397.1, mno-miRn46-5p for XM_010093108.1, and mno-miRn202-1-3p for XM_010093270.1) were selected with a P value ≤ 0.05 as criteria from DL for stem-loop qRT-PCR. Among these, eight target genes were cleaved by eight conserved miRNAs and four target genes were cleaved by four novel miRNAs. With the exception of one miRNA-mRNA pair (mno-miR172a for XM_010106115.1) with P value > 0.05, the remaining 12 pairs had P values lower than 0.05.

As shown in Fig 8A, the expression levels of miR166f, miR166c, miR171a, miR319c, miR156d, miR172a and miRn202-1-3p decreased to the lowest value after 10 d under drought stress and were higher at both 5 d and 15 d. Within the same time period, their corresponding targets XM_010099828.1, XM_010100268.1, XM_010090594.1, XM_010109976.1, XM_010092101.1, XM_010106115.1, and XM_010093270.1 (Fig 8B), exhibited apparently the opposite expression profile. miRn144-5p and its target XM_010102397.1 also exhibited the opposite expression profile. miRn144-5p, however, exhibited down-regulation under drought treatment from 5 d to 15 d, and its target XM_010102397.1 exhibited up-regulation. The expression correlation of these nine miRNAs and their targets illustrated that these miRNAs negatively regulated their targets. miR535 and miR4376 exhibited the same expression profile as their targets, XM_010104394.1 and XM_010106326.1, which declined from 5 d and then increased on 15 d. Two miRNA-mRNA pairs, miRn120-3p for XM_010090594.1 and miRn46-5p for XM_010093108.1, showed non-significant regulation, suggesting that one miRNA may regulate multiple target genes or one gene may be regulated by multiple miRNAs. With increased drought stress time, the expression levels of miRNAs appeared to change significantly, and the expression profiles of many miRNAs increased first, then decreased and increased again. Further studies will be needed to shed light on the regulation network of these miRNAs and their targets in response to drought. Over expression or repressing expression of these miRNAs under drought stress may help to elucidate the regulatory mechanism.

[Figure omitted. See PDF.]

Fig 8. Expression analysis of 13 miRNA-mRNA pairs under drought stress with different processing time.

The expression levels of miRNAs and their mRNA targets were normalized to the level of β-actin. The x axis represents different miRNA (A) or mRNA targets (B), and the y axis represents the relative expression level of miRNAs (A) or mRNA targets (B). The results are mean ± SD of the triplicates of three biological replicates. Letter superscripts above bar indicate the significant difference of the changes between different drought stress time with the level of P > 0.05 (shown as the same or no letter superscripts), 0.01 < P ≤ 0.05 (shown as different small letter superscripts) and P ≤ 0.01 (shown as different capital letter superscripts). (A) Expression profile of miRNAs. (B) Expression profile of mRNA targets.

https://doi.org/10.1371/journal.pone.0172883.g008

Determination of six miRNA targets by RLM-5’ RACE

To further confirm the degradome data, RLM-5’ RACE experiments were successfully employed to validate targets of conserved and novel miRNAs.

Four target genes for conserved miRNAs in category 2, 0, 3, 1, which were XM_010092101.1 for mno-miR156d, XM_010099828.1 for mno-miR166f, XM_010100088.1 for mno-miR399e, and XM_010095399.1 for mno-miR408c, respectively, were successfully validated by RLM 5’-RACE experiments (Fig 9A–9D). Additionally, two target genes for novel miRNAs in category 2 and 1, which were XM_010102397.1 for mno-miRn144-5p and XM_010092821.1 for mno-miRn202-1-3p, respectively, were also amplified and were cleaved at the same position found in the degradome library (Fig 9E and 9F).

[Figure omitted. See PDF.]

Fig 9. Validation of target genes cleaved by four conserved miRNAs and two novel miRNAs in DL.

The nucleotide on the target transcript from 3’ end indicated the cleavage site detected in the degradome. The red line represents the cleaved target tags and is shown in red arrow in T-plots. The cleavage site is shown as a bold letter (nucleotide) and absolute numbers of signature sequences are shown above the nucleotide with the bold letter. The black arrow shows a site verified by gene-specific RLM 5’-RACE and the cleavage frequency as determined at the indicated position is shown in parentheses. The two dots indicate matched RNA base pairs; one dot shows a GU mismatch whereas none dot represent other types of mismatch. (A) XM_010092101.1, a category 2 target for mno-miR156d. (B) XM_010099828.1, a category 0 target for mno-miR166f. (C) XM_010100088.1, a category 3 target for mno-miR399e. (D) XM_010095399.1, a category 1 target for mno-miR408c. (E) XM_010102397.1, a category 2 target for mno-miRn144-5p. (F) XM_010092821.1, a category 1 target for mno-miRn202-1-3p.

https://doi.org/10.1371/journal.pone.0172883.g009

Discussion

High-throughput sequencing of RNA degradome in mulberry

Small RNAs play important and fundamental roles in response to adversity, especially in drought stress response [19,52]. A huge number of miRNAs and siRNAs have been identified by cloning and deep sequencing in higher plants [33,39,53,54]. It is the pivotal problem to find target genes for interpretation of miRNA function. Arabidopsis, rice, maize and soybean miRNA targets have been widely studied by high throughput sequencing [27,30,55,56]. In mulberry, miRNA targets were previously investigated mainly via bioinformatics prediction [3,39]. However, this method yields a high level of mismatches in miRNA:target pairing and is laborious, time-consuming and costly for verification of every single predicted gene. Therefore, only a few conserved miRNA targets have been experimentally validated in mulberry [3]. To overcome these shortages, degradome analysis was used in this study. In order to study the regulation of gene expression during drought exposure, we constructed and sequenced two distinct degradome libraries using mulberry leaves under drought stress conditions. Intriguingly, the clean reads from two libraries match perfectly with the mulberry chromosomes and 862, 439 annotated genomic loci were identified in this study (these data have been deposited in the NCBI’s GEO with accession number GSE84889.). Thus, only a small fraction of the sequences are small RNA targets, implying that the vast majority of them are other types of RNA turnover products.

Importantly, our degradome data verified 50–80% of previously validated targets and predicted 15–30% of potential targets [30,31]. Six targets were validated employing RLM 5’-RACE. These targets were found to be cleaved at the same position with different cleavage frequency. Therefore, degradome sequencing, combining the advantages of high throughput deep sequencing and an effective computational approach, served as an efficient strategy to globally identify small RNA targets in plants. Previous studies have indicated that most miRNA targets are cut in the CDS in plants, which is different from animals [28]. Consistent with this finding, cleavage sites of all of the miRNAs were located in the CDS of the target genes identified by degradome sequencing in mulberry. As described by Md Shamimuzzaman [46], other genes may be potentially regulated by miRNA-guided cleavage in the UTR that have not been detected in our alignment analyses. Our identification of many targets captured by the degradome analysis was consistent with previous reports [5,14,27]. However, target genes of several conserved miRNAs (miR394, miR447, miR482, and miR530) and certain novel miRNAs (24.05%) were not detected through this degradome sequencing in mulberry. This phenomenon may result from a low level of some conserved/known non-conserved target genes of miRNAs. Alternatively, some miRNAs might inhibit target gene expression through translational repression [57–59]. To obtain more information on miRNA targets, degradome libraries from different tissues, organs, and developmental stages should be constructed and integrated with the complete mulberry genome sequence.

Drought-responsive miRNAs and their targets in mulberry

Plant miRNAs play a role in the response to abiotic stress, such as drought, NaCl, and low temperature [11,20,55]. The key to determining the function of a miRNA is identifying its target genes. By deep degradome sequencing plants under drought stress, a total of 225 miRNAs targeting 1323 transcripts were detected, including 373 target genes of 30 conserved miRNA families and 950 target genes of 195 novel miRNAs. Under normal growth conditions, 229 miRNAs targeting 1399 transcripts were identified, including 409 target genes of 30 conserved miRNA families and 990 target genes of 199 novel miRNAs.

Among the conserved miRNA families in DL, we found that mno-miR156, mno-miR172 and mno-miR396 had the highest number of targets, which is consistent with previous research. Most of the mRNA target genes identified by degradome sequenceing were transcription factors, including ARF, SPB, HD-ZIP, MYB, TCP, SCL, GRAS, SPL and so on. These TFs played important roles in the response to abiotic stress in plants [30,34]. We found that miR408c targeted oxysterol-binding (OSBP)-related proteins (ORPs) and dehydration response element B (DREB), which are transcription factors involved in drought tolerance. Hajyzadeh et al. [60] reported that upon drought stress treatment, chickpea plants overexpressing miR408 showed normal growth, while other samples struggled with severe symptoms of stress. Meanwhile, the miR408-overexpressed plants exhibited a shorter plant height compared to vector control plants. Through Nr annotation, we also found that many target genes were important enzymes, such as SK, enzyme E1 catalytic subunit, the type 2 a protein phosphatase, transaldolase and CSP41A. These enzymes also play an important function during stress responses in plants [29]. Our results revealed that miR172a negatively regulates XM_010106115.1, a transcript expressed to SK, a plant hormone signal transduction pathway. Interestingly, different target genes were found in different tissues upon drought stress by Eldem et al. [24]. Of the miRNAs differentially expressed in response to drought stress, 60 miRNAs were found to be tissue-specific. 23 miRNAs were only expressed in leaf and 26 miRNAs were only expressed in roots under drought stress. Eldem et al. also found that in leaves, miR156 expression was 3-fold higher, but they did not detect a significant change in roots compared to control samples. The miR5281 was only expressed in leaf with 4-fold upregulation; in contrast, miR415 was downregulated 14-fold in the root but not expressed in the leaf. Target transcripts (137 for leaf control, 133 for leaf stress, 148 for root control and 153 for root stress) generated significant GO terms related to DNA binding and catalytic activities. In the present study, a number of novel miRNAs were found to target genes that: they were annotated as tubulin gamma-1 chain, pentatricopeptide repeat-containing protein, NEDD8-activating enzyme E1 catalytic subunit, auxilin-related protein 2, arginine/serine-rich-splicing factor RSP31, polyamine transporter At3g13620, AT-rich interactive domain-containing protein, mitochondrial import inner membrane translocase subunit, and ATP binding protein.

According to the expression levels of 12 miRNAs and 11 of their targets by qRT-PCR, we found that miR166f, miR166c, miR171a, miR319c, miR156d, miR172a, miRn202-1-3p and miRn144-5p negatively regulated their targets. miR535 and miR4376, however, showed the same trends in expression as their targets. miR156 and miR172 have been shown to regulate their target genes (SBPs and AP2) predominantly by inhibiting the translation of target genes in many plants [61,62]. miR171 has also been found to negatively regulate its target genes in A. thaliana [63]. We found that with increases in the time of drought stress, the expression levels of miRNAs appeared to change significantly, and the expression of many miRNAs increased first, then decreased and riseincreased again.

Enrichment of transcription factors in miRNA targets

In plants, miRNAs target the genes involved in development and stress response, particularly TFs, metabolic transporters, and signal transduction factors. These genes can unravel a new dimension of the miRNA regulatory network [64]. Our analysis revealed that the majority of these target genes, classified into 51 different annotation categories (Fig 5), were conserved among species [46,65,66]. The target genes are obviously enriched in TFs and transcription regulatory activity, such as NAC, MYB, ARF, SPL, SCL, HD-ZIP, APS, NF-Y subunits, MADS-box, GRAS and TCP (S6 and S7 Tables). This phenomenon was similar to that found in Arabidopsis, rice, soybean, and maize [29,30,34]. A number of the TFs are known to regulate diverse aspects of plant growth, development and responses to environmental changes. For example, MYB family members, such as MYB33 and MYB101, which are targeted by miR159 in rice and soybean, appear to play an important role in the response to abscisic acid (ABA) during seedling stress responses, suggesting that ABA-induced accumulation of miR159 is a homeostatic mechanism to direct MYB33 and MYB101 transcript degradation to desensitize hormone signaling during seedling stress responses. [46,67]. In this study, we detected a number of MYB family TFs regulated by mno-miR159, mno-miR858 and mno-miRn148-5p. miR159 expression is regulated by at least two different hormones, GA and ABA, and regulates the abundance of MYB mRNAs in response to drought stress [67]. Another important family of TFs, the HD-ZIP TF family, determines adaxial/upper cell fate [68]. Our findings indicate that the expression of miR166f increased significantly and its target gene XM_010100268.1, a HD-ZIP TF, was down-regulated during drought stress. According to GO analysis, XM_010100268.1 participates in biological regulation, metabolic process and nucleic acid binding transcription factor activity. NAC domain transcriptional regulator gene was the target of miR394 and miR482 [60], implying that NAC domain transcriptional regulators might also be regulated by other miRNAs. SPL TFs, which are unique to plants, are involved in embryonic development, plastochron length, leaf development, developmental phase transitions, flower and fruit development, fertility, apical dominance, anthocyanin biosynthesis, gibberellin response, light signaling and copper homeostasis. The miR156 family was found to target SPL in this study, consistent with a study involving cotton [60]. Auxin, an important phytohormone, acts as a key player in plant development in higher plants [69]. As transducers of auxin signaling, ARF and TIR play vital roles in plant development under drought stress conditions, including shoot, root and flower formation [70–73]. We also identified miR160, miR172 and some novel miRNAs involved in auxin signaling. NF-Y subunits are known to control a variety of important agronomical traits, including drought tolerance, flowering time, and seed development [74]. In Arabidopsis, overexpression of a miR169-resistant NF-YA5 transgene significantly improves drought resistance through the ABA-dependent pathway by promoting stomatal closure under drought stress [75]. In this study, we discovered two novel miRNAs, miRn71-3p and miRn223-5p, that targeted NF-Y. Four MADS-box genes have previously been validated as targets for miR444 variants in rice and they are regulated by several miR444 variants [29]. In mulberry, we found that miR169 targets the MADS-box.

Regulation networks of miRNAs and targets during drought stress in mulberry

The rapid discovery of miRNAs and experimental evidence for miRNA interactions has ushered in a new era of miRNA research that focuses on networks other than individual miRNA, interactions. The conservation of resistance to adverse environments in plants can be assessed by the summation of gene expression and regulation network connectivity, which can also provide a new avenue for understanding molecular mechanisms and distinguishing functional processes in resistance progression. Here, the consecutive miRNA target networks of each library were constructed. The network structure consists of the basic elements (genes and miRNAs, named nodes) and the connections representing miRNA-target interactions (named edges). After construction of the network, we found that the topological profiles are more similar to a ‘Medusa’model [76], which consists of a regulatory core of hub nodes. This finding indicates that the kernel nodes of the network are determinants in the gene expression levels [40,77].

Generally, miRNA inhibits translation or induces mRNA degradation by binding to the CDS of target mRNAs in plants. As shown in the results, from the overlapped network and the independent network of miRNA-targets in CL and DL groups based on the degradome data (Figs 6 and 7), we identified differential miRNAs and their targets that are involved in drought stress. In degradom, 838 miRNA-mRNA pairs (63.34%), consisting of 272 miRNAs and 732 genes were found only in DL. The cluster results of the network showed that a miRNA can be applied to different mRNAs and a mRNA can also be regulated by a number of different miRNAs. Conversely, the changes in gene expression induced by a miRNA can also act to regulate this that miRNA or other miRNAs. This means that the interactions of miRNAs and mRNAs are not one-to-one, but can form interaction networks with cross regulation. Because of the cross regulation of target genes by miRNAs, small changes in the expression levels of miRNAs or mutations in the pairing sequences will produce serious consequences. This is similar to the cascade amplification effect in the signal transduction process [78]. Butz et al. found that the effect of miRNA on RCC has ‘divergent’ properties whereby the same miRNA targets multiple genes. It is also ‘convergent’ in nature, whereby multiple miRNAs have augmented effects on the same target genes [43]. Hua et al. [79] found that co-transfection of miR-20a and miR-361 had cumulative effects on inhibition of VEGF expression. Wu et al. [80] found that 28 miRNAs can significantly inhibit the expression of p21Cip1/Waf1. Xu et al. [78] constructed the miRNA functional networks (miRFNs) of soybean and found that miRFNs of soybean exhibit a scale-free, small world and modular architecture, with their degrees fit best to power-law and exponential distribution.

Among the miRNA-mRNA interaction pairs, the most important drought responsive target genes were TFs as described above. According to data analysis of miRNA target genes predicted by PicTar, Cui et al. [81] and Boyer et al. [82] also found that miRNAs are more inclined to regulate genes regulated by TFs. Here, we obtained 43 kernel miRNAs and 53 target genes from the network, which are highly correlated with the drought resistance in mulberry. Of the miRNA-mRNA pairs showing significant relationships (P ≤ 0.05), the mno-miR156 family conferred the highest number of targets of the 17 (45.95%) genes (Table 3), and one miR156-mRNA pair was found in CL (Table 2). In rice, the miR164-targeted NAC genes are known to be negative regulators of drought tolerance [83]. Our study revealed that miR156 targeted two target genes (XM_010104767.1 and XM_010092557.1), which belong to the NAC TFs family. Our study also demonstrated that miR156 can target SPL TF genes (miR156 for XM_010090407.1 and XM_010092348.1) and can regulate tolerance to recurring environmental stress as reported by Wang et al. [66] and Stief et al. [84]. These findings indicate that the miR156 family and their target genes may also play important roles in the response to drought stress in mulberry.

Conclusions

This is the first comprehensive identification of drought responsive miRNAs and their targets in mulberry by degradome sequencing. This study provides a framework for further understanding the molecular mechanisms of resistance to drought in mulberry. Further experiments on the identified target genes are necessary to explore potential mechanisms of abiotic stress regulation in mulberry.

Supporting information

[Figure omitted. See PDF.]

S1 Table. The stem-loop RT and qRT-PCR primer sequences for miRNAs expression analysis.

https://doi.org/10.1371/journal.pone.0172883.s001

(DOCX)

S2 Table. The primer sequences for qRT-PCR analysis of target genes.

https://doi.org/10.1371/journal.pone.0172883.s002

(DOCX)

S3 Table. The gene specific primer sequences for RLM-5’RACE of target genes.

https://doi.org/10.1371/journal.pone.0172883.s003

(DOCX)

S4 Table. The information of miRNA-mRNA density in CL.

https://doi.org/10.1371/journal.pone.0172883.s004

(XLSX)

S5 Table. The information of miRNA-mRNA density in DL.

https://doi.org/10.1371/journal.pone.0172883.s005

(XLSX)

S6 Table. Nr annotation of ref genes in mulberry.

https://doi.org/10.1371/journal.pone.0172883.s006

(XLSX)

S7 Table. Target genes blast by KEGG in mulberry.

https://doi.org/10.1371/journal.pone.0172883.s007

(XLSX)

S8 Table. KEGG pathway participated by target genes in mulberry.

https://doi.org/10.1371/journal.pone.0172883.s008

(XLSX)

Author Contributions

1. Conceptualization: Ruixue Li WZ.

2. Data curation: Ruixue Li WZ.

3. Formal analysis: Ruixue Li FH HZ WZ.

4. Funding acquisition: WZ Ruixue Li.

5. Investigation: Ruixue Li DC Rongfang Li Y. Wang.

6. Methodology: Ruixue Li WZ.

7. Resources: Ruixue Li TW RF LL WZ.

8. Validation: Ruixue Li WZ.

9. Visualization: Ruixue Li.

10. Writing – original draft: Ruixue Li.

11. Writing – review & editing: Ruixue Li Y. Wan WZ.

Citation: Li R, Chen D, Wang T, Wan Y, Li R, Fang R, et al. (2017) High throughput deep degradome sequencing reveals microRNAs and their targets in response to drought stress in mulberry (Morus alba). PLoS ONE 12(2): e0172883. https://doi.org/10.1371/journal.pone.0172883

References

1. Bartel DP. MicroRNAs: Target recognition and regulatory functions. Cell. 2009;136(2):215–33. pmid:19167326

2. Kurihara Y, Watanabe Y. Arabidopsis micro-RNA biogenesis through Dicer-like 1 protein functions. Proc Natl Acad Sci USA. 2004;101(34):12753–8. pmid:15314213

3. Wu P, Han SH, Zhao WG, Chen T, Zhou JC, Li L, et al. Genome-wide identification of abiotic stress-regulated and novel microRNAs in mulberry leaf. Plant Physiol Biochem. 2015;95:75–82. pmid:26188501

4. Voinnet O. Origin, Biogenesis, and activity of plant microRNAs. Cell. 2009; 136(4): 669–87. pmid:19239888

5. Xu F, Liu Q, Chen LY, Kuang J, Walk T, Wang JX, et al. Genome-wide identification of soybean microRNAs and their targets reveals their organ-specificity and responses to phosphate starvation. BMC Genomics. 2013;14(1):1–30.

6. Baumberger N, Baulcombe DC. Arabidopsis ARGONAUTE1 is an RNA Slicer that selectively recruits microRNAs and short interfering RNAs. Proc Natl Acad Sci USA. 2015;102(33):11928–33.

7. Chen X. Small RNAs and their roles in plant development. Annu Rev Cell Dev Biol. 2005;25(25): 21–44.

8. Puzey JR, Karger A, Axtell M, Kramer EM. Deep annotation of Populus trichocarpa microRNAs from diverse tissue sets. PLoS One. 2012;7(3): e33034–e33034. pmid:22442676

9. Chen L, Ren YY, Zhang YY, Xu JC, Zhang ZY, Wang YW. Genome-wide profiling of novel and conserved populus microRNAs involved in pathogen stress response by deep sequencing. Planta. 2012;235(5): 873–83. pmid:22101925

10. Qin Y, Duan Z, Xia X, Yin W. Expression profiles of precursor and mature microRNAs under dehydration and high salinity shock in Populus euphratica. Plant Cell Rep. 2011;30(10):1893–907. pmid:21706230

11. Sunkar R, Li YF, Jagadeeswaran G. Functions of microRNAs in plant stress responses. Trends Plant Sci. 2012;17(4):32–41.

12. Fahlgren N, Howell MD, Kasschau KD, Chapman EJ, Sullivan CM, Cumbie JS, et al. High-throughput sequencing of Arabidopsis microRNAs: Evidence for Frequent Birth and Death of MIRNA Genes. PLoS One. 2007;2(2):e219–e219. pmid:17299599

13. Srivastava S, Zheng Y, Kudapa H, Jagadeeswaran G, Hivrale V, Varshney RK, et al. High throughput sequencing of small RNA component of leaves and inflorescence revealed conserved and novel miRNAs as well as phasiRNA loci in chickpea. Plant Sci. 2015;235:46–57. pmid:25900565

14. Zhang JZ, Ai XY, Guo WW, Peng SA, Deng XX, Hu CG. Identification of miRNAs and their target genes using deep sequencing and degradome analysis in trifoliate orange [Poncirus trifoliata L. Raf]. Mol Biotechnol. 2011;51(1):44–57.

15. Budak H, Kantar M, Bulut R, Akpinar BA. Stress responsive miRNAs and isomiRs in cereals. Plant Sci. 2015;235:1–13. pmid:25900561

16. Jagadeeswaran G, Saini A, Sunkar R. Biotic and abiotic stress down-regulate miR398 expression in Arabidopsis. Planta. 2009;229(4):1009–1014. pmid:19148671

17. Li B, Qin Y, Duan H, Yin W, Xia X. Genome-wide characterization of new and drought stress responsive microRNAs in Populus euphratica. J Exp Bot. 2011;62(11):3765–3779(15). pmid:21511902

18. Li L, Tang J, Zhang B, Yang W, Liu M, Wang R, et al. Epigenetic modification of miR-429 promotes liver tumour-initiating cell properties by targeting Rb binding protein 4. Gut 2015;64(1):156–67. pmid:24572141

19. Sunkar R. MicroRNAs with macro-effects on plant stress responses. Semin Cell Dev Biol. 2010;21(8):805–11. pmid:20398781

20. Wang T, Chen L, Zhao M, Tian Q, Zhang WH. Identification of drought-responsive microRNAs in Medicago truncatula by genome-wide high-throughput sequencing. BMC Genomics. 2011;12(12):367.

21. Akdogan G, Tufekci ED, Uranbey S, Unver T. miRNA-based drought regulation in wheat. Funct Integr Genomic. 2016;16(3):221–233.

22. Xie FL, Wang QL, Sun RR, Zhang BH. Deep sequencing reveals important roles of microRNAs in response to drought and salinity stress in cotton. J Exp Bot. 2015;66(3):789–804. pmid:25371507

23. Zhou M, Luo H. Role of microRNA319 in creeping bentgrass salinity and drought stress response. Plant Signal Behav. 2014;9(4):e28700. pmid:25764437

24. Eldem V, Akçay UÇ, Ozhuner E, Bakır Y, Uranbey S, Unver T. Genome-Wide Identification of miRNAs Responsive to Drought in Peach (Prunus persica) by High-Throughput Deep Sequencing. PLoS One. 2012;7(12):e50298. pmid:23227166

25. Millar AA, Waterhouse PM. Plant and animal microRNAs: similarities and differences. Funct Integr Genomics. 2005;5(3):129–135. pmid:15875226

26. Bartel DP. MicroRNAs: genomics, biogenesis, mechanism, and function. Cell. 2004;116(2): 281–297. pmid:14744438

27. Zhou M, Gu LF, Li PC, Song XW, Wei LY, Chen Zhiyu, et al. Degradome sequencing reveals endogenous small RNA targets in rice (Oryza sativa L. ssp. indica). Front Biosci. 2010;5(1):67–90.

28. Jones-Rhoades MW, Bartel DP. Computational identification of plant microRNAs and their targets, including a stress-induced miRNA. Mol Cell. 2004;14(6):787–799. pmid:15200956

29. Li YF, Zheng Y, Addo-Quaye C, Zhang L, Saini A, Jagadeeswaran G, et al. Transcriptome-wide identification of microRNA targets in rice. Plant J. 2010;62(5):742–59. pmid:20202174

30. Addo-Quaye C, Eshoo TW, Bartel DP, Axtell MJ. Endogenous siRNA and miRNA targets identified by sequencing of the Arabidopsis degradome. Curr Biol. 2008;18(10):758–62. pmid:18472421

31. German MA, Pillay M, Jeong DH, Hetawal A, Luo S, Janardhanan P, et al. Global identification of microRNA-target RNA pairs by parallel analysis of RNA ends. Nat Biotechnol. 2008;26(8): 941–6. pmid:18542052

32. Gregory BD, O'Malley RC, Lister R, Urich MA, Tonti-Filippini J, Chen HM, et al. A link between RNA metabolism and silencing affecting Arabidopsis development. Dev Cell. 2008;14(6):854–66. pmid:18486559

33. Zhao M, Tai HH, Sun SZ, Zhang FS, Xu YB, Li WX. Cloning and characterization of maize miRNAs involved in responses to nitrogen deficiency. PloS One. 2012;7(1):e29669–e29669. pmid:22235323

34. Song QX, Liu YF, Hu XY, Zhang WK, Ma B, Chen SY, et al. Identification of miRNAs and their target genes in developing soybean seeds by deep sequencing. BMC Plant Biol. 2011;11(1):5–20.

35. He NJ, Zhao AC, Qin J, Zeng QW, Xiang ZH. Mulberry genome project and mulberry industry. Zhongguo Zhong yao za zhi, Science of Sericulture. 2012;38(1):140–145.

36. Song PH, Zeng QW, Shang JZ, Ma B, Xiang ZH, He JN. Research progress in plant responses to water logging stress. Zhongguo Zhong yao za zhi, Science of Sericulture. 2013;39(1):160–165.

37. German MA, Luo SJ, Schroth G, Meyers BC, Green PJ. Construction of parallel analysis of RNA ends (PARE) libraries for the study of cleaved miRNA targets and the RNA degradome. Nat Protoc. 2009;4(3):356–62. pmid:19247285

38. Addo-Quaye C, Miller W, Axtell M. CleaveLand: a pipeline for using degradome data to find cleaved small RNA targets. Bioinformatics. 2009;25(1):130–1. pmid:19017659

39. Jia L, Zhang DY, Qi XW, Ma B, Xiang ZH, He JN. Identification of the conserved and novel miRNAs in mulberry by High-throughput sequencing. PloS One. 2014;9(8):e104409–e104409. pmid:25118991

40. Allen E, Xie ZX, Gustafson AM, Carrington JC. microRNA-directed phasing during trans-acting siRNA biogenesis in plants. Cell. 2005;121(2):207–221. pmid:15851028

41. Smoot ME, Ono K, Ruscheinski J, Wang PL, Ideker T. Cytoscape 2.8: new features for data integration and network visualization. Bioinformatics. 2011;27(3):431–442. pmid:21149340

42. Chen QL, Lu YY, Peng JH, Dong S, Wei B, Song YN, et al. Dynamical regulation analysis identifies molecular mechanisms of Fuzheng-Huayu formula against hepatitis B-caused liver cirrhosis. Evid Based Complement Alternat Med. 2015;2015(5):238495.

43. Butz H, Szabó PM, Khella HWZ, Nofech-Mozes R, Patocs A, Yousef GM. miRNA-target network reveals miR-124 as a key miRNA contributing to clear cell renal cell carcinoma aggressive behaviour by targeting CAV1 and FLOT1. Oncotarget. 2015;6(14):12543–12557. pmid:26002553

44. Chen C, Ridzon DA, Broomer AJ, Zhou Z, Lee DH, Nguyen JT, et al. Real-time quantification of microRNAs by stem-loop RT-PCR. Nucleic Acids Res. 2005;33(20):e179. pmid:16314309

45. Wu RM, Wood M, Thrush A, Walton EF, Varkonyi-Gasic E. Real-time PCR quantification of plant miRNAs using universal probe library technology. Biochemica. 2007;202(2):12–15

46. Shamimuzzaman M and Vodkin L. Identification of soybean seed developmental stage-specific and tissue-specific miRNA targets by degradome sequencing. BMC Genomics. 2012;13(1):1–14.

47. Boke H, Ozhuner E, Turktas M, Parmaksiz I, Ozcan S, Unver T. Regulation of the alkaloid biosynthesis by miRNA in opium poppy. Plant Biotechnol J. 2015;13(3):409–420. pmid:25735537

48. Schmid R and Blaxter ML. annot8r: GO, EC and KEGG annotation of EST datasets. BMC Bioinformatics. 2008;9(1):1–6.

49. Du Z, Zhou X, Ling Y, Zhang Z, Su Z. AgriGO: a GO analysis toolkit for the agricultural community. Nucleic Acids Res. 2010;38:W64–70. pmid:20435677

50. Nigam D, Kadimi PK, Kumar S, Mishra DC, Rai A. Computational analysis of miRNA-target community network reveals cross talk among different metabolisms. Genom Data. 2015;5:292–296. pmid:26484271

51. Xing Z, Chu C, Chen L, Kong X. The use of Gene Ontology terms and KEGG pathways for analysis and prediction of oncogenes. Biochim Biophys Acta. 2016;In press.

52. Yin FQ, Qin C, Gao J, Liu M, Luo XR, Zhang WY, et al. Genome-wide identification and analysis of drought-responsive genes and microRNAs in tobacco. Int J Mol Sci, 2015;16(3):5714–5740. pmid:25775154

53. Axtell MJ, Bowman JL. Evolution of plant microRNAs and their targets. Trends Plant Sci. 2008;13(7):343–349. pmid:18502167

54. Reinhart BJ, Weinstein EG, Rhoades MW, Bartel B, Bartel DP. MicroRNAs in plants. Genes Dev. 2002;16(13):1616–26. pmid:12101121

55. Zhao YP, Xu ZH, Mo QC, Zou C, Li WX, Xu YB, et al. Combined small RNA and degradome sequencing reveals novel miRNAs and their targets in response to low nitrate availability in maize. Ann Bot. 2013;112(3):633–642. pmid:23788746

56. Rhoades MW, Reinhart BJ, Lim LP, Burge CB, Bartel B, Bartel DP. Prediction of plant microRNA targets. Cell. 2002;110(4):513–520. pmid:12202040

57. Fabian MR, Sonenberg N, Filipowicz W. Regulation of mRNA translation and stability by microRNAs. Annu Rev Biochem. 2010;79:351–79. pmid:20533884

58. Huang J, Liang Z, Yang B, Tian H, Ma J, Zhang H. Derepression of microRNA-mediated protein translation inhibition by apolipoprotein B mRNA-editing enzyme catalytic polypeptide-like 3G (APOBEC3G) and its family members. J Biol Chem. 2007;282(46):33632–40. pmid:17848567

59. Yang XY, Wang LC, Yuan DJ, Lindsey K, Zhang XL (2013) Small RNA and degradome sequencing reveal complex miRNA regulation during cotton somatic embryogenesis. J Exp Bot. 2013;64(6):1521–1536. pmid:23382553

60. Hajyzadeh M, Turktas M, Khawar KM, Unver T. miR408 overexpression causes increased drought tolerance in chickpea. Gene. 2015;555(2):186–193. pmid:25445265

61. Gandikota M, Birkenbihl RP, Hohmann S, Cardon GH, Saedler H, Huijser P. The miRNA156/157 recognition element in the 3’ UTR of the Arabidopsis SBP box gene SPL3 prevents early flowering by translational inhibition in seedlings. Plant J, 2007;49(4):683–693. pmid:17217458

62. Schwab R, Palatnik JF, Riester M, Schommer C, Schmid M, Weigel D. Specific effects of microRNAs on the plant transcriptome. Dev Cell. 2005;8(4):277–284.

63. Liu HH, Tian X, Li YJ, Wu CA, Zheng CC. Microarray-based analysis of stress-regulated microRNAs in Arabidopsis thaliana. RNA. 2008;14(5):836–843. pmid:18356539

64. Chen F, Zhang XF, Zhang N, Wang SS, Yin GH, Dong ZD, et al. Combined small RNA and degradome sequencing reveals novel miRNAs and their targets in the high-yield mutant wheat strain Yunong 3114. PloS One. 2015;10(9):e0137773–e0137773. pmid:26372220

65. Yang JH, Liu XY, Xu BC, Zhao N, Yang XD, Zhang MF. Identification of miRNAs and their targets using high-throughput sequencing and degradome analysis in cytoplasmic male-sterile and its maintainer fertile lines of Brassica juncea. BMC Genomics. 2013;14(1):1–16.

66. Wang YQ, Li L, Tang S, Liu JG, Zhang HS, Zhi H, et al. Combined small RNA and degradome sequencing to identify miRNAs and their targets in response to drought in foxtail millet. BMC Genet. 2016;17(57):1–16.

67. Reyes J, Chua NH. ABA induction of miR159 controls transcript levels of two MYB factors during Arabidopsis seed germination. Plant J. 2007;49(4):592–606. pmid:17217461

68. Nogueira FT, Madi S, Chitwood DH, Juarez MT, Timmermans MC. Two small regulatory RNAs establish opposing fates of a developmental axis. Genes Dev. 2007;21(7):750–755. pmid:17403777

69. Yang XY, Zhang XL, Yuan DJ, Jin FY, Zhang YC, Xu J. Transcript profiling reveals complex auxin signaling pathway and transcription regulation involved in dedifferentiation and redifferentiation during SE in cotton. BMC Plant Biol. 2012;12(4):1–19.

70. Wang D, Pei K, Fu Y, Sun Z, Li S, Liu HQ, et al. Genome-wide analysis of the auxin response factors (ARF) gene family in rice (Oryza sativa). Gene. 2007;394(1–2):13–24. pmid:17408882

71. Si-Ammour A, Windels D, Arn-Bouldoires E, Kutter C, Ailhas J, Meins F, et al. miR393 and secondary siRNAs regulate expression of the TIR1/AFB2 auxin receptor clade and auxin-related development of Arabidopsis leaves. Plant Physiol. 2011;157(2):683–91. pmid:21828251

72. Matsui A, Mizunashi K, Tanaka M, Kaminuma E, Nguyen AH, Nakajima M, et al. tasiRNA-ARF Pathway Moderates Floral Architecture in Arabidopsis Plants Subjected to Drought Stress. Biomed Res Int. 2014;2014(11):303451–303451.

73. Li T, Ma L, Geng YK, Hao CY, Chen XH, Zhang XY. Small RNA and degradome sequencing reveal complex roles of miRNAs and their targets in developing wheat grains. PloS One. 2015;10(10):e0139658–e0139658. pmid:26426440

74. Cao S, Kumimoto RW, Siriwardana CL, Risinger JR, Holt BF. Identification and characterization of NF-Y transcription factor families in the monocot model plant Brachypodium distachyon. PloS One. 2011;6(6):e21805–e21805. pmid:21738795

75. Li WX, Oono YJ, He XJ, Wu JM, Iida K, Lu XY, et al. The Arabidopsis NFYA5 transcription factor is regulated transcriptionally and posttranscriptionally to promote drought resistance. Plant Cell. 2008;20(8):2238–2251. pmid:18682547

76. Guo Y, Feng Y, Trivedi NS, Huang S. Medusa structure of the gene regulatory network: dominance of transcription factors in cancer subtype classification. Exp Biol Med. 2011;236(5):628–636.

77. Zheng S, Tansey WP, Hiebert SW, Zhao Z. Integrative network analysis identifies key genes and pathways in the progression of hepatitis C virus induced hepatocellular carcinoma. Bmc Med Genomics. 2011;4(1):1–10.

78. Xu YG, Guo MZ, Liu X Y, Wang CY, Liu Y. Inferring the soybean (Glycine max) microRNA functional network based on target gene network. Bioinformatics. 2014;30(1):94–103. pmid:24149053

79. Zhong H, Qing L, Ye WB, Chung-Kwun AW, Cai GP, Gu DY, et al. MiRNA-directed regulation of VEGF and other angiogenic factors under hypoxia. PloS One.2006;1(2):e116.

80. Wu S, Huang S, Ding J, Zhao Y, Liang L, Liu T, et al. Multiple microRNAs modulate p21Cip1/Waf1 expression by directly targeting its 3' untranslated region. Oncogene. 2010;29(15):2302–2308. pmid:20190813

81. Cui Q, Yu Z, Pan Y, Purisima EO, Wang E. MicroRNAs preferentially target the genes with high transcriptional regulation complexity. Biochen Bioph Res Co. 2007;352(3):733–738.

82. Boyer LA, Lee TI, Cole MF, Johnstone SE, Levine SS, Zucker JP. et al. Core transcriptional regulatory circuitry in human embryonic stem cells. Cell. 2005;122(6):947–956. pmid:16153702

83. Fang Y, Xie K, Xiong L. Conserved miR164-targeted NAC genes negatively regulate drought resistance in rice. J Exp Bot. 2014;65(8):2119–2135. pmid:24604734

84. Stief A, Altmann S, Hoffmann K, Pant BD, Scheible WR, Baurle L. Arabidopsis miR156 regulates tolerance to recurring environmental stress through SPL transcription factors. Plant Cell. 2014;26(4):1792–1807. pmid:24769482