Comparative Transcriptome Analysis Identifies Key Regulatory Genes Involved in Anthocyanin Metabolism During Flower Development in Lycoris radiata

Lycoris is used as a garden flower due to the colorful and its special flowers. Floral coloration of Lycoris is a vital trait that is mainly regulated via the anthocyanin biosynthetic pathway. In this study, we performed a comparative transcriptome analysis of Lycoris radiata petals at four different flower development stages. A total of 38,798 differentially expressed genes (DEGs) were identified by RNA sequencing, and the correlation between the expression level of the DEGs and the anthocyanin content was explored. The identified DEGs are significantly categorized into ‘flavonoid biosynthesis,’ ‘phenylpropanoid biosynthesis,’ ‘Tropane, piperidine and pyridine alkaloid biosynthesis,’ ‘terpenoid backbone biosynthesis’ and ‘plant hormone signal transduction’ by Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis. The candidate genes involved in anthocyanin accumulation in L. radiata petals during flower development stages were also identified, which included 56 structural genes (especially LrDFR1 and LrFLS) as well as 27 key transcription factor DEGs (such as C3H, GATA, MYB, and NAC). In addition, a key structural gene namely LrDFR1 of anthocyanin biosynthesis pathway was identified as a hub gene in anthocyanin metabolism network. During flower development stages, the expression level of LrDFR1 was positively correlated with the anthocyanin content. Subcellular localization revealed that LrDFR1 is majorly localized in the nucleus, cytoplasm and cell membrane. Overexpression of LrDFR1 increased the anthocyanin accumulation in tobacco leaves and Lycoris petals, suggesting that LrDFR1 acts as a positively regulator of anthocyanin biosynthesis. Our results provide new insights for elucidating the function of anthocyanins in L. radiata petal coloring during flower development.


INTRODUCTION
Plant pigments, such as anthocyanins, carotenoids and chlorophylls, play important roles in affecting the appearance of flower, fruit and seed coloring (Tanaka et al., 2008;Rebecca et al., 2010;Rosas-Saavedra and Stange, 2016;Cui et al., 2021). As an important group of plant pigments, anthocyanins are water soluble and belong to the family of flavonoids. So far, more than 500 different anthocyanins have been isolated from plants (Francavilla and Joye, 2020). They are highly involved in determining flower, seed, fruit and vegetative tissue colors, ranging from pink through scarlet, purple, and blue (Tanaka et al., 2008;Khoo et al., 2017). There are six species of anthocyanins (namely cyanidin, delphinidin, peonidin, malvidin, pelargonidin, and petunidin) in colorful plants (Tanaka et al., 2008;Castaneda-Ovañdo et al., 2009), of which cyanidin is responsible for red-purple coloration, delphinidin contributes to purple or blue-red, and pelargonidin contributes to red and orange (Khoo et al., 2017). Besides, anthocyanins also play various vital functions in plant biological functions, including disease protection, resisting environmental stresses, and promoting pollination (Lev-Yadun and Gould, 2009;Zhang et al., 2014).
The Lycoris species belongs to Amaryllidaceae family, and is a perennial bulb plant native to Northeast Asia, including China, South Korea, and Japan. It consists of about 20 species, of which 15 species and one variety are distributed in China . Among them, Lycoris radiata is considered ornamentally and medicinally valuable, as the colorful and special flowers have been used for decoration and the bulbs are notable to produce alkaloids with various biological activities (Park et al., , 2021. Anthocyanins are abundant in Lycoris flowers and also contribute to their color formation (He et al., 2011;Chun et al., 2013;Yue et al., 2019;Park et al., 2021). For example, four critical anthocyanins, namely cyanidin 3-sophoroside, cyanidin 3-xylosylglucoside, cyanidin 3-sambubioside, and pelargonidin 3-xylosylglucoside in L. longituba tepals of different colors have been well identified (He et al., 2011). In L. radiata flowers, three anthocyanins (cyanidin 3-diglucoside, cyanidin 3-sambubioside, and cyanidin 3-glucoside) were identified (Chun et al., 2013), and their presence during four flower development stages was confirmed more recently (Park et al., 2021). However, the molecular mechanisms of anthocyanins regulating color formation of Lycoris flower remain unclear. Thus, identifying the key genes related to color formation in Lycoris flower would provide a more sufficient genetic resource for manipulation of the related pathways to develop new cultivars with specific flower colors.
In recent years, transcriptome sequencing (RNA-seq) was used as a rapid technique to uncover DEGs, biosynthesis pathways, and TFs related to specific traits in plants (He et al., 2020;Li C. et al., 2020). In this study, we reported the changing profile of anthocyanins and gene expression dynamics in L. radiata petals at four developmental stages by integrated analyses of the physiology and transcriptome. We further identified modules with co-expressed genes and candidate hub genes for anthocyanin accumulation, and revealed LrDFR1 acts as a positive regulator involved in anthocyanin biosynthesis. Our results may serve as a reference for understanding the regulation of key genes and transcription processes in color formation in the flowers of this esthetically important Lycoris.

Plant Materials
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 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.

Measurement of Total Anthocyanins
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 530 -0.25 × A 657 ) × FW −1 , where Q Anthocyanins is the amount of anthocyanins, A 530 and A 657 is the absorption at the indicated wavelengths and FW represents the weight of the fresh sample [g].

Construction of the cDNA Library, Sequencing, and Transcriptome Assembly
Total RNA was extracted with the mirVana miRNA isolation kit (Thermo Fisher Scientific, Waltham, MA, United States) following the manufacturer's protocol. The quality and quantity of the RNA were examined by the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, United States). 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 TM 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).

Functional Annotation
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.

Identification of Differentially Expressed Genes
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.

Transcription Factors Analysis
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 . 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.

Gene Cloning and Construction of Expression Vectors
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 TM 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.
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.

Agrobacterium-Mediated Transient Transformation System of Lycoris Petals
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 600 of 0.6 in a buffer with 10 mM 2-(4-Morpholino) ethanesulfonic acid, 10 mM MgCl 2 , 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.

Validation RNA-Seq by Quantitative Real-Time PCR
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 . qRT-PCR assays were conducted by the SYBR Premix Ex Taq TM 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 Analysis
Statistical analyses were done by SPSS version 10.0 software (IBM Corporation, Armonk, NY, United States). 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).

Anthocyanin Levels in Lycoris radiata Petal During Flower Development Stages
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.

Transcriptome Sequencing and de novo Assembly
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).

Identification of Differentially Expressed Genes in Lycoris radiata Petal During Flower Development Stages
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 2 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 upregulated, 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).

Functional Annotation of Differentially Expressed Genes
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).

Identification of Key Differentially Expressed Genes Responsible for the Anthocyanin Biosynthesis Pathway
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.
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.

Validation of RNA-seq Data by qRT-PCR
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 2 > 0.76), indicating that the RNA-seq data were credible and accurate.

LrDFR1 Is Involved in Anthocyanin Biosynthesis in Lycoris radiata
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).  Figure. (B) Protein-protein interaction network constituted by protein sequences of differentially expressed transcription factors and structural genes involved in anthocyanin synthesis of L. radiata petals. Genes that have the higher weight are depicted in 'yellow and orange,' the 'blue edges' correspond to co-expressed strong links and the 'yellow edges' correspond to co-expressed weak links.  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 . Tobacco leaves were kept in a phytotron at 24 • C under constant lighting for 5 days. DMACA was used to stain proanthocyanidin. Every experiment was performed using 8-10 leaves for each genotype. Experiments were conducted in triplicates, and a representative image is shown. Proanthocyanidin levels of empty vector controls were set as the reference to 1. Asterisks represent significant differences between control and LrDFR1-overexpressing leaves (**p < 0.01). Bars = 1 cm. (C) Phenotypes of anthocyanin accumulation. Arrow indicates the transfected petals. (D,E) Relative anthocyanin levels in transiently transformed Lycoris petals (pBinGFP4: empty vector controls; LrDFR1-OE: LrDFR1-overexpressing petals). Lycoris petals were kept in a phytotron at 24 • C with a constant light for 5 days. Every experiment was performed using 8-10 petals per genotype. Data are shown as mean ± SD. **p < 0.01. Bar = 0.5 cm. (F) Relative expression levels of endogenous anthocyanin biosynthetic genes in pBinGFP4 (empty vector controls) as well as LrDFR1-overexpressing petals. Expression patterns of early biosynthetic genes (CHS, F3H, CHI, and F3 H) as well as late biosynthetic genes (DFR, UFGT, ANS, and 3RT) in petals were investigated. Asterisks represent significant differences between control and LrDFR1-overexpressing petals (**p < 0.01).
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.

Changes of Anthocyanin Contents in the Lycoris radiata Petals During Flower Development Stages
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).

Petals During Flower Development Stages
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 , "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 highquality 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.

Transcriptional Regulation of Color Formation in Lycoris 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 . 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 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., , 2020a. 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.
The LrDFR1 Drives Anthocyanin Accumulation in Lycoris radiata Petals In the anthocyanin biosynthesis pathway, DFR catalyzes dihydroflavonol conversion to leucoanthocyanidins . 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 . 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 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.

CONCLUSION
In this study, we provided a dynamic transcriptome profile of L. radiata petals during flower development stages. Overall, 56 structural genes and 27 key TF DEGs were identified as key genes responsible for L. radiata petal coloration. In the protein-protein interaction network analysis, LrDFR1 was identified as a hub gene in the anthocyanin biosynthesis pathway, and was highly associated with anthocyanin accumulation. Overexpression of LrDFR1 in Lycoris petals and tobacco leaves induced anthocyanin accumulation. In addition, the structural genes and co-expressed TFs reported in this study would serve as useful genetic resources for further functional characterization and molecular breeding programs in L. radiata. Taken together, our results elucidate on the molecular basis of petal development in L. radiata.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://ngdc.cncb.ac. cn/, CRA004779.

AUTHOR CONTRIBUTIONS
NW and ZW designed the research and wrote the manuscript. NW performed most of the experiments and data analysis. XS and FZ collected the experimental materials. TW assisted with the data analysis. WZ provided helpful comments on the manuscript. ZW provided guidance on the whole study and contributed with valuable discussions. All authors read and approved the final manuscript.