Early-flowering variety Paeonia ostii ‘Fengdan’ (FD), an early-flowering mutant line of Paeonia ostii ‘Fengdan’ (MU), and late-flowering variety Paeonia suffruticosa ‘Lianhe’ (LH) were selected as the experiment materials. FD, MU, and LH used in this study were 13-years-old plants with single and white flowers. Fresh petals were collected at 9:00-10:00 am on different days at blooming stage (BS), initial flowering stage (IF), full bloom stage (FB), and decay stage (DE), respectively. Abbreviations of sample and library names were presented in Table 1 . Three biological replicates are different flowers on different stems of the same plant at each developmental stage for each variety respectively. For each flower, all the petals were sampled and pooled prior to freezing by liquid nitrogen. The sampled petals were stored in the freezer (-80°C) for RNA extraction.

In total, 36 libraries were prepared for miRNAome and transcriptome analysis separately ( Table S1 ). Approximately 200 mg petals were used for total RNA extraction. Total RNA integrity was initially assessed by denaturing agarose gel electrophoresis, then confirmed by Bioanalyzer 2100 (Agilent, CA, USA). Total RNA amount and purity quantification were performed on NanoDrop ND-1000 (NanoDrop, Wilmington, DE, USA). Total RNA with a concentration>100 ng/μL, RNA integrity number >7.0, OD260/280>1.8 and amount>50 μg were used for library construction. Petals from four flower development stages of each variety were mixed prior to the degradome library construction. That is, a total of 3 degradome libraries were constructed for degradome sequencing. Kits and reagents used for RNA isolation, purification, quantification, and libraries construction are listed in Table S2 . Transcriptome sequencing was performed by the 2×150 bp paired-end sequencing (PE150) on Illumina Novaseq™ 6000. miRNAome and degradome sequencing were performed by the 1×50 bp single-end sequencing on Illumina Hiseq2500. Libraries construction and sequencing were performed at LC-BIO (Hangzhou, China) according to the vendor’s recommended protocol.

In total, eight miRNA-target pairs associated with floral florescence, development and senescence were randomly selected for quantitative real-time PCR (qRT-PCR) analysis. Fresh petals of FD, MU and LH were collected at developmental stages (BS, IF, FB, DE) respectively. For each developmental stage, petals from three different flowers on different stems of the same plant were sampled individually. Total RNA extraction, miRNA extraction, cDNA synthesis for mRNA and miRNA were performed according to the instructions of manufacturers. Kits information for total RNA/miRNA extraction and cDNA synthesis are shown in Table S2 . The EF1-α and U6 were used as the reference for mRNA and tailing reaction miRNA analysis separately. Primer sequences for qRT-PCR assay are listed in Table S3 . qRT-PCR analysis for miRNA and targets were both performed on a BIORAD CFX96 machine. Three technical replicates per reaction were conducted in the qRT-PCR analyses to ensure statistical validity. The relative quantity was calculated on the basis of 2^−ΔΔCT method (Livak and Schmittgen, 2001).

The date when 80% of flowers reached color exposure (CE), blooming stage (BS), initial flowering (IF), half opening (HO), full blooming (FB), initial decay (ID), and decay (DE) stages were investigated in FD, MU, and LH in 2020 and 2021 respectively ( Figure 1A ). The flower duration time (date from CE to DE) of 80% flowers of FD, MU, LH were also investigated. FD is an early flower variety, while LH is a late flowering variety. MU was a mutant of FD, which blossomed seven days earlier than FD in 2020, and nine days earlier than FD in 2021. Flowering time of FD was 16 days earlier than LH in 2020, and 10 days earlier than LH in 2021. Flowering time of MU was 23 and 19 days earlier than LH in 2020 and 2021, respectively. In addition, floral florescence per plant was 16-17 days for FD, 13-15 days for MU, and 11-12 days for LH, which demonstrated that FD possesses the longest blooming time ( Figures 1B, C ).

Degradome, used for miRNA and siRNA targets characterization by the 5’-ends of uncapped RNAs (German et al., 2008; Gregory et al., 2008). Petals of the four flower developmental stages (BS, IF, FB and DE) for each genotype (FD, MU and LH) were pooled prior the degradome libraries construction (FD, LH, and MU) and used for degradome analysis. Around 10.9 million raw reads were gained from the degradome libraries. Detailed information of total raw reads, unique raw reads, mappable reads, unique mappable reads, mapped reads, unique mapped reads, number of input transcripts, and number of target transcripts are shown in Table S9 . A total of 7,571, 9,250 and 6,457 mRNAs were silenced by miRNAs in FD, MU and LH respectively. The alignment information of miRNAs-mRNA pairs can be found in Table S10 . In total, 11,391 targets showed a differential degradation pattern between FD and LH. In addition, 11,928 targets presented a differential degradation characteristic between MU and FD. Moreover 11,605 targets demonstrated a differential degradation characteristic between MU and LH.

DEGs identified across flower developmental stages (BS, IF, FB, DE) in FD, MU and LH were assembled into a unified set. DEGs identified across varieties (FD, MU, LH) at developmental stages BS, IF, FB and DE were then assembled into another unified set. DEGs from the union sets were then used together for further miRNA-mRNA target pairing confirmation by multi-omics analysis. Subsequently, miRNA-mRNA pairs with opposite regulatory patterns were characterized in terms of the gene-silencing function of miRNAs. With this method, miRNA-mRNA pairs associated with floral florescence, development and senescence were identified.

In total, 32 miRNA targets showed antagonistic regulatory patterns during developmental stages in FD, MU and LH, which might thus be candidates for regulating the floral development and senescence of tree peony ( Table S11 ). Another 191 miRNA targets were differentially expressed across tree peony varieties and developmental stages, which suggests important roles for those in floral florescence regulation (i.e., the timing of flower opening and senescence) ( Table S11 ). Of these, ten miRNA targets showed significantly different expression patterns both across development and across varieties, which suggests possible dual roles in floral florescence regulation in tree peony. Expression patterns of these floral florescence, development and senescence dependent miRNA targets are presented in Figures 6A–E . GO ( Figure 7A ) and KEGG ( Figure 7B ) pathway analysis revealed that these floral florescence-, development- and senescence-dependent miRNA targets were enriched in pathways like plant hormone signal transduction, indole alkaloid biosynthesis, arachidonic acid metabolism, folate biosynthesis, fatty acid elongation, MAPK signaling pathway, etc.

The miRNA-guided floral florescence, development and senescence regulatory networks were complicated, which imply that one specific miRNA might be able to adjust and control many mRNAs, and one specific mRNA also might be targeted by divers miRNAs simultaneously (Liu et al., 2020a). Multiple-to-multiple inter-associations between miRNAs and their target genes which enriched in KEGG pathways and encode TFs were constructed by Cytoscape ( Figure 8 ). The multiple-to-multiple miRNA-mRNA-TF modules were enriched in pyruvate metabolism, carbon fixation in photosynthetic organisms, pentose and glucuronate interconversions, sesquiterpenoid and triterpenoid biosynthesis, aminoacyl-tRNA biosynthesis ( Figure 8A ). The result also displayed that the miRNA-mRNA-TF modules mainly consisted of MYB-related, bHLH, Trihelix, NAC, GRAS and HD-ZIP TF families, demonstrating their potential functions in tree peony floral florescence, development and senescence ( Figure 8B ).

In order to uncover the regulation mechanism of flowering time based on transcriptome data across varieties and flower development stages, weighted gene co-expression network analysis (WGCNA) was performed to detect co-expressed genes to disclose the hub gene which might regulate floral florescence, development and senescence. WGCNA analysis resulted in 43 distinct co-expressed gene modules were exhibited by distinctive colors and shown by a heatmap ( Figure 9A ). Each heatmap represented an expression cluster, which straightly elucidated the relationship between the clusters of three tree peony varieties and four development stages ( Figure 9B ). Then, correlation analysis was performed between modules and samples to find modules with the highest correlation. The candidate hub genes were confirmed by taking the intersection of gene in modules with the highest correlation and the genes used for integrated analysis of miRNA-mRNA-TF. Furthermore, the intersected genes were selected for the hub genes network construction ( Figures 9C, D ). Finally, hub genes psu.G.00014449, psu.G.00003047 and psu.G.00009129 in LH, psu.G.00032165 and psu.G.00007421 in MU may play crucial roles in the floral florescence, development and senescence in tree peony. Regrettably, none intersected hub genes were detected in FD since the lower gene connectivity.

A total of eight DEM and DEG pairs were selected for the qRT-PCR analysis to verify the expression pattern of miRNA and mRNA data obtained from miRNAome, transcriptome, and degradome ( Figure 10 ). The selected miRNA-target pairs were presented as follows: seu-MIR11025-p5 and psu.T.00024044, miR166-5p and psu.T.00024044, PC-5p-564_43386 and psu.T.00034433, mtr-miR396b-5p and psu.T.00010381, PC-3p-602268_25_S and psu.T.00020538, PC-5p-429002_51_S and psu.T.00018467, mtr-MIR2592bj-p3 and psu.T.00015108, PC-5p-143784_277_S and psu.T.00016751. Target t-plots ( Figure S6 ) show the cleavage sites of target genes silenced by miRNAs during developmental stages. In general, the expression levels of DEM-target pairs were consistent with miRNA-guided mRNA cleavage signatures validated by degradome.

As expected, majority of the expression patterns of the examined miRNA-target pairs were similar to the RNA-seq data, with the exception of transcripts psu.T.00010381, psu.T.00016751 and mtr-miR396b-5p, confirming the accuracy and reliability of the sequencing data in general ( Figure 10 ). Our study confirmed an exactly negative correlation for miRNAs and target genes PC-3p-602268_25 and psu.T.00020538 expressed at the four developmental stages in FD and LH, target genes PC-5p-429002_51 and psu.T.00018467 at all four developmental stages in FD, and target genes mtr-miR166g-5p and psu.T.00024044 at the four developmental stages in MU. While the remaining tested miRNA-target pairs did not always show a negative relationship across the four developmental stages, this is consistent with previous research (Liu et al., 2017; Zuluaga et al., 2017).

Previous research revealed that miRNA and their target genes were not always presented a specific one-to-one regulatory relationship, which is due to that a specific miRNA can regulate several target mRNAs simultaneously and a specific mRNA also could be targeted by multiple miRNAs (Liu et al., 2020a; Liu et al., 2020b). Thus, the expression pattern of miRNA and their target genes does not emerge a negative correlation all the time (Liu et al., 2017; Zuluaga et al., 2017; Zuluaga et al., 2018). Our study showed that psu.T.00024044 could be targeted by mtr-miR166g-5p and seu-MIR11025-p5_2ss4CA17CA simultaneously which was consistent with data elucidated previously. Furthermore, changes in expression levels of miRNA and mRNA could have derived from the differences arising from the biological replicates (Pradervand et al., 2009; Xu et al., 2018; Zuluaga et al., 2018; Liu et al., 2020b). Correlation analysis of miRNA-target pairs expression profiles between sequencing and qRT-PCR are shown in Figure S7 .

Floral-associated pathways usually involved the vernalization, autonomous, ambient temperature, photoperiod, gibberellic acid, aging and sugar pathways (Han et al., 2021). Genes involved flowering-time regulation in tree peony (Wang et al., 2019) that were identified including genes involved in floral organ and meristem, vernalization pathway, age pathway, GA pathway, autonomous pathway, photoperiod pathway. High intensity light promotes flowering through the photoperiod pathway which requires the cooperation of chloroplast retrograde signals and silencing transcription of Flowering Locus C (FLC). To response high light induction, transcription factor PTM localized at chloroplast envelope suppresses FLC transcription. It is also known that an intracellular signaling pathway originated from chloroplasts regulates the flowering transition (Feng et al., 2016; Susila et al., 2016). Furthermore, it was discovered that the expression of chloroplast protein CEBP changed during flower development and senescence (Iordachescu et al., 2009). Our results showing that chloroplast and chloroplast envelope genes were expressed in common among both varieties and developmental stages assayed, suggests important roles in floral florescence, development and senescence in tree peony. We also observed a general role for the interaction of turgor pressure and plasmodesmata affecting floral development, perhaps through regulation of plasmodesmata aperture, for which an association with transition to flowering was previously shown (Hernández-Hernández et al., 2019).

Cell division, expansion, differentiation and stress response in organisms were proved to be regulated by hormones in the earlier report (Artur, 2022). ABA signaling presents multiple connections with the photoperiodic pathway (Martignago et al., 2020). Exogenous spraying of ABA result in changes of flowering time, indicating that ABA possibly is an internal factor regulating the floral transition (Conti et al., 2014). GIGANTEA (GI) is a key flowering gene required for photoperiod perception (Mishra and Panigrahi, 2015). ABA signaling integration through GI operates via up-regulation of FT (Riboni et al., 2016). Overexpression of the chrysanthemum R2R3-MYB delays flowering in Arabidopsis (Shan et al., 2012). ABA hypersensitive 1 suppresses frigida-mediated delayed flowering in Arabidopsis (Bezerra et al., 2004). FD and FD-like bZIPs protein complexes play a significant role in modulating ABA signaling (Martignago et al., 2020). ABA activates an intricate regulatory network of signals including TFs that have contrary effects on florescence (Conti et al., 2014). While the role of ABA in flowering in model plants is emerging, the ABA molecular control of flowering still poorly revealed in tree peony. This research found that many of the most highly differentially expressed genes were relevant to plant hormone signal transduction especially ABA. These results might provide a potential basis for further research on mechanisms parsing of ABA regulating floral florescence, development and senescence in tree peony.

MiRNAs have been proved to be of great importance in regulation of gene expression, defense responses, and cell function in plants (Achkar et al., 2016; Choudhary et al., 2021). Recent research have shown that miRNAs play crucial roles in regulating gene expression associated with flowering (Spanudakis and Jackson, 2014). According to the latest report, miR156 and miR172 possibly participate in flowering via an aging pathway (Waheed and Zeng, 2020). MiR167 was reported to be involved in governing floral/fiber-associated agronomic traits in cotton (Arora and Chaudhary, 2021). Inhibition of miR168 in rice could improve yield, prolong flowering and enhance immunity (Wang et al., 2021). While overexpression of miR159 resulted in late flowering, whereas suppression of miR159 leaded to the acceleration of flowering in the ornamental flowering plant gloxinia (Li et al., 2013; Millar et al., 2019).

Researches have demonstrated that TFs play vital roles in floral transition as miRNA targets. MiRNAs and their TF targets regulate gene expression at post-transcriptional level and transcriptional level respectively (Waheed and Zeng, 2020). Previously, the miRNA-mRNA-TF modules like pos-miR319a-3p.2–3p/TCP2, pos-miR159/GAMYB, pos-miR169/nuclear transcription factor Y subunit A, and pos-miR828/WER were identified in variety FD (Han et al., 2020). TF target genes (AP2 and SPL) might have splice sites for PsmiR172a and PsmiR156a, suggesting that miR156 and miR172 probably play important roles during dormancy transition in FD (Zhang et al., 2018). Additionally, studies have revealed that miR156-SPL (Xie et al., 2020; Rao et al., 2021), miR172-AP2 (Ó’Maoiléidigh et al., 2021), miR319-TCP (Zhu et al., 2021), miR159-MYB (Spanudakis and Jackson, 2014) and miR399-PHO2 (Kim et al., 2011) play important roles in floral transition.

Transcription factors MYC2, MYC3, and MYC4 in bHLH family were involved in jasmonate-mediated flowering inhibition in Arabidopsis (Wang et al., 2017). Previous studies showed that the bZIP transcription factors were functionally required for flower development (Strathmann et al., 2001). The TBZF gene encoding bZIP were reported to be abundant in senescing flower buds (Yang et al., 2002). The C2H2 zinc finger family perform functions in pollen development regulation in grapevine (Arrey-Salas et al., 2021). Overexpression of the CcNAC1 gene promotes early flowering in jute (Zhang et al., 2021). Acting as a B3 domain transcription factor, AtREM16 prolongs flowering by coupling on the promoters of SOC1 and FT (Yu et al., 2020). Expression of chrysanthemum Trihelix transcription factors showed that they played important roles in chrysanthemum inflorescences (Song et al., 2016). It has also been shown that expression of chrysanthemum transcription factor ERF can influence flowering time in Arabidopsis (Xing et al., 2019). Additionally, it was shown that CmERF110 interacts with CmFLK to promote flowering by regulating the circadian clock (Huang et al., 2022).

In this study, we identified 16 miRNA-mRNA-TF modules across flower developmental stages ( Table S12 ). The TFs belongs to 16 families. Among these TF families, the bHLH, bZIP, and C2H2 accounted for the largest proportion, followed by B3, Trihelix, ERF, etc. In addition, 71 miRNA-mRNA-TF modules were identified across the tree peony varieties at the flower developmental stages ( Table S12 ). The TFs consisted of 37 TF families, of which the NAC, Trihelix and bHLH families accounted for the largest proportion, followed by ERF, B3, MYB family etc. In this study, the bHLH, NAC, C2H2, bZIP, and Trihelix displayed the most highly differential expression, suggesting that these miRNA-mRNA-TF modules may be crucial factors in floral florescence, development and senescence in tree peony. Their specific functional contributions in tree peony remains to be further explored.

The first draft genome assembly (~13.79 GB) of tree peony variety ‘Luo shen xiao chun’ reported recently represents the largest sequenced genome in dicotyledon to date (Lv et al., 2020). However, due to the unusually large and complex genome, the draft genome assembly of tree peony is still on scaffold level, which hinders miRNA-mRNA pairs identification by transcriptome and sRNAome analysis using the reference genome. A high-quality reference genome is expected to provide substantial fundamental resources for further research in tree peony in the future. Furthermore, lack of a homologous genetic transformation system has hindered functional genomics research in tree peony (Wen et al., 2020). Breakthroughs in virus-induced gene silencing (VIGS) in rose (Cheng et al., 2018; Liang et al., 2020; Cheng et al., 2021; Zhang et al., 2021) and tree peony (Zhao et al., 2017; Xie et al., 2019; Wang et al., 2020; Han et al., 2022), provide hope for the functional characterization and identification of miRNA-mRNA modules for floral florescence, development, and senescence identified in this study in tree peony.