Skip to main content
Download PDF
  • Research article
  • Open access
  • Published:

Analysis of miRNA-mediated regulation of flowering induction in Lilium × formolongi

BMC Plant Biology volume 21, Article number: 190 (2021) Cite this article

Abstract

Background

MicroRNAs play pivotal roles in plant vegetative phase change and flowering induction via integrating into multiple flowering pathways. Lilium × formolongi is an important ornamental lily cultivar that can flower within one year after sowing. However, it remains unresolved how miRNA-mediated regulation networks contribute to the L. × formolongi characteristics of a short vegetative growth period and rapid flowering.

Results

In this study, the small RNA libraries and one degradome library were constructed for L. × formolongi during vegetative growth and flowering initiation, and 366 conserved miRNAs and 32 novel miRNAs were identified. Additionally, 84 miRNAs were significantly differentially expressed during development. A total of 396 targets of 185 miRNAs were identified and validated through degradome sequencing. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses showed that functions of the targets were top enriched in the cold and cadmium ion responses, pentose phosphate pathway and carbon fixation in photosynthetic organisms. Furthermore, among 23 differentially expressed miRNA-target pairs, the miR156s-LfSPL2, miR172a-LfAP2 and miR164a-LfNAC pairs as well as miR159a-LfSPL2 were found to be relevant to flowering based on the correlation analysis of expression profiles in the miRNA libraries, degradome and transcriptome. A coexpression regulatory network focused on differentially expressed pairs was also constructed by WGCNA, and 14 miRNAs were considered putative key miRNAs during vegetative development and flowering induction. miR156a/ d/ e showed particularly strong relationships with other miRNAs in the coexpression network.

Conclusions

This study provides cues for the further exploration of the regulatory mechanisms of short vegetative development and flowering in L. × formolongi.

Background

Lilium (Liliaceae family) is an important genus including crops that are valued for their ornamental characteristics, edible bulbs and medical use. They can be propagated by either sexual or asexual propagation. Sexual propagation through seed reproduction enables breeders to obtain virus-free seedlings and maintain abundant genetic diversity [1]. In most lily cultivars and wild species, at least 1 year of vegetative growth is required for flowering after sowing [2]. However, L. × formolongi, a popular cut flower cultivar with the characteristic of a shortened vegetative phase, can undergo floral transition within one year after germination [3]. Therefore, the investigation of molecular pathways involved in L. × formolongi vegetative development and flowering initiation is necessary to shorten the vegetative growth period and advance the flowering time of lily.

Vegetative growth after germination is necessary for angiosperm plants to achieve competence to flower and respond to external stimuli after the juvenile-adult stage transition, followed by the transition to reproductive growth [4]. As such, the optimal timing of the transition from the vegetative to the reproductive phase and flowering is precisely controlled by a complex genetic network composed of multiple pathways. In the annual plant Arabidopsis thaliana, external factors such as day length and temperature affect the photoperiod and vernalization pathways [5]. They control phase transitions and flowering time in coordination with ageing, gibberellic acid (GA), trehalose 6-phosphate (T6P) and autonomous pathways governed by endogenous factors including age, hormones and carbohydrate status [5, 6]. Currently, a large number of studies suggest that miRNAs are also integrated into flowering pathways, in which they perform orchestrated functions [7].

Mature plant miRNAs exert regulatory effects primarily by directly cleaving mRNAs or repressing the translation of specific mRNAs [8, 9]. In Arabidopsis, several endogenous conserved miRNAs and their targets are integrated into the ageing pathway, in which they influence development phase changes and flowering time [10, 11]. miR156, the orchestrator of the pathway, plays a master regulatory role in vegetative development via transcript cleavage or the translational repression of target SQUAMOSA PROMOTER BINDING PROTEIN-LIKE (SPL/ SBP) genes, resulting in the transition to the adult stage [12, 13]. Conversely, miR172 targets six APETALA 2 (AP2)-like genes that function in flowering repression, and the abundance of miR172 increases with age [14, 15]. In apices, miR172 acts downstream of the miR156/ SPL module, and SPL9 and SPL15 can alter the expression of miR172 by directly binding to the miR172 precursor MIR172b [16]. In addition to the two components of the ageing pathway, the examination of the phenotypes of mutant and transgenic plants has demonstrated that miR159, which is repressed by the DELLA protein in the GA pathway, is involved in the vegetative phase transition, independent of GA, and represses the expression of miR156 [17, 18]. In tomato, miR156-SlSBPs interact with GA and the miR319-LANCEOLATE (LA) module to affect flowering induction and the floral transition [19]. Additionally, overexpressed miR171 can activate miR156 to influence the vegetative-to-reproductive phase change in barley [20].

Many species-specific or less conserved miRNAs have also been identified and further confirmed to be involved in flowering induction in plants. Monocot-specific miRNA miR528 has been found to promote rice flowering induction under long-day conditions [21]. Brassicaceae-specific miRNA miR824 and its target AGAMOUS-LIKE 16 (AGL16) modulate flowering time by repressing the expression of FT under long-day conditions in Arabidopsis [22]. Pooideae-specific miR5200 has been shown to regulate the photoperiod-mediated flowering transition by guiding the sequence-specific cleavage of FT orthologues in Pooideae plants [23].

As a geophyte with underground storage organs, previous studies have confirmed that the flowering time of lily can be affected by multiple factors, including vernalization, photoperiod and bulb size [24, 25]. However, the molecular regulatory factors underlying the vegetative phase transition and responses to various floral induction signals in lily are still unknown. On the other hand, the lily genome is large, complex and lacks a reference genome sequence [26]. Nevertheless, sequencing technologies make it relatively easy to produce sufficient data for the identification of genes related to important agronomic traits, and this research could benefit from computational and degradome analyses that detect the targets of miRNAs rapidly and accurately [27, 28].

Therefore, 15 miRNA libraries and one degradome library were constructed for L. × formolongi at five key developmental stages. Databases, in conjunction with transcriptomic data, were further used to systematically identify and profile miRNAs and their targets during vegetative development and flowering initiation. Potential crucial flowering miRNA-target pairs were also identified. These results may contribute to elucidating the regulatory mechanism of miRNA-mediated short vegetative development in L. × formolongi and accelerating lily flowering.

Results

High-throughput sRNA sequencing in L. × formolongi

To explore comprehensive miRNA data during L. × formolongi development, 15 small RNA libraries were constructed and sequenced from leaf and stem meristem tissues collected at five developmental stages (VJ_I, VJ_II, FI_I, FI_II and FD) (Additional file 1, Fig. S1, S2). In total, 87,720,826 raw reads were generated from these sRNA libraries (Table 1). After removing junk sequences, poor-quality reads and noncoding RNAs, approximately 9.9 million, 10.8 million, 14.7 million, 19.3 million and 5.6 million valid reads that exactly matched the precursors deposited in miRBase 22.0 were obtained from the VJ_I, VJ_II, FI_I, FI_II and FD libraries, respectively. These miRNAs showed unique size distributions. The length distribution of 60,351,351 raw reads and 14,126,406 unique reads was centred on the 21–24 nt range. The most frequent length of valid reads among the libraries was 24 nt, followed by 21 nt for total reads and 22 nt for unique reads (Fig. 1a).

Table 1 High-throughput sequencing statistics of L. × formolongi sRNAs
Full size table
Fig. 1
figure 1

Length distribution of small RNAs sequences, identification of known miRNA families from L. × formolongi. a Length distribution of total reads and unique reads in five libraries. b Distribution of known miRNA family members. c Counts of each known miRNA family

Full size image

Identification of conserved miRNAs in L. × formolongi

The valid reads from all 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 transcriptome localization was obtained for a majority of these miRNAs (Fig. 1b, Additional file 2: Table S1). In addition, the number of homologs was 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 miRNAs 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 all miRNA libraries (Additional file 2: 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 3: Fig. S3.

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 (Additional file 2: 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.

Fig. 2
figure 2

Heat map showing the expression patterns of differentially expressed miRNAs. The miRNAs were clustered by hierarchical clustering at p ≤ 0.05 according to their expression patterns during development in L. × formolongi

Full size image

The majority of miRNAs, including 50 conserved miRNAs and 9 novel miRNAs, presented high expression levels at 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 at the VJ_I or VJ_II stage. Similar to miR156s, most miR166 families exhibited expression peaks at 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 at 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 2: Table S4) and 28 novel miRNAs (Additional file 2: Table S5). According to the cleavage signature abundance at each occupied transcript position, we classified the identified degradome sequences into 5 categories based on the criteria indicated in previous studies [29] (Additional File 2: Table S4). Among these identified targets, 24, 7, 179, 25 and 184 belonged to categories 0, 1, 2, 3 and 4, respectively.

Fig. 3
figure 3

Target plots of miRNA targets identified by degradome sequencing in lily. a-d represent miR5179, miR160a, miR164a and miR172a, respectively. Arrows represent the nucleotide position of cleavage in the target genes

Full size image

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

Fig. 4
figure 4

Gene functional classification of identified target genes. a Gene ontology enrichment and classification of target genes. b KEGG pathway analysis

Full size image

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
Full size table
Fig. 5
figure 5

mRNAs simultaneously targeted by miRNAs belonging to multiple families. a Diagram of LfSPL2 and LfGAMYB. Conserved domains are included, and the boxes represent the miRNA target sites. b Alignment of the miRNA targeting sites of miR156s-LfSPL2, miR159a-LfSPL2 and miR159a/miR319f-LfGAMYB

Full size image

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 differentially expressed genes (DEGs) analysis. There were 15 DEGs corresponding to 23 differentially expressed miRNAs in the transcriptome from the VJ (vegetative growth) to FD stages (Fig. 6). Among them, 13 miRNA-mRNA pairs displayed negative correlations at expression level. The members of the miR156 family were most highly expressed and displayed significantly downregulated expression levels, while the expression of their target LfSPL2 increased during development.

Fig. 6
figure 6

Heat map showing the expression patterns of differentially expressed miRNA-target. These miRNAs targeting differentially expressed mRNAs were clustered by hierarchical clustering at p ≤ 0.05 according to their expression patterns. The black-framed portion represents the expression level of miRNA and mRNA that showed negative correlation relationship at the expression level

Full size image

Based on the results of the correlation and functional analyses of the 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 RT-qPCR analysis. As depicted in Fig. 7a-j, these 10 candidate miRNAs showed substantially similar expression changes to those indicated by the data from the 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, miR5139b-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 at the FI_I stage. Both miR5139b and miR5179 presented expression peaks at the VJ_II stage. All these miRNAs exhibited a remarkable upregulation of transcript abundance from the FD to VB stages, except for miR159f, which exhibited decreased expression at the FD stage.

Fig. 7
figure 7

The expression patterns of miRNAs determined by RT-qPCR during development in L. × formolongi. a-j The expression patterns of miRNAs and their targets. k-l The tissue-specific expression profiles of miR159a and miR160a during lily development processes. The normalized miRNA and target levels at the VJ_I stage were arbitrarily set to 1. Each bar represents the mean ± SE of triplicate assays. * or ** indicates a statistically significant difference relative to the value at VJ_ I for each miRNA at p < 0.05 or 0.01, respectively

Full size image

We also 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 at 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. S4–7). 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 2: 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 2: 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 5: Fig. S8).

Fig. 8
figure 8

The coexpression subnetwork of putative crucial flowering-related miRNAs and their targets. The coexpression subnetwork of putative crucial flowering-related miRNAs and their targets in lily. In the nodes, circles and diamonds represent miRNAs and target transcripts, respectively. Different colours of nodes represent different modules identified by WGCNA

Full size image

Discussion

Potential crucial miRNAs in vegetative development and flowering induction

The role of miRNAs in flowering induction has been confirmed in different plants, such as Arabidopsis, rice, maize, apple and tomato [19, 30, 31]. The mechanism of miRNA functions in the genetic network underlying lily vegetative development and flowering initiation is poorly understood.

In our study, 84 miRNAs showing significant differential expression were identified, and analysis of their expression dynamics during lily development was conducted. The greatest numbers of differentially expressed miRNAs were identified and strongly expressed at the VJ_I and VJ_II stages among the small RNA libraries. These results were consistent with our hypothesis, and we reason that miRNA-mediated modulation is crucial for vegetative development by integrating environmental and endogenous factors to prepare for the induction of flowering.

In the L. × formolongi libraries, miRNAs belonging to the miR156 family were the most frequent among the differentially expressed miRNAs. In Arabidopsis, miR156 is an important orchestrator of the ageing pathway that controls the vegetative-to-reproductive transition and shows very high expression in the vegetative phase [32]. The conserved expression profile of miR156 members has been confirmed in many plants, such as rice [33], maize [34], tobacco [35] and the perennial woody species Acacia confuse, Eucalyptus globulus, and Quercus acutissima [36]. In L. × formolongi, the miR156 families shared similar expression profiles with those observed in other plants during development, demonstrating their evolutionary and functional conservation and their dominant role during vegetative development. Nevertheless, the expression level of miR156k fluctuated and was upregulated at the FI_I stage. In Arabidopsis and rice, Ath-MIR156A/ C and OsmiR156b/ h play master roles among miR156 precursors [37,38,39]. We observed that the precursors of miR156s (MIR156k) were phylogenetically closely related to Arabidopsis sequences and that MIR156k exhibited the closest relationships with Ath-MIR156A/ C, which clustered into the same branch (Additional file 5: Fig. S8), indicating that miR156k precursors may have a greater effect on the abundance of miR156 during lily flowering induction than during the juvenile-adult phase transition.

A prior study in Arabidopsis confirmed that another crucial regulator, miR172, exhibits an inverse expression pattern with miR156 and shows consistently upregulated expression during development to regulate the phase transition as well as floral organ identity [40, 41]. There are five miR172 family members and miR172b could be negatively regulated by miR156e to influence the juvenile-adult phase change [32]. In the perennial species Jatropha curcas, only two miR172s have been identified [42]. However, ten miR172s and only one miR172 (miR172a) were identified among the differentially expressed miRNAs in our libraries. The expression of miR172a was elevated from the VJ to FI_II stages, followed by a decline at the FD stage. Hence, we consider miR172a to exhibit relatively complex or species-specific functions in lily flowering induction or other aspects of development.

In Arabidopsis and rice, miR159a, miR167a and miR160a affect the development of anthers and sepals [43,44,45]. In lily, these highly conserved miRNAs also showed significantly upregulated expression at the VB stage, suggesting their important roles in floral organ formation during early floral development. In addition, miR5179 affects the diversification of perianth organs by cleaving MADS-box genes in orchids [46]. In contrast, miR5179 presented upregulated expression at the VB stage but exhibited strikingly higher expression at the VJ stage, indicating that this less conserved miRNA, which is rarely observed in other plants, may exhibit a more crucial function during the vegetative development of lily than during flower development.

Flowering-related miRNA-mRNA modules in L. × formolongi

miRNAs carry out their biological functions by regulating the expression of specific mRNAs. We identified targets that encoded flowering-related transcription factors or proteins. SPL genes are targeted by miR156 and regulate many aspects of plant vegetative development as important components of multiple pathways [47]. In Arabidopsis, 10 SPL genes are regulated by miR156, and AtSPL3/ 14/ 15/ 9 can promote flowering [4, 48]. SPL7 and SPL8 participate in flowering modulation and are targeted by miR156 in grasses [49]. Among the 11 LfSPL members, two homologs, LfSPL2 and LfSPL3, were targeted by miR156s in L. × formolongi, and only miR156s-LfSPL2 pairs were differentially expressed. Moreover, miR156a/ d/ e were located at the cores of the coexpression networks. Accordingly, miR156a/ d/ e-LfSPL2 modules may be the most critical regulators of the ageing pathway during vegetative development and flowering induction. The mechanism by which the miR156-LfSPL module is involved in the genetic network underlying flowering requires further intensive investigation and verification.

Previous studies in Arabidopsis have confirmed that miR159 regulates the juvenile-to-adult phase transition through a canonical miR159-GAMYB/ MYB33 interaction [50], and MYB33 could also activate miR156 as well as the target SPL9 [17]. Interestingly, we found that miR159a regulated LfSPL2 simultaneously with miR156s, in addition to LfGAMYB in this study. miR159a exhibited high expression levels at the FI_I stage, and the expression level of miR159a and LfSPL2 exhibited an approximately negative relationship from FI_I to VB stage. Hence, the results imply that miR159 may mainly promote the juvenile-to-adult phase change and flowering induction via the miR159-LfGAMYB module. Meanwhile, miR159a might be involved in a new flowering pathway through which it affects flowering initiation by regulating the SPL gene directly, and this regulatory module is worthy of validation in the future.

Crucial factors affecting the vegetative development of lily

A recent study conducted in our laboratory suggested that when the seedlings of L. × formolongi undergo bolting, the seedlings transition into the reproductive stage, and the juvenile-adult transition is accomplished [3]. Therefore, we focused on the miRNA-mediated regulatory mechanism during flowering induction, particularly during vegetative development at the rosette stage of L. × formolongi, and the factors affecting this process.

In the molecular regulatory network, sugar acts upstream of miR156 to alter miR156 expression [37]. T6P influences flowering in the leaves and apical meristem by regulating FT as a representative of carbohydrate conditions [6]. Previous studies have shown that daylength strongly affects the bolting time and quality of L. × formolongi [25]. Temperature also plays an important role in vegetative development. Cold exposure induces flowering by breaking bulb dormancy and influences the growth of bulbous plant species [51]. Our study showed that the functions of the identified targets were enriched in cadmium ion, cold response and sugar metabolism pathways, including the pentose phosphate pathway and carbon fixation in photosynthetic organisms. Thus, vegetative development may be closely associated with carbohydrates, cold or photoperiod in L. × formolongi via miRNAs. LfSPL2 and LfNAM, which are targeted by miR156s and miR164a, respectively, were identified in the sugar metabolism pathway. Hence, the impacts of these factors on the vegetative development of L. × formolongi through the potential miRNA-mRNA modules is also an important topic awaiting future exploration.

Taking the flowering regulation network in Arabidopsis as a reference, the hypothetical schematic diagram was based on the analysis of the small RNA libraries, degradome and transcriptome data (Fig. 9). In Arabidopsis, miR156 and SPLs are master regulators in vegetative development phase transition. SPL proteins can induce floral integrators, SUPPRESSOR OF OVEREXPRESSION OF CONSTANS 1 (SOC1) and AGL42 in leaves and apical meristem [13]. In apical meristem, SPLs can also regulate the floral meristem identity genes LEAFY (LFY) and APETALA 1 (AP1), which are repressed by AP2 transcription factor genes, the target of miR172 [13, 15]. Meanwhile, SPL9 directly promotes the transcription of miR172 to affect vegetative-to-reproductive phase transition [13]. Additionally, ARF 3/4 may modulate the expression of SPL3 via regulating the juvenile-adult phase transition [52].

Fig. 9
figure 9

The putative miRNA-mediated molecular regulation network of flowering induction and development of L. × formolongi. The locations of the genes are determined by their expression patterns. Genes in green background displayed expression peak at the VJ stage, genes in the yellow and blue background displayed up-regulated expression at the FI stage and FD stage respectively. The numbers of differentially expressed genes are also presented. Arrows in dotted lines represents the potential regulatory relationship between two genes

Full size image

Similar to Arabidopsis, the miRNA-target pair miR156-SPL occupies a central position in the putative miRNA-mediated flowering regulation network in L. × formolongi. The decreased expression of miR156 induces LfSPL2 and LfSPL3 expression at the FI_II stage, which may promote the MAD-box gene LfLFY to induce flowering in apical meristem. Meanwhile, miR160s could target LfARF18, and up-regulated LfARF18 at the VJ_II stage might promote the expression of LfSPL2 and LfSPL3. miR159a might be involved in flowering induction and floral meristem formation via repressing the expression of LfGAMYB and LfSPL2, respectively. Additionally, LfGAMYB may also induce the regulators miR156s and LfSPLs in the ageing pathway. LfSPL2 and LfSPL3 may also activate miR172a transcription. miR172a targets LfAP2, and the down-regulated expression of LfAP2 before VB stage may repress LfAP1 and LfSOC1. In addition, miR5179 specifically target MADS-box gene LfAGL15, which may play critical roles at flowering induction stage.

Conclusions

To our knowledge, the present study provides the first integrated transcriptomic, miRNA library and degradome analysis conducted to obtain comprehensive knowledge of miRNAs and their targets during flowering initiation in lily. Crucial miRNAs and miRNA-target pairs involved in vegetative development and floral induction in L. × formolongi were identified via the correlation analysis of sequencing datasets. Factors including carbohydrates, cold and photoperiod and associated enriched target functions may affect vegetative development and flowering induction. These results will help us to gain broader knowledge of the miRNA-mediated mechanism modulating L. × formolongi vegetative development and flowering and will be of great importance for regulating growth and development processes in lily.

Methods

Plant materials and samples collection

The seeds of L. × formolongi (‘Razan No. 2’) (purchased from Murakami Seed, Ibaraki, Japan) were sown after stratification at 4 °C for one month. Seedlings with 2–3 rosette leaves were transplanted into individual pots and grown in a chamber under a 16 h light (25 °C)/ 8 h dark cycle (16 °C) with 70% humidity, illuminated by white light with a 320 μmol·m− 2·s− 1 photosynthetic photon flux density.

A previous study in our laboratory confirmed the key characteristics of vegetative development and flowering initiation in L. × formolongi under long-day conditions through morphological and gene expression pattern investigations [3]. Therefore, fresh leaves and apical stems were collected for small RNA library construction during different developmental stages, including the 2–3 rosette leaf stage (vegetative juvenile_I, VJ_I), 5–6 rosette leaf stage (vegetative juvenile_II, VJ_II), 1–2 internode stage (flowering induction_I, FI_I), 5–6 internode stage (flowering induction_II, FI_II) and 9–10 internode stage (floral differentiation, FD). Similar samples from these stages together with fresh leaves and buds collected from 70 to 71 internodes (visible floral bud stage, VB) were prepared for RT-qPCR analysis. To examine the tissue-specific expression patterns of the differentially expressed miRNAs miR159a and miR160a, fresh leaves, apical meristems, middle stems, bulb scales (in the middle of the bulb), roots and tiny floral bubs (approximately 2 cm) were also collected from one plant from the VJ_I to VB stages. Three biological replicates were collected for each developmental harvesting point. All samples were immediately frozen in liquid nitrogen, followed by storage at − 80 °C for further study.

Small RNA library construction and sequencing

Total RNA was extracted using the EASY spin Plus Plant RNA kit (Aidlab Bio, Beijing, China) according to the manufacturer’s instructions. The high-quality total RNA was used for small RNA libraries generation following the TruSeq Small RNA Sample Prep Kits protocol (Illumina, San Diego, USA). Briefly, after ligating the RNA with 3′ and 5′ adapters, RNAs were used for reverse transcription and PCR amplification. Then gel electrophoresis was performed and gel-purified small RNA libraries of L. × formolongi at different development and flowering initiation stages were constructed followed by quality control analysis. Sequencing was conducted with Illumina HiSeq 2500 platform (LC Bio, Hangzhou, China). Each reaction was performed in triplicate.

Bioinformatic analysis of sRNA data

The sequencing analysis were performed through an in-house program, ACGT101-miR (LC Bio, Hangzhou, China). Briefly, to obtain valid reads without a reference for the L. × formolongi, sequence adapters, low-quality reads and junk reads containing poly-N were filtered out. Subsequently, the filtered reads that matched noncoding RNAs (rRNAs, tRNAs, snRNAs and snoRNA) and repeats in RFam (http://rfam.xfam.org) and Repbase (http://www.girinst.org/education/index.html) databases were discarded. The remaining 18–25 nt valid reads were further processed to summarize the length distribution, and mapping to specific species precursors in miRBase 22.0 via BLAST searches and transcriptome (data not published) analysis of L. × formolongi were performed to identify known miRNAs. The miRNAs are categorized into 3 groups, including known miRNAs (perfectly matched sequences, group 1) and conserved miRNAs (group 2 and group 3), according to whether the reads and precursors could be matched to sequences in miRBase and the transcriptome. Among the conserved miRNAs, the reads could be aligned to precursors in miRbase and the mapped precursors failed to be mapped to the transcriptome. When the reads could also be mapped to the transcriptome, the reads were categorized in group 2, if not, they were categorized in group 3.

After extracting the known miRNAs, the valid reads that were unmapped in miRBase but were subjected to BLAST searches against the transcriptome of L. × formolongi were selected for novel miRNA identification. RNAfold software (http://rna.tbi.univie.ac.at/cgibin/RNAWebSuite/RNAfold.cgi) was employed to predict the secondary structures of potential miRNA sequences. The detailed criteria applied were as follows: (1) numbers of nucleotides in a single bulge in the stem and the mature region were ≤ 12 and 4, respectively; (2) in a predicted hairpin, number of base pairs in the stem region was ≥16; (3) in the mature region, both numbers of biased errors and bulges in one bulge and the number of errors were ≤ 2 [53].

Analysis of differentially expressed miRNAs

To identify the differentially expressed miRNAs in L. × formolongi during the flowering transition phase and development, the deep-sequencing counts were normalized as described in prior reports [54]. In brief, first, a reference data set made up of the copy number median of selected common sequences of all samples was constructed. Second, the copy numbers of total samples from all libraries were transformed into Log2 (copy#) followed by calculation of differences between individuals and reference data set. Third, the linear regressions between the samples and reference set on the subset sequences (|△Log2 (copy#)| < 2), and the arithmetic correction factor were generated to calculate the expected Log2 (copy#). Thus, we calculated the expected Log2 (copy#) values to indicate the original copy numbers of the samples. The differentially expressed miRNAs were identified by one-way ANOVA and p-values were obtained via multiple comparisons between the libraries. The miRNAs with p-values ≤0.05 in each test were defined as significantly differentially expressed miRNAs during lily development.

Degradome library construction and target identification

Total RNAs from the fresh leaves and apical meristem of L. × formolongi at five developmental stages were equally mixed as a single sample to construct the degradome library according to published protocols [55]. Degradome cDNA library sequencing was also processed using Illumina Hiseq 2500 (LC Bio, Hangzhou, China). The raw reads were processed to remove adaptor sequences and low-quality sequences, the clean reads were generated and subjected to potential cleaved target identification by comparison with transcriptome sequences according to the description in previous study [56]. Thus, truncated transcripts resulting from endonucleolytic cleavage guided by miRNAs could be identified [57]. Furthermore, the alignment analysis of comparisons between identified miRNAs and mapped mRNAs was simultaneously performed using computational target prediction algorithms in TargetFinder (https://github.com/carringtonlab/TargetFinder) software to identify miRNA binding sites. Further annotation of targets was conducted in the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway databases for further miRNA-gene regulatory network analysis after all target unigenes were subjected to Blast2GO (https://www.blast2go.com/) analysis and aligned against the enriched KEGG pathways following a previously reported method [58].

Quantitative real-time PCR analysis

Total RNA and miRNA were extracted from leaves and apical meristems at six developmental stages from the seedling to visible floral bud phases under long-day conditions using RNA extraction kits (Aidlab Biotechnology Co., LTD., China; Tiangen Bio, Beijing, China) and miRcute miRNA Isolation kits (Tiangen Bio, Beijing, China) according to the manufacturer’s instructions. The mature miRNA was polyadenylated using the E. coli Poly (A) Poly Polymerase Kit (Invitrogen, Carlsbad, CA, USA) [59]. Then, the obtained RNAs and miRNAs were subjected to quality characterization, poly (A)-enrichment and reverse transcription to obtain cDNA using Quant Script RT Kit (Tiangen Bio, Beijing, China) and ReverTra Ace qPCR RT Master Mix with gDNA Remover (Toyobo, Shanghai, China). In addition, to assess tissue-specific miRNA expression, total miRNAs were separately extracted from different tissues from seedlings in the same developmental phases. RT-qPCR was carried out by using THUNDERBRIRD SYBR qPCR Mix Without Rox (Toyobo, Shanghai, China) on a Bio-Rad CFX96 system (CFX96 Touch, BIO-RAD, USA). Each reaction was performed in triplicate. The 5S RNA and EF genes were used as the internal reference genes for the normalization of the results [60, 61]. The primers for the miRNA and targets were designed using Beacon Designer 7 software, as indicated in Supplementary Table 6. The − 2−∆∆Ct method was employed to statistically analyse these data [62].

Coexpression network construction

To generate a correlation network of miRNAs during lily development, an undirected weighted gene network was constructed with a weighted gene coexpression network analysis (WGCNA) algorithm. A tutorial for the WGCNA package in R is provided at https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA. The normalized miRNA expression data from the libraries were used for correlation network construction after gene screening with a power of 12. We focused the analysis on module eigengenes whose targets were also significantly differentially expressed during lily flowering. The correlation network of interest was selected for visualization in Cytoscape (version 3.6.1). Moreover, the alignment of the sequences of miRNA precursors and the phylogenetic analysis of miR156 families were performed via the ClustalW and NJ methods in MEGA 6.06 software with default parameters. The precursor sequences of other species were obtained from miRBase, Release 21 (http://www.mirbase.org/index.shtml).

Availability of data and materials

The sequencing data have been submitted to NCBI SRA (Sequence Read Archive, http://0-www-ncbi-nlm-nih-gov.brum.beds.ac.uk/sra/), under the accession number PRJNA712953.

Abbreviations

miRNAs:

MicroRNAs

FT :

FLOWERING LOCUS T

T6P:

Trehalose-6-phosphate

nt:

Nucleotide

SPL :

SQUAMOSA PROMOTER BINDING PROTEIN-LIKE

AP2 :

APETALA 2

GA:

Gibberellin

ARF :

AUXIN RESPONSE FACTORS

GO:

Gene Ontology

KEGG:

Kyoto Encyclopedia of Genes and Genomes

NJ:

Neighbor-Joining

RT-qPCR:

Quantitative real time polymerase chain reaction

WGCNA:

weighted gene coexpression network analysis

rRNA:

Ribosomal RNA

tRNA:

Transfer ribonucleic acid

snoRNA:

Small nucleolar RNA

snRNA:

Small nuclear RNA

sRNA:

Small RNA

mRNA:

Messenger RNA

NAC (NAM):

NO APICAL MERISTEM

DEGs:

Differentially expressed genes

LAC :

LACCASE

FAO :

FATTY ALCOHOL OXIDASE

AGL :

AGAMOUS-LIKE

SOC1 :

SUPPRESSOR OF OVEREXPRESSION OF CONSTANS 1

References

  1. Bakhshaie M, Khosravi S, Azadi P, Bagheri H, Tuyl JMV. Biotechnological advances in Lilium. Plant Cell Rep. 2016;35(9):1799–826. https://0-doi-org.brum.beds.ac.uk/10.1007/s00299-016-2017-8.

    Article  CAS  PubMed  Google Scholar 

  2. Fortanier EJ. Reviewing the length of the generation period and its shortening, particularly in tulips. Sci Hort. 1973;1(1):107–16. https://0-doi-org.brum.beds.ac.uk/10.1016/0304-4238(73)90010-1.

    Article  Google Scholar 

  3. Li YF, Zhang MF, Zhang M, Jia GX. Analysis of global gene expression profiles during the flowering initiation process of Lilium × formolongi. Plant Mol Biol. 2017;94(4-5):361–79. https://0-doi-org.brum.beds.ac.uk/10.1007/s11103-017-0612-x.

    Article  CAS  PubMed  Google Scholar 

  4. Huijser P, Schmid M. The control of developmental phase transitions in plants. Development. 2011;138(19):4117–29. https://0-doi-org.brum.beds.ac.uk/10.1242/dev.063511.

    Article  CAS  PubMed  Google Scholar 

  5. Srikanth A, Schmid M. Regulation of flowering time: all roads lead to Rome. Cell Mol Life Sci. 2011;68(12):2013–37. https://0-doi-org.brum.beds.ac.uk/10.1007/s00018-011-0673-y.

    Article  CAS  PubMed  Google Scholar 

  6. Wahl V, Ponnu J, Schlereth A, Arrivault S, Langenecker T, Franke A, et al. Regulation of flowering by trehalose-6-phosphate signaling in Arabidopsis thaliana. Science. 2013;339(6120):704–7. https://0-doi-org.brum.beds.ac.uk/10.1126/science.1230406.

    Article  CAS  PubMed  Google Scholar 

  7. Hong Y, Jackson S. Floral induction and flower formation- the role and potential applications of miRNAs. Plant Biotechnology J. 2015;13(3):282–92. https://0-doi-org.brum.beds.ac.uk/10.1111/pbi.12340.

    Article  CAS  Google Scholar 

  8. Wilczynska A, Bushell M. The complexity of miRNA-mediated repression. Cell Death Differ. 2015;229(1):22–33.

    Article  Google Scholar 

  9. Rogers K, Chen X. Biogenesis, turn over, and mode of action of plant microRNAs. Plant Cell. 2013;25(7):2383–99. https://0-doi-org.brum.beds.ac.uk/10.1105/tpc.113.113159.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Dario M, Griffithsjones S, Kim M. Small RNAs: big impact on plant development. Trends Plant Sci. 2017;22(12):1056–68. https://0-doi-org.brum.beds.ac.uk/10.1016/j.tplants.2017.09.009.

    Article  CAS  Google Scholar 

  11. Luo Y, Guo Z, Li L. Evolutionary conservation of microRNA regulatory programs in plant flower development. Dev Biol. 2013;380(2):133–44. https://0-doi-org.brum.beds.ac.uk/10.1016/j.ydbio.2013.05.009.

    Article  CAS  PubMed  Google Scholar 

  12. He J, Xu M, Willmann MR, McCormick K, Hu T, Yang L, et al. Threshold-dependent repression of SPL gene expression by miR156/miR157 controls vegetative phase change in Arabidopsis thaliana. PLoS Genet. 2018;14(4):e1007337. https://0-doi-org.brum.beds.ac.uk/10.1371/journal.pgen.1007337.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Wang J, Czech B, Weigel D. MiR156- regulated SPL transcription factors define an endogenous flowering pathway in Arabidopsis thaliana. Cell. 2009;138(4):738–49. https://0-doi-org.brum.beds.ac.uk/10.1016/j.cell.2009.06.014.

    Article  CAS  PubMed  Google Scholar 

  14. Mathieu J, Yant LJ, Murdter F, Kuttner F, Schmid M. Repression of flowering by the miR172 target SMZ. PLoS Biol. 2009;7(7):e1000148. https://0-doi-org.brum.beds.ac.uk/10.1371/journal.pbio.1000148.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Aukerman MJ, Sakai H. Regulation of flowering time and floral organ identity by a microRNA and its APETALA2-like target genes. Plant Cell. 2003;15(11):2730–41. https://0-doi-org.brum.beds.ac.uk/10.1105/tpc.016238.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. Wu G, Park MY, Conway SR, Wang JW, Weigel D, Poethig RS. The sequential action of miR156 and miR172 regulates developmental timing in Arabidopsis. Cell. 2009;138(4):750–9. https://0-doi-org.brum.beds.ac.uk/10.1016/j.cell.2009.06.031.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Guo C, Xu Y, Shi M, Lai Y, Wu X, Wang H, et al. Repression of miR156 by miR159 regulates the timing of the juvenile-to-adult transition in Arabidopsis. Plant Cell. 2017;29(6):1293–304. https://0-doi-org.brum.beds.ac.uk/10.1105/tpc.16.00975.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Moon J, Suh SS, Lee H, Choi KR, Hong CB, Paek NC, et al. The SOC1 MADS-box gene integrates vernalization and gibberellin signals for flowering in Arabidopsis. Plant J. 2003;35(5):613–23. https://0-doi-org.brum.beds.ac.uk/10.1046/j.1365-313X.2003.01833.x.

    Article  CAS  PubMed  Google Scholar 

  19. Silva GF, Silva EM, Correa JP, Vicente MH, Jiang N, Notini MM, et al. Tomato floral induction and flower development are orchestrated by the interplay between gibberellin and two unrelated microRNA-controlled modules. New Phytol. 2019;221(3):1328–44. https://0-doi-org.brum.beds.ac.uk/10.1111/nph.15492.

    Article  CAS  PubMed  Google Scholar 

  20. Curaba J, Talbot M, Li Z, Helliwell C. Over-expression of microRNA171 affects phase transitions and floral meristem determinancy in barley. BMC Plant Biol. 2013;13(1):6. https://0-doi-org.brum.beds.ac.uk/10.1186/1471-2229-13-6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Yang R, Li P, Mei H, Wang D, Sun J, Yang C, et al. Fine-tuning of miR528 accumulation modulates flowering time in rice. Mol Plant. 2019;12(8):1103–13. https://0-doi-org.brum.beds.ac.uk/10.1016/j.molp.2019.04.009.

    Article  CAS  PubMed  Google Scholar 

  22. Hu J, Zhou Y, He F, Dong X, Liu L, Coupland G, et al. miR824-regulated AGAMOUS-LIKE16 contributes to flowering time repression in Arabidopsis. Plant Cell. 2014;26:2024–37.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Wu L, Liu D, Wu J, Zhang R, Qin Z, Liu D, et al. Regulation of FLOWERING LOCUS T by a MicroRNA in Brachypodium distachyon. Plant Cell. 2013;25(11):4363–77. https://0-doi-org.brum.beds.ac.uk/10.1105/tpc.113.118620.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Lazare S, Zaccai M. Flowering pathway is regulated by bulb size in Lilium longiflorum (Easter lily). Plant Biol. 2016;18(4):577–84. https://0-doi-org.brum.beds.ac.uk/10.1111/plb.12440.

    Article  CAS  PubMed  Google Scholar 

  25. Li Y, Zhao Y, Zhang M, Jia G. Functional and evolutionary characterization of the CONSTANS-like family in Lilium ×formolongi. Plant Cell Physiol. 2018;59(9):1874–88. https://0-doi-org.brum.beds.ac.uk/10.1093/pcp/pcy105.

    Article  CAS  PubMed  Google Scholar 

  26. Peruzzi L, Leitch IJ, Caparelli KF. Chromosome diversity and evolution in Liliaceae. Ann Bot. 2009;103(3):459–75. https://0-doi-org.brum.beds.ac.uk/10.1093/aob/mcn230.

    Article  CAS  PubMed  Google Scholar 

  27. Lin PC, Lu CW, Shen BN, Lee GZ, Bowman JL, Arteagavazquez MA, et al. Identification of miRNAs and their targets in the liverwort Marchantia polymorpha by integrating RNA-Seq and degradome analyses. Plant Cell Physiol. 2016;57(2):339–58. https://0-doi-org.brum.beds.ac.uk/10.1093/pcp/pcw020.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Edwards D, Batley J. Plant genome sequencing: applications for crop improvement. Plant Biotechnol J. 2010;8(1):2–9. https://0-doi-org.brum.beds.ac.uk/10.1111/j.1467-7652.2009.00459.x.

    Article  CAS  PubMed  Google Scholar 

  29. Han X, Yin H, Song X, Zhang Y, Liu M, Sang J, et al. Integration of small RNAs, degradome and transcriptome sequencing in hyperaccumulator Sedum alfredii uncovers a complex regulatory network and provides insights into cadmium phytoremediation. Plant Biotechnol J. 2016;14(6):1470–83. https://0-doi-org.brum.beds.ac.uk/10.1111/pbi.12512.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Axtell MJ, Bartel DP. Antiquity of microRNAs and their targets in land plants. Plant Cell. 2005;17(6):1658–73. https://0-doi-org.brum.beds.ac.uk/10.1105/tpc.105.032185.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Xing L, Zhang D, Zhao C, Li Y, Ma J, An N, et al. Shoot bending promotes flower bud formation by miRNA-mediated regulation in apple (Malus domestica Borkh.). Plant Biotechnol J. 2016;14:749–70.

    Article  CAS  PubMed  Google Scholar 

  32. Yant L, Mathieu J, Dinh TT, Ott F, Lanz C, Wollmann H, et al. Orchestration of the floral transition and floral development in Arabidopsis by the bifunctional transcription factor APETALA2. Plant Cell. 2010;22(7):2156–70. https://0-doi-org.brum.beds.ac.uk/10.1105/tpc.110.075606.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Xie K, Shen J, Hou X, Yao J, Li X, Xiao J, et al. Gradual increase of miR156 regulates temporal expression changes of numerous genes during leaf development in rice. Plant Physiol. 2012;158(3):1382–94. https://0-doi-org.brum.beds.ac.uk/10.1104/pp.111.190488.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Yang L, Conway SR, Poethig RS. Vegetative phase change is mediated by a leaf-derived signal that represses the transcription of miR156. Development. 2011;138(2):245–9. https://0-doi-org.brum.beds.ac.uk/10.1242/dev.058578.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Feng S, Xu Y, Guo C, Zheng J, Zhou B, Zhang Y, et al. Modulation of miR156 to identify traits associated with vegetative phase change in tobacco (Nicotiana tabacum). J Exp Bot. 2016;67(5):1493–504. https://0-doi-org.brum.beds.ac.uk/10.1093/jxb/erv551.

    Article  CAS  PubMed  Google Scholar 

  36. Wang J, Park MY, Wang L, Koo Y, Chen X, Weigel D, et al. MiRNA control of vegetative phase change in trees. PLoS Genet. 2011;7(2):e1002012. https://0-doi-org.brum.beds.ac.uk/10.1371/journal.pgen.1002012.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Yu S, Cao L, Zhou C, Zhang T, Lian H, Sun Y, et al. Sugar is an endogenous cue for juvenile-to-adult phase transition in plants. eLife. 2013;2:e00269. https://0-doi-org.brum.beds.ac.uk/10.7554/eLife.00269.

    Article  PubMed  PubMed Central  Google Scholar 

  38. Francozorrilla JM, Valli A, Todesco M, Mateos I, Puga MI, Rubiosomoza I, et al. Target mimicry provides a new mechanism for regulation of microRNA activity. Nat Genet. 2007;39(8):1033–7. https://0-doi-org.brum.beds.ac.uk/10.1038/ng2079.

    Article  CAS  Google Scholar 

  39. Xie K, Wu C, Xiong L. Genomic organization, differential expression, and interaction of SQUAMOSA promoter-binding-like transcription factors and microRNA156 in rice. Plant Physio. 2006;142(1):280–93. https://0-doi-org.brum.beds.ac.uk/10.1104/pp.106.084475.

    Article  CAS  Google Scholar 

  40. Zhu Q, Helliwell CA. Regulation of flowering time and floral patterning by miR172. J Exp Bot. 2011;62(2):487–95. https://0-doi-org.brum.beds.ac.uk/10.1093/jxb/erq295.

    Article  CAS  PubMed  Google Scholar 

  41. Jung J, Seo Y, Seo PJ, Reyes JL, Yun J, Chua N, et al. The GIGANTEA-regulated microRNA172 mediates photoperiodic flowering independent of CONSTANS in Arabidopsis. Plant Cell. 2007;19(9):2736–48. https://0-doi-org.brum.beds.ac.uk/10.1105/tpc.107.054528.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Tang M, Bai X, Niu L, Chai X, Chen M, Xu Z. miR172 regulates both vegetative and reproductive development in the perennial woody plant Jatropha curcas. Plant Cell Physiol. 2018;59(12):2549–63. https://0-doi-org.brum.beds.ac.uk/10.1093/pcp/pcy175.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Rubiosomoza I, Weigel D. Coordination of flower maturation by a regulatory circuit of three microRNAs. PLoS Genet. 2013;9(3):e1003374. https://0-doi-org.brum.beds.ac.uk/10.1371/journal.pgen.1003374.

    Article  CAS  Google Scholar 

  44. Tsuji H, Aya K, Ueguchitanaka M, Shimada Y, Nakazono M, Watanabe R, et al. GAMYB controls different sets of genes and is differentially regulated by microRNA in aleurone cells and anthers. Plant J. 2006;47(3):427–44. https://0-doi-org.brum.beds.ac.uk/10.1111/j.1365-313X.2006.02795.x.

    Article  CAS  PubMed  Google Scholar 

  45. Achard P, Herr AJ, Baulcombe DC, Harberd NP. Modulation of floral development by a gibberellin-regulated microRNA. Development. 2004;131(14):3357–65. https://0-doi-org.brum.beds.ac.uk/10.1242/dev.01206.

    Article  CAS  PubMed  Google Scholar 

  46. Aceto S, Sica M, De Paolo S, D’Argenio V, Cantiello P, Salvatore F, et al. The analysis of the inflorescence miRNome of the orchid Orchis italica reveals a DEF-like MADS-box gene as a new miRNA target. PLoS One. 2014;9(5):e97839. https://0-doi-org.brum.beds.ac.uk/10.1371/journal.pone.0097839.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Preston JC, Hileman LC. Functional evolution in the plant SQUAMOSA-PROMOTER BINDING PROTEIN-LIKE (SPL) gene family. Front Plant Sci. 2013;5(4):80.

    Google Scholar 

  48. Wagner D. Making flowers at the right time. Dev Cell. 2016;37(3):208–10. https://0-doi-org.brum.beds.ac.uk/10.1016/j.devcel.2016.04.021.

    Article  CAS  PubMed  Google Scholar 

  49. Gou J, Tang C, Chen N, Wang H, Debnath S, Sun L, et al. SPL7 and SPL8 represent a novel flowering regulation mechanism in switchgrass. New Phytol. 2019;222(3):1610–23. https://0-doi-org.brum.beds.ac.uk/10.1111/nph.15712.

    Article  CAS  PubMed  Google Scholar 

  50. Alonsoperal MM, Li J, Li Y, Allen RS, Schnippenkoetter W, Ohms S, et al. The microRNA159-regulated GAMYB-like genes inhibit growth and promote programmed cell death in Arabidopsis. Plant Physiol. 2010;154(2):757–71. https://0-doi-org.brum.beds.ac.uk/10.1104/pp.110.160630.

    Article  CAS  Google Scholar 

  51. Lucidos JG, Ryu KB, Younis A, Kim C, Hwang Y, Son B, et al. Different day and night temperature responses in Lilium hansonii in relation to growth and flower development. Hortic Environ Biote. 2013;54(5):405–11. https://0-doi-org.brum.beds.ac.uk/10.1007/s13580-013-1241-1.

    Article  CAS  Google Scholar 

  52. RubioSomoza I, Weigel D. MicroRNA networks and developmental plasticity in plants. Trends Plant Sci. 2011;16(5):258–64. https://0-doi-org.brum.beds.ac.uk/10.1016/j.tplants.2011.03.001.

    Article  CAS  Google Scholar 

  53. Jeyaraj A, Zhang X, Hou Y, Shangguan M, Gajjeraman P, Li Y, et al. Genome-wide identification of conserved and novel microRNAs in one bud and two tender leaves of tea plant (Camellia sinensis) by small RNA sequencing, microarray-based hybridization and genome survey scaffold sequences. BMC Plant Biol. 2017;17(1):212. https://0-doi-org.brum.beds.ac.uk/10.1186/s12870-017-1169-1.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. Li X, Shahid MQ, Wu J, Wang L, Liu X, Lu Y. Comparative small RNA analysis of pollen development in autotetraploid and diploid rice. Int J Mol Sc. 2016;17(4):499. https://0-doi-org.brum.beds.ac.uk/10.3390/ijms17040499.

    Article  CAS  Google Scholar 

  55. Addo-Quaye C, Eshoo TW, Bartel DP, Axtell MJ. Endogenous siRNA and miRNA targets identified by sequencing of the Arabidopsis degradome. Curr Biol. 2008;18(10):758–62. https://0-doi-org.brum.beds.ac.uk/10.1016/j.cub.2008.04.042.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  56. Jonesrhoades MW, Bartel DP. Computational identification of plant microRNAs and their targets, including a stress-induced miRNA. Mol Cell. 2004;14(6):787–99. https://0-doi-org.brum.beds.ac.uk/10.1016/j.molcel.2004.05.027.

    Article  CAS  Google Scholar 

  57. Folkes L, Moxon S, Woolfenden HC, Stocks MB, Szittya G, Dalmay T, et al. PAREsnip: a tool for rapid genome-wide discovery of small RNA/target interactions evidenced through degradome sequencing. Nucleic Acids Res. 2012;40(13):e103. https://0-doi-org.brum.beds.ac.uk/10.1093/nar/gks277.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  58. Xie F, Frazier TP, Zhang B. Identification and characterization of microRNAs and their targets in the bioenergy plant switchgrass (Panicum virgatum). Planta. 2010;232(2):417–34. https://0-doi-org.brum.beds.ac.uk/10.1007/s00425-010-1182-1.

    Article  CAS  PubMed  Google Scholar 

  59. Shi R, Chiang VL. Facile means for quantifying microRNA expression by real-time PCR. Bio Techniques. 2018;39:519–25.

    Google Scholar 

  60. Jia XL, Chen YK, Xu XZ, Shen F, Zheng QB, Du Z, et al. miR156 switches on vegetative phase change under the regulation of redox signals in apple seedlings. Sci Rep. 2017;7(1):14223.

    Article  PubMed  PubMed Central  Google Scholar 

  61. Liu Q, Wei C, Zhang MF, Jia GX. Evaluation of putative reference genes for quantitative real-time PCR normalization in Lilium regale during development and under stress. Peer J. 2016;4:e1837. https://0-doi-org.brum.beds.ac.uk/10.7717/peerj.1837.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  62. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2 (−Delta Delta C (T)) method. Methods. 2001;25(4):402–8. https://0-doi-org.brum.beds.ac.uk/10.1006/meth.2001.1262.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgements

Not applicable.

Funding

This work was supported by National Key R&D Program of China (Grant No. 2019YFD1000400) and National Natural Science Foundation of China (Grant No. 31772348). The funders did not have any role in the design of the study, collection, analysis, or interpretation of data or the writing of the manuscript.

Author information

Authors and Affiliations

  1. Beijing 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

    Qian Zhang, Yu-Qian Zhao, Xue Gao & Gui-Xia Jia

Authors
  1. Qian Zhang
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Yu-Qian Zhao
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Xue Gao
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Gui-Xia Jia
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

GXJ designed the experiments. QZ and YQZ conducted the experiments. QZ performed the bioinformatics analysis with the help of XG and drafted the manuscript. GXJ directed and modified the manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Gui-Xia Jia.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Fig. S1

. Pearson correlation between samples. Fig. S2. Principal component analysis based on the expression levels in libraries.

Additional file 2: Table S1

. Conserved miRNA families in lily. Table S2. Novel miRNA families in lily. Table S3. Significant differential expression of miRNAs in lily. Table S4. Identified and annotated target transcripts for the conservative miRNAs in lily. Table S5. Identified and annotated target transcripts for the novel miRNAs in lily. Table S6. Specific RT-qPCR primers of lily. Table S7. The nodes of the gene network weight analysis results generated by WGCNA. Table S8. The edges of the gene network weight analysis results generated by WGCNA.

Additional file 3: Fig. S3

. Predicted hairpin structures of precursors of novel miRNAs. The red-coloured sequences represent mature miRNAs, and the yellow-coloured sequences represent miRNAs

Additional file 4: Fig. S4

. The selection of soft-thresholding power for the coexpression network construction by WGCNA. 12 was selected as soft- thresholding power parameters. Fig. S5. Hierarchical clustering tree. Hierarchical clustering tree showing coexpression modules detected by WGCNA. The major tree branches constitute 5 modules labelled by different colours, including blue, turquoise, yellow, green and brown module. Fig. S6. Cluster tree based of the module eigengenes. Fig. S7. The correlation coefficient heatmap of the coexpression module genes. Each bright spot corresponds to the correlation between each miRNA and other miRNAs. The deeper the colours, the stronger the connectivity between the two miRNAs in the corresponding row and column.

Additional file 5: Fig. S8

. The phylogenetic relationship of miR156 family homologs. Neighbor-joining (NJ) tree constructed using precursor miRNA family sequences from Arabidopsis (ath), Oryza sativa (osa) and L. × formolongi (lfo). Red dots represent precursor sequences from lily, green dots and blue dots represent precursor sequences from the members which were verified important role in regulating flowering time in Arabidopsis and Oryza sativa, respectively.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Zhang, Q., Zhao, YQ., Gao, X. et al. Analysis of miRNA-mediated regulation of flowering induction in Lilium × formolongi. BMC Plant Biol 21, 190 (2021). https://0-doi-org.brum.beds.ac.uk/10.1186/s12870-021-02961-3

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://0-doi-org.brum.beds.ac.uk/10.1186/s12870-021-02961-3

Keywords

BMC Plant Biology

ISSN: 1471-2229

Contact us