The high resistant Lilium oriental hybrid cultivar “Sorbonne” (Gao et al., 2018), was used in this study. Bulbs were stored in a refrigerator (4°C) for 28 days, and then placed at the bottom of 15-cm deep pots filled with turf: vermiculite: perlite substrate (1:1:1, v/v/v). The pots containing the bulbs were placed in a greenhouse for 45 days at 25°C under 12-h light/12-h dark photoperiod, with photosynthetic photon flux density of 240 μmol·m.2·s-1 (Jang et al., 2018).

Botrytis elliptica strain 36423, isolated from symptomatic Lilium plants, was purchased from the Agricultural Cultural Collection of China (http://www.accc.org.cn/). The mycelium was cultured on potato dextrose agar (PDA, pH 5.8; Coolaber, China) medium in Petri dishes (9 cm diameter) at 25°C in the dark for 1 week.

Fully expanded, but not senescent, leaves were collected from Lilium plants at the flower bud stage, and inoculated with B. elliptica, according to the “detached leaves inoculation methods” of Gao et al. (2018). Before inoculation, all utensils and water were autoclaved at 121°C for 20 min, and the detached leaves were wiped clean using cotton soaked with sterile water. To inoculate the detached leaves, B. elliptica mycelium discs (5-mm diameter) were collected from the PDA plates using a sterilizing puncher, and then used to inoculate the abaxial surface of the detached lily leaves in vitro. Each leaf was inoculated with six mycelium discs (inoculated treatment), while those uninoculated with B. elliptica served as the control treatment. The inoculated and uninoculated leaves were placed in Petri dishes (15 cm diameter) lined with moist filter paper. The filter paper and cotton surrounding the leaf petiole were soaked with sterile water to maintain humidity within a range of 90–100%.

The detached leaves were sampled and photographed at 6, 8, 12, 24, 36, 48, and 72 h post inoculation (hpi). The area of each lesion was measured with the ImageJ software (https://imagej.nih.gov/ij/). Superoxide dismutase (SOD) activity was assayed by monitoring the photoreduction of nitro blue tetrazolium (NBT), as describe previously (Li, 2000). Based on the results of the SOD activity assay and the phenotype of disease lesions, leaves in the control, and inoculated treated treatments were sampled at 6, 24, and 48 hpi (hereafter referred to as control_6 h, control_24 h, and control_48 h, and inoculated_6 h, inoculated_24 h, and inoculated_48 h, respectively) for transcriptome and metabolome sequencing. Three biological replicates were performed for each treatment, with each replicate containing three technical repeats.

The leaf samples were freeze-dried under vacuum using the Scientz-100F lyophilizer (Scientz, China), and then crushed at 30 Hz for 1.5 min using a grinder (MM 400; Retsch, Germany) to obtain a fine powder. Then, 100 mg of each powdered sample was extracted in 0.6 ml of 70% methanol. The samples were stored at 4°C overnight, during which time they were vortexed six times to hasten the extraction. Then, each sample was centrifuged at 10 000 × g for 10 min. The supernatant was filtered through a microporous membrane (pore size: 0.22 μm), and stored in sample bottles for ultra-performance liquid chromatography-tandem mass spectrometry (UPLC-MS/MS) analysis.

To analyze the extracts, UPLC (Shim-pack UFLC SHIMADZU CBM30A; https://www.shimadzu.com.cn/) and MS/MS (Applied Biosystems 4500 QTRAP; http://www.appliedbiosystems.com.cn/) were performed using a C18 chromatographic column (Waters ACQUITY UPLC HSS T3 C18; 1.8 μm, 2.1 × 100 mm) with solvent A (0.04% acetic acid in ultrapure water) and solvent B (0.04% acetic acid in acetonitrile) as the mobile phase. A 4-μl aliquot of each sample was injected into the column, and eluted using the following gradient program: 0 min with 95% A and 5% B; 0–10 min with 95–5% A and 5–95% B; 10–11 min with 5% A and 95% B; 11–11.1 min with 5–95% A and 95–5% B; and 11.1–14 min, 95% A and 5% B. The flow rate and column temperature were maintained at 0.35 mL/min and 40°C, respectively. The mass spectrometer parameters were set as follows: temperature of electrospray ionization: 550°C; voltage: 5.5 kV; curtain gas: 30 psi; collision-activated dissociation: high. Each ion pair was scanned and detected based on optimized declustering potential and collision energy in the triple quadrupole system.

Qualitative and quantitative analyses of the metabolites were performed using the multiple reaction monitoring (MRM), KEGG compound database, and MetWare database. Metabolites were identified based on their molecular weight, Mass Spectrometry (MS²) fragments, MS² fragments isotope distribution, and retention time (RT). Through the MetWare self-developed intelligent secondary spectrum matching method, the secondary spectrum and RT of the metabolites in the project samples are intelligently matched one by one with the MetWare database. The MS tolerance and MS² tolerance are set to 2 and 5 ppm, respectively. The peak area integral of all the mass spectrum peaks was derived after obtaining the metabolic substance spectrum analysis data of different samples, followed by an integral correction performed for the mass spectrum peak of the same metabolite occurring in different samples (Fraga et al., 2010). Quality control (QC) samples, i.e., samples prepared from a mixture of sample extracts, were used to analyze the reproducibility of the instrument under the same treatment method.

The “MetaboAnalystR 1.0.1” package in the R computing platform v 3.5.0 (R Core Team, 2018, Austria) was used to statistically analyze the metabolomics data and to generate plots. Orthogonal partial least squares discriminant analysis (OPLS-DA) was conducted using MetaboAnalystR to identify differentially accumulated metabolites (DAMs), with Variable Importance in Projection (VIP) score ≥ 1 and absolute Log₂fold change (FC) ≥ 1. Pathways with significantly regulated metabolites mapped to it were then subjected to metabolite set enrichment analysis; their respective statistical significance was determined using the hypergeometric test and its p-value.

A total of 18 cDNA libraries were sequenced using the NEBNext® Ultra™ RNA Library Prep Kit for Illumina® (NEB; USA). Clean reads obtained from each cDNA library were assembled de novo into unigenes using the Trinity (2.6.6.) platform (Haas et al., 2013). Functional annotations of the unigenes were determined using NCBI non-redundant (NR) protein, UniProt, KEGG, Gene Ontology (GO), and Clusters of Orthologous Groups of proteins (COG) databases (Kanehisa et al., 2008). Transcription factors and gene coding sequence (CDS) were used ITAK web online (http://itak.feilab.net/cgi-bin/itak/index.cgi) and TransDecoder (5.3.0, https://github.com/TransDecoder) software, respectively. Differentially expressed genes (DEGs) were analyzed using the DESeq2 (v1.22.2) software package in R, and subjected to the Benjamini–Hochberg method for multiple hypothesis testing (i.e., |log2FC| ≥ 1, FDR [false discovery rate] <0.05) (Love et al., 2014, Varet et al., 2016). Heat maps were constructed using the R package “pheatmap” (v1.0.12) and TBtools software v0.66836 (Chen et al., 2020). Venn diagrams were generated using Venn v1.6. Weighted gene coexpression network analysis (WGCNA) was performed using the R package “WGCNA” (Zhang and Horvath, 2005), and visualized using Cytoscape v1.7.251 (https://cytoscape.org/index.html).

To validate the RNA-seq data, the expression of 10 defense-related DEGs, including four DEGs involved in plant–pathogen interactions (Cui et al., 2018a) and six DEGs involved in the phenylpropanoid and flavonoid biosynthesis pathways, was analyzed by qRT-PCR. The qRT-PCR was performed using the SsoFast EvaGreen Supermix (Bio-Red, USA) on the CFX96 Real-Time PCR Detection System (Bio-Rad, USA) under the following conditions: 95°C for 30 s, followed by 39 cycles of 95°C for 5 s and 57°C for 5 s, melt curve 65–95°C, increment 0.5°C for 5s. The specificity of the primers was verified based on the unimodality of the melt curve. Actin (ACT) and elongation factor 1 (EF1) served as reference genes (Cui et al., 2018a). Three biological replicates were performed for each treatment, with each replicate containing three technical repeats. Data were analyzed using the 2^−ΔΔCt method (Livak and Schmittgen, 2001). Primers used for qRT-PCR are listed in Supplementary Table 1.

Water-soaked lesions of the same diameter (0.5 cm) as the plug used for inoculation were observed at 12 hpi (Figure 1A). Over time, the lesions first became rotten and brown (36 hpi) and then turned necrotic (48–72 hpi). The area of lesions expanded with time after inoculation (Figure 1B). Changes in SOD activity were also observed in inoculated leaves (Figure 1C). Significant differences were detected in SOD activity between the control_6 h and inoculated_6 h samples. Superoxide dismutase activity declined before 8 hpi but then increased over time. The sampling time points for metabolomics and transcriptomics were determined based on both the disease symptoms and SOD activity. Given that transcriptional changes usually precede the associated physiological and phenotypic changes, here we considered 6 hpi as the early stage of infection. As the regular ellipse-shaped lesions formed, the 24- and 48-hpi time points were considered as the middle and late stages stage of infection, respectively, for sampling.

To determine the metabolic changes induced by B. elliptica infection, we performed non-targeted metabolome analysis using control_6 h, control_24 h, control_48 h, inoculated_6 h, inoculated_24 h, and inoculated_48 h samples. A total of 524 metabolites common to all samples were identified based on their chromatographic and mass spectrometric parameters (Supplementary Table 2). Orthogonal partial least squares discriminant analysis and principal components analysis (PCA) uncovered differences in the metabolites of all samples (Supplementary Figures 1, 2). High correlation was observed among the QC samples (Supplementary Figure 3). Pearson correlation coefficients were consistently high (r > 0.819) across all three biological replicates (Supplementary Figure 4). Heat map, based on the hierarchical clustering analysis of metabolite levels, revealed significant differences in metabolites levels between inoculated and control treatments, and these differences became more pronounced over time (Supplementary Figure 5).

The DAMs were similarly identified by comparing the different temporal stages of pathogen infection, with cutoff values of log₂ FC ≥ 1 and VIP score ≥ 1. Venn diagram shows the overlap of DAMs among the three stages of infection (Supplementary Figure 6; Supplementary Table 3). A total of 115 DAMs were detected, of which 8, 12, and 62 DAMs were expressed only at the 6-, 24- hpi, and 48-hpi time points. Only eight DAMs, including benzyl salicylate, eriodictyol, tryptamine, esculetin, butin, caffeic acid, dihydrokaempferol, and N'-feruloyl putrescine, were common to all three infection stages (Supplementary Figure 7). Among these eight DAMs, five were involved in phenolic acid metabolism and are known as important secondary metabolites for pathogen resistance (Supplementary Figure 7). Moreover, the accumulation of 15 metabolites was up-regulated and that of 6 metabolites was down-regulated in inoculated_6 h compared with control_6 h (Supplementary Figure 8). Similarly, 39 and 4 metabolites were up-regulated and down-regulated, respectively, in inoculated_24 h compared with control_24 h, and 85 and 5 metabolites were up-regulated and down-regulated, respectively, in inoculated_48 h compared with control_48 h (Supplementary Figure 8). Furthermore, according to KEGG enrichment analysis, the DAMs were enriched in phenylpropanoid biosynthesis, flavonoid biosynthesis, metabolic pathway, and biosynthesis of secondary metabolites at all three infection stages. Phenylpropanoid biosynthesis was significantly enriched at 6 hpi (Supplementary Figure 9A); tryptophan and indole alkaloid biosynthesis were significantly enriched at 24 hpi (Supplementary Figure 9B); and isoflavonoid biosynthesis, fructose and mannose metabolism, and galactose metabolism were significantly enriched at 48 hpi (Supplementary Figure 9C).

Samples were used for non-targeted metabolome analysis were subjected to RNA-seq to profile the genome-wide changes in gene expression upon the inoculation of leaves with B. elliptica. After removing low-quality reads, 435.16 Gb clean reads were obtained, with an average GC content of 49.87% (Supplementary Table 4). A total of 430,835 transcripts were obtained, averaging 701 bp in length, with N50 and N90 values of 1,213 and 266 bp, respectively (Supplementary Table 5). Using KEGG, NR, Swiss-Prot, GO, and Trembl databases, a total of 283,213 unigenes, with an average length of 534 bp, were functionally annotated (Supplementary Tables 5, 6). In the NCBI NR database, Asparagus officinalis (8.94%), Elaeis guineensis (8.72%), Phoenix dactylifera (7.39%), Vitis vinifera (3.63%), and Cajanus cajan (3.45%) gave the top BLASTx hits (Supplementary Figure 10). Hierarchal clustering analysis was performed to examine the significant changes in unigene expression. The results showed a stage-specific transcriptome profile after inoculation with B. elliptica (Supplementary Figure 11). To validate the RNA-seq data, the expression profiles of 10 defense-related DEGs were evaluated in inoculated and control leaves at 6, 12, 24, 36, and 48 hpi by qRT-PCR (Supplementary Figure 12). All RNA-seq data can be downloaded from NCBI (BioProjects: PRJNA742853).

The DEGs were identified by comparing the RNA-seq data of inoculated leaf samples (inoculated_6 h, inoculated_24 h, and inoculated_48 h) with those of control leaf samples (control_6 h, control_24 h, and control_48 h) using cutoff values of log₂(FC) ≥ 1 and Padj ≤ 0.05. All DEGs identified at the three temporal stages were analyzed using the K-means clustering algorithm. The DEGs could be grouped into six clusters and classified into four types (Figure 2A): up-regulated (clusters I, II); down-regulated (clusters IV, V); first up-regulated, then down-regulated (clusters VI); and first down-regulated, then up-regulated (cluster III) (Figure 2A). Gene ontology enrichment analysis showed that one cytomembrane-related term, two chloroplast-related terms, two cytoderm-related terms, three signal transduction receptor-related terms, and six photosynthesis-related terms were enriched among the down-regulated genes (Figure 2B). This suggests that certain biological processes, such as photosynthesis, chlorophyll synthesis, cell wall biogenesis, anchored component of membrane, and signal transduction, might be inhibited by the destruction of leaf tissue by the necrotrophic pathogen. The up-regulated gene clusters were found to be related to chitinase, oxidoreductase activity, cell recognition, secondary metabolites, and the phenylpropanoid biosynthetic process, pointing to their potential positive role in the production of resistant metabolites and proteins in B. elliptica-inoculated leaves. The results of KEGG enrichment analysis were consistent with those of GO enrichment analysis: Clusters I and II were mainly enriched in the phenylpropanoid biosynthesis, flavonoid biosynthesis, plant hormone signal transduction, and plant–pathogen interaction pathways; clusters IV and V were mainly enriched in the pathways of photosynthesis and starch and sucrose metabolism (Figure 2C). Taken together, the up-regulated clusters resolved by the K-means clustering algorithm potentially play a pivotal part in the resistance to B. elliptica, while the down-regulated clusters related to photosynthesis, cell wall biosynthesis, and other metabolic pathways seem to be negatively affected by B. elliptica.

A Venn diagram was used to demonstrate the DEGs identified at all three stages. Overall, 5,295 DEGs were detected, of which 125, 434, and 3,089 were uniquely expressed at 6, 24, and 48 hpi, respectively, and 72 DEGs were common to all three stages (Supplementary Figure 13; Supplementary Table 7). Kyoto Encyclopedia of Genes and Genomes enrichment analysis was performed to further characterize these DEGs. At 6 hpi, 401 DEGs (86 up-regulated, 315 down-regulated) were significantly enriched in “biosynthesis of secondary metabolites”, “biosynthesis of amino acid”, and “carbon metabolism” pathways (Supplementary Figures 14, 15A). At 24 hpi, 1,887 DEGs (309 up-regulated, 1,578 down-regulated) were significantly enriched in six pathways: “biosynthesis of secondary metabolites”, “plant hormone signal transduction”, “phenylpropanoid biosynthesis”, “flavonoid biosynthesis”, “isoquinoline alkaloid biosynthesis”, and “MAPK signal pathway” (Supplementary Figures 14, 15B). At 48 hpi, 4,726 DEGs (2,772 up-regulated, 1,954 down-regulated) were significantly enriched in five pathways: “biosynthesis of secondary metabolites”, “plant hormone signal transduction”, “plant–pathogen interaction”, “flavonoid biosynthesis”, and “carbon metabolism” (Supplementary Figures 14, 15C). Taken together, these results indicate that “biosynthesis of secondary metabolites”, “flavonoid biosynthesis”, “phenylpropanoid biosynthesis”, “plant hormone transduction signal”, “plant-pathogen interaction”, and “MAPK signal pathway” were the chief pathways underpinning the responses of “Sorbonne” to B. elliptica infection at the transcript level.

The phenylpropanoid and flavonoid biosynthesis pathways could be integrated into a transcriptional cascade and metabolic network, which together play a critical role in how the Lilium hybrid “Sorbonne” responds to B. elliptica infection. All the DEGs and DAMs participating in the phenylpropanoid and flavonoid biosynthesis pathways that exhibited up- or down-regulated trends in the different infection stages were selected and mapped to this network (Figure 3; Supplementary Figure 16).

Phenylalanine is the first metabolite in the phenylpropanoid pathway. Although our non-targeted metabolome analysis indicated a significant difference in its accumulation between the inoculated and control treatments at 48 hpi, we identified significant expression of phenylalanine ammonia-lyase (PAL), a structural gene that catalyzes phenylalanine, clustering at 24 hpi. The expression of PAL (DN129152_c1_g2) in the inoculated treatment was 3.5- and 6.5-fold higher at 24 and 48 hpi, respectively, compared with the control (Supplementary Figure 16; Supplementary Table 2), implying that PAL participates in the defense response. Moreover, the levels of phenylalanine derivatives, such as phenethylamine and 4-hydroxybenzoic acid, were significantly higher in the inoculated leaves at 24 and 48 hpi compared with the control, indicating that the phenylpropanoid pathway was activated by the infection. Besides, the same expression pattern was exhibited by the structural gene (DN134629_c3_g2) that catalyzes 4-hydroxybenzoic acid. The accumulated level of caffeic acid, a precursor of the lignin pathway, was 1.6-, 3.8-, and 11-fold higher in inoculated leaves at 6, 24, and 48 hpi, respectively, compared with the control (Supplementary Table 2). The level of esculetin, a downstream metabolite of caffeic acid, increased rapidly in leaves after inoculation (Supplementary Table 2). Among the downstream metabolites of caffeic acid, ferulic acid and coniferyl alcohol, which act as precursors of guaiacyl, were down-regulated at 6 hpi but up-regulated at 24 hpi and 48 hpi. Sinapinaldehyde and sinapic acid, precursors of syringyl lignin, were up-regulated in all three stages. The accumulation of syringaresinol-hex and syringaresinol, downstream metabolites of S-lignin, increased significantly in leaves post inoculation (Supplementary Table 2).

The accumulation of naringenin chalcone, naringenin, and dihydrokaempferol, which represent metabolites in the flavonoid biosynthesis pathway, increased significantly at different stages of infection. Additionally, the expression of chalcone synthase (CHS), a structural gene that promotes the synthesis of the three abovementioned metabolites, was induced at all three stages (Supplementary Figure 16; Supplementary Table 2). Among all the DAMs detected in the flavonoid biosynthesis pathway, the greatest difference in accumulation occurred in eriodictyol, which was undetectable in all the control treatments and was 5- and 37-fold higher at 24 and 48 hpi, respectively, than at 6 hpi in the inoculated treatments (Supplementary Table 2). The accumulation level of kaempferol in the inoculated treatment was 190- and 14-fold higher than that in the control treatment at 6 and 24 hpi, respectively, but its level was similar between the inoculated and control treatments at 48 hpi (Supplementary Table 2). The expression of UDP-glycosyltransferase (DN130728_c1_g2) in the inoculated treatment increased significantly at 6 hpi, which could potentially explain the accumulation of kaempferol derivatives at 24 and 48 hpi. The accumulation of four downstream metabolites of kaempferol, namely kaempferol 3-O-β-d-glucopyranoside, kaempferol 3-O-β-glucosyl(1 → 2)(6′-O-acetyl)-β-D-galactoside, kaempferol 7-O-rhamnoside, and kaempferol 3-O-β-(2″-O-acetyl-β-D-glucuronide), was significantly increased in the inoculated treatment at 48 hpi compared with the corresponding control (Figure 3). In addition, the accumulation of the derivatives of luteolin (6-hydroxyluteolin 5-glucoside) and hesperetin (hesperetin 5-O-glucoside) in the inoculated treatment increased significantly at 24 and 48 hpi. The level of diosmetin, a downstream product of luteolin, was five-fold higher in the inoculated treatment than in the control treatment.

To investigate the gene regulatory network of Lilium Sorbonne after its inoculation with B. elliptica, we performed WGCNA of all DEGs identified from the RNA-seq analysis of the 18 cDNA libraries. After preprocessing the RNA-seq data, 14 gene co-expression modules, each comprising 80–2,659 genes, were discovered (Figure 4A). Module-trait and sample relationship analyses showed the eigengenes of these modules were correlated with the different infection stages. Unlike most modules in which the trend of genes was that of almost no difference between the inoculated and control treatments, a few modules did differ significantly during the defense response (Figure 4B). Next, we identified key genes (i.e., hub genes) that played crucial roles during the infection, based on the WGCNA for these notable modules. The levels of gene expression in the black, dark-green, and green modules showed peak up-regulation at 6, 24, and 48 hpi, respectively, which correspond to the early, middle, and late stages of infection, respectively (Figure 5). Gene expression in the dark-gray module was higher in the inoculated treatment than in the control treatment at all-time points, indicating that genes in this module play a prominent role in the defense response against B. elliptica (Figure 6).

In the black module (Figure 5A), seven hub genes related to plant defense were identified, such as two LURP homologs (DN126030_c0_g1, DN126030_c0_g2), one MYB30 homolog (DN134583_c1_g1), and one AS1 homolog (DN125418_c0_g1). Among these, DOF5.5 (DN129414_c0_g2) showed a regulatory relationship with the three defense response-related hub genes. In addition, one WRKY70 gene (DN139965_C0_G2) showed significant difference in expression between inoculated and control treatments at 6 hpi. Additionally, an ERF homolog (DN131819_c0_g10) and a HCT homolog (DN141355_c0_g2) were up-regulated at the early stage of infection. In the dark-green module (Figure 5B), we identified three hub genes related to plant–pathogen interactions: a PBL19 homolog (DN136771_c2_g1), a CRK2 homolog (DN142167_c0_g1), and a WAK5 homolog (DN136340_c0_g1). Notably, PBL19 showed a strong regulatory relationship with other two hub genes, CRK2 and WAK5. A CKX9 homolog (DN129547_c0_g1), a I2′H homolog (DN140336_c1_g1), and a CAD homolog (DN127300_c1_g1) were significantly up-regulated in the dark-green module. At the late stage of the infection (Figure 5C), 11 hub genes were identified, including a WRKY22 homolog (DN139438_c8_g1), a WRKY30 homolog (DN126506_c0_g1), three JAZ homologs (DN139462_c0_g2, DN143599_c2_g1, DN138678_c2_g5,), a DELLA homolog (DN126417_c0_g11), a MAPK17/18 homolog (DN137625_c0_g), and a CML27 homolog (DN125403_c0_g1). Furthermore, both an RPS2 homolog (DN131875_c0_g1) and an RPM1 homolog (DN143855_c2_g1) displayed strong regulatory relationships with other genes in the green module. Genes in the dark-gray module (2,659, i.e., more than half of all identified DEGs) were up-regulated at all three infection stages (Figure 6A). Analysis of DEGs revealed 41 hub genes, which included a MYB61 homolog (DN137495_c4_g1), a WRKY33 homolog (DN123783_c0_g1), MYC2 homolog (DN143930_c9_g2), and a BRI1 homolog (DN141312_c0_g1), among others (Figures 6B,C). A number of hub genes were structural genes involved in the phenylpropanoid and flavonoid biosynthesis pathways (e.g., two PAL homologs [DN129152_c1_g2, DN143696_c1_g2], an HCT homolog [DN139546_c0_g1], two CHS homologs [DN143809_c2_g2, DN122144_c0_g1], an F3′H homolog [DN130296_c0_g1], and a FLS homolog [DN128107_c0_g1], and plant–pathogen interactions (e.g., three FLS2 homologs [DN135524_c1_g1, DN135020_c4_g1, DN143767_c1_g3] and an RIN4 homolog [DN137830_c2_g1], in addition to several PR genes. Lastly, a gene called Pti5 (DN139615_c5_g1), which encoded an AP2-EREBP family TF, also exerted a great influence on plant–pathogen interactions.

Phenolic compounds, such as metabolites of the phenylpropanoid and flavonoid pathways, act as fungitoxic and antimicrobial defense compounds in host–pathogen interactions (Martinez et al., 2017; Xu et al., 2018, 2019). Caffeic acid is a key metabolite in the phenylpropanoid pathway. In our study, significantly more caffeic acid accumulated over time in the inoculated treatment than in the control treatment. Consistent with this result, the structural genes involved in the synthesis of caffeic acid (HCT and CSE) were significantly up-regulated in the inoculated treatment (Figure 3; Supplementary Figures 12, 16). In apple, caffeic acid has been shown to effectively promote lignin accumulation and inhibit gray mold infection (Zhang et al., 2020). Although the accumulation of lignin was not detected in our study, we detected significant accumulation of the precursors precursor of G-lignin (coniferyl alcohol) and S-lignin (sinapic acid and sinapinaldehyde) (Figure 3). The structural genes for the synthesis of lignin, COMT, CCR, and CAD, were also significantly up-regulated in the inoculated treatment (Supplementary Figure 16). Additionally, two transcripts, DN141508_c1_g3 and DN141508_c1_g4, encoding the CCoAOMT enzyme, were significantly up-regulated at all stages of infection (Supplementary Table 7). CCoAOMT is reportedly a lignin synthase that can enhance the resistance of Lilium to B. cinerea infection (Fu et al., 2020). Our results showed that caffeic acid may enhance the disease resistance of Lilium by modulating the lignin biosynthetic pathway.

In our network, the flavonoid pathway begins with naringenin chalcone to produce naringenin, which is then converted into apigenin, dihydrokaempferol, and eriodictyol (Figure 3). Subsequently, apigenin and dihydrokaempferol are converted into luteolin and kaempferol, respectively (Figure 3). The entire pathway was specifically activated in Lilium after infection with B. elliptica, as evident from the up-regulation of several structural genes, including CHS, CHI, F3′H, and FLS (Supplementary Figures 12, 16). Eriodictyol was recently identified as a novel antibacterial compound; however, its role in plant–pathogen interactions has not yet been determined (Ho et al., 2018). In the current study, eriodictyol accumulation showed significant differences between the inoculated and control treatments at all stages of infection. Interestingly, evidence shows that eriodictyol synthase F3′H (DN130296_c0_g1) enhances host plant resistance against invading pathogens (Mizuno et al., 2014; Hutabarat et al., 2016). Further analysis of eriodictyol may thus open new opportunities for the cultivation of resistant Lilium. Despite no significant difference in the accumulation of kaempferol at the late stage of infection, the expression of its biosynthetic gene FLS (DN128107_c0_g1) and conjugated compounds, which can be stored in vacuoles for long periods, increased with time post-inoculation. FLS and F3′H are responsible for the synthesis of kaempferol and quercetin, and were recently shown to be critical for disease resistance in plants (Zhang et al., 2016). Together, these findings suggest the possibility that kaempferol and eriodictyol are the main metabolites of the flavonoid pathway that shape the response of Lilium to B. elliptica infection, and their levels are affected by a series of enzymes such as CHS, F3H, and FLS.

Analysis of the time-course RNA-seq data and WGCNA results also support the central role of phenylpropanoid and flavonoid biosynthesis pathways in the defense response of Lilium. At the early stage of infection, the MYB30 gene was significantly up-regulated in inoculated leaves (Figure 5A). The MYB TFs have been shown to confer resistance to rice plants by regulating the phenylpropanoid pathway (He et al., 2020). Consistently, the accumulation of phenylpropanoid pathway-related metabolites and the expression of MYB30 significantly increased at the early stage of pathogen infection. At the middle stage, CAD and I2′H, which represent structural genes in the lignin and flavonoid pathways, respectively, were significantly up-regulated in the dark-green module (Figure 5B). At all three stages, MYB61 was highly up-regulated, which affected the expression of structural genes involved in the phenylpropanoid and flavonoid pathways, such as PAL, CSE, CHS, CHI, and POD (Figure 6). This finding is consistent with the significant accumulation of caffeic acid, dihydrokaempferol, and kaempferol in the inoculated treatments at the three stage of infection.

Analysis of the transcriptome time series and WGCNA results let us further elucidate important signaling pathways underlying the defense responses at transcriptional level. In B. elliptica-inoculated leaves with no lesions, the defense regulator gene DOF5.5 was significantly up-regulated, indicating its involvement in the plant response to infection, similar to that of MYB30 and AS1 (Figure 5A). The DOF TFs function as upstream regulators of the salicylic acid (SA) signaling pathway, and enhance the resistance to pathogens in plants (Kang et al., 2003, 2016; Yu et al., 2019). In Arabidopsis, changes in MYB30 expression levels modulated the SA content and SA-associated gene expression levels (Raffaele et al., 2006). Interestingly, the WGCNA results showed \ a regulatory relationship between DOF5.5 and AS1. The AS1 gene acts a positive regulator of SA-independent extracellular defenses but a negative regulator of the defense response by selectively binding to the promoters of genes controlled by the immune activator, jasmonic acid (JA) (Nurmberg et al., 2007). Additionally, WRKY70 differed significantly between the inoculated and control treatments at the early stage of infection. The AtWRKY70 is an important node of convergence for the SA- and JA-mediated defense signaling pathways (Li et al., 2017).

With the expansion of necrotic lesions on B. elliptica-inoculated “Sorbonne” leaves (middle stage of infection), genes including PBL19, CRK2, and WAK5, which are all members of the Interleukin-1 Receptor-Associated Kinase (IRAK) gene family, were instrumental in the defense against B. elliptica (Figure 5B). The IRAK family proteins are conserved upstream signaling molecules that can regulate various stress adaptation programs in plants (Srideepthi et al., 2020). Our results demonstrated a strong correlation between the expression of PBL19 and CKX9. OsPBL1 exhibits 67% amino acid sequence identity to a positive regulator of ETI, and the expression of OsPBL1 in transgenic Arabidopsis increased after an exogenous treatment of cytokinin and SA (Lee and Kim, 2015). CKX5, the only known gene involved in cytokinin catabolism, was recently proven to respond to B. cinerea infection in various ways that are differently modulated by JA and ethylene biosynthesis pathways in Arabidopsis (Li et al., 2021).

At the late stage of B. elliptica infection (green module), necrotic lesions on leaves expanded, and JAZ, WRKY30, and WRKY22 genes played central roles in plant defense (Figure 5C). The JAZ genes encode transcriptional repressors of JA-responsive genes and major components of the JA receptor complex (Thatcher et al., 2016). WRKY30 and WRKY22 are known to enhance disease resistance in plants via the SA and JA defense systems (Peng et al., 2012; Han et al., 2013; Jiang et al., 2015; Kloth et al., 2016).

Furthermore, several differentially expressed TF-encoding genes were identified in the dark-gray module, indicating that these TFs may be significantly involved in the defense response at all three stages (Figure 6). The results revealed a major TF involved in SA and JA signaling, WRKY33 (Figure 5C). WRKY33 showed a strong regulatory relationship with MYC2 and NPR1, and was highly expressed in the inoculated leaves. AtWRKY33 plays a major role in the crosstalk between JA and SA signaling pathways and metabolic responses in response to B. cinerea infection (Birkenbihl et al., 2012; Liu et al., 2017). MYC2 positively regulates JA-mediated flavonoid biosynthesis, while NPR1 modulates SA and JA antagonism in plants (Dombrecht et al., 2007; Knoth et al., 2009). In addition, the FLS2, BAK1, BRI1, and MYB61 genes were significantly up-regulated across all infection stages in our results. FLS2 is a phylogenetically related cell surface pattern recognition receptor and a coreceptor for BAK1. BAK1 is a coreceptor for the brassinolide (BR) receptor BRI1 (Tian et al., 2014), which is responsible for initiating the events of BR signal transduction. While the BR signaling pathway contributes to the growth–defense tradeoff by suppressing the expression of defensin and glucosinolate biosynthesis genes (Liu et al., 2018; Liao et al., 2020). MYB61 was also positively regulated by BRI1 in our results, but whether MYB61 is involved in BR signaling remains unknown.

Taken together, our analysis indicates that SA and JA signaling pathways play pivotal roles in B. elliptica– Sorbonne” interaction. This supports the proposal that a large number of transcripts related to B. elliptica resistance in L. regale were involved in the JA and phenylpropanoid pathways (Cui et al., 2018a). Moreover, as reported in the interactions between other hosts and Botrytis spp. (e.g., Arabidopsis–B. cinerea and Arabidopsis–Alternaria brassicicola), the SA and JA signaling pathways play a crucial role in the response of “Sorbonne” to fungal infection (Zheng et al., 2006; Ederli et al., 2015, 2021; Liao et al., 2020). Besides, the BR-related genes were significantly up-regulated during across all infection stages in our results. The crosstalk between BR and JA signaling affects the growth-defense tradeoff in Arabidopsis–B. cinerea interaction (Liu et al., 2018; Liao et al., 2020), while exogenous BR before the inoculation of B. cinerea enhances the defense response in rose petals (Liu et al., 2018). The role of BR pathway in Lilium defense response to the gray mold remains to be further elucidated.

In addition to phenolics and hormone signaling pathways, our results showed that other signal transduction pathways play a prominent role in the defense response. At the early stage of stage, LURP was significantly up-regulated in the plant response to infection, and regulated the transcript levels of POD and HSP70 (Figure 5A). POD and HSP70 are important for plant resistance (Baig, 2018). AtLURP1 shows an unusually pronounced transcriptional up-regulation in response to infection (Knoth and Eulgem, 2008). Together, these data suggest that LURP plays a pivotal role in the response of “Sorbonne” to B. elliptica infection at the early stage.

At the middle stage, the PBL19, CRK2, and WAK5 genes were instrumental in the defense response of “Sorbonne” against B. elliptica (Figure 5B). The pbl13-2 knockout mutant shows higher level of reactive oxygen species (ROS) and greater flagellin-induced activation of mitogen-activated protein kinases (MAPKs) than the wild type (Lin et al., 2015). Moreover, overexpression of CRKs enhances PTI in transgenic Arabidopsis (Yeh et al., 2015). Wall-associated kinases (WAKs) localize to the cell wall and participate in pathogen recognition and signal transduction in plants (Kurt et al., 2020). Collectively, the transduction of various signals contributed to the defense response of Lilium against B. elliptica at the middle stage of infection.

At the late stage of infection, CML27, DELLA, and JAZ played central roles in plant defense. CML27 encodes a Ca²⁺-binding protein involved in the Ca²⁺ signaling pathway. Ca²⁺ regulates diverse cellular processes and functions as a secondary messenger, enabling plants to sense and quickly respond to extracellular stimuli (Edel et al., 2017). AtRPS2 and AtRPM1 activate Ca²⁺-dependent protein kinases (CPKs) for mediating bifurcate immune responses (Gao et al., 2013). In our results, RPS2 and RPM1 were significantly up-regulated among the hub genes, indicating the important role of the Ca²⁺ signaling pathway in the defense response of “Sorbonne”.

In all three stages (Figure 6), MYB61 was highly expressed, implying that it is central to the defense response. B. elliptica infects Lilium via the stomata, and MYB61 regulates stomatal aperture (Hsieh et al., 2001; Gao et al., 2018; Romero-Romero et al., 2018), which may explain its role in the defense response. In the current study, the expression levels of PR10 and PR4 were strongly correlated with that of Pti5 (Figure 6). Overexpression of the tomato TF genes, Pti4/5/6, in Arabidopsis showed that Pti4/5/6 activate the expression of a wide array of PR genes in vivo, resulting in enhanced defense against certain fungal pathogens (Gu et al., 2002; Gonzalez-Lamothe et al., 2008). The expression levels of RIN4, RPM1, RPS2, NPR1, and TGA were significantly up-regulated at all infection stages. In Arabidopsis, the phosphorylation and cleavage of RIN4 activates RPM1 and RPS2, which encode R proteins involved in ETI (Li et al., 2019; Ray et al., 2019; Yang et al., 2021; Zhao et al., 2021). NPR1-mediated DNA binding of TGA2 is critical for the activation of defense related genes. Both NPR1 and TGA1 act as master redox-sensitive transcriptional regulators of PR genes in plants (Fu and Dong, 2013). Altogether, our results suggest the transduction of diverse signals jointly shape the defense response of Lilium hybrid “Sorbonne” to B. elliptica.