Skip to main content


Front. Plant Sci., 18 May 2017
Sec. Plant Pathogen Interactions

Transcriptome-Wide Analysis of Botrytis elliptica Responsive microRNAs and Their Targets in Lilium Regale Wilson by High-Throughput Sequencing and Degradome Analysis

Xue Gao1, Qi Cui1, Qin-Zheng Cao1, Qiang Liu1, Heng-Bin He1, Dong-Mei Zhang2,3 and Gui-Xia Jia1*
  • 1Beijing Key Laboratory of Ornamental Plants Germplasm Innovation and Molecular Breeding, National Engineering Research Center for Floriculture, Beijing Laboratory of Urban and Rural Ecological Environment and College of Landscape Architecture, Beijing Forestry University, Beijing, China
  • 2Shanghai Academy of Landscape Architecture Science and Planning, Shanghai, China
  • 3Shanghai Engineering Research Center of Landscaping on Challenging Urban Sites, Shanghai, China

MicroRNAs, as master regulators of gene expression, have been widely identified and play crucial roles in plant-pathogen interactions. A fatal pathogen, Botrytis elliptica, causes the serious folia disease of lily, which reduces production because of the high susceptibility of most cultivated species. However, the miRNAs related to Botrytis infection of lily, and the miRNA-mediated gene regulatory networks providing resistance to B. elliptica in lily remain largely unexplored. To systematically dissect B. elliptica-responsive miRNAs and their target genes, three small RNA libraries were constructed from the leaves of Lilium regale, a promising Chinese wild Lilium species, which had been subjected to mock B. elliptica treatment or B. elliptica infection for 6 and 24 h. By high-throughput sequencing, 71 known miRNAs belonging to 47 conserved families and 24 novel miRNA were identified, of which 18 miRNAs were downreguleted and 13 were upregulated in response to B. elliptica. Moreover, based on the lily mRNA transcriptome, 22 targets for 9 known and 1 novel miRNAs were identified by the degradome sequencing approach. Most target genes for elliptica-responsive miRNAs were involved in metabolic processes, few encoding different transcription factors, including ELONGATION FACTOR 1 ALPHA (EF1a) and TEOSINTE BRANCHED1/CYCLOIDEA/PROLIFERATING CELL FACTOR 2 (TCP2). Furthermore, the expression patterns of a set of elliptica-responsive miRNAs and their targets were validated by quantitative real-time PCR. This study represents the first transcriptome-based analysis of miRNAs responsive to B. elliptica and their targets in lily. The results reveal the possible regulatory roles of miRNAs and their targets in B. elliptica interaction, which will extend our understanding of the mechanisms of this disease in lily.


Lily (Lilium spp.) is one of the most valuable commercial market flower bulbs in the world, mainly owing to its various ornamental functions as a cut flower or a potted plant. However, the production and reproduction of lily have been impeded by Botrytis species, an aggressive fungal pathogen. It can infect more than 200 plant species, of which Botrytis elliptica specializes in lily (Van Baarlen et al., 2007; Weiberg et al., 2013). This pathogen infects leaves, stems and flowers of lily during plant growth as well as survival. To limit yield losses caused by gray mold, agronomic, genetic, and biological approaches have been proposed, however, cultivars are still threatened (Angelini et al., 2014). Although certain cultivars show some level of resistance to B. elliptica, few Lilium species or lily cultivars are absolutely immune to B. elliptica infection (Kohl et al., 1995; Lu et al., 2007). L. regale Wilson is native to western Sichuan Province in China, which is specifically distributed in the river basin of Minjiang. It was found to possess extremely high resistance to viruses and fungal pathogens, such as F. oxysporum and B. elliptica (Lim et al., 2003; Du et al., 2014). Therefore, it has been widely regarded as a main plant to investigate the disease-resistant mechanism of lily. For example, Rao et al. (2013) selected L. regale as materials to isolate the genes differentially expressed in a resistant to F. oxysporum reaction using the suppression subtractive hybridization. Liu et al. (2016) fully investigated the reference genes of L. regale under different abiotic and biotic stresses because of the advantages in breeding for further research. Thus, L. regale was utilized as experimental material in our study to explore the biochemical and genetic bases of resistance to B. elliptica.

MicroRNAs (miRNAs), ~21–24 nucleotides (nt) in length, are a subset of small, endogenous non-coding RNAs that mediate gene expression at the post-transcriptional level by targeting mRNAs for degradation or translational repression (Carrington and Ambros, 2003; Jover-Gil et al., 2005). Increasing evidences indicates that miRNAs are intimately involved in developmental processes (Wu et al., 2009; Kang et al., 2012; Zhang et al., 2012), flowering (Wang et al., 2012, 2014; Luo et al., 2013), and also many adaptive responses to both abiotic and biotic stresses in plants (Jones-Rhoades et al., 2006; Zhao C. et al., 2015; Niu et al., 2016; Shu et al., 2016). Roles of miRNAs in plant defense responses have been identified in Solanum lycopersicum (Jin et al., 2012; Jin and Wu, 2015), Arabidopsis thaliana (Kurubanjerdjit et al., 2014), Paulownia fortune (Niu et al., 2016), Arachis hypogaea L. (Zhao C. et al., 2015), and Carica papaya (Abreu et al., 2014), showing that some miRNAs can silence genes involved in immunity. For example, by using microarray technology, three known miRNAs (miR160, miR169, and miR171) were identified to exhibit differential expression profiles under Botrytis cinerea infection, which regulated adaptions of tomato seedlings at the post-transcriptional level (Jin et al., 2012). Recently, high-throughput sequencing technologies have provided unprecedented opportunities to generate both known and novel miRNAs by constructing small RNA libraries. Thus, it has been extensively applied to discovering miRNA in plant disease. For instance, via high-throughput sequencing, 41 significantly differentially expressed known miRNAs and 39 novel miRNAs were identified in response to Rhizoctonia solani (Luo et al., 2015). Among these miRNAs, miR171, miR408, and miR398 were found to be upregulated at a high level after R. solani treated 6, 12, and 24 h (Luo et al., 2015). Besides, Luan et al. (2015) found a total 70 miRNAs which were manifested to change significantly in samples treated with P. infestans using high-throughput sequencing, including 50 downregulated miRNAs and 20 upregulated miRNAs. Moreover, miRNAs that may be involved in plant disease can be further studied combined with the analysis of target genes. In peanut, the accumulation of some miRNAs was altered after infection with bacterial wilt disease, and more than 10% of their targets were shown to be defense response genes (Zhao C. et al., 2015). Functional identification and analysis of miRNAs and their targets in response to pathogen infection could provide valuable information for understanding the mechanisms underlying disease resistance. Concerning gray mold disease, Solanum lycopersicum (Jin and Wu, 2015) and Paeonia lactiflora (Zhao D. et al., 2015) have been reported to have miRNAs with functional roles in resistance to B. cinerea. However, the diversity of miRNAs and their roles in response to B. elliptica-infection in lily remain largely unexplored.

Given their huge and complex genomes, research on lily miRNA has been delayed relative to that on model plants, such as Arabidopsis (Llave et al., 2002). However, better opportunities to characterize miRNAs in plants have arisen via high-throughput sequencing platforms, especially in non-model organisms for which the full genome sequence has yet to be acquired, such as Cunninghamia lanceolata (Wan et al., 2012), Rehmannia glutinosa (Yang et al., 2011a), Citrusre ticulata (Guo et al., 2016), and Asparagus officinalis (Chen et al., 2016). Moreover, cleaved miRNA targets can also be identified by degradome sequencing at the transcriptome level, overcoming the problem that predicted target genes must otherwise be independently validated (Addo-Quaye et al., 2008; Li et al., 2010). To examine the differential accumulation of miRNAs that respond to B. elliptica, three libraries of sRNAs from L. regale under B. elliptica-mock treatment and elliptica-infection for 6 and 24 h were constructed, which were sequenced using the Solexa/Illumina platform with the transcriptome of L. regale. In addition, the corresponding target genes of these miRNAs were also identified by degradome sequencing in lily. The expression patterns of miRNAs under B. elliptica infection and their targets were identified and further verified by quantitative real-time PCR (qRT-PCR). The results suggest the presence of regulatory miRNAs in the economically important ornamental plant lily, and also shed new light on the miRNA-mediated regulatory networks that respond to fungal infection in lily.

Materials and Methods

Plants and B. elliptica Inoculation

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 × 104 conidia·mL−1 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.

RNA Extraction, Small RNA Library Construction and Sequencing

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).

Analysis of sRNA Sequencing Data

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 ( and the NCBI GenBank database ( 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 ( by BlastN to identify conserved miRNAs. Furthermore, the unannotated small RNAs were analyzed for the prediction of novel miRNAs.

Identification and Validation of Novel miRNAs of L. regale

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.

Identification of B. elliptica-Responsive miRNAs

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 = log2 (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).

Degradome Library Construction and Target Identification

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.

Quantitative Real-Time PCR (qPT-PCR) Analysis of miRNA Expression

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 ddH2O. 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.


High-Throughput sRNAs Sequencing in L. regale Wilson

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).


Table 1. Distribution of small RNAs among different categories in three libraries.

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).


Figure 1. Length distribution of small RNAs sequences in three libraries.

Identification of Known miRNAs in L. regale Wilson

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).


Figure 2. Conserved miRNA families in L. regale and across species. All miRNAs of L. regale were identified based on sRNA sequencing data, and those of other plants were obtained from miRBase 20.0.


Figure 3. Number and abundance of identified known miRNA families from lily (A: distribution of known miRNA family members in lily; B: counts of each known miRNA family in lily).

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).

Prediction of Novel miRNAs in L. regale Wilson

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−1 (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.


Figure 4. The hairpin structures for precursors of novel miRNAs. The red colored sequences represent mature miRNAs, and the yellow colored sequences represent the miRNA* (A represents miRn3, B represents miRn4, and C represents miRn5).

miRNAs Differentially Expressed in Response to B. elliptica

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, log2Ratio > 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.


Figure 5. Heatmaps show the differential expression of 24 conserved and 7 novel miRNAs of L. regale under B. elliptica infection 0, 6, and 24 h with p ≤ 0.01. Red indicates higher levels of miRNAs and green indicates lower levels. The absolute signal intensity ranges from –3.0 to +3.0, with corresponding color changes from green to red.

Identification of the B. elliptica-Responsive miRNA Targets by Using Degradome Analysis

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.


Figure 6. Target plots (t-plots) of miRNA targets identified by degradome sequencing in lily. Red arrows represent the nucleotide position of cleavage on the target genes. (A–D) Represents miR168a, miR5059, miR5142, and miR5072, respectively.


Table 2. Identified target transcripts for the known and novel miRNAs in lily.

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.

Validation of B. elliptica-Responsive miRNAs and Their Targets by qRT-PCR

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.


Figure 7. Analysis of the differential expression of miRNAs with high abundance by comparing the deep sequencing (right y-axis) and qRT-PCR (left y-axis) in L. regale. The fold-change of miRNAs was transformed as Log2 (treat/CK) from date of the deep sequencing and qRT-PCR.

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.


Figure 8. qRT-PCR validation of B. elliptica-responsive miRNAs in lily. The blue bar represents the resistant L. regale, and the red bar represents the susceptible “Tresor.” The amount of expression was normalized to the level of 18S rRNA. The normalized miRNA levels at 0 h were arbitrarily set to 1. Each bar shows the mean ± SE of triplicate assays. * or ** indicates a statistically significant difference as relative to the value at 0 h for each miRNAs at p < 0.05 or 0.01, respectively.

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.


Figure 9. qRT-PCR validation of B. elliptica-responsive miRNA target genes in L. regale (Conting 1715, Conting 39906, Conting 5555, Conting 6476, Conting 4725, and Conting 5180 represent genes encoding EF1, SODIUM/HYDROGEN EXCHANGER 6, HSP, HYPOTHETICAL PROTEIN 1, TCP2, and VESICLE-ASSOCIATED MEMBRANE PROTEIN). The amount of expression was normalized to the expression level of the CLATHRIN gene. The normalized mRNA levels at 0 h were arbitrarily set to 1. Each bar shows the mean ± SE of triplicate assays. * or ** indicates a statistically significant difference as relative to the value at 0 h for each miRNAs at p < 0.05 or 0.01, respectively.


miRNA Library Construction of L. regale Wilson

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.

Characteristics of B. elliptica-Responsive miRNAs 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.

Role of miRNAs-Target Gene Mediated Gene Expression Involved in Plant Defense

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.


In conclusion, to the best of our knowledge, this is the first study to perform transcriptome-wide identification of elliptica-responsive miRNAs and their targets using small RNA sequencing and degradome analysis in lily. A total of 71 miRNAs, belonging to 47 miRNA families, were identified in L. regale leaves subjected to mock B. elliptica treatment and B. elliptica infection. Among these, 31 miRNAs were differentially expressed in the elliptica-infected leaves and responsive to B. elliptica stress. Using degradome sequencing, 22 targets cleaved by miRNAs were confirmed in lily. These targets were functionally predicted to transcription factors, biosynthesis metabolic enzymes, biotic, and abiotic stress-responsive proteins. The expression patterns of these differentially regulated miRNAs and their targets were shown to be regulated by B. elliptica infection. These results could provide new information for the further identification and characterization of miRNAs in lily, and advance our understanding of the function of miRNAs and their targets in regulating plant responses to B. elliptica.

Author Contributions

XG, QL, and GJ designed experiments; XG, QC, and DZ carried out experiments and analyzed experimental results. XG and QzC wrote the manuscript. XG and HH modified the manuscript.


This work was supported by National High Technology Research and Development Program of China (863 Program) (Grant No. 2013AA102706) and the National Natural Science Foundation of China (Grant No. 31470106).

Conflict of Interest Statement

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.

Supplementary Material

The Supplementary Material for this article can be found online at:

Supplementary Figure S1. Stem-loop RT-PCR electrophoresis of novel miRNAs.

Supplementary Table S1. Abundance of conserved miRNAs from three libraries in lily.

Supplementary Table S2. Identification and expression analysis of novel miRNAs in lily.

Supplementary Table S3. Differential expression of B. elliptica-responsive miRNAs in lily.

Supplementary Table S4. Stem-loop primers of lily.

Supplementary Table S5. Specific qRT-PCR primers of lily.


Abreu, P. M. V., Gaspar, C. G., Buss, D. S., Ventura, J. A., Ferreira, P. C. G., and Fernandes, P. M. B. (2014). Carica papaya microRNAs are responsive to Papaya meleira virus infection. PLoS ONE 9:e103401. doi: 10.1371/journal.pone.0103401

PubMed Abstract | CrossRef Full Text | Google Scholar

Addo-Quaye, C., Eshoo, T. W., Bartel, D. P., and Axtell, M. J. (2008). Endogenous siRNA and miRNA targets identified by sequencing of the Arabidopsis degradome. Curr. Biol. 18, 758–762. doi: 10.1016/j.cub.2008.04.042

PubMed Abstract | CrossRef Full Text | Google Scholar

An, F., Liang, Y., Li, J., Chen, X., Han, H., and Li, F. (2014). Construction and significance analysis of the MicroRNA expression profile of Hemerocallis fulva at low temperature. Biosci. Biotech. Bioch. 78, 378–383. doi: 10.1080/09168451.2014.878214

PubMed Abstract | CrossRef Full Text | Google Scholar

An, F. M., and Tsair, C. M. (2012). Transcriptome-wide characterization of miRNA-directed and non-miRNA-directed endonucleolytic cleavage using Degradome analysis under low ambient temperature in Phalaenopsis aphrodite subsp. formosana. Plant Cell Physiol. 53, 1737–1750. doi: 10.1093/pcp/pcs118

PubMed Abstract | CrossRef Full Text | Google Scholar

Angelini, R. M. D. M., Rotolo, C., Masiello, M., Gerin, D., Pollastro, S., and Faretra, F. (2014). Occurrence of fungicide resistance in populations of Botryotinia fuckeliana (Botrytis cinerea) on table grape and strawberry in southern Italy. Pest Manag. Sci. 70, 1785–1796. doi: 10.1002/ps.3711

CrossRef Full Text | Google Scholar

Bari, R., Pant, B. D., Stitt, M., and Scheible, W. (2006). PHO2, microRNA399, and PHR1 define a phosphate-signaling pathway in plants. Plant Physiol. 141, 988–999. doi: 10.1104/pp.106.079707

PubMed Abstract | CrossRef Full Text | Google Scholar

Beauclair, L., Yu, A., and Bouché, N. (2010). MicroRNA-directed cleavage and translational repression of the copper chaperone for superoxide dismutase mRNA in Arabidopsis. Plant J. 62, 454–462. doi: 10.1111/j.1365-313X.2010.04162.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Carrington, J. C., and Ambros, V. (2003). Role of microRNAs in plant and animal development. Science 301, 336–338. doi: 10.1126/science.1085242

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, H., Jones, A. D., and Howe, G. A. (2006). Constitutive activation of the jasmonate signaling pathway enhances the production of secondary metabolites in tomato. FEBS Lett. 580, 2540–2546. doi: 10.1016/j.febslet.2006.03.070

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, J., Zheng, Y., Qin, L., Wang, Y., Chen, L., He, Y., et al. (2016). Identification of miRNAs and their targets through high-throughput sequencing and degradome analysis in male and female Asparagus officinalis. BMC Plant Biol. 16:80. doi: 10.1186/s12870-016-0770-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, L., Ren, Y., Zhang, Y., Xu, J., Zhang, Z., and Wang, Y. (2012). Genome-wide profiling of novel and conserved Populus microRNAs involved in pathogen stress response by deep sequencing. Planta 235, 873–883. doi: 10.1007/s00425-011-1548-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Debigaré, R., and Price, S. R. (2003). Proteolysis, the ubiquitin-proteasome system, and renal diseases. Am. J. Physiol. Renal Physiol. 285:F1. doi: 10.1152/ajprenal.00244.2002

PubMed Abstract | CrossRef Full Text | Google Scholar

Du, Y.-P., He, H.-B., Wang, Z.-X., Wei, C., Li, S., and Jia, G.-X. (2014). Investigation and evaluation of the genus Lilium resources native to China. Genet. Resour. Crop Evol. 61, 395–412. doi: 10.1007/s10722-013-0045-6

CrossRef Full Text | Google Scholar

Duan, H., Lu, X., Lian, C., An, Y., Xia, X., and Yin, W. (2016). Genome-wide analysis of microRNA responses to the phytohormone abscisic acid in Populus euphratica. Front. Plant Sci. 7:1184. doi: 10.3389/fpls.2016.01184

PubMed Abstract | CrossRef Full Text | Google Scholar

Ginzinger, D. G. (2002). Gene quantification using real-time quantitative PCR: an emerging technology hits the mainstream. Exp. Hematol. 30, 503–512. doi: 10.1016/S0301-472X(02)00806-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Glazov, E. A., Cottee, P. A., Barris, W. C., Moore, R. J., Dalrymple, B. P., and Tizard, M. L. (2008). A microRNA catalog of the developing chicken embryo identified by a deep sequencing approach. Genome Res. 18, 957–964. doi: 10.1101/gr.074740.107

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, R., Chen, X., Lin, Y., Xu, X., Thu, M. K., and Lai, Z. (2016). Identification of novel and conserved miRNAs in leaves of in vitro grown Citrus reticulata “Lugan” plantlets by solexa sequencing. Front. Plant Sci. 6:1212. doi: 10.3389/fpls.2015.01212

PubMed Abstract | CrossRef Full Text | Google Scholar

Jia, H., Zhang, C., Pervaiz, T., Zhao, P., Liu, Z., Wang, B., et al. (2016). Jasmonic acid involves in grape fruit ripening and resistant against Botrytis cinerea. Funct. Integr. Genomics 16, 79–94. doi: 10.1007/s10142-015-0468-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, D. H., Yin, C. S., and Yu, A. P. (2006). Duplication and expression analysis of multicopy miRNA gene family members in Arabidopsis and rice. Cell Res. 16, 507–518. doi: 10.1038/

PubMed Abstract | CrossRef Full Text | Google Scholar

Jin, W., and Wu, F. (2015). Characterization of miRNAs associated with Botrytis cinerea infection of tomato leaves. BMC Plant Biol. 15:1. doi: 10.1186/s12870-014-0410-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Jin, W., Wu, F., Xiao, L., Liang, G., Zhen, Y., Guo, Z., et al. (2012). Microarray-based analysis of tomato miRNA regulated by Botrytis cinerea. J. Plant Growth Regul. 31, 38–46. doi: 10.1007/s00344-011-9217-9

CrossRef Full Text | Google Scholar

Jones-Rhoades, M. W., Bartel, D. P., and Bartel, B. (2006). MicroRNAs and their regulatory roles in plants. Plant Biol. 57, 19–53. doi: 10.1146/annurev.arplant.57.032905.105218

PubMed Abstract | CrossRef Full Text | Google Scholar

Jover-Gil, S., Candela, H., and Ponce, M. (2005). Plant microRNAs and development. Int. J. Dev. Biol. 49, 733–744. doi: 10.1387/ijdb.052015sj

PubMed Abstract | CrossRef Full Text | Google Scholar

Kang, M., Zhao, Q., Zhu, D., and Yu, J. (2012). Characterization of microRNAs expression during maize seed development. Acta Trop. 13, 1–11. doi: 10.1186/1471-2164-13-360

PubMed Abstract | CrossRef Full Text | Google Scholar

Kohl, J., Molhoek, W. M. L., Van Der Plas, C. H., and Fokkema, N. J. (1995). Effect of Ulocladium atrum and other antagonists on sporulation of Botrytis cinerea on dead lily leaves exposed to field conditions. Phytopathology 85, 393–401. doi: 10.1094/Phyto-85-393

CrossRef Full Text

Kurubanjerdjit, N., Tsai, J. J. P., Huang, C., and Ng, K. (2014). Disturbance of Arabidopsis thaliana microRNA-regulated pathways by Xcc bacterial effector proteins. Amino Acids 46, 953–961. doi: 10.1007/s00726-013-1646-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Y., Zheng, Y., Addo-Quaye, C., Zhang, L., Saini, A., Jagadeeswaran, G., et al. (2010). Transcriptome-wide identification of microRNA targets in rice. Plant J. 62, 742–759. doi: 10.1111/j.1365-313X.2010.04187.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Lim, J. H., Rhee, H. K., Kim, Y. J., and Van Tuyl, J. M. (2003). Resistance to Fusarium oxysporum f.sp. lilii in Lilium. Acta Hortic. 620, 315–318. doi: 10.17660/actahortic.2003.620.38

CrossRef Full Text | Google Scholar

Liu, Q., Wei, C., Zhang, M., and Jia, G. (2016). Evaluation of putative reference genes for quantitative real-time PCR normalization in Lilium regale during development and under stress. PeerJ 4:e1837. doi: 10.7717/peerj.1837

PubMed Abstract | CrossRef Full Text | Google Scholar

Llave, C., Kasschau, K. D., Rector, M. A., and Carrington, J. C. (2002). Endogenous and silencing-associated small RNAs in plants. Plant Cell 14, 1605–1619. doi: 10.1105/tpc.003210

PubMed Abstract | CrossRef Full Text | Google Scholar

Long, W., Mai, Y. X., Zhang, Y. C., Qian, L., and Yang, H. Q. (2010). MicroRNA171c-targeted SCL6-II, SCL6-III, and SCL6-IV genes regulate shoot branching in Arabidopsis. Mol. Plant 3, 794–806. doi: 10.1093/mp/ssq042

CrossRef Full Text | Google Scholar

Lu, Y., Liu, Y., and Chen, C. (2007). Stomatal closure, callose deposition, and increase of LsGRP1-corresponding transcript in probenazole-induced resistance against Botrytis elliptica in lily. Plant Sci. 172, 913–919. doi: 10.1016/j.plantsci.2006.12.020

CrossRef Full Text | Google Scholar

Luan, Y., Cui, J., Zhai, J., Li, J., Han, L., and Meng, J. (2015). High-throughput sequencing reveals differential expression of miRNAs in tomato inoculated with Phytophthora infestans. Planta 241, 1405–1416. doi: 10.1007/s00425-015-2267-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Luo, M., Lin, H., Gao, J., Li, W., Shen, Y., Zhao, M., et al. (2015). Identification and characterization of differentially expressed microRNAs in response to Rhizoctonia solani in maize. Acta Physiol. Plant. 37, 1–17. doi: 10.1007/s11738-015-2001-x

CrossRef Full Text | Google Scholar

Luo, Y., Guo, Z., and Li, L. (2013). Evolutionary conservation of microRNA regulatory programs in plant flower development. Dev. Biol. 380, 133–144. doi: 10.1016/j.ydbio.2013.05.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Maksymiec, W., Wianowska, D., Dawidowicz, A. L., Radkiewicz, S., Mardarowicz, M., and Krupa, Z. (2005). The level of jasmonic acid in Arabidopsis thaliana and Phaseolus coccineus plants under heavy metal stress. J. Plant Physiol. 162, 1338–1346. doi: 10.1016/j.jplph.2005.01.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Man, M. Z., Wang, X., and Wang, Y. (2000). POWER_SAGE: comparing statistical tests for SAGE experiments. Bioinformatics 16, 953–959. doi: 10.1093/bioinformatics/16.11.953

PubMed Abstract | CrossRef Full Text | Google Scholar

McConn, M., Creelman, R. A., Bell, B., Mullet, J. E., and Browse, J. (1997). Jasmonate is essential for insect defense in Arabidopsis. Proc. Natl. Acad. Sci. U.S.A. 94, 5473–5477. doi: 10.1073/pnas.94.10.5473

PubMed Abstract | CrossRef Full Text | Google Scholar

Meyers, B. C., Axtell, M. J., Bartel, B., Bartel, D. P., Baulcombe, D., Bowman, J. L., et al. (2008). Criteria for annotation of plant MicroRNAs. Plant Cell 20, 3186–3190. doi: 10.1105/tpc.108.064311

PubMed Abstract | CrossRef Full Text | Google Scholar

Naqvi, A. R., Haq, Q. M. R., and Mukherjee, S. K. (2010). MicroRNA profiling of tomato leaf curl new delhi virus (tolcndv) infected tomato leaves indicates that deregulation of miR159/319 and miR172 might be linked with leaf curl disease. Virology J. 7, 1–16. doi: 10.1186/1743-422X-7-281

PubMed Abstract | CrossRef Full Text | Google Scholar

Niu, S., Fan, G., Deng, M., Zhao, Z., Xu, E., and Cao, L. (2016). Discovery of microRNAs and transcript targets related to witches' broom disease in Paulownia fortunei by high-throughput sequencing and degradome approach. Mol. Genet. Genomics 291, 181–191. doi: 10.1007/s00438-015-1102-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Njemini, R., Demanet, C., and Mets, T. (2004). Inflammatory status as an important determinant of heat shock protein 70 serum concentrations during aging. Biogerontology 5, 31–38. doi: 10.1023/B:BGEN.0000017684.15626.29

PubMed Abstract | CrossRef Full Text | Google Scholar

Rao, J., Liu, D., Zhang, N., He, H., Ge, F., and Chen, C. (2013). Identification of genes differentially expressed in a resistant reaction to Fusarium oxysporum in Lilium regale by SSH. Ieri Procedia 5, 95–101. doi: 10.1016/j.ieri.2013.11.076

CrossRef Full Text | Google Scholar

Schommer, C., Palatnik, J. F., Aggarwal, P., Chtelat, A., and Cubas, P. (2008). Control of jasmonate biosynthesis and senescence by miR319 targets. PLoS Biol. 6:e230. doi: 10.1371/journal.pbio.0060230

PubMed Abstract | CrossRef Full Text | Google Scholar

Schreiber, A. W., Shi, B., Huang, C., Langridge, P., and Baumann, U. (2011). Discovery of barley miRNAs through deep sequencing of short reads. BMC Genomics 12:129. doi: 10.1186/1471-2164-12-129

PubMed Abstract | CrossRef Full Text | Google Scholar

Shen, D., Suhrkamp, I., Wang, Y., Liu, S., Menkhaus, J., Verreet, J., et al. (2014). Identification and characterization of microRNAs in oilseed rape (Brassica napus) responsive to infection with the pathogenic fungus Verticillium longisporum using Brassica AA (Brassica rapa) and CC (Brassica oleracea) as reference genomes. New Phytol. 204, 577–594. doi: 10.1111/nph.12934

PubMed Abstract | CrossRef Full Text | Google Scholar

Shu, Y. J., Liu, Y., and Li, W. (2016). Genome-wide investigation of microRNAs and their targets in response to freezing stress in Medicago sativa L., based on high-throughput sequencing. G3 6, 755–765. doi: 10.1534/g3.115.025981

PubMed Abstract | CrossRef Full Text | Google Scholar

Sunkar, R., and Zhu, J. (2004). Novel and stress-regulated microRNAs and other small RNAs from Arabidopsis. Plant Cell 16, 2001–2019. doi: 10.1105/tpc.104.022830

PubMed Abstract | CrossRef Full Text | Google Scholar

Turner, M., Yu, O., and Subramanian, S. (2011). Genome organization and characteristics of soybean microRNAs. BMC Genomics 13:169. doi: 10.1186/1471-2164-13-169

PubMed Abstract | CrossRef Full Text | Google Scholar

Van Baarlen, P., Woltering, E. J., Staats, M., and Van Kan, J. A. (2007). Histochemical and genetic analysis of host and non-host interactions of Arabidopsis with three Botrytis species: an important role for cell death control. Mol. Plant Pathol. 8, 41–54. doi: 10.1111/j.1364-3703.2006.00367.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Wan, L., Wang, F., Guo, X., Lu, S., Qiu, Z., Zhao, Y., et al. (2012). Identification and characterization of small non-coding RNAs from Chinese fir by high throughput sequencing. BMC Plant Biol. 12:146. doi: 10.1186/1471-2229-12-146

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, R., Xu, L., Zhu, X., Zhai, L., Wang, Y., Yu, R., et al. (2015). Transcriptome-wide characterization of novel and heat-stress-responsive microRNAs in radish (Raphanus Sativus L.) using next-generation sequencing. Plant Mol. Biol. Rep. 33, 867–880. doi: 10.1007/s11105-014-0786-1

CrossRef Full Text | Google Scholar

Wang, T., Pan, H., Wang, J., Yang, W., Cheng, T., and Zhang, Q. (2014). Identification and profiling of novel and conserved microRNAs during the flower opening process in Prunus mume via deep sequencing. Mol. Genet. Genomics 289, 169–183. doi: 10.1007/s00438-013-0800-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, W., Lin, F., Chang, W., Lin, K., Huang, H., and Lin, N. (2009). MiRExpress: analyzing high-throughput sequencing data for profiling microRNA expression. BMC Bioinform. 10:328. doi: 10.1186/1471-2105-10-328

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, W., and Luan, Y. (2015). The advance of tomato disease-related microRNAs. Plant Cell Rep. 34, 1089–1097. doi: 10.1007/s00299-015-1782-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Y., Liu, W., Shen, H., Zhu, X., Zhai, L., Xu, L., et al. (2015). Identification of radish (Raphanus sativus L.) miRNAs and their target genes to explore miRNA-mediated regulatory networks in lead (Pb) stress responses by high-throughput sequencing and degradome analysis. Plant Mol. Biol. Rep. 33, 358–376. doi: 10.1007/s11105-014-0752-y

CrossRef Full Text | Google Scholar

Wang, Z. J., Huang, J. Q., Huang, Y. J., Li, Z., and Zheng, B. S. (2012). Discovery and profiling of novel and conserved microRNAs during flower development in Carya cathayensis via deep sequencing. Planta 236, 613–621. doi: 10.1007/s00425-012-1634-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Weiberg, A., Wang, M., Lin, F. M., Zhao, H., Zhang, Z., Kaloshian, I., et al. (2013). Fungal small RNAs suppress plant immunity by hijacking host RNA interference pathways. Science 342, 118–123. doi: 10.1126/science.1239705

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, G., Park, M. Y., Conway, S. R., Wang, J., Weigel, D., and Poethig, R. S. (2009). The sequential action of miR156 and miR172 regulates developmental timing in Arabidopsis. Cell 138, 750–759. doi: 10.1016/j.cell.2009.06.031

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, C., Li, D., Mao, D., Liu, X., Ji, C., Li, X., et al. (2013). Overexpression of microRNA319 impacts leaf morphogenesis and leads to enhanced cold tolerance in rice (Oryza sativa L.). Plant Cell Environ. 36, 2207–2218. doi: 10.1111/pce.12130

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, Y., Chen, X., Chen, J., Xu, H., Li, J., and Zhang, Z. (2011a). Identification of novel and conserved microRNAs in Rehmannia glutinosa L. by solexa sequencing. Plant Mol. Biol. Rep. 29, 986–996. doi: 10.1007/s11105-011-0293-6

CrossRef Full Text | Google Scholar

Yang, Y., Chen, X., Chen, J., Xu, H., Li, J., and Zhang, Z. (2011b). Differential miRNA expression in Rehmannia glutinosa plants subjected to continuous cropping. BMC Plant Biol. 11:53. doi: 10.1186/1471-2229-11-53

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, J., Zhang, S., Han, S., Wu, T., Li, X., Li, W., et al. (2012). Genome-wide identification of microRNAs in larch and stage-specific modulation of 11 conserved microRNAs and their targets during somatic embryogenesis. Planta 236, 647–657. doi: 10.1007/s00425-012-1643-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, W., Xie, Y., Xu, L., Wang, Y., Zhu, X., Wang, R., et al. (2016). Identification of microRNAs and their target genes explores miRNA-mediated regulatory network of cytoplasmic male sterility occurrence during anther development in radish (Raphanus sativus L.). Front. Plant Sci. 7:1054. doi: 10.3389/fpls.2016.01054

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, C., Xia, H., Cao, T., Yang, Y., Zhao, S., Hou, L., et al. (2015). Small RNA and degradome deep sequencing reveals peanut microRNA roles in response to pathogen infection. Plant Mol. Biol. Rep. 33, 1013–1029. doi: 10.1007/s11105-014-0806-1

CrossRef Full Text | Google Scholar

Zhao, D., Gong, S., Hao, Z., and Tao, J. (2015). Identification of miRNAs responsive to Botrytis cinerea in herbaceous peony (Paeonia lactiflora Pall.) by high-throughput sequencing. Genes (Basel). 6, 918–934. doi: 10.3390/genes6030918

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, Y., Thilmony, R., Bender, C. L., Schaller, A., He, S. Y., and Howe, G. A. (2003). Virulence systems of Pseudomonas syringae pv. tomato promote bacterial speck disease in tomato by targeting the jasmonate signaling pathway. Plant J. 36, 485–499. doi: 10.1046/j.1365-313X.2003.01895.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Keywords: lily (Lilium regale Wilson), Botrytis elliptica, miRNA, high-throughput sequencing, degradome

Citation: Gao X, Cui Q, Cao Q-Z, Liu Q, He H-B, Zhang D-M and Jia G-X (2017) Transcriptome-Wide Analysis of Botrytis elliptica Responsive microRNAs and Their Targets in Lilium Regale Wilson by High-Throughput Sequencing and Degradome Analysis. Front. Plant Sci. 8:753. doi: 10.3389/fpls.2017.00753

Received: 09 December 2016; Accepted: 21 April 2017;
Published: 18 May 2017.

Edited by:

Pierre-Emmanuel Courty, Institut National de la Recherche Agronomique (INRA), France

Reviewed by:

Zhongxiong Lai, Fujian Agriculture and Forestry University, China
Yao-Cheng Lin, Academia Sinica, Taiwan

Copyright © 2017 Gao, Cui, Cao, Liu, He, Zhang and Jia. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Gui-Xia Jia,