Transcriptomic Profiling of Heat-Treated Oriental Lily Reveals LhERF109 as a Positive Regulator in Anthocyanin Accumulation

1. Introduction

The suppression of anthocyanin biosynthesis under high temperature is a widespread phenomenon observed across the plant kingdom, affecting both fruits and ornamental flowers [1]. This temperature-sensitive response has gained increasing research attention due to global climate warming and its significant impact on horticultural industries [2].

Mounting evidence suggests that high temperature can directly impacts MYB transcription factors and influence anthocyanin accumulation through complex upstream regulatory pathways. In Oriental lilies, elevated temperatures disrupt MYB transcription factor balance by suppressing R2R3-MYB activators (LhMYB12) while inducing MYB repressors (LhMYBC2), leading to reduced expression of anthocyanin biosynthetic genes including CHS, DFR, and ANS [3,4]. Similar direct effects on MYB regulation have been documented in other species: in grape berries, high temperature significantly reduces the expression of VvMYBA1 and VvMYBA2 activators [5]; in apples, heat stress downregulates MdMYB10 [6]; and in Arabidopsis, elevated temperatures suppress AtPAP1/AtMYB75 activator expression while enhancing the repressive activity of AtMYBL2 [7]. In Malus profusion, high temperatures suppress R2R3-MYB activators (MpMYB10) while inducing MYB repressors (MpMYB15), reducing expression of key anthocyanin biosynthetic genes MpCHS, MpDFR, MpLDOX, and MpUFGT [8].These examples collectively demonstrate that temperature-induced modulation of the balance between MYB activators and repressors represents a conserved regulatory mechanism affecting anthocyanin biosynthesis across diverse plant species.

In fruits and Arabidopsis thaliana, the mechanisms by which high temperature inhibits anthocyanin accumulation have been extensively characterized. In apple (Malus domestica), thermal stress modulates the expression and activity of MYB10 through the COP1-HY5 photothermal integration module [9]. The interaction between high temperature and light signaling involves B-box (BBX) transcription factors, where MdCOL4 (BBX24) interacts with MdHY5 to inhibit MdMYB1 expression [10]. Similar upstream regulatory networks have been identified in Arabidopsis, and the same COP1-HY5-MYBL2 signaling module has been confirmed to regulate anthocyanin biosynthesis in response to high temperature [11,12]. In contrast to fruits and model plants, ornamental flowering crops show diverse regulatory mechanisms for anthocyanin biosynthesis under high temperatures. Some heat-tolerant cultivars maintain stable pigmentation at elevated temperatures, suggesting the HY5-COP1 module may not be the primary regulatory pathway in these species [6,13]. Instead, ornamental flowering crops exhibit more complex regulatory mechanisms involving multiple external cues (like light and hormones) that may dilute the influence of the HY5-COP1 module in controlling floral pigmentation [14,15,16]. Although research on upstream mechanisms of high-temperature suppression of anthocyanin synthesis is limited, studies show that MYB transcription factors, key regulators of anthocyanin biosynthesis, are regulated by various upstream factors [17,18,19,20]. Among these, APETALA2/ETHYLENE RESPONSIVE FACTOR (AP2/ERF) family members may play important roles in integrating temperature and hormone signals to control floral pigmentation. For instance, in ‘Viviana’ lilies (Lilium spp.), LvERF113, induced by ethylene, promotes anthocyanin synthesis in tepals by repressing the negative regulator LvMYB1, forming the LvMYB5-LvERF113-LvMYB1 module under LvMYB5’s positive regulation[21]. Similarly, in apple (Malus domestica), MdERF109 enhances fruit anthocyanin biosynthesis via methyl jasmonate (MeJA) signals, interacting with MdWER[22]. These examples underscore ERF’s role in coordinating hormonal and environmental cues for pigmentation.

Oriental lily (Lilium spp., Oriental Group) is renowned for its exceptionally large, fragrant flowers with vibrant colors and distinctive recurved petals [23]. These outstanding horticultural traits establish its significant position in both cut flower industry and landscape applications [24]. As an important breeding parent, Oriental lily has been widely used in interspecific hybridization, resulting in series like LLO (Longiflorum × Oriental lily) [25] and AOA (Asiatic × Oriental lily) [26,27]. However, pink-flowered cultivars commonly exhibit heat sensitivity, with flower color fading under high temperature stress [28]. This undesirable trait is often inherited by progeny, seriously affecting product quality and economic value. Research has shown that flower color fading in Oriental lilies under high temperature is primarily regulated by MYB transcription factors, however, the upstream regulatory mechanisms remain unclear. Previous studies have identified the crucial roles of MYB transcription factors such as LhMYB12 and LhMYBC2 in response to high temperature, but little is known about the upstream molecular networks that regulate these MYB factors.

In this study, we performed transcriptome sequencing analysis of ‘Souvenir’ Oriental lily under control and high-temperature conditions, and conducted weighted gene co-expression network analysis (WGCNA) with known MYB transcription factors, anthocyanin biosynthesis structural genes, and related transcription factors to explore the upstream regulatory mechanisms of high temperature-induced color fading.

2. Materials and Methods

2.1. Plant Materials and Temperature Treatments

Oriental lily ‘Souvenir’ plants at developmental St₂-St₃ were selected for uniformity and subjected to controlled temperature treatments of 20℃ (control) and 35℃ (high temperature) in growth chambers. Environmental conditions were standardized across both treatments (16h light/8h dark photoperiod, 70% relative humidity, 12,000 lx light intensity). Tepal samples were systematically collected at 12h, 24h, 48h, and 72h post-treatment (designated as T1, T2, T3, and T4, respectively), with three biological replicates per treatment combination to enable robust color measurement and physiological analyses.

2.2. Determination of Anthocyanin Content

Total anthocyanin content was measured using the pH differential method, with modifications based on Song et al. [29]. 0.5 g of finely ground tepal tissue was extracted with 5 mL of acidified methanol (0.1% HCl, v/v) for 24 h at 4°C in darkness. The extract was centrifuged at 12,000 g for 10 min at 4°C, and the supernatant was collected. Absorbance was measured at 530 nm and 657 nm using a spectrophotometer (UV-1800, Shi madzu, Japan). Total anthocyanin content was calculated as X = (A₅₃₀ - 0.25A₆₅₇) × . where X is the total anthocyanin content (relative units), A530 and A645 are the absorbance values at 530 nm and 645 nm, respectively, and is the sample weight (g). Measurements were performed in triplicate for each treatment.

For specific anthocyanin analysis, the method proposed by Xia et al. [30] was adopted, and the determination was performed using a high-performance liquid chromatography (HPLC) instrument (Shimadzu Corporation, Japan). Anthocyanins were separated on a Shim-pack Scepter C18 column (4.6 mm × 250 mm, 5 μm) at 30°C. The mobile phase consisted of (A) 0.1% formic acid and (B) acetonitrile. The gradient elution program was as follows: 0-8 min, 10-22.2% B; 18-25 min, 22.2-22.8% B; 25-28 min, 22.8-44% B; 28-45 min, 44-50% B; 45-50 min, 50-60% B; 50-70 min, 60-85% B. The flow rate was 1.0 mL/min, and the injection volume was 10μL. Cyanidin-3-O-glucoside was identified by comparing retention times and spectral characteristics with authentic standards (Shanghai yuanye), and quantified using standard curves. Each sample run for HPLC was repeated three times under the same conditions.

2.3. Determination of Sugar Content

Samples preserved in a -80 °C freezer were retrieved for the determination of sugar content. The contents of sucrose, reducing sugars, and soluble sugars were determined according to the kit protocol (G0506W, G0502W, G0501W Grace Biotechnology, Suzhou, China).

2.4. RNA Extraction and Transcriptome Sequencing

Total RNA was extracted from ‘Souvenir’ lily tepal samples using the Universal Plant Total RNA Extraction Kit (RC112-01, Vazyme, Nanjing, China) according to the manufacturer's instructions. RNA quality and integrity were assessed via agarose gel electrophoresis, with purity verified using a Nanodrop 2000 spectrophotometer (Thermo Fisher Scientific Inc., USA) [31]. High-quality RNA samples (OD260/280 ratio of 1.8 to 2.0) from each time point (12, 24, 48, and 72 h) with three biological replicates were used for library construction. Library preparation and RNA-Seq, performed by Biomarker Technologies (Beijing, China), involved poly(A) mRNA enrichment using oligo (dT) magnetic beads, followed by fragmentation, first and second-strand cDNA synthesis, purification (QiaQuick PCR Extraction Kit, Qiagen), end-repair, poly(A) addition, and adapter ligation. The resulting libraries were sequenced on the Illumina HiSeqTM 2000 platform using a PE150 strategy. Raw sequencing data were processed using Trimmomatic to remove adapters, reads with >10% unknown nucleotides, and low-quality sequences (base quality ≤10). Due to the lack of a reference genome for Oriental lily, de novo assembly was performed using Trinity software (v2.8.5) with default parameters (K-mer = 31, K-mer cover = 6 for specific runs). Candidate coding regions were identified with TransDecoder (v5.5.0), and ESTscan was used for sequences that could not be aligned. The assembled unigenes were functionally annotated via BLASTx searches (E-value ≤ 1×e^-5, identity > 70%, query coverage ≥ 80%) against multiple databases including Nr, Swiss-Prot, KEGG, GO, and COG. Additionally, protein domains were detected using InterProScan (v5.33-72.0) based on the InterPro database. Read mapping and expression quantification were performed using Bowtie2 and RSEM, respectively, with gene expression levels calculated as FPKM (fragments per kilobase of transcript per million mapped reads).

2.5. Weighted Gene Co-Expression Network Analysis (WGCNA)

Weighted gene co-expression network analysis was performed using the WGCNA package (v1.69) in R. Genes with low expression (FPKM < 1 in more than half of the samples) were filtered out before analysis. A signed network was constructed with a soft threshold power of β = 12, determined by the scale-free topology criterion. Module identification was performed using dynamic tree cutting with a minimum module size of 30 genes and a module similarity threshold of 0.25. Module-trait relationships were calculated by correlating module eigengenes with anthocyanin content data. The gene co-expression network was visualized using Cytoscape (v3.8.0) with a weight threshold of 0.95.

2.6. Phylogenetic Analysis

The full-length coding sequence of LhERF109 was identified from the transcriptome data. Multiple sequence alignment (MSA) of ERF proteins from Arabidopsis thaliana, retrieved from The Arabidopsis Information Resource (TAIR) database (https://www.arabidopsis.org/) [32], was performed using MUSCLE (v3.8.31) with default parameters. The phylogenetic tree was constructed using the maximum likelihood method in MEGA 7.0 with 1,000 bootstrap replicates. The evolutionary distances were computed using the JTT matrix-based method. The tree was visualized and annotated using iTOL (Interactive Tree of Life).

2.7. Transient Transformation of Lily Tepals

For transient transformation experiments, the full-length coding sequence of LhERF109 was amplified by PCR and cloned into the pS1300 expression vector under the control of the CaMV 35S promoter using KpnI and SalI restriction sites. For virus-induced gene silencing (VIGS), a 300-bp fragment of LhERF109 was cloned into the pTRV2 vector. All constructs were verified by sequencing before use.

For transient expression and VIGS experiments, St₂ flower buds of Oriental lily (Lilium ‘Souvenir’) were selected based on their developmental uniformity. Whole flower branches, each bearing a single St₂ bud, were used for infiltration. Approximately 100 μL of the Agrobacterium suspension was injected into the abaxial side of each tepal using a 1-mL needleless syringe, ensuring even distribution of the inoculum. Post-infiltration, the entire flower branch was placed in a container with its stem inserted into sterile distilled water to maintain hydration and support the flower during the experiment. The setup was then transferred to a growth chamber under controlled conditions. The samples were initially subjected to a 12-hour dark treatment at 20°C to minimize stress, followed by culture at 20°C under a 16-hour light/8-hour dark photoperiod [33]. Tepal samples were collected 72 hours post-infiltration for phenotypic observation, anthocyanin content measurement (via spectrophotometric analysis), and gene expression analysis, including quantitative real-time PCR (qRT-PCR) to assess LhERF109 expression levels.

2.8. Quantitative Real-Time PCR (qRT-PCR)

The quantitative expressions of crucial genes involved in anthocyanin synthesis was analyzed by qRT-PCR referring to previous description [34]. The mRNA expression level was assessed with 2^-△△Ct method using an internal reference gene LhActin.

2.9. Statistical Analysis

All experiments were performed with at least three biological replicates. Statistical analysis was performed using SPSS software (version 25.0). Data were analyzed using Student’s t-test or one-way analysis of variance (ANOVA) followed by Tukey’s HSD test. Differences were considered statistically significant at p < 0.05. Data are presented as means ± standard deviation (SD).

3. Results

3.1. Effect of High Temperature on the Physiology and Biochemistry of Oriental Lily ‘Souvenir’

Under the 20°C treatment, the St₂ (Figure 1A) of the lily ‘Souvenir’ exhibited a progressive petal - coloring phenomenon. In sharp contrast, upon initiating the 35℃ treatment, the red coloration of the lily petal began to decline (Figure 1B). Subsequently, this color continued to fade, eventually approaching an almost white state in the lily petals (Figure 1B). Consistent with these visual observations, the total anthocyanin content and Cyanidin 3-O-β-rutinoside content in tepals under 20°C treatment showed a continuous increase, while both decreased steadily in tepals under 35°C treatment (Figure 1C).

The determination of the changes in the sugar content of Oriental lily ‘Souvenir’ flowers under high temperature showed that over time, the sucrose content in the petals generally increased under suitable temperature conditions, while it generally decreased under high temperature conditions (Figure 1D).

To validate the accuracy of RNA-seq detection, this study conducted an investigation into the expression patterns of structural genes that have been firmly established to be involved in the regulation of anthocyanin synthesis. The results (Figure 1E) demonstrated congruence with those derived from the transcriptomic analysis, thereby indicating that the sequencing results obtained in this study exhibit a high level of reliability.

Figure 1. (A) Diagram of the developmental stages of lilies. (B) Comparison of the color of the lily 'Souvenir' petals under different durations of high-temperature treatments T1, T2, T3, and T4. (C) The total anthocyanin content of the lily ‘Souvenir’ petals under high-temperature treatments at T1, T2, T3, and T4. (D) The Cyanidin 3-O-β-rutinoside content of the lily ‘Souvenir’ petals under high-temperature treatments at T1, T2, T3, and T4 (* p < 0.1, ** p < 0.01, *** p < 0.001).

Figure 1. (A) Diagram of the developmental stages of lilies. (B) Comparison of the color of the lily 'Souvenir' petals under different durations of high-temperature treatments T1, T2, T3, and T4. (C) The total anthocyanin content of the lily ‘Souvenir’ petals under high-temperature treatments at T1, T2, T3, and T4. (D) The Cyanidin 3-O-β-rutinoside content of the lily ‘Souvenir’ petals under high-temperature treatments at T1, T2, T3, and T4 (* p < 0.1, ** p < 0.01, *** p < 0.001).

3.2. Gene Expression Differences Between NT- and HT-Treated Oriental Lily

RNA-seq analysis of Oriental lily ‘Souvenir’ flowers subjected to 20°C and 35°C treatments (collected at four time points) generated 143.9 Gb of high-quality clean data with Q30 > 94.25%. In the absence of a reference genome for Lilium Oriental, clean reads were de novo assembled using Trinity software, yielding 63,394 Unigenes with an N50 length of 1.58 kb. Among these, 18,504 Unigenes exceeded 1 kb in length. Comprehensive functional annotation against multiple databases successfully classified 31,101 Unigenes (49.06%), with transcript abundance quantified using FPKM values. Differential expression analysis identified 8,354 differentially expressed genes (DEGs) across four pairwise comparisons (|log₂FC| > 1, adjusted p < 0.01). As illustrated in Figure 2A, the 35°C vs. 20°C comparison revealed distinct temporal expression patterns: at T1, 1,665 genes were up-regulated and 1,055 down-regulated; at T2, 1,581 genes showed up-regulation while 1,107 exhibited down-regulation; at T3, 1,885 genes were up-regulated and 1,854 down-regulated; and at T4, 1,946 genes displayed increased expression while 1,957 showed decreased expression. The Venn diagram analysis (Figure 2B) revealed 738 DEGs consistently differentially expressed across all four time points (T1, T2, T3, and T4; corresponding to 12, 24, 48, and 72 hours post-temperature treatment). Notable time-specific expression patterns were observed, with 931, 569, 684, and 892 unique DEGs identified at T1, T2, T3, and T4, respectively. This temporal diversity in transcriptional profiles indicates a sophisticated and dynamic gene expression reprogramming in response to elevated temperature stress over time.

Figure 2. Histogram and Venn diagram showing numbers of differentially expressed genes (DEGs) in the comparisons 35 °C vs. 20 °C at T1, T2, T3, and T4. (A) Histogram showing numbers of up- and down-regulated DEGs in the comparisons 35 °C vs. 20 °C at T1, T2, T3, and T4. (B) Venn diagram presenting the overlap of DEGs in the comparisons 35 °C vs. 20 °C at T1, T2, T3, and T4. T1, T2, T3, and T4 represent oriental lily after 12, 24, 48, and 72 hours of 20 °C and 35 °C treatment, respectively.

Figure 2. Histogram and Venn diagram showing numbers of differentially expressed genes (DEGs) in the comparisons 35 °C vs. 20 °C at T1, T2, T3, and T4. (A) Histogram showing numbers of up- and down-regulated DEGs in the comparisons 35 °C vs. 20 °C at T1, T2, T3, and T4. (B) Venn diagram presenting the overlap of DEGs in the comparisons 35 °C vs. 20 °C at T1, T2, T3, and T4. T1, T2, T3, and T4 represent oriental lily after 12, 24, 48, and 72 hours of 20 °C and 35 °C treatment, respectively.

3.3. Functional Classification of DEGs

GO and KEGG pathway analyses were performed to understand the biological functions of DEGs across different time points. GO analysis revealed a dynamic response pattern to high-temperature stress (Figure 3). At T1, biological processes were predominantly enriched in stress response terms, including ‘response to hydrogen peroxide’, ‘response to reactive oxygen species’, and ‘response to heat’ (-log₁₀(p-value) > 12.5), indicating the activation of immediate stress response mechanisms. The cellular component analysis showed significant enrichment in membrane system terms, particularly ‘integral component of membrane’. At T2, RNA modification emerged as the most significantly enriched term (-log₁₀(p-value) = 25), while stress response pathways remained active. Notably, cellular components showed substantial changes in ‘intracellular membrane-bounded organelle’ and organelle-related terms. By T3, the enriched terms shifted towards light response and protection mechanisms, with significant enrichment in ‘photoprotection’ and ‘cellular response to red light’, suggesting adaptation of photosynthetic machinery. At T4, ‘RNA modification’ and ‘mitochondrial mRNA processing’ became prominently enriched, accompanied by metabolic regulation terms.

Based on the results of KEGG enrichment analysis, we can clearly observe distinct temporal patterns in metabolic and signaling pathway responses (Figure 4). In the early phase (T1), protein processing in the endoplasmic reticulum pathway was significantly enriched, reflecting how plant cells initially respond to heat damage through protein quality control mechanisms. As stress continued to 24 hours (T2), secondary metabolism-related pathways, particularly flavonoid biosynthesis, began to show significant responses, indicating a transition from immediate protective mechanisms to metabolic adjustment phases. By 48-72 hours (T3-T4), plant-pathogen interaction pathways and defense-related mechanisms became more prominent, while anthocyanin biosynthesis showed specific enrichment at T3, directly correlating with the observed decrease in pigment content. Notably, the plant hormone signal transduction pathway maintained a consistent moderate level of enrichment throughout the stress process, suggesting that hormone signaling plays a continuous regulatory role in coordinating heat stress responses. This temporal metabolic network reorganization reveals the complex adaptation mechanisms of lilies to high temperature stress, involving a complete process from protein homeostasis maintenance to metabolic adjustment and defense system activation. These findings provide important clues for understanding the regulatory network of pigment synthesis, supporting our identification of LhERF109 as a key regulatory factor connecting heat stress response and anthocyanin biosynthesis.

Figure 3. Gene Ontology (GO) enrichment analysis of differentially expressed genes (DEGs) in Oriental lily tepals under high temperature stress. The DEGs were categorized into three main types: biological process (BP), cellular component (CC), and molecular function (MF). T1, T2, T3, and T4 represent samples collected at 3, 6, 9, and 12 days after 20°C and 35°C treatment, respectively. The color scale indicates -log10(p-value) of enriched GO terms, and the size of dots represents the count of genes in each term.

Figure 3. Gene Ontology (GO) enrichment analysis of differentially expressed genes (DEGs) in Oriental lily tepals under high temperature stress. The DEGs were categorized into three main types: biological process (BP), cellular component (CC), and molecular function (MF). T1, T2, T3, and T4 represent samples collected at 3, 6, 9, and 12 days after 20°C and 35°C treatment, respectively. The color scale indicates -log10(p-value) of enriched GO terms, and the size of dots represents the count of genes in each term.

Figure 4. KEGG pathway enrichment analysis of differentially expressed genes (DEGs) in Oriental lily tepals at 35 °C vs. 20 °C comparisons at T1, T2, T3, and T4. T1, T2, T3, and T4 represent tepals after 12h, 24h, 48h, and 72h of temperature treatment, respectively.

Figure 4. KEGG pathway enrichment analysis of differentially expressed genes (DEGs) in Oriental lily tepals at 35 °C vs. 20 °C comparisons at T1, T2, T3, and T4. T1, T2, T3, and T4 represent tepals after 12h, 24h, 48h, and 72h of temperature treatment, respectively.

3.4. Anthocyanin-Related DEGs Revealed by Analysis of Co-Expression Networks

The weighted gene co-expression network analysis (WGCNA) resulted in the identification of eight distinct co-expression modules (Figure 5A and 5B). Among these eight modules, the black module (MEblack) exhibited the strongest positive correlations with cyanidin-3-rutinoside (cy-3) and total anthocyanin content (TCA) (correlation coefficients of 0.9 and 0.69, respectively, P<0.01), suggesting that genes within this module likely play crucial roles in anthocyanin accumulation under temperature stress. The gene co-expression network was constructed using known anthocyanin-related genes, including the key transcription factor MYB12 and structural genes involved in the anthocyanin biosynthetic pathway (ANS, DFR, UGT78, and F3’H), along with potential regulatory transcription factors identified in the black module (Figure 5C). Network topology analysis revealed that MYB12 exhibited high connectivity with several structural genes and other transcription factors, serving as a central hub gene in the network. The structural genes, particularly ANS and DFR, showed strong connections with multiple transcription factors, including members of the ERF, NAC, and BHLH families, suggesting complex regulatory mechanisms controlling anthocyanin biosynthesis under temperature stress. Network visualization was conducted using Cytoscape with a weight threshold of 0.95.

The gene expression heatmap (Figure 5D) further corroborated these findings, showing significant differences in gene expression patterns between 20°C and 35°C treatments. Most genes in the black module, including MYB12, ANS, DFR, and F3’H, exhibited higher expression levels at 20°C and significantly decreased expression at 35°C, consistent with the trends observed in anthocyanin content. Notably, several transcription factors from the ERF, NAC, and WRKY families displayed temperature-sensitive expression patterns, indicating their potential roles as key regulators of temperature-responsive anthocyanin synthesis.

Overall, our WGCNA analysis established a comprehensive co-expression network centered around known anthocyan in biosynthetic genes and identified key regulatory transcription factors associated with temperature-dependent anthocyanin accumulation, providing important candidate genes for further functional validation.

Figure 5. Weighted gene co-expression network analysis (WGCNA) of DEGs in apple at 20°C and 35°C treatment. (A) Hierarchical cluster tree presenting modules of co-expressed genes with color annotation. (B) Module–anthocyanin correlations and corresponding p-values. (C) Connection network between 5 anthocyanin structural genes and 21 TFs (transcript factors) in the ‘black’ module. (D) The heat map for the expression of 12 anthocyanin structural genes and 21 TFs (transcript factors) in the ‘black’ module.

Figure 5. Weighted gene co-expression network analysis (WGCNA) of DEGs in apple at 20°C and 35°C treatment. (A) Hierarchical cluster tree presenting modules of co-expressed genes with color annotation. (B) Module–anthocyanin correlations and corresponding p-values. (C) Connection network between 5 anthocyanin structural genes and 21 TFs (transcript factors) in the ‘black’ module. (D) The heat map for the expression of 12 anthocyanin structural genes and 21 TFs (transcript factors) in the ‘black’ module.

3.5. Phylogenetic Tree Construction of LhERF109

Phylogenetic analysis positioned LhERF109 within ERF subfamily Group X, clustering with AtERF109 (At4g34410) with high bootstrap support (92%) (Figure 6A).

LhERF109 plays a critical role in regulating anthocyanin biosynthesis in Oriental lily ‘Souvenir’. Compared with the control (S1300), transient overexpression of LhERF109 (LhERF109-S1300) significantly enhanced anthocyanin accumulation in tepals, resulting in deeper purple-red pigmentation (Figure 7A). Conversely, silencing of LhERF109 (LhERF109-TRV) markedly reduced anthocyanin accumulation, leading to lighter coloration compared to the TRV control (Figure 7A). Quantitative RT-PCR confirmed the successful overexpression and silencing of LhERF109 in the respective treatments (Figure 7B).

Consistent with the phenotypic changes, anthocyanin biosynthetic genes showed differential expression patterns across treatments. Early and middle-stage genes including LhMYB12, LhWD40a, LhCHSb, and LhANS were significantly upregulated in LhERF109-S1300 compared to the control, while they were downregulated in the LhERF109-TRV treatment (Figure 7C). Notably, LhBHLH2 expression remained relatively stable in the overexpression treatment but was significantly reduced in the silenced group.

For late-stage and transport-related genes, most followed a similar pattern, with LhUGT2, LhUGT3, and LhGSTF10 showing strong upregulation in LhERF109-S1300 and downregulation in LhERF109-TRV (Figure 7D). Interestingly, LhF3H displayed an opposite regulation pattern, being downregulated in LhERF109-S1300 and upregulated in LhERF109-TRV. LhDFR exhibited unique behavior with upregulation in both overexpression and silencing treatments, though through potentially different mechanisms. LhF3’H showed minimal response to LhERF109 manipulation, suggesting alternative regulatory control.

These results demonstrate that LhERF109 functions as a positive regulator of anthocyanin biosynthesis in Oriental lily by modulating the expression of key structural and regulatory genes in the anthocyanin pathway, with gene-specific regulatory patterns that suggest a complex control network.

Figure 6. Phylogenetic analysis of ERF transcription factors from Lilium and Arabidopsis thaliana.

Figure 6. Phylogenetic analysis of ERF transcription factors from Lilium and Arabidopsis thaliana.

Figure 7. LhERF109-regulated anthocyanin biosynthesis in Oriental lily ‘Souvenir’. (A) The phenotypes of S1300, LhERF109-S1300, TRV, and LhERF109-TRV lily tepals. (B) The expression levels of LhERF109 in S1300, LhERF109-S1300, TRV, and LhERF109-TRV treatments. (C) The expression levels of early and middle-stage anthocyanin biosynthesis genes in S1300, LhERF109-S1300, TRV, and LhERF109-TRV treatments. (D) The expression levels of late-stage anthocyanin biosynthesis and transport-related genes in S1300, LhERF109-S1300, TRV, and LhERF109-TRV treatments. Values represent means ± SD of three independent biological replicates. Statistical significance was computed with the use of Student’s t-test (* p < 0.1, ** p < 0.01, *** p < 0.001).

Figure 7. LhERF109-regulated anthocyanin biosynthesis in Oriental lily ‘Souvenir’. (A) The phenotypes of S1300, LhERF109-S1300, TRV, and LhERF109-TRV lily tepals. (B) The expression levels of LhERF109 in S1300, LhERF109-S1300, TRV, and LhERF109-TRV treatments. (C) The expression levels of early and middle-stage anthocyanin biosynthesis genes in S1300, LhERF109-S1300, TRV, and LhERF109-TRV treatments. (D) The expression levels of late-stage anthocyanin biosynthesis and transport-related genes in S1300, LhERF109-S1300, TRV, and LhERF109-TRV treatments. Values represent means ± SD of three independent biological replicates. Statistical significance was computed with the use of Student’s t-test (* p < 0.1, ** p < 0.01, *** p < 0.001).

4. Discussion

High temperature significantly inhibits anthocyanin accumulation in plant tissues, but the regulatory mechanisms differ among species and tissues [13]. Our transcriptome analysis of Oriental lily tepals under heat stress revealed a complex metabolic network reconstruction. KEGG analysis showed that protein processing in endoplasmic reticulum was most significant in the T1 stage, while flavonoid biosynthesis, anthocyanin biosynthesis, and phenylpropanoid metabolism were continuously enriched from T2 to T4 stages. In addition, circadian rhythm, plant hormone signal transduction, and chlorophyll metabolism pathways also showed significant enrichment. The GO enrichment analysis further revealed a progressive transition from early stress response (ROS-related terms) to later metabolic adaptation (membrane system modifications) in lily tepals. This temporal pattern differs from what has been observed in apple fruits, where carbon metabolism was consistently the most enriched pathway [8,35]. Our findings also contrast with observations in grape berries, where heat stress primarily affects flavonoid pathway genes through direct transcriptional suppression [36]. In the chrysanthemum study, functional annotation and enrichment analysis revealed significant enrichment of heat shock response and flavonoid biosynthesis pathways, with 63 differentially expressed genes associated with ‘flavonoid biosynthesis’ highlighting the temperature sensitivity of this pathway [37]. These multi-species comparisons suggest that the regulation of heat-induced anthocyanin accumulation varies significantly among different species.

High temperature stress significantly affects sugar metabolism and anthocyanin synthesis in the petals of Oriental lily ‘Souvenir’. Our research demonstrates that under high temperature conditions, sucrose content consistently decreases while soluble sugars show a dynamic pattern of initial decrease followed by increase, which aligns with the significant enrichment of heat response and abiotic stress response pathways in GO analysis. Sucrose, serving as both a precursor and signaling molecule for anthocyanin synthesis, directly influences anthocyanin accumulation [38,39]. Heat stress disrupts the normal relationship between sucrose and anthocyanin synthesis, leading to significant inhibition of anthocyanin biosynthesis. The accumulation of soluble sugars serves as a protective mechanism to regulate osmotic pressure and lower metabolic rates under high temperature [40,41]. These findings provide insights for enhancing plant heat tolerance through exogenous sugar applications [42,43].

Recent research by Yang et al. [4] suggests that heat stress redirects metabolic resources from anthocyanin biosynthesis toward isoflavone synthesis in lily petals. However, our results challenge this interpretation. Our analysis revealed that total anthocyanin content significantly decreases under high temperature (35°C), with cyanidin, a major anthocyanin in lily tepals, showing particularly severe suppression compared to normal conditions (20°C). This overall reduction in anthocyanin content indicates that multiple mechanisms beyond simple metabolic redirection are involved in the heat response. This pattern resembles findings in cereal crops, where heat stress triggers coordinated adjustments in primary metabolism, hormone signaling, and secondary metabolite production [44]. These findings highlight the importance of species-specific approaches in heat tolerance breeding.

Through WGCNA analysis and co-expression network construction within the black module, we successfully identified LhERF109 as a key regulator associated with anthocyanin-related genes. Based on comprehensive evidence across multiple species and our experimental findings, we propose that LhERF109 functions as a critical molecular integrator connecting heat stress responses with anthocyanin biosynthesis in Oriental lily. Phylogenetic analysis positions LhERF109 within the Group X ERF subfamily, clustering closely with Arabidopsis AtERF109 with high bootstrap support (92%), indicating evolutionary conservation with Group IX ERFs, which are known to mediate stress responses and hormone signaling pathways such as jasmonic acid (JA), indole-3-acetic acid (IAA), and ethylene (ETH) [45,46,47,48,49]. Transient transformation experiments demonstrated that LhERF109 overexpression (LhERF109-S1300) significantly enhanced tepal anthocyanin accumulation, evidenced by deeper purple-red pigmentation (Figure 7A), with robust upregulation of the regulatory gene LhMYB12 and structural genes including LhCHSb, LhDFR, LhUGT2, LhUGT3, and LhGSTF10, as well as LhANS (Figure 7C, D). In contrast, silencing (LhERF109-TRV) resulted in lighter coloration and reduced anthocyanin content, accompanied by significant downregulation of LhMYB12, LhBHLH2, LhWD40a, LhUGT2, LhUGT3, LhGSTF10, and LhANS (Figure 7A, C, D). The expression of LhBHLH2 exhibited an asymmetric pattern, remaining stable during overexpression but significantly downregulated upon silencing (Figure 7C), suggesting a complex regulatory mechanism beyond simple transcriptional activation.

The pronounced upregulation of LhANS during overexpression suggests that LhERF109 may directly target anthocyanin biosynthesis genes, enhancing the conversion of leucoanthocyanidins to anthocyanidins, a rate-limiting step in anthocyanin biosynthesis. Notably, the significant induction of late-stage anthocyanin biosynthesis genes, particularly transport-related genes such as LhUGT2, LhUGT3, and LhGSTF10, suggests that LhERF109 may play a critical role in facilitating anthocyanin transport and stabilization under heat stress conditions. This regulatory mechanism is supported by similar observations in other species, such as MdERF1B in apple, which binds MdANS and MdUFGT promoters [47], and MdERF38, which interacts with MdDFR and MdUFGT under drought stress [50].

The asymmetric regulation of LhBHLH2, with minimal upregulation during overexpression and significant downregulation during silencing, mirrors patterns in other species. In apple, MdMYB10 overexpression strongly activates MdUFGT with little effect on MdBHLH33, whereas RNAi silencing significantly downregulates MdBHLH33 [50]. In grape, VvMYBA1 overexpression promotes VvUFGT with stable VvBHLH1, while silencing markedly reduces VvBHLH1 [51]. Similarly, in petunia, PhMYB27 overexpression enhances anthocyanin with limited PhBHLH3 upregulation, but silencing leads to a pronounced PhBHLH3 reduction [52], and in Arabidopsis, AtMYB75 overexpression induces DFR with minimal AtBHLH42 (GL3) change, while silencing significantly downregulates AtBHLH42 [53]. These parallels suggest that LhERF109 may prioritize MYB and ANS-mediated activation under heat stress, while indirectly maintaining BHLH, a strategy potentially conserved across species.

In conclusion, our study demonstrates that ERF109, as a heat-responsive transcription factor, influences lily flower coloration by regulating multiple stages of the anthocyanin biosynthesis pathway. This finding not only enriches our understanding of plant color regulation but also provides potential molecular targets for breeding heat-resistant varieties with stable flower coloration. Future research should focus on the potential interaction mechanisms between ERF109 and the MBW complex, as well as how high temperatures affect ERF109 expression and activity, to further understand the temperature-transcription factor-color formation regulatory network.

5. Conclusions

This study reveals the molecular mechanisms of flower color changes in Oriental lily under heat stress through transcriptome analysis and Weighted Gene Co-expression Network Analysis (WGCNA). We found that high temperature inhibits anthocyanin accumulation, leading to color fading. Using WGCNA, we identified LhERF109 as a key transcription factor regulating this process for the first time and verified its regulatory effects on critical anthocyanin biosynthesis genes. Transient transformation experiments demonstrated that LhERF109 overexpression significantly enhanced tepal anthocyanin accumulation, while silencing resulted in lighter coloration and reduced anthocyanin content. Future research will focus on the precise mechanisms of LhERF109, providing support for breeding heat-tolerant ornamental crops with stable flower coloration.