Three developmental stages of the flower of M. candidum including closed buds with white petals (McI), closed buds with pink petals (McII), and opened buds with pink petals (McIII) were collected and immediately frozen in liquid nitrogen. Three replicates of each stage sample were taken from three seedlings. Total RNA was extracted by using an OminiPlant RNA Kit (DNase I) (CW2598, CWBIO, Taizhou, China) according to the manufacturer’s protocol. An Agilent 2100 Bioanalyzer was utilized to assess the quality of the RNA (Agilent Technologies, Palo Alto, CA, USA). To generate a cDNA library, we fragmented the mRNA using a fragmentation buffer and reverse-transcribed the resulting small fragments into cDNA using random primers. The second strand cDNA was synthesized by employing DNA polymerase I, RNaseH, dNTP, and buffer. After purification with the poly(A) and PCR extraction kit (Qiagen), the synthesized products were further purified using the QiaQuick PCR extraction kit (Qiagen, Venlo, The Netherlands) and then ligated to Illumina sequencing adapters. Afterward, the second-strand cDNA was digestedusing the enzyme UNG (Uracil-N-Glycosylase), and the resulting products were separated by size on an agarose gel before being amplified by PCR. The PCR products were then sequenced on Illumina HiSeq TM 4000 platforms.

After obtaining sequence data, raw reads consisting of adapters or low-quality bases were filtered by the fastp software (version 0.18.0) (Chen et al., 2018b) with steps: 1) removing reads containing adapters; 2) removing reads containing more than 10% of unknown nucleotides (N); 3) removing low quality reads containing more than 50% of low quality (Q-value ≤ 20) bases. After filtering, the clean short reads were aligned to the ribosome RNA (rRNA) database using Bowtie2 (version 2.2.8) (Langmead and Salzberg, 2012) to eliminate the rRNA-mapped reads. After building an index of the M. candidum genome, the paired-end clean reads were mapped to the reference genome by using HISAT2 (version 2.1.0) (Kim et al., 2015) with a default parameter “-rna-strandedness RF”.

The reconstruction of transcripts was carried out with the software Stringtie (version 1.3.4) (Pertea et al., 2015). To identify the novel transcripts among the reconstructed transcripts, all the reconstructed transcripts were aligned to the reference genome and were then categorized into 12 categories using the Cuffcompare program (Trapnell et al., 2010). A novel transcript was defined if it had one of the following class codes - u, i, j, x, c, e, or o. To further identify new genes, we used the following parameters: length of transcript > 200 bp, number of exons > 1 (Lu et al., 2019). All the novel transcripts were then aligned to the Nr, KEGG, and GO databases to obtain protein functional annotations.

Software including Coding-Non-Coding Index (CNCI, version2) (Sun et al., 2013) and Coding Potential Calculator 2 (CPC2, version 0.9-r2) (Kong et al., 2007) (http://cpc.cbi.pku.edu.cn/) were used to assess the protein-coding potential of novel transcripts. Only transcripts which meet the protein-coding-score criteria (CNCI sequence-score < 0 and CPC2 coding probability < 0.5) were considered as lncRNAs. Based on the location relative to protein-encoding genes (PCG), lncRNAs are categorized into five classes: intergenic lncRNAs (located between two PCGs), bidirectional lncRNAs (situated on the opposite strand but within 1 kb of the promoter on the sense strand), intronic lncRNAs (located within an intron of a PCG on the sense strand), antisense lncRNAs (transcribed from the opposite strand of a PCG), and sense lncRNA (spanning multiple introns or exons within a PCG) (Ahmad et al., 2021).

The abundance of the transcript was quantified using StringTie (Pertea et al., 2015). An FPKM (fragment per kilobase of transcript per million mapped reads) value was calculated for each transcription region to mitigate the effect of varying transcript lengths and sequencing data amounts on expression. DEGs of coding RNAs and lncRNAs were analyzed separately by using DESeq2 software (Li and Dewey, 2011) between two versus groups by the following steps: 1) normalization of read counts; 2) calculation of p-value; 3) correction of value to get a false discovery rate (FDR) value. Differentially expressed lncRNAs were screened with the threshold: Fold change ≥ 2, FDR ≤ 0.05.

For miRNAs, expression level was calculated and normalized by transcripts per million (TPM) methods. EdgeR (Robinson et al., 2010) were used to conducted differential analysis. DEmiRNAs were identified with the absolute threshold log₂ (fold change) ≥ 0.585, p-value ≤ 0.05. For mRNAs, HISAT2 (Kim et al., 2015) were used to align to the genome of M. candidum (http://evolution.sysu.edu.cn/Sequences.html). The FPKM value was used to quantify expression abundance and variations. DEmRNAs were identified by DESeq2 (Love et al., 2014) package with the absolute threshold log₂ (fold change) ≥ 1 and p-value ≤ 0.05. Metabolites were analyzed by using an LC-ESI-MS/MS system (UPLC, Shim-pack UFLCSHIMADZU CBM30A, http://www.shimadzu.com.cn/; and MS/MS (Applied Biosystems 6500 QTRAP, http://www.appliedbiosystems.com.cn/) (Chen et al., 2013) and variable importance in projection (VIP) score of (O)PLS model was applied to rank the metabolites and distinguish two versus groups. DSmets were identified with the absolute threshold log₂ (fold change) ≥ 0.585, vip ≤ 1.

In order to identify three direct-acting models of lncRNAs and mRNAs including antisense-regulation, cis-regulation, and trans-regulation, the software RNAplex (version 0.2) (Shen et al., 2014) (http://www.tbi.univie.ac.at/RNA/RNAplex.1.html) was used to predict the complementary correlation and detect the interaction between antisense lncRNAs and mRNAs. Antisense lncRNAs can form complementary base pairing with mRNAs. The ViennaRNA package (Lorenz et al., 2011) in R software was used to predict the best base pairing based on the calculation of minimum free energy through thermodynamics structure. Cis lncRNAs could regulate neighboring genes on the same allele. LncRNAs with the unknown region were annotated again to identify cis-regulators. Trans lncRNAs could regulate co-expressed genes far from them. The correlation of expression between lncRNAs and protein-coding genes was used to identify a trans-regulation relationship with a Pearson correlation of more than 0.999.

In order to verify the lncRNA-seq result, we chose 12 lncRNAs to perform qRT-PCR. The primers of lncRNAs used in qRT-PCR were listed in Table S1 . 0.5 µg total RNA of three stages flowers was reverse-transcribed into first-strand cDNA using the PrimeScript RT reagent Kit gDNA Eraser (Takara, Dalian, China). Then, the SYBR @Premix Ex Taq TMII (Takara, Dalian, China) was used according to the manufacturer’s instruction for qRT-PCR of lncRNAs on the Illumina Eco real-time PCR system (Illumina, USA). The α -tubulin gene of M. candidum was used as the internal reference gene. Ct values were then calculated by the 2^−ΔΔCt algorithm. Primers sequences, length of PCR products and PCR amplification efficiency for each pair of primers are listed in Table S1 . Relative expression levels of selected lncRNAs at different stages of flowers were analyzed using One-way ANOVA with multiple comparison by GraphPad Prism 9. The graphs were visualized by GraphPad Prism 9 and Adobe Illustrator 2020.

The ceRNA network was constructed based on the following rules: 1) negative correlation between miRNA and lncRNAs, as well as miRNA and mRNAs with calculation of Spearman correlation coefficient (SCC); 2) positive correlation between lncRNAs and mRNAs with calculation of Spearman correlation coefficient (SCC); 3) enrich degree of lncRNA binding to miRNAs and miRNAs binding to mRNAs. Relationships among ceRNA members and metabolites were visualized by the ggalluvial package in R3.6.1 software.

To verify the lncRNA-miRNA-mRNA regulatory chains, the transient coexpression experiments were performed in leaves of Nicotiana benthamiana according to Lu et al. (Lu et al., 2019). Precursors of miRNA (pre-miRNA), lncRNAs, and mRNAs were cloned from mixture of cDNA of three-stages flowers. For miRNAs that could not obtain pre-miRNAs in the cDNA, we performed overlap PCR to ligate these mature miRNAs into the vector. Pre-miRNAs and lncRNAs were ligated into pCAMBIA2300 at SmaI and XbaI sites. mRNAs were ligated into pCAMBIA1300 at XbaI and KpnI sites. The vector harboring pre-miRNAs, lncRNAs, and mRNAs were then transformed into Agrobacterium tumefaciens GV3101 strain. Equal amounts of agrobacterial cell cultures containing miRNAs, lncRNAs, and mRNAs were mixed respectively. Then, the mixture was infiltrated into leaves of N. benthamiana. After being incubated in dark for two days, the infiltrated tobacco leaves were observed under laser scanning confocal microscopy (LSCM) with EGFP = 780, and collected for total RNA extractions (CW2598, CWBIO, Taizhou, China) and qRT-PCR (Takara, Dalian, China). Tobacco L23 was used as an internal reference (Liu et al., 2012). The primers were listed in Table S1 .

Prior to conducting differentially expressed analyses, we performed basic statistical analyses on the sequencing data, such as examining correlations between samples, evaluating expression distributions, and identifying types of lncRNAs. Sample correlation was determined based on the expression level of lncRNAs ( Figure 1A ). Correlations among samples of each stage had a relatively high coefficient, indicating these samples meet our analysis requirements. A relatively high average expression level was observed at the late stage ( Figure 1B ). In order to increase the reliability of the results, we used CNCI and CPC2 to identify lncRNAs in sequenced samples. A total of 4,508 lncRNAs were identified by CNCI and a total of 6,567 lncRNAs were identified by CPC2. To ensure greater accuracy, we took the intersection of the results from both software and obtained 3,955 lncRNAs ( Figure 1C ). Five types of lncRNA were identified, including sense lncRNAs, antisense lncRNAs, intronic lncRNAs, bidirectional lncRNAs, and intergenic lncRNAs. It should be noted that some lncRNAs in the graph cannot be assigned a precise category. Among the identified lncRNA, the proportion of antisense lncRNAs was the highest, while that of intronic lncRNAs was the lowest ( Figure 1D ). Overall, the lncRNAs were suitable for further analysis.

For all lncRNAs identified, we performed differential expression analyses among different stages ( Figure 2 ). There were 190 up-regulated and 228 down-regulated lncRNAs in the McII vs McI group ( Figure 2A ; Table S2 ). In the comparison between the McIII and McI group, we identified 552 up-regulated and 602 down-regulated lncRNAs ( Figure 2B ; Table S3 ). 426 lncRNAs were up-regulated and 454 lncRNAs were down-regulated in the McIII vs McII group ( Figure 2C ; Table S4 ). Differential expression analyses results indicated huge differences of lncRNA expression levels between McIII and the other two stages. Then, we gave a statistics among the three comparison groups using Venn method which yieled 67 common lncRNAs ( Figure 2D ). As lncRNAs play a variety of roles during the different stages of development, we used the union of three comparison groups to conduct further analysis and identify the roles of lncRNAs during the different stages of flower development.

It has been reported that lncRNAs interact with mRNAs in three distinct ways including cis, trans, and antisense interactions (Nie et al., 2012). Cis model means that lncRNAs could induce reconstruction and histone modification to influence the combination of transcript factors with promoter and enhancer or combine with transcription elements directly to affect the expression of the protein-encoding genes. In the antisense model, it is assumed that lncRNAs could directly bind to mRNAs to influence alternative splicing of mRNAs or to decay mRNAs into siRNAs, which is closely related to the stability of mRNAs. The trans model is often able to be judged by the opposite expression trends between lncRNAs and their mRNA counterparts. LncRNAs could sometimes bind to proteins translated by mRNAs, altering their activity, structure, or position to influence their corresponding mRNAs. Here we also confirmed three functional models for differentially expressed lncRNAs and then GO and KEGG analyses were performed for their corresponding mRNAs ( Figures 3A, B ; Tables S5-S10 ). In Figures 3A, B , some GO terms associated with flower development and color formation were chosen for visualization, including alkaloids, pigments, flavonoids, hormones, and flower development-related terms or path. As a result of the three acting models, the trans model enriched more genes than the other two functional models. However, we could not find the antisense models of these DElncRNAs in the flavone biosynthetic process of the result of GO analysis as well as the isoflavonoid biosynthetic process, the folate biosynthesis, and the brassinosteroid biosynthetic process of the results of KEGG analysis ( Figures 3A, B ). Similarly, there was also lack of cis models of these DElncRNAs in the Tropane Alkaloid biosynthetic process and the Alkaloid metabolic process of the results of GO analysis, and the isoflavonoid biosynthesis and Anthocyanin biosynthesis of the results of KEGG analysis ( Figures 3A, B ).

To further investigate the potential roles of lncRNAs in regulating gene expression by antisense and cis models, we built lncRNAs-target genes networks ( Figures 3C, D ). In the cis model, we identified 26 lncRNAs that function on genes involved in flavonoid biosynthesis. Of these, three lncRNAs, lnc-MSTRG.23628.1, lnc-MSTRG.5535.1, and lnc-MSTRG.23626.1, were related to anthocyanin biosynthesis, while five lncRNAs, lnc-MSTRG.14710.1, lnc-MSTRG.26825.1, lnc-MSTRG.28377.1, lnc-MSTRG.6553.1, and lnc-MSTRG.30768.1 were related to flavone biosynthesis. In addition to flavonoid biosynthesis, we found one lncRNA in carotene biosynthesis and three lncRNAs in carotenoid biosynthesis ( Figure 3C ). In contrast to the cis model, we observed fewer lncRNAs involved in flavonoid biosynthesis in the antisense model, with only three lncRNAs, lnc-MSTRG.15809.1, lnc-MSTRG.12803.1, and lnc-MSTRG.30533.1, were detected. No lncRNAs involved in flavone biosynthesis process were identified ( Figure 3D ). Overall, our findings suggest that there are more cis models than antisense models involved in regulating flower development and color formation in M. candidum. The antisense model involves lncRNAs forming complementary relationships with genes. To investigate this further, we examined several genes, McMADS4a, 4b, McMYB1,17, McCHSJ, and McCHI1, along with their corresponding lncRNAs to assess their complementary relationship. Our results showed a high degree of complementarity between the lncRNAs and their target genes ( Figure 3E ).

In order to confirm the accuracy of our lncRNA-seq data, we performed qRT-PCR assays on 12 selected lncRNAs. The expression patterns of most of these lncRNAs were found to be consistent with the results obtained from the sequencing analysis. Specifically, lnc-MSTRG.17619.1, lnc-MSTRG.29477.1, and lnc-MSTRG.34245.1 were highly expressed during the McI stage, while lnc-MSTRG.10215.1, lnc-MSTRG.11402.2, lnc-MSTRG.223785.1, lnc-MSTRG.24719.1, lnc-MSTRG.29558.1, and lnc-MSTRG.29930.1 were highly expressed during the McII stage Figure 4A . Additionally, lnc-MSTRG.28008.2 and lnc-MSTRG.11415.1 were highly expressed during the McIII stage, while lnc-MSTRG.24133.1 was highly expressed in both the McII and McIII stages Figure 4A . We found that the results obtained from the lncRNA-seq analysis and the qRT-PCR experiments were positively correlated, with a slope of 1.147 and an R value of 0.727 for the McII vs McI group, and a slope of 0.968 and an R value of 0.645 for the McIII vs McII group (as shown in Figure 4B ). Based on these findings, we concluded that the lncRNA-seq analysis was valid and accurate.

We investigated the expression pattern of these lncRNAs during different stages of the flower development process. In the McII vs McI group, a total of six DElncRNAs were confirmed, of which four were down-regulated and two were up-regulated ( Figures 5A, D ). Ten DElncRNAs were differentially expressed between the McIII and McI groups, of which three were downregulated and seven were upregulated ( Figures 5B, E ). Five DElncRNAs related to flower color formation were found between McIII and McII, including one that was down-regulated and four that were up-regulated ( Figures 5C, F ). During the transition from the second to the third stage of flower development, lnc-MSTRG.10215.1 and lnc-MSTRG.24133.1 were up-regulated and exhibited an upward trend in expression, indicating their corresponding competing endonous mRNAs were key regulator in regulating this period of development. Conversely, lnc-MSTRG.34245.1, lnc-MSTRG.29477.1, and lnc-MSTRG.17619.1 showed a decline in expression from the second stage onwards ( Figures 5A, B ), indicating their corresponding competing endonous mRNAs were less involved in this process. It was observed that lnc-MSTRG.11415.1 and lnc-MSTRG.29558.1 began to be highly expressed in the third period (McIII), suggesting that these two lncRNAs will play a critical role in flower development after color formation ( Figures 5B, C, E, F ). There is only a low expression of lnc-MSTRG.11402.2 in the second stage, showing that its target gene plays a negative role in flower coloration process.

To further clarify the regulatory mechanism of the aforementioned lncRNAs, we first examined the supplementary relation between lncRNAs and miRNAs, miRNAs and target genes ( Figure S1 ). Subsequently, we constructed a ceRNA regulatory network, which included lncRNAs, miRNAs, mRNAs, and metabolites for flower formation-related pathways ( Figure 6A ). As shown in the graph, lnc-MSTRG.10215.1 and lnc-MSTRG.11402.2 were found to regulate McF3H through interacting with miR5207, ultimately influencing dihydromyricetin (DHM), dihydroquercetin (DHQ), and pinobanksin (PBA). Dihydroquercetin (DHQ), also known as 3, 5, 7, 3, 4 - pentahydroxy flavanone or taxifolin, is a bioactive flavonoid that is considered one of the rarest and most effective natural antioxidants (Yu et al., 2021). Both dihydroquercetin and dihydromyricetin are colorless dihydroflavones that can be reduced to leucoanthocyanidins under the action of the dihydroflavonol-4-reductase (DFR) enzyme. Subsequently, under the catalysis of downstream enzymes, the leucoanthocyanidins undergo a process of the transformation into orange pelargonidin, reddish-purple cyanidin, and violet-blue delphinidin (Tanaka and Brugliera, 2013; Liu et al., 2019). Pinobanksin, one of the most common phenolic constituents of pine heartwood, was first isolated from Pinus Banksiana (Erdtman, 1944). When combined with zinc ion, this compound is specific to 3-hydroxyflavanones to produce much deeper colors during the reducing process. Pinobanksin usually produces orange-red when reduced with magnesium or zinc and hydrochloric acid (Lindstedt, 1950). P-coumaroyl Shikimic acid is a by-product of the anthocyanin pathway, derived from 4-coumaroyl-CoA by the formation of an ester bond with shikimic acid catalyzed by Shikimate O-hydroxycinnamoyl transferase (HCT) (Paliyath et al., 2009; Li et al., 2021). In the presence of the enzymes C4H1 and C3H3, this compound can be catalyzed again to transform into caffeoyl shikimic acid (Chen et al., 2011). Both p-coumaroyl Shikimic acid and caffeoyl shikimic acid are key compounds in Phenylpropanoid biosynthesis during lignin formation. The Sankey graph showed that miR3704, miR480, and miR1220 had cleavage effects on McHST, McAHT1, and McPHT3, respectively. Additionally, we discovered that lnc-MSTRG.24133.1 could adsorb miR3704 to maintain a balance between miRNA and target gene. Three lncRNAs, lnc-MSTRG.17619.1, lnc-MSTRG.29477.5, and lnc-MSTRG.34245.1, were found to bind to miR480 while another set of three lncRNAs, lnc-MSTRG.23785.1, lnc-MSTRG.29558.1, and lnc-MSTRG.29930.1 were found to bind to miR1220. These results indicate that the three lncRNAs that ceRNA regulatory mechanism play an important role in the formation of p-coumaroyl Shikimic acid, which ultimately determine flower color. Three specific genes, namely McUGT88B1, McUGT88F3, and McCYP98A2, were found to have no associated with metabolites. This may indicate that the three genes related to ceRNAs do not significant affect metabolism. It is also possible that annotation of the metabolite database is incomplete.

We conducted a correlation analysis to determine the relationship between ceRNA members and metabolites in specific pathways ( Figure 6B ). Each ceRNA regulatory routine showed a negative correlation between miRNAs and both lncRNAs and mRNAs, whereas a positive correlation was observed between lncRNAs and mRNAs, indicating the presence of a true competitive mechanism between lncRNAs and mRNAs. Lnc-MSTRG.28008.2 and McASP2 were negatively correlated with L-phenylalanine, however, miR8131-3p was positively correlated with L-phenylalanine, indicating that this ceRNA regulatory chain is responsible for L-phenylalanine synthesis ( Figure 6B ). It is known that both the flavonoid and anthocyanin pathways originate from the amino acid L-phenylalanine, which is deaminated by the action of phenylalanine ammonia-lyase (PAL) to produce trans-cinnamic acid and ammonia (Fasoula et al., 1995). Lnc-MSTRG.29477.1, lnc-MSTRG.29558.1, McAHT1, and McPHT3 were positively correlated with p-coumaroyl Shikimic acid, whereas, miR480 and miR1220 were negatively correlated with p-coumaroyl Shikimic acid, indicating that these two ceRNAs regulatory chains were key regulator in the synthesis of p-coumaroyl Shikimic acid. However, the rest of ceRNA regulatory networks did not display obvious putative correlations, which could be attibuted to more complex regulatory mechanisms like multiple ceRNA chains, histone modification, RNA methylation modification process, etc. Further studies are needed to gain a better understanding of these complex regulatory mechanisms.

To further validate the identified lncRNA-mediated ceRNA regulatory network, we employed pCAMBIA1300 harboring target genes, and pCAMBIA2300 harboring lncRNAs and miRNAs to performed tobacco injection and qRT-PCR experiment. Ordinarily, miRNAs have cleavage effects on the target genes, when lncRNAs bind with miRNAs, the cleavage effects will be weakened. We selected 7 lncRNA-mediated ceRNA regulatory routines which included 6 genes, 5 miRNAs, and 6 lncRNAs to validate the ceRNA regulatory mechanism in tobacco ( Figure 7 ). In accordance with our anticipation, when only injecting the miRNAs and target genes into leaves of tobacco, the fluorescence signal became weak ( Figure 7 ), and the expression level of target genes begun to decrease. When we added the lncRNAs in corresponding ceRNA chains, the fluorescence signal was regained, and the expression level of the target gene also began to increase ( Figure 7 ). These results provided evidence for the involvement of lncRNA-mediated ceRNA regulation in M. candidum.

Based on the aforementioned results, we summarized some key lncRNA-mediated ceRNA regulatory relationship during flower development ( Figure 8 ; Figure S2 ). LncRNAs could bind with a post-transcriptional product of genes directly to degrade them into siRNAs (small interference RNAs). Another common path of lncRNAs is to compete with mRNAs to bind to miRNAs, thereby easing the cleavage effect of miRNA and facilitating more efficient translation of mRNAs ( Figure 8A ). We observed some ceRNA regulatory mechanisms operating in the flavonoid pathway. Specifically, lnc-MSTRG.11402.2-miR5207-McF3H formed a regulatory chain that acts as a catalyst for the transformation from pinocembrin to pinobanksin ( Figure 8B ). Lnc-MSTRG.28008.2-miR8131-3p-McASP2 catalyzes the conversion of phenylpyruvate into L-phenylalanine. Lnc-MSTRG.29930.1, lnc-MSTRG.23785.1, and lnc-MSTRG.29558.1 could bind with miR1220 to regulate McPHT3 genes; lnc-MSTRG.24133.1 could affect McHST by interacting with miR3704; Lnc-MSTRG.34245.1, lnc-MSTRG.29477.1, and, lnc-MSTRG.17619.1 completed with miR480 to bind to McAHT1. All these three combinations together influenced the metabolic process from p-Coumaroyl-CoA to p-Coumaroyl shikimic acid ( Figure 8B ). Out of all the observed lncRNA-miRNA-mRNA-metabolites regulatory chains, lnc-MSTRG.29930.1, McPHT3, and p-Coumaroyl shikimic acid were lowly expressed in the second stage, whereas miR1220 was highly expressed in the second stage, indicating the presence of a ceRNA regulatory mechanism. It was found that lnc-MSTRG.24133.1, McHST, and p-Coumaroyl shikimic acid were highly expressed in the third stage, but miR3704 was lowly expressed, suggesting that this ceRNA regulatory chain also plays an important role in the biosynthesis of p-Coumaroyl shikimic acid ( Figure S4 ). We also found some irregular ceRNA regulatory chains, for example, lnc-MSTRG.11402.2-miR5207-McF3H, which was thought to be involved in the synthesis of dihydroquercetin, dihydromyricetin, and pinobanksin. However, we did not observe any corresponding expression trends for these metabolites; lnc-MSTRG.17619.1, lnc-MSTRG.29477.1, lnc-MSTRG.34245.1-miR480-McAHT1-p-Coumaroyl shikimic acid displayed putative trends in the first stage, whereas, in the third stage, they were inconsistent with the putative results ( Figures S3 , S4 ). These irregular phenomena may be caused by other post-transcription regulating or modifications mechanisms.