Table 1
High-throughput sequencing statistics of L. × formolongi sRNAs
Category
|
VJ_I
|
VJ_II
|
FI_I
|
FI_II
|
FD
|
Total
|
unique
|
Total
|
unique
|
Total
|
unique
|
Total
|
unique
|
Total
|
unique
|
Raw reads
|
14,240,103
|
4,165,620
|
13,012,678
|
5,070,705
|
25,100,402
|
4,417,437
|
24,430,380
|
5,493,407
|
10,937,263
|
2,363,859
|
3ADT & length filter
|
3,253,881
|
1,696,297
|
1,726,013
|
1,370,377
|
8,960,668
|
1,903,704
|
4,503,056
|
1,379,069
|
5,065,486
|
855,271
|
Junk reads
|
43,445
|
18,415
|
44,301
|
25,086
|
55,442c
|
15,965
|
59,104
|
21,952
|
31,103
|
10,933
|
Rfam
|
rRNA
|
625,844
|
17,010
|
103,024
|
7,038
|
830,461
|
18,275
|
152,864
|
8,762
|
135,278
|
6,417
|
tRNA
|
334,519
|
3,726
|
271,546
|
3,448
|
429,462
|
3,870
|
364,860
|
3,249
|
74,368
|
1,669c
|
snoRNA
|
4,261
|
380
|
5,215
|
436
|
20,391
|
926
|
19,562
|
890
|
17,995
|
843
|
snRNA
|
4,524
|
430
|
4,826
|
523
|
12,468
|
951
|
8,310
|
761
|
5,221
|
489
|
others
|
35,921
|
1,242
|
8,421
|
678
|
97,281
|
1,808
|
11,959
|
1,010
|
14,357
|
683
|
Total
|
1,005,069
|
22,788
|
393,032
|
12,123
|
1,390,063
|
25,830
|
557,554
|
14,672
|
247,219
|
10,101
|
Repeats
|
8,622
|
451
|
8,427
|
458
|
19,758
|
680
|
10,301
|
708
|
3,934
|
323
|
valid reads
|
9,932,250
|
2,427,801
|
10,843,854
|
3,662,791
|
14,681,011
|
2,471,404
|
19,303,311
|
4,077,107
|
5,590,925
|
1,487,303
|
Identification of conserved miRNAs in L. × formolongi
The valid reads from the five libraries were subjected to mapping with the precursors deposited in miRBase 22.1 and the transcriptome of L. × formolongi to identify conserved miRNAs. A total of 345 conserved miRNA precursors and 366 unique miRNAs belonging to 116 miRNA families were discovered, and the genomic localization was obtained for a majority of these miRNAs (Fig. 1b, Additional file 1: Table S1). In addition, the numbers of homologs were counted and varied among the different miRNA families. The largest family was the miR156 family, with 22 members, followed by the miR166 and miR169 families, which had 17 and 14 members, respectively. However, miR162, miR391, miR828 and miR5072 were each represented by only one member. Furthermore, two miR2916 members, one miR5523 member and one miR11602 member were identified as known miRNAs (group 1) whose reads and precursors could be perfectly mapped to the reported Lilium miRNAs/ pre-miRNAs sequences in miRBase and the transcriptome.
The relative abundance of these conserved miRNAs during lily development was characterized. As shown in Fig. 1c, the miRNAs exhibited remarkably different levels of abundance. The top read counts of miRNAs were obtained for the miR166 family, followed by the miR159 and miR408 families during developmental stages.
Identification of novel miRANs in L. × formolongi
The unannotated reads that were mapped to the transcriptome of L. × formolongi were filtered for novel miRNA prediction according to the criteria for the hairpin structure of the precursor. Collectively, 32 novel miRNAs were found among the five miRNA libraries (Additional file 1: Table S2). The lengths of these mature novel miRNAs and the hairpin lengths of the novel pre-miRNAs were between 19 and 25 nt and 55 and 191 nt, respectively. The majority of the miRNAs in these intervals were 21 nt and 164 nt long, respectively. The free energy of these secondary structures required to maintain stable hairpin structures ranged from -125.5 to -20.9 kCal/ mol. Precursors forming stem loop structures among the novel miRNAs, such as miRn5, miRn8 and miRn10, were predicted and are presented in Additional file 2: Fig. S1.
Differentially expressed miRNAs during development
To identify the key miRNAs playing regulatory roles during vegetative growth and flowering initiation in lily, 75 conserved miRNAs and 9 novel miRNAs were identified as showing significant differences in expression (Additional file 1: Table S3). Cluster analysis of their expression patterns was performed, and a heat map showing their distinct expression profiles was generated (Fig. 2). These miRNAs were characterized by three expression profiles with different expression changes from the VJ_ I to FI_II stages before FD.
The majority of miRNAs, including 50 conserved miRNAs and 9 novel miRNAs, presented high expression levels in the VJ (VJ_I and VJ_II) stage and were downregulated during later stages of development. In the miR156 family, most miR156s displayed consistently decreased expression, with the highest expression levels in the VJ_I or VJ_II stage. Similar to miR156s, most miR166 families exhibited expression peaks in the VJ_II stage. By contrast, 20 conserved miRNAs exhibited elevated expression from the VJ to FI_II stage or FD stage. miR159d, miR159f, miR172a and miR171 homologs exhibited upregulated expression levels from VJ_I to FI_II, and several less conserved miRNAs, such as miR5540, miR8708 and miR3630a/ b, maintained increased expression from the VJ to FD stages. In addition, many differentially expressed miRNAs with fluctuating transcription levels were found in the libraries. The expression of miR156k, miR2592bj and miR167e/ a was downregulated in the VJ_II and FI_II stages.
Degradome sequencing
By combining the results of degradome and computational target prediction for the annotated corresponding mRNAs, 396 target transcripts of 235 miRNAs were globally identified in the degradome and analysed by detecting mRNA fragments that were diagnostic of miRNA-directed cleavage, including 157 conserved miRNAs (Fig. 3, Additional file 1: Table S4) and 28 novel miRNAs (Additional file 1: Table S5). By distinguishing occupied transcript positions on the basis of abundance in the degradome, we classified the identified targets into 5 categories according to the criteria indicated in previous studies [28]. Among these identified targets, 24, 7, 179, 25 and 184 belonged to categories 0, 1, 2, 3 and 4, respectively.
In total, 6473 of the target sites were confirmed by degradome sequencing, and 4076 targets of 367 miRNA reads were simultaneously detected by computational target prediction and experimental sequencing. The enriched GO items were analysed to elucidate the functions of the identified targeted genes. Among biological processes, ‘oxidation-reduction process’ (GO: 0055114), ‘response to cold’ (GO: 0009409) and ‘response to cadmium ion’ (GO: 0046686) were the most abundant categories. ‘Nucleus’ (GO: 0005634) and ‘ATP binding’ (GO: 0005524) were the most highly represented cellular component and molecular function categories, respectively (Fig. 4a). The enriched KEGG pathway analyses showed that the pentose phosphate pathway (ko00030), carbon fixation in photosynthetic organisms (ko00710) and flavonoid biosynthesis (ko00941) were the most enriched pathways (Fig. 4b).
According to the functional annotations, putative flowering-related targets were identified. They consisted of squamosal promoter-binding, floral homeotic protein APETELA2-like, MADS-box, NAC domain-containing, R2R3-MYB, GAMYB and hormone signalling auxin response factor 18-like proteins (Table 2). We also found that the functions of LfSPL2 and LfNAM, targeted by miR156s and miR164a, were enriched in starch and sucrose metabolism, pentose and glucuronate interconversion pathways. Furthermore, both miR156s and miR159a were indicated to regulate the common target LfSPL2 according to the identification of different target sites (Fig. 5a). miR159a and miR319f targeted LfGAMYB sharing a common targeting site (Fig. 5b). These transcripts targeted by miRNAs were considered candidate key regulators during vegetative development and flowering induction in lily.
Table 2
Identification of putative flowering-related target transcripts in L. × formolongi
miR_name
|
Transcript
|
Transcript Annotation
|
Symbol
|
Term
|
lfo-miR156a,k,e,ac,d,g,y,r,h
|
Unigene18229_All
|
squamosa promoter-binding-like protein 12-like
|
LfSPL2
|
Starch and sucrose metabolism-
|
lfo-miR156a,e,ac,d,g,r,y
|
Unigene9969_All
|
predicted protein
|
LfSPL3
|
-
|
lfo-miR159a,f,c
|
Unigene2305_All
|
transcription factor GAMYB
|
LfGAMYB
|
developmental process
|
lfo-miR159a
|
Unigene18229_All
|
squamosa promoter-binding-like protein 12-like
|
LfSPL2
|
Starch and sucrose metabolism-
|
lfo-miR159e
|
CL730.Contig4_All
|
transcription factor TCP2-like
|
|
-
|
lfo-miR160a,f,g,b,d
|
CL6847.Contig1_All
|
auxin response factor 18-like
|
LfARF
|
root cap development; auxin mediated signaling
|
lfo-miR164a
|
CL6145.Contig1_All
|
NAC domain-containing protein
|
LfNAC
|
-
|
|
Unigene9281_All
|
no apical meristem-related protein
|
LfNAM
|
Starch and sucrose metabolism; Pentose and glucuronate interconversions;
|
lfo-miR172a
|
CL1418.Contig1_All
|
vesicle-associated membrane protein
|
|
response to salt stress
|
|
Unigene18548_All
|
floral homeotic protein APETALA 2-like
|
LfAP2
|
-
|
lfo-miR319f
|
Unigene2305_All
|
transcription factor GAMYB
|
LfGAMYB
|
developmental process
|
lfo-miR394
|
Unigene14281_All
|
F-box family protein
|
|
regulation of auxin mediated signaling
|
lfo-miR828
|
Unigene3673_All
|
transcription factor R2R3-MYB
|
|
|
lfo-miR5179
|
Unigene7245_All
|
MADS-box transcription factor
|
LfAGL15
|
regulation of transcription, DNA-dependent
|
Correlation analysis of miRNA expression profiles and miRNA target genes
To investigate the crucial participants and their variable trends during lily vegetative development and flowering initiation, the 28 detected target transcripts of significantly differentially expressed miRNAs obtained from the degradome were applied for the analysis of differentially expressed genes (DEGs). There were 15 DEGs corresponding to 23 differentially expressed miRNAs in the transcriptome. A heatmap was generated to show the reverse expression profiles of 23 miRNA-target pairs in the transcriptome and miRNA libraries from the VJ (vegetative growth) to FD stages (Fig. 6). Among these pairs, the members of the miR156 family were most highly expressed and displayed significantly downregulated expression levels, while the expression of their target LfSPL2 was increased during development. Putative flowering-related miRNA-target pairs such as miR172a-LfAP2, miR156s/ miR159a-LfSPL2 and miR164a-LfNAC were also identified in the correlation analysis.
Based on the results of the correlation analysis and the functional analysis of targets, 12 miRNA-target pairs that showed distinct expression patterns and may take part in flowering induction were selected to verify the accuracy of the sequencing results by performing qRT-PCR analysis. As depicted in Fig. 7a-j, these 10 candidate miRNAs were indicated to show substantially similar expression changes to those indicated by the data from miRNA libraries. A portion of these miRNAs were inversely correlated with their corresponding targets, such as miR156a-LfSPL2/ LfSPL3, miR159a-LfGAMYB, miR397a-LfLACs, miR160a-LfARF18, miR5139-LfFAO2, and miR5179-LfAGL15, which confirmed the negative relationship in the expression profile of miRNA-mRNA modules. Moreover, the expression profiles of the miRNAs from the VJ_I to FD stages were diverse. miR160a was gradually upregulated after the VJ_I stage, while four miRNAs, including miR159a, miR172a, miR167a and miR171b, exhibited upregulated expression in the FI_I stage. Both miR5139b and miR5179 presented expression peaks in the VJ_II stage. Moreover, all of these miRNAs exhibited a remarkable upregulation of transcript abundance from the FD to VB stages, with the exception of miR159f, which exhibited decreased expression at the FD stage.
In addition, we assayed the tissue-specific expression profiles of two significantly differentially expressed miRNAs, miR159a and miR160a (Fig. 7k, l). They showed similar expression patterns in different organs and tissues from the VJ to VB stages. Both presented high expression in roots and stems in the flowering induction (FI_I, FI_II) stages. Expression peaks of miR159a and miR160a were found in leaves at the FI_I stage. However, the level of miR159a increased from the VJ to FD stages in the apical meristem and bulb scales. In visible floral buds, the strikingly highest expression of mR159a was observed. Likewise, miR160a showed the highest expression level at the FD stage in the apical meristem.
Gene coexpression network analysis
To mine the highly correlated miRNA clusters that were relevant to flowering induction, we used miRNA data obtained from high-throughput sequencing to conduct correlation network analysis by WGCNA. After testing for abnormal samples or genes, correlation networks were constructed by setting the soft threshold to 12 (Additional file 4: Fig. S3-6). As a result, 5 miRNA coexpression modules comprising 70 (blue module), 40 (brown module), 53 (green module), 82 (turquoise module) and 33 (yellow module) interconnected miRNAs were found, as shown in Additional file 1: Table S7. We also recognized 14 miRNAs belonging to differentially expressed miRNA-target pairs in the blue module (miR169b, miR156a/ h/ e/ y), brown module (miR164a, miR172a, miR2592bj), green module (miR160a), turquoise module (miR156d/ g/ ac) and yellow module (miR398a, miR397a). The greatest number of these miRNAs belonging to differentially expressed pairs was found most in the blue module; in this regard, the blue module could be considered the key module related to critical stages of vegetative development and flowering initiation in lily.
We also extracted a subnetwork containing 795 edges connected with the differentially expressed miRNAs, which also included their targets (Additional file 1: Table S8). In the regulatory network, miR156a/ d/ e exhibited the greatest numbers of interactions with other sequences, suggesting their core roles among the coexpressed genes (Fig. 8). Moreover, phylogenetic analysis demonstrated that miR156s were most closely related to Arabidopsis. They are likely candidate miRNAs involved in lily flowering induction, and the specific roles of the members of this family require further validation and investigation (Additional file 3: Fig. S2).