Molecular and Metabolic Insights into Anthocyanin Biosynthesis for Spot Formation on Lilium leichtlinii var. maximowiczii Flower Petals

Plants exhibit remarkable diversity in their petal colors through biosynthesis and the accumulation of various pigments. Lilium, an important cut and potted flower, has many coloring pattern variations, including bicolors and spots. To elucidate the mechanisms regulating spot formation in Lilium leichtlinii var. maximowiczii petals, we used multiple approaches to investigate the changes in petal carotenoids, spot anthocyanins, and gene expression dynamics. This included green petals without spots (D1-Pe and D1-Sp), yellow–green petals with purple spots (D2-Pe and D2-Sp), light-orange petals with dark-purple spots (D3-Pe and D3-Sp), and orange petals with dark-purple spots (D4-Pe and D4-Sp). D3-Pe and D4-Pe contained large amounts of capsanthin and capsorubin and small amounts of zeaxanthin and violaxanthin, which contributed to the orange color. In addition to cyanidin-3-O-glucoside, pelargonidin-3-O-rutinoside, cyanidin-3-O-rutinoside, and peonidin-3-O-rutinoside may also contribute to L. leichtlinii var. maximowiczii‘s petal spot colors. KEGs involved in flavonoid biosyntheses, such as CHS, DFR, and MYB12, were significantly upregulated in D2-Sp and D3-Sp, compared with D1-Sp, as well as in spots, compared with petals. Upregulated anthocyanin concentrations and biosynthesis-related genes promoted spot formation and color transition. Our results provide global insight into pigment accumulation and the regulatory mechanisms underlying spot formation during flower development in L. leichtlinii var. maximowiczii.


Introduction
Color is a visible trait of plants and is greatly significant to plant growth and development. Flower color formation relies on the type and content of pigment. In most plants, flower colors are mainly determined by flavonoids, especially anthocyanins, carotenoids, and betalains [1]. Flavones, flavanones, chalcones, and other compounds in flavonoids form yellow colors. Anthocyanins form red, purple, and blue colors. Carotenoids are yellow, orange, and red. Betalains form yellow, orange, and red colors, but they only exist in some Caryophyllaceae plants and cannot coexist with anthocyanins [2].
The genus Lilium (family Liliaceae) is a bulbous monocot comprising more than 100 species [13]. As important cut and potted flowers, they exhibit many horticultural features, such as flower color, shape, and fragrance. The flower colors are principally red, pink, orange, yellow, white, and some intermediate colors. There are also some special colors, such as green lily (L. fargesii) and purple lily (L. Souliei). After years of breeding, many new commercial lilies have been cultivated, showing extensive variation in flower color [14]. The Royal Horticulture Society in the UK divided lily varieties into nine hybrids according to their characteristics. Among them, the typical characteristics of Asiatic hybrids include a wide variation in flower colors, from yellow to orange, pink, red, or white. The main flower colors of Oriental hybrids are pink, reddish violet, and white, while L. longiflorum hybrids are mainly white [15].
Lilies also have many variations in flower coloring patterns, such as bicolors and spots [15]. In bicolor lilies, different anthocyanin types and contents accumulate in different tepals parts; in L. regale and some trumpet lilies, anthocyanins accumulate in the throat. Some Asiatic hybrids have different colors in the upper and lower tepals. In Oriental hybrids and some of their parents, the anthocyanin coloring in the vascular bundle region of tepals is star-shaped, which is inconsistent with the other tepal parts.
In many wild lilies and lily hybrids, the inner surfaces of the tepals are blotched with red or dark spots. These spots are generally divided into three types. The first type is raised spots, located on the adaxial surface of tepals. The increase in epidermal and parenchymal cells can cause the surface of the tepals to uplift. Anthocyanins accumulate in both epidermal and parenchymal cells, and cell division may increase before pigment synthesis [15,16]. The second type is splatters: many small spots are splattered in the lower part of the tepals. The spot area is not uplifted and the surface is smooth. Anthocyanins only accumulate in epidermal cells, and there is no difference in morphology between pigmented cells and those without pigments [16,17]. The third type is brush marks: there include brush-like nut-brown stripes along the vascular bundles at the base of the tepals [15].
L. leichtlinii var. maximowiczii (section Sinomartagon), a wild species known for its tolerance of cold native to Japan, North Korea, the Far East of Russia, and north and northeastern China, is an important species used for breeding Asiatic hybrids [18]. Meanwhile, the bulbs of L. leichtlinii var. maximowiczii have been used as food and medicine in China, Japan, and Korea as typical edible lily species [19]. The flower is recurred and pendulous, and its petals are orange-red, with raised black spots. However, anthocyanins in the spots on the petals have not been sufficiently identified or quantified. In this study, transcriptomic and metabolomic technologies were used to identify differentially expressed genes (DEGs) and to reveal anthocyanin changes in petals that determine spot formation in L. leichtlinii var. maximowiczii. Through RNA-seq abundance and high-performance liquid chromatography (HPLC), we investigated the transcript profiles of the petal spots during multiple floral development stages to elucidate the relationship between structural genes and TFs for anthocyanin biosynthesis in L. leichtlinii var. maximowiczii petal spot. Our results provide new insights for the identification of functional genes and metabolites in spot formation, which lays a foundation for spot improvement using advanced breeding technology in lily varieties.

Morphology Analysis of Petal Color Transitions and Spots Formation
The petal color of every flower transformed continuously from green to yellowish green, and then from light orange to orange-red during floral development in L. leichtlinii var. maximowiczii. According to the petal color transition and spot appearance processes, the floral developmental stages were described as follows: no obvious anthocyanin or carotenoid pigmentation was observed at stage 1 (D1), anthocyanin coloring at spots started to appear at stage 2 (D2), carotenoid coloring of the petal background was obvious at stage 3 (D3), full pigmentation of the petal background occurred at stage 4 (1 d before anthesis, D4), and flowers blossomed at stage 5 (D5) ( Figure 1A). Morphological analysis showed no orange coloration in the petals before floral development stage D3. The raised spot formation at D2 ( Figure 1C) was ahead of the petal color transitions at D3. Based on morphological anatomy analysis, we collected petals (Pe) and spots (Sp) at D1, D2, D3, and D4 for transcriptome sequencing and metabolome identification ( Figure 1B). (green petals without spots), D2 (yellowish green petals with purple spots), D3 (light-orange petals with dark-purple spots), D4 (orange petals with dark-purple spots, one day before anthesis), and D5 (the day of anthesis). Bar = 1 cm. (B) Sampling strategy: petals and spots were collected at D1, D2, D3, and D4 for transcriptome sequencing and metabolome identification. (C) Scanning electron microscopy appearance of raised spots on D2 petals. Bar = 500 µm.

Carotenoid Identification and Quantification in Petals
To obtain an accurate understanding of carotenoid accumulation in petals, carotenoid profiling was performed on L. leichtlinii var. maximowiczii petals using LC-ESI-MS/MS during petal color transitions. Fourteen carotenoids were detected in D1-Pe, D2-Pe, D3-Pe, and D4-Pe (Table 1). The major carotenoids in D3-Pe and D4-Pe were capsanthin, capsorubin, zeaxanthin, and violaxanthin. The levels of these four carotenoids showed an upregulation trend from D1 to D4. The capsanthin content in D1-Pe and D2-Pe was 1.22 µg/g and 2.27 µg/g, respectively. Then, the content drastically increased hundredsfold to 406.27 µg/g in D3-Pe and 940.63 µg/g in D4-Pe. The capsorubin content significantly increased from 0.25 µg/g in D1-Pe and 1.65 µg/g in D2-Pe to 215.28 µg/g in D3-Pe and 493.43 µg/g in D4-Pe. The zeaxanthin and violaxanthin contents were 31.24 µg/g and 29.07 µg/g in D4-Pe, respectively, which were both much lower than the contents of capsanthin (940.63 µg/g) and capsorubin (493.43 µg/g) in D4-Pe.
Lutein had the highest carotenoid content in D1-Pe, which showed a decreasing trend. The lutein content slightly decreased from 158.63 µg/g in D1-Pe to 114.62 µg/g in D2-Pe and then significantly decreased to 32.10 µg/g in D3-Pe and 1.27 µg/g in D4-Pe.
The 17 common anthocyanins included the top 10 most abundant anthocyanins, including cyanidin, pelargonidin, and peonidin, especially cyanidin-3-O-glucoside, pelargonidin-3-O-rutinoside, cyanidin-3-O-rutinoside, and peonidin-3-O-rutinoside, which play a significant role in spot formation in L. leichtlinii var. maximowiczii (Table 2). Anthocyanins identified in the three comparison groups showed similar accumulation, which corresponded to the spot color intensity of L. leichtlinii var. maximowiczii. We propose that the DAMs found in anthocyanin biosynthesis, especially the top ten highest anthocyanidin contents, may contribute to spot coloring in L. leichtlinii var. maximowiczii petals.   indicates the metabolite contents from low to high. Identification of DAMs between three comparison groups was performed using variable importance in projection values ≥1 and fold change ≥2 or ≤0.5; in addition, the content must be ≥ 1 µg/g DW in spots in at least one stage.

De Novo Assembly and Gene Function Annotation
RNA-seq analysis was performed to further study the molecular mechanism of spot formation. Twenty-four libraries were established using Pe and Sp at D1, D2, D3, and D4 stages (three biological replicates for each sample), and a total of 8.41 to 11.05 G clean base for each sample were obtained through sequencing. A total of 1,598,209,652 raw reads and 1,541,848,276 clean reads were obtained from twenty-four samples. The Q20 and Q30 values of each library were ≥97.93 and 94.17%, respectively. The GC content of each sample ranged from 49.37-50.79%, with an average of 50.22% (Table S3).
A total of 277,166 unigenes were obtained, with an average length of 743 bp and an N50 of 1088 bp, using the Trinity (v2.11.0) method (Table S4). The total number of unigenes was much higher than the other eukaryotic; we thought that there were quite a number of lncRNAs in the assembled unigenes. Unigenes without coding potential, as predicted by CNCI (version 2), CPC (version 0.9-r2), CPAT (1.2.4), and the intersection of the nonprotein-coding potential results were chosen as novel long non-coding RNAs (lncRNAs). As a result, 151,202 (54.55%) unigenes were identified as lncRNAs, the remaining 125,964 (45.45%) unigenes were identified as coding unigenes, and our subsequently RNA-seq analysis mainly focused on the protein-coding unigenes.

Identification of DEGs and Co-Expression Network Analysis
For each unigene, the fragment per kilobase of transcript per million mapped reads value was calculated to quantify the expression abundance and variation using RSEM software. Unigenes with a corrected p value < 0.05 and |log2 (fold change) | > 1 were considered DEGs. A total of 36,405 differentially expressed unigenes were identified from 16 comparisons, and the number of differentially expressed unigenes in each comparison is shown in Table S6.
To obtain a comprehensive understanding of genes expressed in the petals and spots across the period of flowering and to identify the specific genes highly associated with spot formation, weighted gene co-expression network analysis was performed using the differentially expressed unigenes from 16 previously identified comparisons. Co-expression networks were constructed based on pairwise correlations of gene expression across all samples. Modules were defined as clusters of highly interconnected genes, and genes within the same cluster had high correlation coefficients. This analysis identified 22 distinct modules (labeled with different colors), as shown in the dendrogram in Figure 5A. Notably, three out of 22 co-expression modules, light yellow, dark green, and dark magenta, were composed of genes highly correlated with spot formation and anthocyanin accumulation at D2 and D3, with an upregulation trend. GO and KEGG enrichment analyses were subsequently performed for each module. The light yellow module consisted of 1618 genes. Their expression significantly increased in D2 and D3 and then decreased in D4. KEGG enrichment analysis showed significant enrichment of the flavonoid biosynthesis (ko00941) pathway, including 28 flavonoid-related unigenes, such as CHS, CHI, F3H, F3 H, and MT ( Figure 5C; Table S7). Among the 76 identified TFs in the light yellow module, 24 candidates were involved in anthocyanin accumulation, including MYB, bHLH, WD40, WRKY, bZIP, and ERF ( Figure 5B; Table S8).
A total of 318 genes were identified in the dark magenta module, displaying a modest increase from D1-Sp to D2-Sp, compared with the light yellow module, and the flavonoid biosynthesis (ko00941) pathway was the most enriched ( Figure 5D

Discussion
The main aim of this study was to understand the transcriptional control of spot formation and anthocyanin distribution in L. leichtlinii var. maximowiczii flowers. In Asiatic hybrids, carotenoids are the main color components of yellow and orange flowers, anthocyanin is predominant in pink and brown flowers, and red flowers contain both anthocyanins and carotenoids [20,21]. A previous study reported that lutein, violaxanthin, and β-carotene are major carotenoids in green tissues, which are conserved among plants [22]. In contrast, carotenoid compositions in flowers and fruits are different from those in green tissues, and large variations have been observed among species [15,23].
Previous studies have shown that the color components of orange petals, such as L. amabile, L. davidii var. willmottiae, and L. leichtlinii var. maximowiczii include capsanthin, capsorubin, zeaxanthin, violaxanthin, and antheraxanthin [24]. The orange color of tiger lily (L. lancifolium 'Splendens') flowers is primarily due to the accumulation of two k-xanthophylls: capsanthin and capsorubin [25]. Conversely, no anthocyanins have been detected in white Asiatic and Oriental hybrids and L. longiflorum [26][27][28]. Generally, the difference in the flower color of lily hybrids is caused by the difference in anthocyanin content; a large amount of cyanidin 3-O-β-rutinoside and a small amount of cyanidin 3-O-β-rutinoside-7-O-β-glucoside accumulate in some red varieties of Oriental and Asiatic hybrids [28]. The main component of spots on tepals is also cyanidin 3-O-β-rutinoside [29].
In our study, in the green petal at D1-Pe, the carotenoid with the highest content was lutein (158.63 µg/g), with small amounts of neoxanthin (13.47 µg/g), β-carotene (12.85 µg/g), and violaxanthin (11.66 µg/g) were also present. When the petal background turned light orange at D3-Pe and orange at D4-Pe, all carotenoids, except violaxanthin, showed a downward trend (Figure 1, Table 1). Capsanthin and capsorubin drastically increased at D3-Pe and D4-Pe during petal color transitions. A large amount of capsanthin and capsorubin, with small amounts of zeaxanthin and violaxanthin, may contribute to the orange color of L. leichtlinii var. maximowiczii petals.
Anthocyanins are synthesized through a series of enzymatic reactions along the phenylalanine pathway and are further modified by glycosylation, methylation, and acylation until they are finally stored in vacuoles [32]. Anthocyanin biosynthesis gene expression levels determine the anthocyanin content in plants [33]. In a recent study, eight structural genes (CHSa, CHSb, CHIa, CHIb, F3H, F3 H, DFR, and ANS) and three TF genes (MYB12, MYB15, and bHLH2) were isolated, and anthocyanin accumulation was examined during flower development in four cultivars of two Lilium species. High transcript accumulations were observed in the red tepals of 'Gran Tourismo,' followed by the pink tepals of 'Perth,' and both tepals contained cyanidin. The white tepals of 'Rialto' and 'Lincoln' did not contain anthocyanin and showed the lowest transcript accumulations [34]. The expression of these structural genes is strongly correlated with MYB12 and MYB15 expression [34]. Three CHS genes and one DFR gene were cloned in the Asiatic hybrid 'Montreux.' LhCHSA and LhCHSB were expressed in the tepals, filaments, and pistils. LhCHSC is specifically expressed in anthers [35,36]. DFR is only expressed in the organs where anthocyanin accumulates, whereas CHS is also expressed in some sites where anthocyanin does not accumulate, indicating that DFR is more closely related to anthocyanin coloring than CHS in Asiatic lilies. Coloring analysis of the filial generation showed that CHS and DFR independently control tepal and spot coloring [35].
In our study, seven key structural genes were activated during anthocyanidin biosynthesis in L. leichtlinii var. maximowiczii. Among them, the expression of DEGs encoding CHS, CHI, F3H, F3 H, DFR, and ANS was higher in D2-Sp and D3-Sp than in D2-Pe and D3-Pe, which may result in higher anthocyanidin levels during spot formation ( Figure 6). Twenty CHSs, including CHSA, CHSB, and CHSC, were significantly upregulated in spots (D2-Sp, D3-Sp, and D4-Sp), but could not be detected in petals (D1, D2, D3, D4-Pe), as the RNA-seq results showed ( Figure 6B; Table S11). CHS is a key enzyme in anthocyanin biosynthesis, and the loss of its activity typically results in albino flowers [37,38]. Similarly, in the tree peony cultivar 'Qing Hai Hu Yin Bo', PsCHS was strongly and specifically expressed in petal blotches that appeared at the base of the white petals [39], which was most likely the direct cause of anthocyanin blotch formation. Three CHI, three DFR, three ANS, one F3H, and two F3 H genes were identified in RNA-seq, showing the same expression profiles as CHS.
Although our study and other studies have clarified some key genes associated with flower coloration in Lilium, the relationship between gene expression and flower color development needs to be more thoroughly investigated. In many plants, additional genes have been identified as important factors that affect anthocyanin synthesis in flowers and fruits. In Petunia hybrida petals, CHS silencing inhibited the accumulation of anthocyanins and the mRNA levels of the corresponding endogenous targets, such as CHI and DFR were unaffected [40]. Sequence-specific degradation of CHS mRNA in petal sectors along the central veins is the known cause of the P. hybrida 'Red Star' pigmentation pattern [41]. In Lycium, CHS, CHI, F3 5 H, DFR, ANS, and UFGT expression levels were highly consistent with anthocyanin accumulation in black and red fruit [42,43]. In purple-red Padus virginiana L. leaves, the flavonoid biosynthesis genes (PAL, CHS, and CHI) and their transcriptional regulators (MYB, HD-Zip, and bHLH) exhibited increased expression during purple-red periods [44]. These findings indicate the importance of CHS, CHI, ANS, F3H, F3 H, DFR, and UFGT in anthocyanin biosynthesis.
In addition to key structural genes, TFs, such as MYB, bHLH, WD, and WRKY, also play important roles in anthocyanin biosynthesis [7,31]. In many Lilium species, MYB12 usually interacts with bHLH2 to form a MYB12/bHLH2 complex that upregulates structural gene transcription [28,45]. By inhibiting the transcription of MYB12, for example, in high heat (35 • C), the transcription accumulation of synthetic genes decreases, and the anthocyanin accumulation is thus reduced [9]. Two TFs are involved in regulating anthocyanin synthesis in the Asiatic hybrid 'Montreux' [46]. In the white tepals of 'Rialto' and L. speciosum, the W-to-L substitution in the R2 repeat of MYB12 was responsible for the absent transcriptional activation of anthocyanin structural genes, resulting in a lack of anthocyanin accumulation [27,47]. In Clarkia gracilis, a species with petal spots/blotches similar to tree peony, allelic variation in the cis-regulatory region of the R2R3-MYB gene, CgMYB1, restricts its expression and subsequent anthocyanin accumulation to either a basal or central petal spot/blotch in different subspecies [48]. Further, an R2R3-MYB transcription factor (CgsMYB12) was identified to be responsible for anthocyanin pigmentation of the basal region ('cup') in the petal of C. gracilis ssp. sonomensis. Additionally, two R2R3-MYB genes, CgsMYB6 and CgsMYB11, were involved in petal background pigmentation [49]. In this study, two MYB (cluster-12430.165019 and cluster-12430.132482 MYB12) and one bHLH (cluster-12430.66365) genes were specifically expressed in spots between D2-Sp and D2-Pe (Table S10), and they were also enriched in the light yellow module (Table S8). They may positively regulate the expression of most structural genes involved in flavonoid biosynthesis with the MBW complex during spot formation. In the Oriental hybrid lily 'Sorbonne', MYB12 regulates both whole tepal and raised spot pigmentation [50]. However, MYB12 (cluster-12430.132482) was specifically expressed in spots between D2-Sp and D2-Pe and D3-Sp and D3-Pe. At the D3 stage, the petals appeared light orange in color, indicating that more MYB was involved in petal coloration in L. leichtlinii var. maximowiczii flowers. Some TFs, such as AP2/ERF, bZip, and WRKY, display differential expression during spot formation, suggesting that they may synergistically regulate the expression of genes involved in flavonoid biosynthesis, as mentioned in a recent report, where the MYB and HD-Zip TFs exhibited similar expression, possibly acting as key regulators of flavonoid biosynthesis.

Plant Materials and Sampling
Wild lily L. leichtlinii var. maximowiczii species were used in this study. Bulbs were planted in pots on 3 March, and grown in a greenhouse (unheated) at the Shanghai Academy of Agricultural Sciences, Shanghai, China under natural light with daytime and nighttime temperatures of 17-26 • C. On 22 April, the flower buds were obvious. The floral developmental stages are described as follows, D1 (green petals without spots, bud length <3 cm), D2 (yellowish green petals with purple spots, bud length <3-4.5 cm), D3 (light-orange petals with dark-purple spots), D4 (orange petals with dark-purple spots, one day before anthesis), and D5 (the day of anthesis) ( Figure 1A). Petals and spots were collected from flowers at D1, D2, D3, and D4 ( Figure 1B). All materials were sampled, immediately frozen in liquid nitrogen, and stored at −80 • C.

Sample Preparation and Extraction
For Flavonoids: the spots samples were freeze-dried, ground into powder (30 Hz,1.5 min), and stored at −80 • C until needed. Fifty milligrams of powder were weighed and extracted with 0.5 mL methanol/water/hydrochloric acid (500:500:1, v/v/v). Then, the extract was vortexed for 5 min, ultrasonicated for 5 min, and centrifuged at 12,000× g at 4 • C for 3 min. The residue was re-extracted by repeating the aforementioned steps under the same conditions. The supernatants were collected and filtered through a membrane filter (0.22 µm, Anpel) before LC-MS/MS analysis.
For carotenoids: the petals samples were freeze-dried, ground into powder (30 Hz, 1.5 min), and stored at −80 • C until needed. Fifty milligrams of powder was weighed and extracted with a 0.5 mL mixed solution of n-hexane: acetone: ethanol (1:1:1, v/v/v), 10 µL of internal standard mixed solution (20 µg/mL) was added into the extract as internal standards (IS). The extract was vortexed for 20 min at room temperature. The supernatants were collected after centrifuging at 12,000 r/min for 5 min at 4 • C. The residue was reextracted by repeating the above steps again. To the supernatants, 0.5 mL saturated sodium chloride solution was added and vortexed, and the upper layer was collected, this step was repeated two times more. Then, the supernatant was evaporated to dryness and dissolved in 0.5 mL MTBE, then 0.5 mL 10% KOH-MeOH was added, the mixture was vortexed, and the reaction was allowed to take place at room temperature overnight. After the reaction, 1 mL saturated sodium chloride solution and 0.5 mL MTBE were added and vortexed, and the upper layer was collected, this step was repeated two times and the supernatant was evaporated to dryness and reconstituted in 100 µL mixed solution of MeOH/MTBE (1:1, v/v). The solution was filtered through a 0.22 µm membrane filter (0.22 µm, Anpel) for further LC-MS/MS analysis.

Metabolite Profiling
Flavonoids in the spots (D1-Sp, D2-Sp, D3-Sp, and D4-Sp) and carotenoids in the petals (D1-Pe, D2-Pe, D3-Pe, and D4-Pe) of flowers at the four stages were identified and quantified using an LC-ESI-MS/MS system (HPLC, UFLC SHIMADZU CBM20A system, www.shimadzu.com.cn/, accessed on 23 August 2021; MS, Applied Biosystems 6500 QTrap, www.appliedbiosystems.com.cn/, accessed on 23 August 2021). Three biological samples were evaluated for each part of the spots or petals (i.e., samples of the same segment from at least five individuals were mixed for one replicate) for a total of twenty-four samples.
Metabolites were identified and quantified by Wuhan MetWare Biotechnology Co., Ltd. (Wuhan, China) based on the self-built MWDB database and a public database. Carotenoids and flavonoids were analyzed using scheduled multiple reaction monitoring (MRM). Data acquisitions were performed using Analyst 1.6.3 software (Sciex). All metabolites were quantified using Multiquant 3.0.3 software (Sciex). The identified metabolites with significant differences in content were set with 0.5 ≥ fold change or ≥ 2, p value < 0.05, and VIP ≥ 1. To study the accumulation of specific metabolites, principal component analysis and orthogonal partial least squares-discriminant analysis were performed using R (www.r-project.org/, accessed on 23 August 2021).

RNA Extraction, Quantification, and Sequencing
Petals and spots from flowers at D1, D2, D3, and D4 were collected, and three independent biological replicates were used. The amount of 0.1 g of flower tissue was used for RNA extraction. Purity, quantification, and integrity of total RNA were estimated using the NanoPhotometer ® spectrophotometer (IMPLEN, Westlake Village, CA, USA), Qubit ® RNA Assay Kit in Qubit ® 2.0 Fluorometer (Life Technologies, Carlsbad, CA, USA), and the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA). mRNA was purified from a total of 1 µg of RNA per sample using poly-T oligo-attached magnetic beads. A NEB fragmentation buffer was added to break the RNA into short segments, and the first-strand cDNA was synthesized by reverse transcriptase. The second-strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. In order to select cDNA fragments of preferentially 250~300 bp in length, the library fragments were purified with the AMPure XP system (Beckman Coulter, Beverly, USA). Finally, 24 libraries representing the collected petal and spot samples were prepared following standard procedures of the Illumina HiSeq TM 4000 system (San Diego, CA, USA) according to the manufacturer's instructions. Next, cDNA libraries were sequenced on an Illumina sequencing platform and 150 bp paired-end reads were generated by Metware Biotechnology Co., Ltd. (Wuhan, China).
The resulting unigenes were evaluated by identifying the complete BUSCO hits [51]. Following gene annotation, the potential transcript sequences encoded by each gene were subjected to DIAMOND or HMMER analysis against the NCBI non-redundant protein, eggNOG, Swissprot, Trembl, and Pfam databases to identify homologous proteins in other species. Possible GO terms in the protein sequences were identified using SwissProt and Trembl database ID mapping. The genes were annotated with KEGG terms using KOBAS (version 3.0), and TF analysis was performed using iTAK. TransDecoder software (https://github.com/TransDecoder/TransDecoder/wiki, accessed on 30 August 2021) was used to predict the coding regions within the transcript sequences. Unigenes without functional annotation were used to perform lncRNA prediction using CPC (v0.9-r2), CPAT (v1.2.4), and CNCI (version Feb 2 28, 2014). Unigenes predicted as lncRNAs using all three software programs were identified as lncRNAs, and a total of 151,202 lncRNAs were identified.
For each sample, paired-end reads were mapped back to the unigenes using bowtie (v2.3.4.1). The rsem-calculate-expression function within the software package RSEM (v1.3.1) was used to estimate the expression levels for each gene using the expectationmaximization algorithm. The expression profiles generated for each sample were combined and used to detect DEGs using DESeq2 (v1.22.2), and unigenes with fold change values > 2 and an adjusted p value < 0.05, according to the False Discovery rate method of Benjamini and Hochberg, were selected as DEG candidates. Subsequently, the WGCNA package (v1.71) was used for co-expression network analysis. All software parameters used in this study are listed in Table S14.

Petal Surface Analysis Using Scanning Electron Microscopy
The morphological characteristics of L. leichtlinii var. maximowiczii petals were examined using scanning electron microscopy. Samples were handled according to Salehi et al. [52] with minor modifications. The fresh petals and spots from the four stages were dissected at approximately 5 mm, fixed in 2.5% glutaraldehyde solution at room temperature for 2 h, and then washed five times for 10 min each in phosphate buffer (0.1 M), pH 7.4. Afterward, the samples were dehydrated by washing in a series of ethanol solutions [2 × 50% (30 min), 75% (30 min), 90% (30 min), 95% (30 min), and 2 × 100% ethanol (30 min)]. Then, the samples were put in isoamyl acetate for 20 min and critical point dried in an HCP-2 Hitachi (Tokyo, Japan) equipment, wherein they were placed first in a 50-80% liquid carbon dioxide (L-CO 2 ) at 10 • C for 20 min, and then at 40 • C for 5 min. For coating with 10 nm of gold, the dried samples were placed in metal stubs and placed in an ion sputter (MC1000, Hitachi, Tokyo, Japan). A scanning electron microscope (SU8100, Hitachi, Tokyo, Japan) was used for observations.

RT-qPCR Analysis
To validate the transcriptome data, the relative expression of nine DEGs identified in the transcriptome analysis was evaluated via RT-qPCR using three biological and three technical replicates. Total RNA was extracted using TRIzol reagent according to the manufacturer's instructions (Aidlab Biotechnology Co., Ltd., Beijing, China). The firststrand cDNA was synthesized using 2 µg of total RNA as a template with a first-strand cDNA Synthesis Kit (Toyobo, Japan) in a 20 µL reaction volume. Gene-specific primers for RT-qPCR were designed according to the selected sequences derived from RNA-seq (Table S15). The 18 S rRNA gene was used as an internal control. RT-qPCR was performed using the TaKaRa SYBR Green Mix kit (TaKaRa, Beijing, China) and the ABI 7500 fast real-time detection system. The RT-qPCR program was initiated with a denaturation step (95 • C for 30 s), followed by 40 cycles of PCR (95 • C for 15 s, 58 • C for 30 s), melting (95 • C for 15 s, 60 • C for 1 min, 95 • C, 15 s). Relative expression analyses of quantitative data were performed using the 2 −∆∆Ct method [53].

Conclusions
In this study, metabolomic and transcriptome analyses were used to identify the key anthocyanins and genes responsible for L. leichtlinii var. maximowiczii spot coloration in the petals. Overall, 40 anthocyanin metabolites were detected in L. leichtlinii var. maximowiczii petal spots. Cyanidin-3-O-glucoside, pelargonidin-3-O-rutinoside, and cyanidin-3-O-rutinoside are the main anthocyanin components in L. leichtlinii var. maximowiczii spots. Moreover, nine enzymes (especially CHS and DFR) in the anthocyanin biosynthesis pathway and seven differentially expressed TFs (especially MYB12 and bHLH) were identified as candidate regulators contributing to the color diversity. The results of this study provide valuable information and new insights for the further evaluation of L. leichtlinii var. maximowiczii genetic diversity.