Lycoris radiata (L’Her.) Herb. plants were grown in Experimental Plantation of Institute of Botany, Jiangsu Province and Chinese Academy of Sciences, Nanjing, China. According to the studies reported previously (Yue et al., 2019; Park et al., 2021), three biological replicates of L. radiata flowers were sampled at four development stages, which were FB (floral bud stage), FL1 (partially opening flower stage), FL2 (fully opened flower stage) and R (senescent flower stage), as shown in Figure 1A. Each biological replicate was taken from petals of five flowers and pooled together. For gene expression analysis, different L. radiata tissues, including scape, stamen, pistil, flower stalk, and petal samples were obtained during flowering time, while leaf, root, as well as bulb samples from the same plants were obtained during the vigorous vegetative growth stage. The fresh samples were harvested and instantly frozen in liquid nitrogen, then kept at –80°C until use.

Extraction and determination of anthocyanins of L. radiata flowers was performed following the protocol of Mehrtens et al. (2005) with minor modifications. Briefly, approximately 0.1 g fresh petals were ground in 1 mL of acidic methanol (0.1 mol L^–1 HCl) and then incubated overnight in the dark at 4°C with gentle shaking. After centrifugation for 10 min at 12,000 rpm, the supernatant was diluted four times with acidic methanol and the absorbance was measured at 530 and 657 nm using a UV-1600 spectrophotometer (SHIMADZU, Kyoto, Japan). The concentration of anthocyanins was calculated using the following formula: Q_Anthocyanins = (A₅₃₀ – 0.25 × A₆₅₇) × FW^–1, where Q_Anthocyanins is the amount of anthocyanins, A₅₃₀ and A₆₅₇ is the absorption at the indicated wavelengths and FW represents the weight of the fresh sample [g].

Total RNA was extracted with the mirVana miRNA isolation kit (Thermo Fisher Scientific, Waltham, MA, USA) following the manufacturer’s protocol. The quality and quantity of the RNA were examined by the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Samples with RNA Integrity Number (RIN) ≥ 7 were subjected to cDNA library construction using the TruSeq Stranded mRNA LTSample Prep Kit (Illumina, San Diego, CA, United States). Sequencing of the cDNA libraries was done on the Illumina sequencing platform (Illumina HiSeq™ 2500) by Shanghai OE Biotech. Co. Ltd. (Shanghai, China). Reads were cleaned by removing adapters, as well as low-quality and ambiguous regions, then subjected to de novo assembly using the Trinity software (Grabherr et al., 2011).

Alignment of the assembled unigenes was done against public databases including National Center for Biotechnology Information (NCBI) non-redundant protein (Nr) and nucleotide (Nt) database, the Swiss-Prot protein database, Gene Ontology (GO) database, Protein Family (Pfam) database, Kyoto Encyclopedia of Genes and Genomes (KEGG) database, Eukaryotic Ortholog Groups (KOG) database, and eggNOG (evolutionary gene genealogy: Non-supervised Orthologous Groups) database.

The expression level of unigenes was calculated using fragments per kilobase per million fragments mapped (FPKM) method (Mortazavi et al., 2008). Identification of DEGs among samples at four development stage was done using the DESeq2 package implemented in R software, with cutoff values of |log2 (fold change)| > 1 and p-value < 0.05 algorithms (Young et al., 2010). To visualize the differential expression profiles, we generated a heatmap for the Trimmed Mean of M-values (TMM) normalized against FPKM via the pheatmap package in R.

To predict TFs involved in color formation of L. radiata, we utilized the getorf database (mini-size 150) to find the open reading frame (ORF) (Rice et al., 2000) and then used the HMM search database (version 3.0) to align the ORFs to the TF protein domain. The aligned sequences were described according to the TF families available on the PlantTF database version 3.0 (Zhang et al., 2011). Moreover, the Pearson’s correlation coefficient (PCC) between these differentially expressed TFs, structure genes and total anthocyanin content of samples was calculated. The TFs with |PCC| > 0.8 were selected for subsequent analysis. The TF expression data, which included expression levels for MYB, bHLH, WD40, and the DEGs identified in the flavonoid biosynthetic pathway, was screened using blastx software, with an e-value of 1e-10. The target gene sequence was aligned to the protein sequence of the reference species contained in the string database, and the protein interaction relationship of the reference species was used to construct an interaction network. Network visualization for the interaction network related to DFR and DEGs was performed using Cytoscape version 3.6.1.

Cloning of LrDFR1 was based on putative ORFs of unigenes from the RNA-seq database. Primers (Supplementary Table 1) were synthesized for ORF sequence amplification using Tks Gflex™ DNA Polymerase (Takara, Dalian, China) from L. radiata petal cDNA. Reaction conditions were: 5 min of 95°C, 35 cycles for 30 s at 94°C, 30 s at 60°C, 1 min at 72°C, with extension at 72°C for 10 min. PCR products were cloned into pMD19-T simple vectors (Takara, Dalian, China). Afterward, those T-vectors were transferred into DH5α competent cells (Takara, Dalian, China) for amplification. The overexpression vectors of LrDFR1 were established by linking their ORFs into a linear plant transformation vector, pBinGFP4, using the One Step Cloning Kit (Vazyme, Nanjing, China). Then the 35S:LrDFR1 recombinant vectors were transformed into Agrobacterium tumefaciens EHA105 competent cells.

The pBinGFP4 vector with LrDFR1-GFP was transformed into the Agrobacterium tumefaciens strain EHA105, and transferred into Nicotiana benthamiana epidermal cells (Sheludko et al., 2007). Cultivation of the transformed N. benthamiana leaves was done for 2–6 days. For co-localization with membrane-localized marker, 35S: PIP2;1-mCherry construction was used (Huang et al., 2019). Assessment of transformed N. benthamiana epidermal cells was observed with confocal laser scanning microscopy (Zeiss LSM780 META, Jena, Germany). For staining of the nuclei, 10 mg/mL 4′6-diamidino-2-phenylindole (DAPI) was infiltrated into N. benthamiana leaves 6 h before observation.

Staining of proanthocyanidin was conducted as described by An et al. (2015). Briefly, light-treated N. benthamiana leaves were decolorized in a solution of ethanol: glacial acetic acid (3:1). A dimethylaminocinnamaldehyde (DMACA) reagent staining solution (Sigma-Aldrich, St. Louis, MO, United States) was then added for staining.

The A. tumefaciens harboring 35S:LrDFR1-GFP construct and the control pBinGFP4 vector were prepared for injecting into Lycoris petals, respectively. The recombinant Agrobacterium strains were cultured in YEB broth containing 50 μg mL^–1 kanamycin and incubated at 28 °C. Then, the collected recombinant Agrobacterium strains were resuspended to OD₆₀₀ of 0.6 in a buffer with 10 mM 2-(4-Morpholino) ethanesulfonic acid, 10 mM MgCl₂, and 120 μM acetosyringone. Transformed Lycoris petals were stored for 48 h in the dark after which they were transferred to a phytotron at a constant photon flux density of 100 μmol m^–2 s^–1. With 5 days cultivation, Lycoris petals were obtained for anthocyanin level assessment and RNA extraction.

For validating gene expression using qRT-PCR, 32 unigenes associated with anthocyanin biosynthesis and phytohormone metabolism were randomly selected (Supplementary Table 2). Total RNA isolation was conducted by using the RNAprep Pure Plant Kit (Tiangen, Beijing, China). First-strand cDNA was synthesized with TransScript One-Step gDNA Removal and cDNA Synthesis SuperMix kit (Takara, Dalian, China), and the extracted RNA was used as template according to manufacturer’s instructions. A list of gene-specific primers is provided in Supplementary Table 1. The quantified expression levels of the tested genes were normalized against the house keeping genes TIP41-like protein (TIP41) according to previous study on L. aurea (Ma et al., 2016). qRT-PCR assays were conducted by the SYBR Premix Ex Taq™ II kit (Tli RNaseH Plus) (Takara, Dalian, China) in a Bio-Rad iQ5 Gradient RT-PCR system. Reaction conditions were: 30 s of denaturation at 95°C and 40 amplification cycles (5 s at 95°C, 30 s at 60°C). Calculation of relative target gene expression levels was done using the 2^–ΔΔCt method (Livak and Schmittgen, 2001). Experiments were conducted using three independent biological and three technical replicates.

Statistical analyses were done by SPSS version 10.0 software (IBM Corporation, Armonk, NY, USA). The significant difference among sets of data was determined by one-way analysis of variance (ANOVA) with Duncan’s multiple range test (p < 0.05) or a significant t-test (^**p < 0.01, *p < 0.05). All the results are presented as the mean ± standard deviation (SD).

During the red flower development of L. radiata, petals underwent a rapid color change from slight red to brilliant red (Figure 1A). At the flower bud (FB) stage, a slight red color was observed, then the color intensity was significantly increased with rapid elongation of petals in FL1. Subsequently, the intensity of L. radiata decreased at FL2 and R stages (Figure 1A). We thus investigated the changes of anthocyanin contents in L. radiata at four different petal development stages. Notably, anthocyanin content at FL1 stage was significantly higher than that of FB, FL2 and R stages (Figure 1B), suggesting that changes in anthocyanin levels could be the main reason for red color formation of L. radiata.

To further study the molecular mechanism of L. radiata petal coloring during flower development, twelve libraries were established using samples at four flower development (FB, FL1, FL2, and R) stages (three biological replicates for samples at each development stage), and a total of 644.93 million raw reads as well as 96.73 Gb raw bases were obtained. After eliminating the adaptor, poor-quality sequences, and ambiguous reads, 634.09 million clean reads and 89.86 Gb clean bases were retrieved from 12 samples (Supplementary Table 3). The quality score above 30 (Q30) of each library was 93.75–94.91%, and GC percentages ranged from 44.99–46.51% (Supplementary Table 3). By using Trinity software, the de novo assembly of 12 petal transcriptomes totally generated 87,584 unigenes with an average length of 942 bp (Supplementary Table 4). Sequence length distribution showed that 27,073 (30.91%) unigenes had a mean length ≥ 1000 bp (Supplementary Figure 1 and Supplementary Table 4). The N50 was determined to be 1,334 bp, which indicated that the quality of sequence assembly was good.

FPKM values were used to estimate the transcription levels of unigenes. More than 50.0% of unigenes had FPKM values above 1 (Supplementary Figure 2). In addition, the use of relative unigene expression obtained from FPKM for principal component analysis (PCA) showed 52.10% variability among the samples (Supplementary Figure 3). Moreover, heatmap coefficient matrix analysis of the samples based on the FPKM values showed that most biological replicates (except FB3 sample) exhibited similar expression patterns, indicating relatively high reliability of our sequencing data (Supplementary Figure 4).

All of the unigenes were annotated by BLAST search against the public databases. The results revealed that 40,974 (46.78%), 29,476 (33.65%), 37,487 (42.8%), and 22,318 (25.48%) unigenes were annotated to the Nr, Swiss-Prot, eggNOG, and Pfam databases, respectively. Taken the entire public databases together, a total of 41,534 (47.42%) unigenes could be successfully annotated (Supplementary Table 5). To elucidate their main biological functions, GO, KOG, and KEGG pathway assessments were also performed (Supplementary Figures 5–7). Consequently, 27,296 (31.17%) unigenes were assigned into three main categories including “biological process” (BP), “cellular component” (CC), and “molecular function (MF),” which could be further distributed under 50 GO terms (Supplementary Figure 5). In addition. 15,122 (17.27%) unigenes were associated with 126 KEGG pathways, and category ‘Metabolism’ (6187 unigenes) was the most abundant (Supplementary Figure 6). Moreover, the KOG analysis showed that 23,858 (27.24%) annotated unigenes were assigned into 25 classes (Supplementary Figure 7).

To identify the key DEGs involved in L. radiata petal color transitions, six pair-wise comparison groups (FL1 vs. R, FB vs. R, FL2 vs. R, FL2 vs. FL1, FL1 vs. FB, and FL2 vs. FB) were conducted (Figure 2). A total of 38,798 DEGs were identified among all samples based on a |log₂ fold change| >1 at p < 0.05. Among these comparison groups, the largest abundance of DEGs (23,202) was found between FB and R libraries, of which 10,958 and 12,244 genes were down-regulated and up-regulated, respectively (Figure 2A). Conversely, the smallest abundance of DEGs (9,057) was recorded between FL2 and FL1 libraries, with 5,033 and 4,024 of them down-regulated and up-regulated, respectively (Figure 2A). Furthermore, the overlap DEGs among the six comparison groups were screened. The results indicated that 38 genes were differentially expressed among all the comparisons, which indicated that these DEGs might have key functions in the color expression of different petals (Figure 2B and Supplementary Table 6).

To elaborate the functions of DEGs and identify genes involved in regulating anthocyanin accumulation in L. radiata, all the DEGs were firstly subjected to GO analyses, and 14,555 of the 38,798 DEGs were assigned to GO annotations (Supplementary Table 7 and Figure 3A). In the biological process category, most of the DEGs were mapped to ‘cellular process’ (9,476, 20.23%), ‘metabolic process’ (8,073, 17.23%), and ‘response to stimulus’ (4,181, 8.92%) terms. In the cellular component category, more than 63.08% of DEGs were enriched in ‘cell,’ ‘cell part’ and ‘organelle’ terms, but for molecular function, nearly 86.13% of DEGs were mapped to ‘catalytic activity’ and ‘binding’ terms (Figure 3A). For the KEGG annotation results, 7,631 DEGs among all samples were also mapped to 126 KEGG pathways (Supplementary Table 7). Comparisons across the samples at four petal development stages revealed significant enrichment of DEGs in ‘flavonoid biosynthesis,’ ‘phenylpropanoid biosynthesis,’ ‘Tropane, piperidine and pyridine alkaloid biosynthesis,’ ‘terpenoid backbone biosynthesis’ as well as ‘plant hormone signal transduction’ pathways (Figure 3B and Supplementary Figure 8). For example, the significantly enriched KEGG pathway term ‘Tropane, piperidine and pyridine alkaloid biosynthesis’ was shared in all the comparisons. The ‘flavonoid biosynthesis’ pathway was enriched in FL1 vs. R, FB vs. R, FL2 vs. R, FL2 vs. FL1, and FL2 vs. FB, but not in FL1 vs. FB. In addition, the ‘plant hormone signal transduction’ pathway was enriched in FB vs. R, FL2 vs. FB, and FL2 vs. FL1 (Supplementary Figure 8).

To elucidate the molecular basis underlying difference in anthocyanin biosynthesis among the four flower development stages in L. radiata, DEGs involved in the anthocyanin synthesis pathway were identified. The results revealed that 56 DEGs were enriched in the anthocyanin synthesis pathway, including PAL, C4H, 4CL, CHS, CHI, F3H, F3′H, DFR, ANS, UFGT, FLS, ANR, and LAR (Figure 4A). Moreover, the Pearson’s correlation coefficient between the expression level of these DEGs and the total anthocyanins content was further calculated (Figure 4B). The results showed that 23 DEGs negatively regulated anthocyanin synthesis, whereas 33 DEGs positively regulated the anthocyanin synthesis. Among them, the expression level of two DEGs, namely LrDFR1 (DN43960) and LrDFR2 (DN42380) indicated a significant positive correlation with the total anthocyanins content in petals during the flower development stages, while LrFLS (DN37334) indicated a significant negative correlation with the total anthocyanins content (|PCC| > 0.8, Figure 4B and Table 1), suggesting that these three DEGs may have an essential role in anthocyanin accumulation.

Transcription factors were subsequently predicted to whether modulate anthocyanin accumulation and biosynthesis in L. radiata petals during flowering development stages. In this study, a total of 1,631 TFs were identified by searching the TF database. The classified results indicated that most of these TFs belonged to the MYB, C2C2, AP2/ERF, C2H2, and bHLH family (Supplementary Figure 9). Furthermore, the differentially expressed TFs (721) were characterized by analyzing their FPKM values (Supplementary Figure 9). Importantly, co-expression modules of these 721 TF DEGs were analyzed with Short Time-series Expression Miner (STEM) software. In all, six clusters of 272 TF DEGs were detected (Figure 5A and Supplementary Table 8).

On the other hand, by calculating the PCC between the expression level of 721 TF DEGs and the total anthocyanins content, 27 TFs genes (|PCC| > 0.8) involved in the accumulation of anthocyanins were identified, including 10 positive regulators and 17 negative regulators (Table 1). These 10 positive regulators, including MYB (1), AP2/ERF (1), bHLH (1), bZIP (1), NAC (1), NF-X1 (1), and Trihelixs (4) genes, likely act to improve anthocyanin synthesis during L. radiata petal development stages. However, the 17 negative regulators, including Alfin-like (1), AP2/ERF (1), GATA (1), GRFs (2), bHLHs (2), MYBs (2), C2H2s (2), C3Hs (4), MADS (1), and NAC (1), might act as repressors in L. radiata anthocyanin biosynthesis (Table 1). Notably, 11 of 17 negative TF regulators (subclass 4, 5, and 6) and one positive regulator (subclass 2) were enriched in TF co-expression modules (Supplementary Table 8).

Previous studies have reported that bHLH, MYB and WD40 TFs regulate anthocyanin biosynthesis thereby activating or repressing transcription of anthocyanin structural genes. We then performed unigenes regarding to MYB, bHLH and WD40, as well as 56 DEGs involved in anthocyanin biosynthesis (Figure 5) to analyze their interaction network and hope to identify the hub TF genes that could affect anthocyanin biosynthesis pathway. The results showed that four DFRs, four MYBs, two WD40s, two 4CLs, one F3′H, one UFGT, one CHS, one ANS, one FLS, and one CHI were selected as hub genes based on their connection position in the network modules, expression pattern and functional annotation (Supplementary Table 9a and Supplementary Figure 10). Furthermore, those genes (shown in Supplementary Table 9a) and 27 key TF genes (Table 1) were selected to build the interaction network for further analysis. Among them, LrDFR1 (DN43960) and LrFLS (DN37334) could be regarded as key hub genes for participating anthocyanin biosynthesis. Two MYBs (DN45447 and DN33872), two NACs (DN39353 and DN39797), one C3H (DN42291), and one GATA (DN30983) TF genes were identified as hub genes in regulating anthocyanin biosynthesis (Figure 5B and Supplementary Table 9b). The above results indicate that these eight genes may play essential roles in anthocyanin synthesis in L. radiata during petal development.

To validate the accuracy and transcription profiles revealed by the RNA-seq data, 32 unigenes were selected for qRT-PCR assays. The relative expression levels of these 32 genes were normalized to the expression of LrTIP41, and compared with the RNA-Seq data, as shown in Figure 6A. Further linear regression analysis revealed that the expression levels of these genes were well correlated with the RNA-Seq results (Figure 6B, R² > 0.76), indicating that the RNA-seq data were credible and accurate.

In this study, we cloned LrDFR1 gene (DN43960) from L. radiata. The full-length cDNA of LrDFR1 is 1113 bp in length and it encodes a 370 amino acid protein with a molecular weight of 41.67 kDa (Supplementary Table 10). The deduced amino acid sequence of LrDFR1 revealed a high similarity with DFR proteins from Agapanthus praecox (75.33%), Muscari armeniacum (74.74%), and Hyacinthus orientalis (72.72%) (Figure 7A). Multiple amino acid sequence alignments showed the highly preserved NADPH-binding motif (VTGAAGFIGSWLIMRLLERGY) (Gang, 2005) and the substrate-binding domain (T128–K154) (Johnson et al., 2001) in the LrDFR1 sequence (Supplementary Figure 11). qRT-PCR was then performed to assess whether expression patterns of LrDFR1 in different tissues and flower development stages were coincided with anthocyanin accumulation in L. radiata. LrDFR1 was found to be expressed in all tissues, with the highest expression levels in petals (Figure 7B). Moreover, expression levels of LrDFR1 were significantly increased from stage FB to stage R, peaking at stage FL1 (Figure 7C). These findings imply tissue-specific expression levels for LrDFR1, which is associated with anthocyanin accumulation in L. radiata petals.

Moreover, we transiently expressed LrDFR1 in tobacco epidermal cells to assess subcellular localization of LrDFR1. As shown in Figure 7D, the fluorescent signal of LrDFR1-GFP was localized into the nucleus, cytoplasm and cell membrane, while GFP was evenly distributed in the cell (Figure 7D and Supplementary Figure 12). To determine the roles of LrDFR1 in regulating anthocyanin as well as proanthocyanidin biosynthesis in L. radiata, an LrDFR1-overexpressing plasmid was transfected into Lycoris petals and tobacco epidermal cells (Figure 8A). Overexpression of LrDFR1 in tobacco and Lycoris petals markedly enhanced proanthocyanidin and anthocyanin accumulation (Figures 8B–E). To assess the effects of LrDFR1 on endogenous Lycoris petals genes that are involved in anthocyanin synthesis, the expression levels of CHS, CHI, F3H, F3′H, DFR, ANS, UFGT, and 3RT were determined (Figure 8F). Among them, the expressions of LrDFR1, ANS, UFGT and 3RT were significantly higher in LrDFR1-overexpressing plants than in control plants (Figure 8F). These results suggest that LrDFR1 may play important roles in anthocyanins biosynthesis of Lycoris petals.

Flowering plants exhibit a wide variation in their flora, foliage, and fruit colors, as a result of genetic factors and variations in environments. Flavonoids/anthocyanins, betalains and carotenoids are the major metabolites for coloration in plant reproductive organs (Griesbach, 2005; Tanaka et al., 2008). Most of the red, purple, and blue-colored flowers (such as red rose, lavender, and blue chicory) as well as fruits (such as berries, currants, and grapes) contained high anthocyanins content (Khoo et al., 2017). The genus Lycoris is used as a garden flower due to the colorful and special flowers, and the flower colors of Lycoris are diverse. For example, the flower color of L. radiata and L. rosea was red, that of L. aurea and L. chinensis was yellow. L. sprengeri and L. haywardii showed red and blue color, while L. longituba displays an exceptionally wide diversity of flower colors from purple, red, orange, to yellow (He et al., 2011). Similar to the flowers of other species, the petals of Lycoris are rich in anthocyanins, and their color formation are largely related to anthocyanins (He et al., 2011; Chun et al., 2013; Yue et al., 2019; Park et al., 2021). In this study, we determined the content of anthocyanins in the petals of L. radiata, and the results showed that the color intensity of the L. radiata petals was changed with the different anthocyanin contents. The anthocyanins increased then decreased during the flower development stages (Figure 1), which are similar to the results recently reported by Park et al. (2021).

To data, transcriptome sequencing is highly employed for predicting novel genes, gene function, and genome evolution for plant breeding and horticulture research (Rameneni et al., 2020). For example, transcriptome analysis has revealed the role of anthocyanin in flower color formation in several horticultural crops, such as Camellia sinensis (Zhou et al., 2020), “Tiny Padhye” (Lilium spp.) (Xu et al., 2017), lilies (Lilium spp.) (Suzuki et al., 2016), Magnolia sprengeri (Shi et al., 2014), Paeonia lactiflora (Zhao et al., 2014), Paeonia delavayi (Shi et al., 2015), and Silene littorea (Casimiro-Soriguer et al., 2016). For better understanding of petals color formation during flower development stages in L. radiata, a comparative transcriptomics analysis was carried out. The results showed that approximately 70.27 GB of high-quality data, and 87,584 unigenes were obtained. Further analyses, based on NR, Swiss-Prot, KEGG, KOG, GO, Pfam, and eggNOG databases, predicted 38,798 DEGs associated with a specific or general function (Supplementary Tables 4, 5).

The variations in floral coloration emanates from different processes, such as pathways competition, expression levels of structural genes involved in pigment formation, and mutations of structural or regulatory genes (Grotewold, 2006; Cui et al., 2021). In plants, phenylpropanoids represent a vital group of physiologically active secondary metabolites derived from phenylalanine, and anthocyanins, flavonols, isoflaconoids and flavonols have a similar metabolism pathway during their biosynthesis (Ferrer et al., 2008). KEGG pathway analysis showed that the ‘phenylpropanoid biosynthesis,’ ‘flavonoid biosynthesis,’ as well as ‘flavone and flavonol biosynthesis’ pathways were enriched between each two transcriptomes of L. radiata petals during flower development stages (Figure 3B and Supplementary Figure 8). Given that anthocyanin biosynthesis pathway is well known to modulate color formation in plants, we mainly focused on them as the candidate pathways to elucidate their involvement in petal/flower color formation in L. radiata. Subsequently, we identified the main functional genes participated in the anthocyanin biosynthetic pathway, and found that most of structural genes such as F3′H, UFGT, DFR, and FLS were elevated in L. radiata petals at FL1 and FL2 stages (Figure 4A). Therefore, these genes might have contributed to the increasing anthocyanin content in petals from the FB stage to the FL1 and FL2 stage, as evidenced in Figure 1. For example, three F3′H genes (DN41001, DN43758, and DN46768) were highly expressed in petals at FL1 and FL2 stages (Figure 4A). Another prominent gene, UFGT (DN44965), which glycolyzes anthocyanidin into anthocyanin (Xie et al., 2003), was also highly expressed in petals at FL1, FL2 and R stages, as compared to that of the samples at FB stage (Figure 4A). All of these genes were positively correlated with the biosynthesis of anthocyanins (Niu et al., 2010). Notably, two DFR genes (DN42380, DN43960) and one FLS (DN37334) (Table 1) were found to be highly associated with the total anthocyanins content (|PCC| > 0.8), suggesting they may have an essential function in the phenotypic expression of petal color (Figure 4B). In anthocyanin biosynthesis, DFR catalyze the reduction of dihydroquercetin to leucoanthocyanidins, and the level of DFR expression have been associated with flower color changes (Nakatsuka et al., 2003; Zhao et al., 2012). qRT-PCR also indicated that the hub gene LrDFR1 were mostly expressed the most in the FL1 samples (Figure 7C). Our results suggest that these enzymes may be the most important enzymes to catalyze anthocyanin biosynthesis in L. radiata petals.

Transcription factors play critical functions in flavonoid biosynthesis, by regulating expression of structural genes. For example, the class of TFs identified were previously implicated in regulation of petal color formation in roses (Li D. et al., 2020). Particularly, MYB-bHLH-WD40 complexes have been implicated in multi-level regulation of flavonoid biosynthesis (Gallego et al., 2018), whereas the R2R3-MYB family was shown to play a vital role in regulation of spatiotemporal expressions of genes involved in anthocyanin biosynthetic in plants (Gonzalez et al., 2008; Zhao and Tao, 2015). Besides, MYB-domain TFs are important mediators of anthocyanin accumulation and participate in colorations of various organs in horticultural as well as ornamental plants (Tang et al., 2017; Liu Y. et al., 2019; Xiang et al., 2019; Zhai et al., 2019; He et al., 2020; Sun et al., 2020; Wang et al., 2020; Zhong et al., 2020).

In this study, the most abundant TFs including AP2/ERF, bHLH, bZIP, C2C2, HSF, MYB, NAC, TIFY, and WRKY families were predicted (Supplementary Table 8). In addition, we employed a K-means clustering, as proposed earlier by Handhayani and Hiryanto (2015), which permitted the clustering of 272 TF unique genes among the samples (FB, FL1, FL2, and R) into six sub-clusters with some members in Cluster 2 associated with genes from the MYB and bHLH TFs (Figure 5A). Based on the expression level of TFs obtained from the transcriptome data, 27 TFs (Table 1) were found to highly associate with the total anthocyanin content (|PCC| > 0.8), and these TFs may have an essential function in the phenotypic expression of L. radiata petal color. Interestingly, among these TFs, three MYBs showed two different expression patterns. The expression level of two MYBs (DN45447, DN51496) was highest in FB, followed by FL1, FL2, and R, which was contrary to the total anthocyanin content trend. Conversely, the expression of LrMYB1 (DN33872) exhibited a similar trend to the total anthocyanin content in the L. radiata petals (Supplementary Figure 13), indicating that MYBs (DN45447 and DN51496) negatively regulated anthocyanin accumulation, whereas LrMYB1 (DN33872) was identified as one of the eight hub genes may positively regulate anthocyanin accumulation in L. radiata (Table 1 and Figure 5B).

Subsequently, two negatively correlated bHLHs (DN36174 and DN41224) and one positively correlated LrbHLH1 (DN48856) were identified (Table 1). In plants, MYB often forms protein complexes with bHLH and WD40 to participate in anthocyanin biosynthesis rather than regulate anthocyanin biosynthesis directly (Feng et al., 2020). In apple, MdMYB1, MdMYB9, MdMYB10, and MdMYBA act as positive modulators of anthocyanin biosynthesis, by activating the expressions of MdDFR and MdUF3GT (Takos et al., 2006; Ban et al., 2007; Espley et al., 2007; An et al., 2015). On the contrary, downregulation of MdMYB1 inhibits anthocyanin accumulation mediated by ethylene, abscisic acid (ABA), wounding, drought, and different light intensities (An et al., 2018, 2019, 2020a,b). Notably, our results also revealed a significant upregulation of LrMYB1 (DN33872) and LrbHLH1 (DN48856), of which the expression was positively correlated with LrDFR1, LrCHS, LrCHI, F3′H, LrUFGT and LrANS genes during petal development stages (Supplementary Figure 13). This is similar to that of LhMYB12-Lat, which has previously been associated with activation of accumulation of anthocyanin in lily petals (Yamagishi et al., 2014). In our co-expression networks, the module that was positively correlated with anthocyanin contents and modules negatively correlated with anthocyanin content were identified. Overall, whether these MYB TFs interact with bHLH TFs to regulate anthocyanin biosynthesis in L. radiata remains to be further investigated.

In the anthocyanin biosynthesis pathway, DFR catalyzes dihydroflavonol conversion to leucoanthocyanidins (Zhang et al., 2014). DFR belongs to the superfamily of short chain dehydrogenase reductase (SDR), which has a highly preserved NADPH-binding domain “VTGAAGFIGSWLIMRLLERGY” as well as a substrate-binding domain in plants (Martens et al., 2002; Haselmair-Gosch et al., 2018). In this study, based on the expression level of the anthocyanin structure genes obtained from the transcriptome data, LrDFR1 and LrDFR2 (Table 1) were found to highly associate with the total anthocyanin content (PCC > 0.8), suggesting DFR may have an essential function in the phenotypic expression of L. radiata petal color. LrDFR1 was then identified as one of the hub genes (Figure 5B) and important to positively regulate anthocyanin production in L. radiata petals. Multiple amino acid alignments showed that LrDFR1 contains the NADPH-binding domains and substrate-binding domains. Phylogenic tree analysis revealed a high similarity between LrDFR1 and other characterized DFRs, implying that LrDFR1 belongs to the monocot DFR family and exhibits catalytic characteristics.

The DFR genes of Iris and Gentiana have been reported to be associated with the absence of brick-red flowers (Noda et al., 2017). Moreover, heterologous MaDFR expressions in N. tabacum has been associated with enhanced anthocyanin accumulation, which leads to darker flower colors, suggesting that MaDFR is involved in flower color development (Liu H. et al., 2019). After the introduction of maize (Zea mays) DFR into white-flowered petunia varieties, transgenic plant flowers accumulate non-native pelargonidin, which results in novel brick red-flower varieties (Meyer et al., 1987). In this study, the expression patterns of LrDFR1 was first temporally and spatially tested in various tissues and petal development stages of L. radiata. It showed that the expression levels of LrDFR1 were correlated with total anthocyanin accumulation. These findings imply that LrDFR1 is associated with petal color development in L. radiata (Figures 1, 7C,D). The spatial and temporal expression characteristics of LrDFR1 gene were found similarly in several other species (Liu H. et al., 2019; Lim et al., 2020). In order to investigate the functional divergence of LrDFR1 gene in the flavonoid biosynthesis, we performed transient expression analyses using Lycoris petals and tobacco leaves. Overexpressed LrDFR1 was associated with significantly elevated anthocyanin content and proanthocyanidin content in Lycoris petals and tobacco leaves. Interestingly, overexpression of LrDFR1 also enhanced the expression of downstream genes (LrANS and LrUFGT) involved in anthocyanins biosynthesis in transgenic Lycoris petals (Figures 8B–F). In addition, for plant breeders, a single DFR gene maybe ideal for determining flower colors. DFR is vital for pigmentation, when compared to other anthocyanin biosynthetic genes, which only regulate plant flower color hue. Thus, whether LrDFR1 has a high preference for dihydromyricetin, and is accountable for the limited flower colors in L. radiata needs further study.