The B. elliptica-resistant Lilium regale Wilson and the B. elliptica-susceptible Asian hybrid “Tresor” as host plants were planted on culture medium under a 12 h day/night cycle at 25/22°C and placed in a greenhouse at Beijing Forestry University. The plants for treatment were conducted during the flower bud stage. Botrytis elliptica isolated from diseased lily leaves was grown on potato dextrose agar medium under near-UV light for 7 days. For biotic stress treatment, leaves detached from the fifth to tenth sites below the apex were first carefully washed under distilled water and then surface-sprayed with 5 × 10⁴ conidia·mL⁻¹ of B. elliptica (Liu et al., 2016). Subsequently, all of the treated leaves were placed in trays covered with preservative film to ensure a relative humidity of 90–100%. When leaf tissues were harvested at 0, 6, 12, and 24 h postinoculation (hpi), sample pools were established from at least five independent plants for each sampling time, The sampling time set from early experiments (Liu et al., 2016) and gene expression variation in the transcriptome sequencing in our laboratory (data not published). All of the collected samples were immediately frozen in liquid nitrogen and stored at −80°C until RNA extraction.

Total RNA from lily leaves of the control and treated samples was extracted using the EASYspin Plus Plant RNA kit (Aidlab Bio, Beijing, China) according to the manufacturer's protocol. Three small RNA libraries of L. regale (B. elliptica-mock: “BE0h,” B. elliptica-infected for 6 h: “BE6h,” B. elliptica-infected for 24 h: “BE24h”) were constructed by Solexa/Illumina sequencing (LC Bio, Hangzhou, China). Briefly, the population of small RNAs was initially isolated, purified, and subsequently ligated to Solexa adaptors at each end. Then, the assembled small RNAs were reverse-transcribed to cDNA, followed by PCR, and the purified PCR products were eventually used for cluster generation and sequencing analysis with the ACGT101-miR program (version 4.2; LC Bio). The results were deposited in the Sequence Read Archive (SRA) at NCBI database (accession number: SRP103981).

After the sequencing reactions, the raw reads were first subjected to filtering out adaptor dimers, junk reads and low-quality tags, and then mapped onto transcriptome of L. regale. Subsequently, perfectly matched sequences were compared against non-coding RNAs from Rfam database (http://rfam.sanger.ac.uk/) and the NCBI GenBank database (http://www.ncbi.nlm.nih.gov/genbank) to classify the degradation fragments of non-coding RNAs, which were excluded from further analysis. The remaining sequences were aligned with the mature and precursors of miRNAs in miRBase 20.0 (http://mirbase.org) by BlastN to identify conserved miRNAs. Furthermore, the unannotated small RNAs were analyzed for the prediction of novel miRNAs.

For novel miRNAs identification, the hairpin RNA structures of novel miRNAs were predicted by RNAfold software and the criteria were accorded as follows: (1) the miRNA and miRNA^* are located in opposite stem-arms and form a duplex with two nucleotide 3′ overhangs; (2) the number of mismatched bases between miRNA and miRNA^* <4; (3) asymmetric bulges are minimal in size (typically two or less); (4) the potential miRNA precursor with minimal folding energy indexes (MFEI) >0.6 (Jones-Rhoades et al., 2006; Meyers et al., 2008).

Stem-loop RT-PCR was performed to validate the identification novel miRNAs. Total miRNAs were extracted from samples and reverse-transcribed to cDNA using QuantScript RT Kit (Tiangen Bio, Beijing, China) with specific RT primer (Supplementary Table S4). The RT-PCR reaction was continued at 94°C for 3 min, followed by 30 cycles of 95°C for 30 s, 55°C for 30 s, and 72°C for 15 s, ending with 72°C for 5 min with 2 × Taq PCR MasterMix (Aidlbab Bio, Beijing, China). Amplification products were separated with 2.5% agarose gel electrophoresis.

To identify B. elliptica-responsive miRNAs of L. regale, the miRNA abundance in the three libraries was normalized by the total number of miRNAs per sample (Zhang et al., 2016). If the normalized read count of a given miRNA was zero, the expression value was modified to 0.01 for further analysis. The expression levels of miRNAs between stress and control groups (BE6h/BE0h, BE24h/BE0h) were calculated as follows: fold-change = log₂ (BE6h/BE0h or BE24h/BE0h). When the fold-change with a p ≤ 0.01 was more than 2.0 or less than −0.5, the miRNA was considered to be significantly differentially expressed in response to B. elliptica. The p-value was calculated according to previously established methods (Man et al., 2000).

To investigate the potential target mRNAs, a degradome library was constructed as previously described (Addo-Quaye et al., 2008). Similar to the small RNA libraries, the degradome cDNA library was sequenced on Illumina GAIIx (LC Bio, Hangzhou, China). After removing adaptor sequences and low quality sequencing reads, the clean reads were used to identify potentially cleaved targets based on L. regale assembled transcriptome sequences. Meanwhile, the degradome analysis and the classification of target categories were performed with CleaveLand 3. Moreover, Cytoscape 3.2 software was used to generate and visualize the miRNA-mediated gene regulatory network.

The expression of 12 selected miRNAs and 6 targets from the different times postinoculation was assayed using qRT-PCR to validate results obtained from high-throughput sequencing of miRNAs and their targets. Total RNAs and miRNAs were extracted from samples and reverse-transcribed to cDNA using ReverTra Ace® qPCR RT Master Mix with gDNA Remover (Toyobo, Shanghai, China) following the manufacturer's instructions (Chen et al., 2012; Liu et al., 2016). qRT-PCR was performed using a total reaction volume of 20 μL, which contained 2 μL of diluted cDNA, 1 μL primer mix, 10 μL of 2 × qPCR mix, and 7 μL ddH₂O. Each reaction was performed on a MiniOpticon Real-time PCR Detection System (Bio-Rad, Hercules, CA, USA), using SYBR Premix Ex Taq (Takara, Otsu, Japan) in triplicate. Amplification reactions were performed as follows: 95°5 for 5 s, 69°9 for 20 s, and 72°2 for 20 s. All reactions were performed in triplicate. The 18S rRNA and CLATHRIN genes were used as the reference genes for normalization and the relative expression levels were calculated using the 2^−heCT method (Ginzinger, 2002). The primers for amplification were designed using the Beacon Designer 7 software, as listed in Supplementary Table S5.

To reveal Botrytis elliptica-responsive microRNAs in lily, three sRNA libraries (BE0h, BE6h, and BE24h) were constructed and sequenced using Illumina Genome Analyzer IIx (LC Bio, Hangzhou, China). Deep sequencing of the libraries generated a total of 26.44 million raw reads representing 3.52 million unique sequences from the three sRNA libraries (Table 1). Based on removal of the low-quality tags and adaptor contaminants, 2,556,296, 6,981,340, and 5,049,471 clean sequences were obtained for these three libraries, respectively, with a size of 17–25 nt. After further filtering out the Rfam (rRNA, tRNA, snRNA, snoRNA, and other Rfam RNAs) and Repbase sequences, the remaining reads were queried against miRBase 20.0 database. From that, totals of 2,065,391 (representing 431,479 unique sequences), 5,981,827 (representing 724,767 unique sequences), and 4,162,571 (representing 587,701 unique sequences) sRNA sequences, respectively, were obtained and used for the identification of conserved and novel miRNAs (Table 1).

The length distributions of the sRNA sequences, which ranged from 17 to 25 nt, were similar among the three libraries (Figure 1), with 24 nt representing the most common length, and followed by 21 and 22 nt. These results are consistent with the observed length distribution of mature plant miRNAs, such as in Arabidopsis thaliana (Sunkar and Zhu, 2004), Solanum lycopersicum (Jin and Wu, 2015), and Raphanus sativus L. (Wang Y. et al., 2015; Figure 1). Specifically, the relative abundances of 24-nt sRNAs in the BE6h and BE24h groups were markedly lower than those in the elliptica-mock library, suggesting that the 24-nt miRNA classes might be repressed under B. elliptica infection. Conversely, both 21- and 22-nt sequences in the BE6h and BE24h were notably more abundant than those in BE0h library. Thus, the sRNAs with different sizes probably perform different functions: 21-nt sRNAs usually mediate gene silencing at the post-transcriptional level, while 24-nt sRNAs typically perform gene silencing mediated by RNA-dependent DNA methylation and heterochromatin maintenance (Duan et al., 2016).

To identify known miRNAs in the three libraries, the clean reads were compared with known miRNA precursor or mature miRNA sequences in miRBase 20.0, allowing no more than two mismatches. Ultimately, 71 non-redundant miRNAs belonging to 47 conserved families were identified in the three libraries (Supplementary Table S1), which showed high sequence similarity to those of other plant species in miRBase (Figure 2). Among the 47 miRNA families, the numbers of the conserved miRNA family members were dramatically different (Supplementary Table S1 and Figure 3A). For example, the miR166 family was the most highly represented with eight members, followed by miR159 and miR396 with five members each. Of the remaining miRNA families, most possessed a single member, although these families are known to contain more members in other plants (e.g., miR158, miR171, miR172, and miR319). Moreover, unique sequences, found in few other plant species, were detected in the lily libraries. For instance, miR5072 has only been registered in Oryza sativa, and miR5054 in Brachypodium distachyon (Schreiber et al., 2011).

The relative abundance of different conserved miRNA families can be estimated by the number of reads, which exhibited large divergences. A few conserved miRNA families, such as miR159, miR166, miR167, miR396, miR894, and miR5368, showed extraordinarily high expression abundance among the three libraries, of which the miR166 family exhibited the most abundant expression with reads accounting for over 50% of all conserved miRNA reads (Figure 3B). Moreover, miR162, miR164, miR827, miR5072, miR6173, miR6300, and miR6478 were moderately abundant in the three libraries, although they contained only a single member each. However, a number of miRNA families, such as miR395, miR482, miR858, and miR1883 showed relatively low expression levels, with no more than five reads in the three libraries (Figure 3B). The abundance of miRNA families can dramatically vary, suggesting functional variation within each family. In addition, distinct members of the same miRNA family also showed dramatically different levels of expression. For instance, the abundance of miR166 members varied from 0 to 8,911 reads in the BE0h library (Supplementary Table S1), likely due to function-specific expression of different miRNA members (Jiang et al., 2006).

The remaining reads that had not been mapped to selected miRNAs in miRBase, but perfectly matched to the transcriptome of L. regale, were identified as potentially novel miRNAs, of which the precursors with canonical stem-loop structures were further analyzed using a series of stringent filter strategies to ensure that they satisfied the common criteria for annotation of plant miRNAs established by the research community (Jones-Rhoades et al., 2006; Meyers et al., 2008). Based on these criteria, 24 novel miRNA sequences were identified, of which only five were determined to contain complementary miRNA^*s (Supplementary Table S2). Because of the instability of miRNA^* in cells, miRNAs without miRNA^*s were also regarded as novel miRNAs. The length of mature novel miRNAs varied from 17 to 25 nt, with the majority being 21 nt, and their precursors ranged from 55 to 282 nt with a mean length 151 nt (Supplementary Table S2). The hairpin structures of the precursors of miRn3, miRn4, and miRn5 are presented in Figure 4. The values of minimum folding free energy (MFE) of these predicted pre-miRNAs ranged from −185.8 to −19 kcal mol⁻¹ (Supplementary Table S2). Consistent with the characteristics of miRNAs, the MFEI ranged from 0.6 to 3. The expression levels of the novel miRNAs varied significantly in the three libraries (Supplementary Table S2). Although most of the novel miRNAs had relatively low expression, these miRNAs might be specific to lily. To validate these novel miRNAs, stem-loop RT-PCR was performed in miRn3, miRn5, miRn22, and miRn24 as examples. As shown in Supplementary Figure S1, these novel miRNAs had signal on agarose gels, indicating that they can be considered as novel miRNAs of lily.

To systematically identify B. elliptica-responsive miRNAs in L. regale, a normalized differential expression pattern analysis was performed between the B. elliptica-mock library and the two B. elliptica-infected libraries. Among all of the identified conserved and novel miRNAs, 24 conserved and 7 novel miRNAs were identified to be significantly differentially expressed upon comparing these three libraries (p ≤ 0.01, log₂Ratio > 2 or < −0.5; Figure 5, Supplementary Table S3). Of these, 8 known miRNAs, namely, miR160, miR166a, miR166c, miR319, miR393a, miR399b, miR408, miR5368b, and 3 novel miRNAs were upregulated continuously in elliptica-infected leaves (BE6h and BE24h) compared with their levels in elliptica-mock (BE0h); in contrast, 13 known and 4 novel miRNAs were downregulated in BE6h and BE24h (Figure 5). Intriguingly, two miRNAs (miR168a and miR6478) firstly declined in BE6h and then peaked in BE24h; while miR396c exhibited the opposite trend. Among all of the markedly differentially expressed miRNAs, miRNA319 was considered to be the most significantly upregulated known miRNA, with a 3.85-fold increase of expression (BE6h/BE0h), while miR164 was the most downregulated. We found that the most differentially expressed miRNAs were downregulated under B. elliptica treatment. These significant variations in miRNA expression implied that these miRNAs could play crucial roles in the response to B. elliptica of lily.

Transcriptome-wide analysis of the miRNA-cleaved mRNAs using degradome sequencing technology was performed to understand the biological function of miRNAs in lily. After discarding the low-quality sequences, a total of 13,026,025 raw reads with 3,146,014 unique raw reads were remained. Then, 5,550,212 of the sequences were mapped to the L. regale transcriptome to identify the fragments of degraded mRNAs. In total, 22 target transcript sequences that could potentially be cleaved by 10 miRNAs (including 9 known and 1 novel ones) were identified through the Cleveland pipeline (Figure 6, Table 2). Generally, several targets were regulated by only a single miRNA. For instance, miR166b and miR5059 had the most target genes (miR166b targeted five genes and miR5059 targeted seven genes), indicating that these two families might play diverse roles in resistance to B. elliptica stress. Furthermore, six known miRNAs (miR319, miR408, miR1862, miR2275, miR5072, and miR5142) had only one target each.

The target genes of miRNAs were classified in a centralized functional distribution. Most genes were involved in metabolism, including diverse metabolic enzymes such as 12-OXOPHYTODIENOATE REDUCTASE 3 (OPR3) and AGMATINE DEIMINASE (ADI), which were cleaved by miR5059 and miR1862, respectively (Table 2). Besides, two target genes, encoding the transcription factors EF1a and TCP2, were targeted by miR166b and miR319 (Table 2). These transcription factors could participate in various aspects of plant development and stress responses and also act as the main nodes in gene expression networks in plants. HEAT SHOCK PROTEIN (HSP), which could be involved in response to biotic and abiotic stresses, was a target gene of miR5059. Generally, most of the targets were involved in various types of signal sensing and transduction, such as reactive oxidative species (ROS) signaling, the biosynthesis of jasmonic acid (JA), and cell wall synthesis, which were considered to be involved in the response to B. elliptica. Although sequencing revealed that the novel miRNAs were present at a lower expression level than the known miRNAs, their targets were also identified by degradome sequencing analysis. However, unknown, hypothetical, or predicted proteins and even no targets were identified for other miRNAs in the degradome sequencing data, partly due to the limited number of accessible lily reference sequences.

To confirm the results obtained from small RNA deep sequencing, quantitative real-time PCR was performed to confirm the expression levels of 12 miRNAs among the different times of B. elliptica infection of L. regale. As shown in Figure 7, most of the expression patterns of B. elliptica-responsive miRNAs from qRT-PCR in L. regale were similar to those obtained by deep sequencing. However, the fold-change of expression obtained by qRT-PCR was not completely consistent with the results of bioinformatic analysis, which partially due to the differences in the sensitivity, specificity, and algorithm between the two techniques.

To further identify the function of conserved and novel miRNAs in lily, the expression levels were detected by qRT-PCR in both the B. elliptica-resistant L. regale and the B. elliptica-susceptible Asian hybrid “Tresor.” As shown in Figure 8, the expression of most miRNAs seems to share the same trend in the two species during pathogen attack and much lower expression levels were detected in the susceptible cultivar “Tresor.” Of these, five miRNAs (miR319, miR399b, miR408, miR5059, and miRn22) were upregulated and peaked at 6 or 12 h in both the two species. Another three miRNAs (miR156a, miR167a, and miR5072) had a downregulated expression pattern. Meanwhile, three miRNAs seem to show different expression patterns of the two species with different resistant to B. elliptica. For example, the level of miR164 transcripts sharply dropped at 6 h, and then continued at a low level in L. regale. However, the expression of miR164 was raised at 24 h in “Tresor.” Besides, miR166b and miR171a were exhibited an opposite expression profile in the two species. In the resistant L. regale, miR166b was downregulated after B. elliptica infection; to the contrary, miR166b was gradually upregulated at 24 h in “Tresor.” These differently expressed miRNAs might play key roles in resistant of B. elliptica.

To further confirm the dynamic correlation between the miRNAs and their corresponding targets, the expression patterns of six target genes were examined by qRT-PCR of L. regale at different time points of B. elliptica infection (Figure 9). As expected, an approximately negative correlation was observed between miRNAs and their targets. For example, EF1a and NHE-6 (Conting 1715 and Conting 39906), two targets of miR166b, were upregulated at 6 h; whereas miR166b sharply declined at that time. Moreover, the expression level of TCP2 was significantly reduced at 6 h, the time at which the expression of miR319 was relatively enhanced. It could be found that the relative expression trends between the B. elliptica-responsive miRNAs and their corresponding targets were generally inversed.

miRNAs, as key factors in numerous cellular events, play multiple roles in regulating many biological processes in plants. With the development of high-throughput sequencing and bioinformatic approaches, miRNAs have been identified from various plant species without fully sequenced genomes, including Hemerocallis fulva (An et al., 2014), Phalaenopsis Aphrodite (An and Tsair, 2012), Asparagus officinalis (Chen et al., 2016), and Rehmannia glutinosa (Yang et al., 2011b). Little research has been conducted to identify miRNAs from L. regale, because of the limited genomic data which restricts sRNAs annotation and identification of the many sources of Lilium sRNA production. In this study, we conducted the first screen for L. regale miRNAs by deep sequencing with three B. elliptica-infection libraries. To increase the accuracy, an mRNA transcriptome from B. elliptica-infected L. regale provided by our laboratory was also used as a reference sequence for miRNA high-throughput sequencing and could provide more valuable information for degradome sequences analyses. Moreover, the expression pattern of transcriptome mRNAs could also be a reference of miRNA-target gene interaction.

This research resulted in the characterization of 71 known miRNAs across 47 conserved families in lily. Most miRNA families identified in lily have also been found in other model plants. For example, the sequence of lre-miR166 family had a perfect match with miR166 in other species, such as Glycine max, Populus trichocarpa and Oryza sativa (Figure 2). This implies that the miRNAs are well-developed in both monocotyledons and dicotyledons, and might play the similar roles in their growth and development. Moreover, some miRNAs, that had been identified in just a few species, even in only one other species, were found in L. regale. For example, miR5059 and miR6300, which have to date only been registered in Brachypodium distachyum (Schreiber et al., 2011) and Glycine max (Turner et al., 2011), also appeared in lily. Furthermore, 24 novel miRNAs with low abundance were discovered in Lilium, and these miRNAs with limited expression might be specifically involved in the response to B. elliptica in lily.

Solexa sequencing not only provides the foundation for identifying known and novel miRNAs, but is also a way of estimating the expression profiles of miRNA (Wang et al., 2009). Here, via high-throughput sequencing, we found that 24 known and 7 novel miRNAs were significantly responsive to B. elliptica infection compared with their levels in the uninfected control. Among the differentially expressed miRNAs, the majority were downregulated under diseases stress, which is in accordance with findings in other species suffering biotic stress (Luo et al., 2015; Wang and Luan, 2015). For instance, including miR164, miR166, and miR168, the expression of these miRNAs were induced under Verticillium longisporum infection in Brassica napus (Shen et al., 2014), which were exhibit the same expression pattern in lily under B. elliptica-infection. Moreover, other miRNAs could be specifically regulated by B. elliptica, mainly Lilium-specific ones, including miR5059, miR5139, miR6300, and novel miRNAs, which have rarely been reported in other species. On the other hand, several previously reported disease-responsive miRNAs, such as miR394 in tomato (Jin et al., 2012; Jin and Wu, 2015), were not detected in the current study, which suggested that these miRNAs might be species-specific and that their expression was not specifically altered in lily, or perhaps these miRNAs did not exhibit altered expression during the duration of the study.

In order to further investigate the roles of miRNAs in response to B. elliptica, L. regale, and “Tresor” with different levels of resistance to this pathogen were processed for qRT-PCR. It was observed that the higher levels of pathogen-induced expression in resistant plants compared with the pathogen infected susceptible plants, suggesting overexpression of miRNAs might enhance the resistance to B. elliptica. For example, miR319 was upregulated eight times higher at 6 h in the resistant L. regale than in the susceptible “Tresor.” Previous reports revealed that miR319 was involved in plant responses to drought, salinity, and biotic stresses (Glazov et al., 2008; Naqvi et al., 2010; Jin and Wu, 2015). The different expression level of miR319 between the resistant and susceptible lily suggested that this miRNA might be responsive to B. elliptica. Furthermore, miR166 and miR171a were detected to show different expression trend after inoculation between the B. elliptica-resistant L. regale and the B. elliptica-susceptible “Tresor.” It has been reported that miR171 was regulated leaf formations and development (Long et al., 2010). In tomato, miR171a was also detected to be downregulated in response to B. cinerea through microarray analysis (Jin et al., 2012). The expression level of miR171 was low possibly due to stress from B. elliptica infection to leaves. These results implied that the expression of miRNAs was significantly different because of the plant adaptation and the resistance to B. elliptica. However, there is still a need for further studies involving deep profiling of the differential expression patterns and to confirm the precise regulatory roles of these B. elliptica-responsive miRNAs in lily.

To further explore the function of miRNAs in defense against B. elliptica, a transcriptome-wide analysis of the degradome was performed, and the target transcripts for the known and novel miRNAs were determined. Through this approach, some miRNAs were easily detectable, but their targets had not been identified yet. Thus, we applied strict selection criteria and identified only 22 targets of these miRNAs by degradome analysis (Table 2). These targets were shown to be involved in signal transduction and plant responses to disease, including transcription factors, enzymes, or proteins in signal pathways, which were regulated by miRNAs through transcript cleavage or translational repression (Sunkar and Zhu, 2004). In our study, it was shown that TCP2 was likely targets of miR319, which was upregulated after B. elliptica inoculation as indicated by decrease of its targets TCP transcript, the same as in other reports (Jin and Wu, 2015). As a transcription factor, TCP2 could most likely downregulated the expression of the ubiquitin proteolytic pathway, which is extensively involved in plant growth, development, and signal transduction of many physiological and biochemical responses (Debigaré and Price, 2003; Yang et al., 2013). Moreover, several studies have shown that TCP TFs could bind the TCP-recognized motif (GGACCAC) in the promoter of LIPOXYGENASE (LOX), which is an enzyme involved in JA biosynthetic pathway (Maksymiec et al., 2005; Schommer et al., 2008). Generally, JAs as endogenous hormones can respond to various biotic and abiotic stresses, which have already been reported in Arabidopsis (McConn et al., 1997) and Vitis (Jia et al., 2016). Furthermore, OPR3, which also controls one of the most important steps of JA biosynthesis (Chen et al., 2006), was identified as a target of miR5059 and upregulated upon 6 h B. elliptica infection. These results implied that disease could activate the biosynthesis and accumulation of JA (Zhao et al., 2003), and the JA signaling pathway might increase plant resistance against Botrytis (Jia et al., 2016). Another main target of miR5059 was COPPER CHAPERONE (CCH), a scavenger enzyme of ROS detoxifying superoxide radicals (Beauclair et al., 2010). Downregulation of this miRNA resulted in an increase in CCH expression, which could enhance the tolerance to oxidative stress and reduce the damage to plants. E2, targeted by miR2275 in lily, which was also targeted by miR399 in rice, could participate in inorganic phosphate (Pi)-signaling pathways, and Pi is a common limiting factor for plant growth (Bari et al., 2006). Several miRNA targets could encode proteins involved in responses to environmental stress. For instance, the HSP family is one of the largest groups of proteins induced by heat shock stress; these proteins are present in almost all organisms and have significant functions in cellular homeostasis under adverse environmental conditions (Njemini et al., 2004). In the present study, miR5059 could regulate the expression of HSPs, which were also previously reported to be targeted by miR2616 in Medicago sativa (Shu et al., 2016) and miR2592 in Raphanus sativus (Wang R. et al., 2015), showing that the expression of HSPs could be regulated by B. elliptica infection. In summary, target genes mediated by miRNAs, which could participate in JA biosynthetic and ROS pathway, enhance the level of defense or tolerance to B. elliptica in L. regale.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.