Identification and Characterization of Salvia miltiorrhizain miRNAs in Response to Replanting Disease

Replanting disease is a major factor limiting the artificial cultivation of the traditional Chinese medicinal herb Salvia miltiorrhiza. At present, little information is available regarding the role of miRNAs in response to replanting disease. In this study, two small RNA libraries obtained from first-year (FPR) and second-year plant (SPR) roots were subjected to a high-throughput sequencing method. Bioinformatics analysis revealed that 110 known and 7 novel miRNAs were annotated in the roots of S. miltiorrhiza. Moreover, 39 known and 2 novel miRNAs were identified and validated for differential expression in FPR compared with SPR. Thirty-one of these miRNAs were further analyzed by qRT-PCR, which revealed that 5 miRNAs negatively regulated the expression levels of 7 target genes involved in root development or stress responses. This study not only provides novel insights into the miRNA content of S. miltiorrhiza in response to replanting disease but also demonstrates that 5 miRNAs may be involved in these responses. Interactions among the differentially expressed miRNAs with their targets may form an important component of the molecular basis of replanting disease in S. miltiorrhiza.


Introduction
Salvia miltiorrhiza Bunge is a very popular traditional Chinese medicinal plant. The economic importance of S. miltiorrhiza results from the medicinal activity of extracts of its tuberous roots, which include tanshinones and salvianolic acid, and the herb has been used extensively to treat a variety of conditions, especially cardiovascular and cerebrovascular diseases [1][2][3]. Unfortunately, while several studies have described how to improve the major medicinal bioactive constituents of S. miltiorrhiza [4][5][6][7][8], little attention has been paid to S. miltiorrhiza yields and quality.
The cultivation area of S. miltiorrhiza in China has expanded over the years because of increased demands for the plant in the domestic and international markets. Environmental been grown in the previous year. For convenience, we named members of the former group as first-year plants and members of the latter group as second-year plants. Root samples were obtained from five independent plants at the tuberous root expansion stage (Aug. 10,2015), and their RNA content was extracted using an improved CTAB-LiCl method [18].

Identification of novel miRNA in S. miltiorrhiza
Raw read sequences in the two libraries were combined into one small RNA library for novel miRNA prediction, and all reads matching the tRNA, rRNA, or conserved miRNA sequences with two mismatches were removed using Bowtie [19]. The remaining reads were mapped to the assembled transcript sequences (downloaded from http://bi.sky.zstu.edu.cn/Bio111/ DsTRD/home.php) [20] using Bowtie with perfect matching [19]. With one end anchored 20 bp away from the mapped small RNA location, 120-360 bp sequences, each with an extension of 20 bp covering the small RNA region, were collected using a Perl script. Secondary structures of each sequence were predicted using RNAfold from the Vienna package (version 1.8.2) [21]. Under conditions similar to those proposed by Thakur et al. [22], a stem-loop structure with 3 gaps involving 8 bases at the small RNA location was considered a candidate miRNA precursor.

Identification of replant-responsive miRNAs
The frequency of the miRNAs from two libraries was normalized to 1 million by total clean reads of miRNAs (RPM) in each sample. If the normalized read count of a given miRNA was zero, the expression value was modified to 0.01 for further analysis. The fold change between the SPR and FPR libraries was calculated as follows: fold change = log 2 (SPR/FPR). miRNAs with fold changes > 1 or < −1 and with P 0.001 were respectively considered to be upregulated or downregulated in response to replanting disease. The P-value was calculated according to previously established methods [23].

Quantitative RT-PCR analysis
Total RNA was treated with RNase-free DNase I (TaKaRa, Dalian, China) to remove genomic DNA. Forward primers for 5 selected miRNAs were designed based on the sequence of the miRNAs and are listed in S1 Table. The reverse transcription reaction was performed with the One Step PrimeScript miRNA cDNA Synthesis Kit (TaKaRa, Dalian, China) according to the manufacturer's protocol.
SYBR Green PCR was performed following the manufacturer's instructions (Takara, Japan). In brief, 2 μl of cDNA template was added to 10 μl of 2× SYBR Green PCR master mix (Takara), 1 μM each primer, and ddH 2 O to a final volume of 20 μl. The reactions were amplified for 1min at 95°C, followed by 40 cycles of 95°C for 5 s and 60°C for 34 s. All reactions were performed in triplicate, and the controls (no template and no RT) were included for each gene. The threshold cycle (C T ) values were automatically determined by the instrument. The fold-changes were calculated using 2 -Delta DeltaCt method, where DeltaDeltaCt = (Ct target − Ct inner ) treatment − (Ct target − Ct inner ) control [26].

Results
Deep sequencing of sRNAs in S.miltiorrhiza root Two sRNA libraries were constructed from the root tissues of FPR and SPR to identify miRNAs responding to replanting disease in S. miltiorrhiza roots. The sequencing results were submitted to the GenBank SRA database (accession numbers: SRR3056582 and SRR3056449). A total of 14.73 million raw reads were generated from the two sRNA libraries. Low-quality tags and adaptor contaminations were removed to obtain 7,374,092 (representing 2,050,638 unique sequences) and 3,581,868 (representing 1,512,866 unique sequences) clean reads from the FPR and SPR libraries, respectively; these reads ranged from 18 nt to 30 nt in both groups. Sequences of about 21 and 24 nt in length were dominant in both libraries (Fig 1). This result is consistent with those presented in previous studies on other plant genera, such as Arabidopsis [27], Oryza [28], Medicago [29], and Populus [30].

Identification of known miRNAs in S. miltiorrhiza
Unique sRNA sequences were mapped to miRBase 21.0 with perfect matches using Perl script [31]. A total of 110 unique sequences belonging to 31 conserved miRNA families were identified in the FPR and SPR libraries (S2 Table). Among these sequences, 87miRNAs belonging to 27 families (except for miR408, miRNA827, miR5139 and miR8155) were detected in FPR, while 96 miRNAs belonging to 28 families (except for miR397, miR399 and miR2111) were detected in SPR. Between the two libraries, the more abundant sequences were found in the miR156, miR396, miR 319, and miR166 families. The miR319 family, which was composed of 14 members, dominated the sequences (Fig 2).
Identification and validation of novel miRNAs in S. miltiorrhiza sRNA reads that were homologous to known miRNAs and other non-coding sRNAs (Rfam 10) were excluded, and the secondary structures of the precursors of the remaining 19-24 nt sRNAs were analyzed using RNAfold software to determine novel miRNAs in S. miltiorrhiza. Precursors with canonical stem-loop structures were further analyzed using more stringent filters to ensure that they satisfied the criteria established by the research community [32,33]. Thirty-two miR-NAs candidates derived from 71 loci satisfied the screening criteria. The 32 putative miRNA precursors were then used to extract miRNA Ã s, which are considered strong evidence of DICER-LIKE-1 (DCL1)-derived products. Thirteen of the 32 miRNA candidates contained miRNA star (miRNA Ã ) sequences identified from the same libraries; the remaining 19 candidates did not contain any miRNA Ã sequences (S3 Table). The 13 candidates with miRNA Ã sequences were considered to be novel S. miltiorrhiza miRNAs, whereas the19 remaining candidates without miRNA Ã sequences were considered potential S. miltiorrhiza miRNAs. The stem-loop structures and miRNA Ã sequences of the13 novel miRNAs are shown in S1 Fig. RT-PCR was performed to validate the 13 novel miRNAs and determine their expression patterns in the root, stem, leaf, and flower of S. miltiorrhiza. The electrophoresis results indicated that the 7 novel miRNAs all showed amplification in the root, stem, leaf, and flower of S. miltiorrhiza (S2 Fig). Seven novel miRNAs were analyzed for tissue-specificity by qRT-PCR, and the results showed obvious tissue specificity. The expressions of all of the novel miRNA tissues, except that of miR017, were higher in the stem and leaf than in the root and flower. MiR021a, miR028, and miR031 expression levels were highest in the leaf, while miR025a expression levels were highest in the stem. The expression levels of miR031 in the root were lower by about 6 times than those in the stem, leaf, and flower of S. miltiorrhiza. miRNA017 gene expression was lower in the stem and leaf than in the root and flower; the highest expression of this gene was observed in the flower (Fig 3).The 7 novel miRNAs miR001a, miR008a, miR012a, miR017, miR021a, miR028 and miR031 were renamed as smi-miR35829, smi-miR35830, smi-miR35831, smi-miR35832, smi-miR35833, smi-miR35834 and smi-miR35835, respectively.

Differential expression identification of known and novel miRNAs involved in response to replanting disease
To identify the miRNAs involved in the response of S. miltiorrhiza to replanting disease, differential expressions in the two libraries were estimated from the read counts determined through high-throughput sequencing. We selected RPMs that exceeded 1 and presented a P value of less than 0.001 from both libraries. miRNAs that met this condition and exhibited a log2 (SPR RPM/FPR RPM) fold-change higher than 1 were designated as upregulated, whereas those with a log2 (SPR RPM/FPR RPM) fold-change of less than −1 were designated as downregulated (S2 and S3 Tables).
The read counts of 110 known unique sequences from both libraries were retrieved to determine miRNAs responding to replanting disease. A total of 39 known miRNAs from 15 families and 2 novel miRNAs were differentially expressed in response to replanting disease (Table 1). Among these differentially expressed miRNAs, 38 were upregulated and only 3 were downregulated in the SPR library compared with those in the FPR library. Among these 41 sequences, the miR166, miR319 and miR396 families, which possessed 8, 5 and 5 members, respectively, were in a dominant position. Only pab-miR160a-like was absent from the FPR library.

Identification of S. miltiorrhiza miRNAs by high-throughput sequencing
Identification of miRNAs in medicinal plants, such as R. glutinosa [13], Panax ginseng [34], and Pinellia pedatisecta [35], using high-throughput sequencing or miRNA arrays has previously been reported. High-throughput sequencing is particularly useful for identifying involved abiotic and biotic stress responses [36,37]. In this study, miRNA libraries constructed from the FPR and SPR were used to identify novel and replant-responsive miRNAs in S. miltiorrhiza, and their mechanisms of action were subsequently investigated. Among the miRNAs identified using high-throughput sequencing, 65% of those already known were expressed at low levels (less than 10 raw reads; S1 Table). This finding suggests that high-throughput sequencing is a powerful tool for identifying poorly expressed miRNAs in plants. The miR160, miR165/166, and miR396 families, which are believed to target auxin response factor (ARF), homeodomain leucine zipper (HD-ZIP III), and growth-regulating factor (GRF) transcription factors, respectively, were abundantly represented in both libraries. By comparing the expression levels of all members of a miRNA family, dominant members, such as miR166 and miR396, may be found. These dominant members may perform key regulatory roles in response to stress. Some family members, such as Smi-miR166a/b/c/e/f/g/h/I in the miR166 family, exhibited comparable expression levels, thereby indicating that several members of a family may exert synergistic effects in the relevant regulatory network.
With the discovery of miRNAs as ubiquitous regulators of plant growth and development, including nearly all biological, metabolic, and stress processes, researchers have focused on miRNAs in S. miltiorrhiza. However, given that the complete genome of S. miltiorrhiza is currently unavailable, the discovery of miRNAs from S. miltiorrhiza is relatively limited.
Xu et al. [24] studied tissue-specific miRNAs in S.miltiorrhiza and identified 164 miRNAs after redundancy elimination; of these miRNAs, 28 were known, 22 were novel, and another 114 miRNAs were considered homologous with known miRNAs. Some of these 114 miRNAs were 1-4 bases to the left or right of known miRNA sequences or differed from known miRNA sequences in terms of one base. In this study, we identified 117 miRNAs, 110 of which belonged to 31 known families and 7 of which were novel miRNAs. While twenty-three known miRNAs were consistently identified between the present study and Xu et al.'s research [24] (S5 Table), no novel miRNA were identified between both works, likely because assembled transcriptome sequences were used as reference RNAs for the miRNA precursors analysis in our study whereas ESTs were used in Xu et al.'s work.

Medicinal plant miRNAs responding to replanting
Yang et al. [13] identified 16 miRNAs from 13 miRNA families, including 2 novel miRNAs, responding to R. glutinosa replanting disease. In this study, we identified 31 miRNAs from 14 miRNA families, including 1 novel miRNA responding to S. miltiorrhiza replanting disease. In addition, three replant-responsive miRNAs, namely sim-miR156a and sim-miR160a, were identified both in S. miltiorrhiza and R. glutinosa (S6 Table), suggesting that these 2 miRNAs might have important role in respons to replanting disease.

Conclusion
Our study provides a global perspective of differential miRNA expression in S. miltiorrhiza in response to replanting disease for the first time. We identified 110 known and 7 novel miRNAs, of which 5 known miRNAs were expressed differentially in the first-and second-year crops. Further investigation revealed that the target genes of differentially expressed miRNAs were related to root growth and development. Combining with the phenotype of poor growth in second-year crop, our results revealed that replanting disease upregulates the expressions of smi-miR156a-1, pab-miR160a-like, smi-miR164a-1, ath-miR165a-3p-like, and gma-miR396e-like. Correspondingly, their target genes including SPL13, ARF18-like, NAC100-like, athb-14-like and GRF3-like which were related to root growth were downregulated and then might be presented a growth inhibition.