Unraveling verticillium wilt resistance: insight from the integration of transcriptome and metabolome in wild eggplant

Verticillium wilt, caused by Verticillium dahliae, is a soil-borne disease affecting eggplant. Wild eggplant, recognized as an excellent disease-resistant resource against verticillium wilt, plays a pivotal role in grafting and breeding for disease resistance. However, the underlying resistance mechanisms of wild eggplant remain poorly understood. This study compared two wild eggplant varieties, LC-2 (high resistance) and LC-7 (sensitive) at the phenotypic, transcriptomic, and metabolomic levels to determine the molecular basis of their resistance to verticillium wilt. These two varieties exhibit substantial phenotypic differences in petal color, leaf spines, and fruit traits. Following inoculation with V. dahliae, LC-2 demonstrated significantly higher activities of polyphenol oxidase, superoxide dismutase, peroxidase, phenylalanine ammonia lyase, β-1,3 glucanase, and chitinase than did LC-7. RNA sequencing revealed 4,017 differentially expressed genes (DEGs), with a significant portion implicated in processes associated with disease resistance and growth. These processes encompassed defense responses, cell wall biogenesis, developmental processes, and biosynthesis of spermidine, cinnamic acid, and cutin. A gene co-expression analysis identified 13 transcription factors as hub genes in modules related to plant defense response. Some genes exhibited distinct expression patterns between LC-2 and LC-7, suggesting their crucial roles in responding to infection. Further, metabolome analysis identified 549 differentially accumulated metabolites (DAMs) between LC-2 and LC-7, primarily consisting of compounds such as flavonoids, phenolic acids, lipids, and other metabolites. Integrated transcriptome and metabolome analyses revealed the association of 35 gene–metabolite pairs in modules related to the plant defense response, highlighting the interconnected processes underlying the plant defense response. These findings characterize the molecular basis of LC-2 resistance to verticillium wilt and thus have potential value for future breeding of wilt-resistant eggplant varieties.


Introduction
Verticillium dahliae (V.dahliae), a worldwide soil-borne fungus causing verticillium wilt, is one of the major diseases that harm agricultural production (Klosterman et al., 2009).Studies have shown that the pathogenic fungi overwinter through microsclerotia that are stress resistant, surviving in the soil for long periods and thus playing an important role in the disease cycle (Tzima et al., 2010).The hyphae enter host plants by formation of an infection structure that develops into conidia and expands upward from the roots during transpiration, entering the stems and leaves of plants and resulting in leaf curl, necrosis, defoliation, and discoloration, ultimately leading to plant withering and death (Sink and Grey, 1999).Verticillium wilt can occur throughout the entire growth period of eggplant, usually occurring during the fruiting period.In China 20%-40% of the eggplant fields are affected by V. dahliae, causing large-scale reductions in yield (Fradin and Thomma, 2006).V. dahliae contains multiple transposons and repetitive genes as well as genes involved in metabolic regulation and signal transmission (Zhang et al., 2020;Feng et al., 2023).This demonstration of genetic variation suggested that the mechanisms of pathogenicity may be diverse, and thus further analysis of disease resistance mechanisms is needed.
Researchers have isolated and cloned many genes related to the resistance to verticillium wilt from different plant species (Song et al., 2020).For example, the transient silencing of GhNDR1 and GhMKK2 by Agrobacterium-mediated VIGS in cotton (Gossypium hirsutum) led to increased sensitivity of cotton to V. dahliae, while silencing of GhNPR1 had no significant effect on antagonism toward the pathogen (Gao et al., 2011).In addition, GbNRX1 can regulate cotton resistance by clearing extracellular ROS (Li et al., 2016).VdMKK1 is significantly upregulated during the invasion of host cotton by V. dahliae, and it plays a crucial role in the invasion process by regulating the expression of cell wall synthesis genes and thereby maintaining the integrity of cell walls (Li et al., 2023).GbHyPRP1 negatively regulates the resistance of cotton plants to V. dahliae by regulating cell wall structures and ROS levels (Yang et al., 2018).In addition to defense-related proteins, further studies have been conducted on the role of transcription factors in Verticillium resistance (Song et al., 2020).A full-length cDNA library and expressed sequence tag (EST) sequencing in cotton inoculated with V. dahliae screened two ethylene-related genes, GbERF1 and GbERF2; overexpression of GbERF2 induced the expression of related genes and PR proteins by regulating the ethylene pathway and enhancing the resistance of tobacco (Zuo et al., 2007).GhWRKY1 enhanced resistance to verticillium wilt by increasing lignification (Hu et al., 2021).VdMRTF1 can mediate cotton resistance to V. dahliae by regulating various processes, including melanin biosynthesis, development of microsclerotia, and resistance to elevated Ca 2+ levels (Lai et al., 2022).An NAC transcription factor, GhATAF1, and GhMYB108 were both induced by V. dahliae infection and promoted defense responses (Cheng et al., 2016;He et al., 2016).Some oxidoreductases and hydrolases, including peroxidase and cell wall proteolytic enzymes, play important roles in the resistance response to verticillium wilt.Several studies of cotton resistance to verticillium wilt have isolated these protease genes (Zuo et al., 2005;Xu et al., 2011a), and related studies have shown that after inoculation with pathogens, the peroxidase activity of resistant varieties is significantly higher than that of susceptible varieties (Gayoso et al., 2010;Xu et al., 2011a).In addition, hydrolytic enzymes such as chitinase and glucanase are involved in the resistance to verticillium wilt by degrading fungal cell walls (Dubery and Slater, 1997).
The mechanism of resistance to verticillium wilt has been extensively studied in the Solanaceae (Kawchuk et al., 2001;Kruijt et al., 2005).The Ve genes, including Ve1 and Ve2, were early discovered genes for resistance to verticillium wilt.The genes encode leucine-rich repetitive (LRR) proteins (Kruijt et al., 2005).When these two Ve genes were expressed separately in susceptible potato varieties, both produced resistance to isolated strains of V. dahliae (Kawchuk et al., 2001).Similar resistance mechanisms have also been found in homologous gene phenotypes of Ve in other Solanaceae crops, including tomato (Chai et al., 2003) and eggplant (Liu et al., 2012).Jue et al. (2014) cloned two genes conferring strong resistance to verticillium wilt, StoNPR1 and StoWRKY6, from wild eggplant Torubam.The genes enhanced the resistance of eggplant to verticillium wilt by synthesizing proteins that mediated the salicylic acid signaling pathway.The interaction between miR482e and nucleotide binding sites (NBS) and leucine-rich repeat (LRR) motifs in potatoes led to increased severity of V. dahliae infection (Yang et al., 2015).
Although various genes have been identified in the resistance system to verticillium wilt, little is known about the complex molecular mechanisms underlying defense responses in eggplant.In recent years, the rise of genome, transcriptome, metabolome, and proteome analyses has resulted in an in-depth understanding of the resistance mechanisms toward verticillium wilt (Zhang et al., 2020;Feng et al., 2023).Xu et al. (2011a) identified 3,442 defense-responsive genes from the transcriptomic profiles of V. dahliae-infected cotton using RNA-seq technology.Further research on these differentially expressed genes revealed the key role of lignin metabolism in the resistance of cotton to verticillium wilt (Xu et al., 2011b).A comprehensive analysis utilizing transcriptome and metabolome revealed the intrinsic components and DEGs in Arabidopsis inoculated with V. dahliae (Su et al., 2018).In eggplant, Tomah et al. ( 2023) identified genes encoding enzymes involved in microsclerotia degradation based on the analysis of transcriptomics.Differentially expressed genes between the V. dahliae infected and uninfected samples identified 10 different gene network modules and 22 hub genes with potential roles in regulating cotton defense against V. dahliae infection (Xiong et al., 2021).
To better understand the molecular networks of plantpathogen interactions in eggplant, we investigated two wild eggplant varieties with different resistance phenotypes, LC-2 (high resistance) and LC-7 (sensitive) and analyzed the differences in enzyme activity associated with resistance.We identified transcriptome changes occurring during the process of V. dahliae infection in order to characterize the biological processes involved based on Gene Ontology (GO) terms of the differentially expressed genes using RNA-seq technology.A co-expression analysis identified key hub genes involved in regulating the resistance pathways of eggplant to verticillium wilt.Additionally, a metabolomic analysis identified key differential metabolites.A correlation analysis between transcriptome and metabolome was used to explore the relationship between differentially expressed genes and metabolic components, thereby providing a theoretical basis for the in-depth study of the mechanism of action of eggplant against verticillium wilt.

Plant material and Verticillium dahliae infection treatment
Eggplant (Solanum melongena L.) local varieties LC-2 and LC-7 were collected from Langcang County, Pu'er City, Yunnan Province.V. dahliae was provided by the Institute of Vegetables and Flowers, Chinese Academy of Agricultural Sciences.Seeds of each variety were soaked in 600 mg•L -1 gibberellic acid for germination and sown in a plug containing a substrate (perlite: peat soil = 1:2 V/V).The plug was placed in a constant temperature incubator with a temperature of 30°C during the day and 22°C at night, with a light/dark cycle of 14 h/10 h.After 5-6 true leaves had grown, the V. dahliae was artificially inoculated at a concentration of 1×10 8 cfu•mL -1 using the root-injury-root immersion and root irrigation methods.The fibrous roots of the seedlings were cut off at 0.2 cm, immersed in the V. dahliae solution for 20 min, and planted in plastic pots.Each pot was transferred back to the constant temperature incubator under the same conditions as before and irrigated with 20 mL V. dahliae solution every two days.Uninoculated individuals were used as controls, the leaves of three plants were randomly selected and sampled, each plant was taken as a biological replicate.In the treatment group, two time points were selected at 20 days and 40 days after inoculation, three biological replicates were set at each inoculation time point for LC-2 and LC-7, with one biological replicate investigating 20 plants.

Identification of the disease resistance level
The classification of disease level was as follows: Level 0: asymptomatic; Level 1: less than 25% of the leaves turning yellow; Level 2: 25%-50% of the leaves turning yellow; Level 3: more than 50% of the leaves turning yellow with some diseased leaves fall off; Level 4: the leaves of the whole plant are yellow, and even the whole plant dies.

Measurement of enzyme activities
The third and fourth leaves from the top were selected as experimental materials.Experiments were conducted with three biological replicates by spectrophotometry.The kits were purchased from Beijing Boxbio Science & Technology Co., Ltd., and the determination of each index was performed according to the manufacturer's instructions.A total of six physiological indicators were measured, including polyphenol oxidase (PPO, kit code: AKAO004C), superoxide dismutase (SOD, kit code: AKAO001C), peroxidase (POD, kit code: AKAO005C), phenylalanine ammonia lyase (PAL, kit code: AKAM012U), b-1,3 glucanase (GUN, kit code: AKSU038C), and chitinase (CHT, kit code: AKSU045C).The activities were normalized by the weight of leaf.The change of light absorption value of A410nm of the reaction system by 0.005 per minute per g tissue was defined as a unit of enzyme activity of PPO.When the inhibition rate of xanthine oxidase coupling reaction system in each g tissue was 50%, the activity of SOD in the reaction system was defined as an enzyme activity unit at 560 nm.The change of A470nm per minute by 0.01 per g tissue per mL system was one unit of enzyme activity of POD.One enzyme activity of PAL was defined as a 0.1 change in absorbance value at 290 nm per minute per g of tissue per mL of the reaction system unit.The production of 1 mg of reducing sugar per g of tissue sample per hour was defined as a unit of enzyme activity of GUN at 540 nm.The production of 1 mg of N-acetylglucosamine per g of tissue per hour was defined as a unit of enzyme activity of CHT at 585 nm.

Metabolomic analysis
The 18 LC-2 and LC-7 leaf samples (two varieties × three biological replicates × three time points) of treatment and control groups at 0, 20, and 40 days were freeze-dried and crushed into powder using a mixer mill (MM 400, Retsch).Three biological replicates were set for each sample.Metabolite profiling was performed using a widely targeted metabolome method by Wuhan Metware Biotechnology Co., Ltd.(Wuhan, China) (http:// www.metware.cn/).The samples were obtained from a total of 100 mg of powder extracted at 4°C overnight with 1.0 ml of 70% aqueous methanol and then centrifuged at 10,000 × g for 10 min.The extracts were detected by UPLC-ESI-MS/MS (UPLC: SHIMADZU Nexera X2; MS: Applied Biosystems 4500 QTRAP).The following method was employed for metabolite quantification and analysis.Scans using both Linear ion trap (LIT) and triple quadrupole (QQQ) were acquired on a triple quadrupole-linear ion trap mass spectrometer (API 4500 Q TRAP LC/MS/MS System) equipped with an ESI Turbo Ion-Spray interface.Multiple-reaction monitoring (MRM) was utilized in conjunction with a selfcompiled database (Metware database, MWDB) to perform metabolite quantification.Differences between LC-2 and LC-7 samples at the same time point were analyzed using orthogonal partial least-squares discrimination analysis (OPLS-DA) and variable importance in projection (VIP).Identification of differential metabolites was performed within and between the local varieties.Within each local variety, samples from the 0-day of each variety (representing the initial state) were used as controls to identify differentially accumulated metabolites (DAMs) at 20 and 40 days.Between local varieties, LC-7 was used as a control to identify DAMs between LC-2 and LC-7 at the same treatment time point.Metabolites with significant differences in content were set with thresholds of variable importance in projection (VIP) ≥ 1 and fold change ≥ 2 or ≤ 0.5 (Guo et al., 2022).The Kyoto Encyclopedia of Genes and Genomes (KEGG) compound database (http://www.kegg.jp/kegg/compound/)was utilized for metabolite annotation.

Transcriptome sequencing and data analysis
Transcriptome sequencing was performed by using the leaf samples as for the metabolomic analysis, and three biological replicates in each sample were set.Total RNA was extracted from the samples using Trizol reagent (Invitrogen, USA), and RNA quality was assessed through 1% agarose gel electrophoresis.A NanoDrop 2000 (NanoDrop Technologies, Wilmington, DE, United States) and an Agilent 2100 Bioanalyzer (Agilent Technologies, USA) were used to determine the RNA concentration and purity, respectively.Highquality RNAs were then used for RNA sequencing by Wuhan Metware Biotechnology Co., Ltd.(Wuhan, China) (http:// www.metware.cn/)using the Illumina HiSeq platform (Illumina Inc., San Diego, CA, USA).After a standard library construction process (Li et al., 2022a), clean reads were aligned to the eggplant reference genome (http://eggplant-hq.cn/Eggplant/home/index) using the Hisat2 program.The mapped reads of each sample were assembled by StringTie, and featureCounts v1.5.0-p3 was used to count the number of reads mapped to each gene to calculate FPKM.Differentially expressed genes (DEGs) were identified within and between local varieties by using the same comparison strategy as used in the DAMs identification.DESeq2 software was used to identify the differentially expressed genes using the criteria absolute value of the log2 fold change ≥ 1 and false discovery rate (FDR) < 0.05.The Gene ontology enrichment for DEGs was visualized using WEGO (https:// wego.genomics.cn/)and REVIGO (http://revigo.irb.hr/).The Kyoto Encyclopedia of Genes and Genomes (KEGG) database (https:// www.genome.jp/kegg/) was used to analyze the pathways of DEGs.

Gene co-expression network analysis
A gene co-expression network was constructed by using the WGCNA R package.The FPKM values of DEGs were normalized by using the limma R package.According to the results of multiple soft thresholding powers calculation (Supplementary Figure S1), the network was constructed by setting a power of 18, a minModuleSize of 30, and a mergeCutHeight of 0.25.Module eigengenes (MEs) were defined as the first principal component of the expression matrix of the corresponding module and were calculated by using the moduleEigengenes function with default parameters for each module.After network construction, the top 30% of weighted edges were analyzed, as millions of edges were generated.The module hub genes were identified as the nodes with the top 10% degree in each module.

Correlation analysis of metabolome and transcriptome data
The potential associations between the transcriptome and metabolome data were identified by assessing the correlations between metabolite content and the expression pattern of each gene and the ME of each gene co-expression module.Briefly, Pearson's correlation coefficient (PCC) for the content of DAMs and the expression of DEGs was calculated using the "cor" function in R, and significant correlations were screened using a criterion of an absolute value greater than 0.80 and a p-value less than 0.05.The PCC was also calculated between metabolite content and ME using the "moduleEigengenes" function in the WGCNA R package.Considering the capacity of ME to partially represent the average expression of module genes, it is important to note that the average per se may mask some of the correlations between genes and metabolites.Therefore, a lower criterion than DAMs-DEGs correlation was used, with absolute coefficients greater than 0.6 and p-values less than 0.05 considered as indicating a significant correlation between metabolites and gene modules.Finally, metabolites that demonstrated significant correlations with both the modules and the genes within the modules were selectively filtered, thereby establishing a comprehensive association between the metabolome and transcriptome.

Quantitative real-time PCR (qRT-PCR) validation
To validate gene expression patterns discovered by transcriptome sequencing, 13 transcription factors were randomly selected, and their expression patterns were validated by qRT-PCR.Total RNA of each LC-2 and LC-7 sample was reverse transcribed into cDNA using a PrimeScriptTM RT reagent kit (Takara, Dalian, China).The qRT-PCR reaction was performed using TB GreenTM Premix Ex TaqTM II (Takara, Dalian, China) on a Roche LightCycler 96 system (Roche) following the manufacturers' instructions.Relative quantitative analysis of data was performed by the 2 −DDCT method, with GAPDH as the reference gene (Zhao et al., 2013;Zhou et al., 2014).The primers used for RT-qPCR are listed in Supplementary Table S1.Three independent biological samples of eggplant plants of LC-2 and LC-7 under the same treatment were used, which were different from the transcriptome samples.

Results
LC-2 exhibits greater resistance and physiology shifts in response to V. dahliae than LC-7 LC-2 and LC-7 were grown in a facility greenhouse to ensure consistent environmental conditions.Agronomic traits were investigated when the respective first panicle fruits were ripe.
LC-22 and LC-7 showed significant differences in main stem color, main stem prickle, leaf spines, corolla color, single fruit weight, fruit shape, and fruit surface furrows (Figure 1; Supplementary Table S2).
To clarify the differences in resistance between LC-2 and LC-7, the resistance level was identified at 20 and 40 days after V. dahliae inoculation.The results showed that the overall resistance of LC-2 was stronger than that of LC-7.Specifically, 20 days after inoculation, LC-2 showed high resistance and LC-7 showed moderate resistance to V. dahliae, and 40 days after inoculation, LC-2 showed susceptibility, while LC-7 showed high susceptibility (Figure 2).
To analyze the characteristics of changes in the enzymes related to disease resistance in LC-2 and LC-7, the six indices PPO, SOD, POD, PAL, GUN, and CHT were measured in leaves at 20 and 40 days after V. dahliae inoculation.The results showed that the six activities demonstrated gradually increasing trends with time after inoculation with V. dahliae.The values of the six indexes in LC-2 were higher than those in LC-7 under normal growth conditions.All six indices were positively correlated with resistance to V. dahliae (Table 1).
The DAMs of LC-2 and LC-7 in response to V. dahliae infection To gain a better understanding of the differences in metabolites between LC-2 and LC-7 in response to V. dahliae infection, metabolic profiling was conducted using samples from the treatment and control groups of LC-2 and LC-7 at 0, 20, and 40 days.A total of 840 metabolites were detected; these included flavonoids (175), phenolic acids (142), lipids (130), alkaloids (85), amino acids and derivatives (60), organic acids (58), nucleotides and derivatives (38), terpenoids (31), lignans and coumarins (27), steroids (16), quinones (3), and other metabolites (75).A principal component analysis (PCA) revealed that the three biological replicates of each sample tended to group together.The 0 d and 20 d samples had a closer distance between cultivars in the 2D PCA analysis, while the 40 d sample had a greater distance between LC-2 and LC-7 (Figure 3A).
Within each variety, differential accumulation of metabolites was compared between 20-day and 0-day samples as well as between 40day and 0-day time points.Totals of 593 and 600 DAMs were identified in LC-2 and LC-7, respectively (Supplementary Table S3).Based on the differential expression patterns at the two time points, we classified the expression patterns of the DAMs into eight categories: 20up_40NA, 20down_40NA, 20NA_40up, 20NA_40down, 20up_40up, 20down_40up, 20up_40down, and 20down_40down ("Up," "down," and "NA" represent upregulation, downregulation, and no differential expression at that time point, respectively).The results showed that LC-2 had the most 20up_40uptype metabolites, followed by 20up_40NA-type metabolites (Figure 3B), indicating that many metabolites in LC-2 responded to V. dahliae infection and increased in content at 20 days.In contrast, LC-7 had the largest number of differential metabolites in the 20NA_40 up category (Figure 3B), suggesting that metabolites in LC-7 responded to V. dahliae infection and began to accumulate at 40 days.A comparison of LC-2 and LC-7 identified 549 DAMs between local varieties.A total of 346 DAMs (176 upregulated and 170 downregulated) were identified at the 20-day time point, with the most common metabolites being flavonoids (96), phenolic acids (63), and lipids (49) (Supplementary Figure S2A).These DAMs were significantly enriched in pathways such as flavone and flavonol biosynthesis, lysine degradation, and tryptophan metabolism (Figure 3C).At the 40-day time point, there were 399 DAMs (180 upregulated and 219 downregulated) between LC-2 and LC-7, with the most common being phenolic acids (85), flavonoids (77), and lipids (65) (Supplementary Figure S2B).These DAMs were significantly enriched in pathways such as biosynthesis of Performance of typical agronomic traits of LC-2 and LC-7 in the field, including petal color, leaf spines, and fruit traits.Li et al. 10.3389/fpls.2024.1378748Frontiers in Plant Science frontiersin.orgsecondary metabolites, flavone and flavonol biosynthesis, and cysteine and methionine metabolism (Figure 3D).
Transcriptomic signatures of LC-2 and LC-7 in response to V. dahliae infection After raw data filtering and examination of the sequencing error rate and the distribution of GC content, an average of 6.74 Gb of clean data was generated for each library.The average error rate was 0.03; the average Q30 was 92.59%, and the GC content ranged from 42.13% to 46.54% (Supplementary Table S4).The difference in genetic background led to distinct average mapping rates for LC-2 and LC-7: LC-2 achieved 81.28%, while LC-7 attained 70.26%, with fluctuations within 3% for each library.The lower mapping rate in LC-7 is not attributed to contamination by V. dahliae, as samples without infection also displayed mapping rates similar to other samples.For the sequences mapped to the reference genome,  85.42%-89.78% of the clean reads were mapped to exon regions, and approximately 3.90%-6.31% of the clean reads were mapped to intron regions (Supplementary Table S5).Finally, 36568 transcripts were identified.
The temporal expression pattern of each gene was determined by comparing 20-day vs. 0-day and 40-day vs. 0-day levels within each variety.Using the criteria for identifying DEGs, 6,878 and 7,814 DEGs were identified in LC-2 and LC-7, respectively.The DEGs in LC-2 were predominantly distributed among the 20down_40NA, 20up_40NA, 20NA_40up, and 20up_40up categories, with a similar number of DEGs in each category (Figure 4A).In contrast, the DEGs in LC-7 were largely distributed in the 20NA_40up and 20NA_40down categories (Figure 4A).These findings suggest that the response of LC-7 to V. dahliae infection occurred later than that of LC-2.
The comparative analysis of biological processes associated with DEGs revealed significant differences between LC-2 and LC-7.The percentage of genes enriched in developmental processes, cell wall organization or biogenesis, and regulation of biosynthetic processes was substantially higher in LC-2 than in LC-7.Conversely, the percentage of genes enriched in processes related to response to fungus was similar between the two varieties (Figure 4B).
The comparison of DEGs between LC-2 and LC-7 revealed that LC-2 and LC-7 had 4,017 and 4,953 unique DEGs, respectively (Figure 4C).These genes exhibited variety-specific differential expression patterns (Figures 4D, E).Among these, a total of 2,861 DEGs were common to both LC-2 and LC-7 (Figure 4C).Further analysis of their expression patterns at the same time points between varieties (LC-2 vs. LC-7) revealed that 351 genes did not meet the criteria for differential expression between varieties and exhibited similar expression patterns over time (Figure 4F).The remaining 2,510 genes showed completely different expression patterns between LC-2 and LC-7 (Figure 4G).

Diverse types of DEGs were involved in various biological processes
The DEGs of each type were subjected to Gene Ontology (GO) enrichment analysis.The analysis showed that among the 4017 LC-2-specific DEGs, the genes were involved in defense response, biological synthesis processes (such as cell wall biogenesis, spermidine biosynthetic process, cinnamic acid biosynthetic process, and regulation of cutin biosynthetic process), and developmental processes (Supplementary Figure S3A).In contrast, the 4953 LC-7-specific DEGs were primarily involved in defense response to other organisms, plant hormone biosynthetic process, and metabolic processes (such as abscisic acid and salicylic acid) (Supplementary Figure S3B).The 351 DEGs with similar expression patterns between varieties were involved in defense response, plant hormone response processes (such as salicylic acid and abscisic acid), and compound biosynthesis processes (Supplementary Figure S3C).Moreover, the 2510 genes that were shared between varieties but with different expression patterns were involved in defense response and antibiotic metabolic and response processes as well as other processes such as sterol biosynthetic process and cell wall organization (Supplementary Figure S3D).Moreover, a total of 408 DEGs were involved in the plant-pathogen interaction pathway (ko04626, Figure 5A) and the disease resistance-related pathway in plant hormone signal transduction (ko04075, Figure 5B).Most of these genes in these two pathways showed different expression patterns between LC-2 and LC-7 (Figures 5A, B).Supplementary Table S6 provides detailed information on these genes.

Gene co-expression modules related to defense response against V. dahliae infection
The co-expression network analysis resulted in the identification of 29 modules that encompassed a total of 6,258 nodes.To confirm the genes associated with defense response, we analyzed the distribution of genes related to the plant-pathogen interaction pathway (ko04626) and the disease resistance-related pathway in plant hormone signal transduction (ko04075) across different modules.The findings revealed that 72% of the disease resistance-related genes were recruited in the turquoise, green, tan, blue, and yellow modules, which collectively accounted for 49.1% of the total nodes in the network.Thus, these five modules were considered to be associated with the defense response against V. dahliae infection.A closer examination of the network's topological structure revealed that the blue, turquoise, tan, and green modules exhibited a transitional connection structure, and the blue module was linked to the yellow module through several co-expression genes (Figure 6A).The functional analysis demonstrated that the genes in the blue module were primarily involved in biological processes such as growth and development, epigenetic regulation, and compound biosynthetic processes (Figure 6B).In contrast, the turquoise, tan, green, and yellow modules were enriched with genes related to defense response, hormone signaling pathway, and stress response processes (Figures 6C-F).The green module was notably enriched with a large number of genes related to cell wall biogenesis and thickening (Figure 6F).
Module hub genes were identified by selecting nodes with the top 10% degree in each module.A total of 332 hub genes were screened in the five key modules (Supplementary Table S7), including 13 transcription factors (TFs) from the AP2, Bhlh, GRAS, NAC, NF-YA, RWP-RK, HB-HD-ZIP, GARP-G2-like, and other TF families (Table 2).KEGG pathway annotation revealed that some of these TFs were involved in plant-pathogen interactions, plant hormone transduction, and other pathways related to pathogen resistance.Additionally, there were some unannotated transcription factors whose functions and roles in The expression and functional characteristics of V. dahliae-responsive genes in LC-2 and LC-7.biological stress required further research and verification.These 13 key transcription factors were co-expressed with genes of hormone response pathways, defense response pathways, and defense-related metabolite synthesis.Furthermore, their neighboring genes in the network were co-expressed with cell wall thickening genes (Figure 7A).The analysis of the expression patterns of these key TFs revealed that some had different response patterns during V. dahliae infection between LC-2 and LC-7.For example, the transcription factors Smechr0101299 and Smechr0802554 of the HB-HD-ZIP family, and Smechr0303661 of the GRAS family were induced in expression in the LC-2 after 20 days of V. dahliae infection treatment, while they were continuously downregulated in expression in LC-7 (Figure 7B).The differences in these key transcription factors between LC-2 and LC-7 may be one of the important reasons for the difference in disease resistance.

Correlations of module genes and metabolites
The correlation analysis revealed that 236065 DEG-DAM pairs were significantly correlated, whereas only 460 metabolites were significantly correlated with the ME of the five key modules.
By combining these two results, we identified 261 metabolites significantly associated with 1190 genes in key modules, as well as with the MEs of the key modules (Figure 8A).We further screened the metabolites and related genes to check whether they participated in the same biological pathways.This led to the identification of 35 pairs of gene-metabolite combinations in the key modules, where metabolites and related genes were involved in the same processes related to defense response.These processes included flavonoid biosynthesis, tropane, piperidine and pyridine alkaloid biosynthesis, biosynthesis of various antibiotics, cutin, suberine and wax biosynthesis, plant hormone signal transduction, neomycin, kanamycin, and gentamicin biosynthesis (Figure 8B).Additional information regarding these 35 genemetabolite combinations in the key modules is provided in Supplementary Table S8.
Confirmation of the transcriptome data using RT−qPCR Expression levels of 13 genes were validated using RT-qPCR assays.The results of RT-qPCR assays were consistent with the transcriptome analysis, and demonstrated the accuracy of RNAseq in this study.Notably, three of the 13 transcription factors (Smechr0303661 Smechr0101299, and Smechr0802554) were induced in LC-2 and suppressed in LC-7 (Figure 9).These results indicated that these three genes play potentially important roles in disease resistance to V. dahliae.

Discussion
In recent years, eggplant production in China has become increasingly industrialized, and continuous cropping has led to increasingly severe eggplant wilt disease, causing serious losses in eggplant production (Pisuttu et al., 2020).The research on the resistance mechanism of verticillium wilt has made breakthroughs in crops such as cotton (Gao et al., 2011), tomato (Chai et al., 2003;Fu et al., 2016), and potato (Yan et al., 2022), but little is known about the mechanism of resistance to verticillium wilt in eggplant.
The rapid development of molecular biology technology has provided a foundation for the in-depth exploration of the mechanism of resistance to wilt disease in eggplant.Wild eggplant resources show excellent resistance to verticillium wilt.There are abundant species of wild eggplant in Yunnan, China.Our research group had collected more than 40 kinds of wild eggplant in Yunnan Province.The genetic relationship of these resources and their resistance to vegetable soil-borne diseases were identified (unpublished).Two eggplant species LC-2 and LC-7 from Lancang County, Puer City, Yunnan Province showed different resistance to verticillium wilt and closer genetic relationship compared with other eggplant resources.Therefore, LC-2 and LC-7 were selected as experimental materials for the combined transcriptome and metabolome analysis.
Plants activate defense mechanisms to protect themselves when attacked by pathogens, responding quickly with either direct or indirect reactions.The degree of curling of plant leaf epidermis    Li et al. 10.3389/fpls.2024.1378748Frontiers in Plant Science frontiersin.org(Singkaravanit-Ogawa et al., 2021), thickness of stratum corneum (Kempel et al., 2011) and structure of cell wall (Bacete et al., 2018), can serve as physical barriers for defense.Plants may attempt to limit the invasion and spread of pathogens by curling.In our study, a significant difference in the degree of leaf curling between LC-2 and LC-7 was observed (Figure 1), indicating potential differences in their defense mechanisms against diseases.The defense of eggplant not only rely on physical defense but also on innate defense mechanisms like antioxidative systems.When plants are infected with pathogens, they produce reactive oxygen species.
Previous studies have shown that the accumulation of reactive oxygen species is closely related to plant disease resistance (Gayoso et al., 2010;Xu et al., 2011a).However, high levels of reactive oxygen species can have adverse effects on plant growth, and thus the expression levels of various antioxidant enzymes such as SOD, POD, and CAT will also change accordingly to protect plants from oxidative damage (Asadi et al., 2017;Huang et al., 2019).The enhanced activities of antioxidant enzymes can effectively scavenge harmful reactive oxygen species to maintain the stability of plant defense systems (Mittler et al., 2022).In tomato seedlings, the enzyme activities of the defense-related antioxidants CAT, SOD, POD, and PAL significantly increased upon inoculation with V. dahliae compared to controls (Pei et al., 2022).The activities of POD and PAL in cotton increased due to the infection of the V991 strain (Xu et al., 2011a).There were positive correlations between POD, PAL, and PPO activities and resistance to verticillium wilt in eggplant caused by V. dahliae (Zhou, 2012).In this study, we found that the six indexes of PPO, SOD, POD, PAL, GUN, and CHT in LC-2 and LC-7 showed gradually increasing trends with time after inoculation with V. dahliae.Furthermore, the values of the six indexes in LC-2 were higher than those in LC-7 under normal growth conditions and inoculation with V. dahliae (Table 1).This was similar to a previous research report (Zhou, 2012).Cultivated crop wild relatives provide valuable genetic resources for crop improvement (Bohra et al., 2022).LC-2, as a wild relative of cultivated eggplant, may possess verticillium wilt resistance as a result of the combined action of multiple resistance-associated genes, given that plant disease resistance itself is a complex phenotype.Identifying genes that potentially play a role in its resistance is one of the key objectives of this study.Genetic background differences can manifest as differences in gene content, including variations in copy number and sequence.However, regardless of the form, the action of disease-resistant genes must be realized through expression in response to V. dahliae infection.
From a transcriptomic standpoint, LC-2 and LC-7 exhibit notable differences in gene expression in response to V. dahliae infection.While both strains possess a similar number of genes associated with fungal/pathogen defense mechanisms, LC-2 demonstrates a significant portion of DEGs responding to V. dahliae infection at an early stage (20d), whereas LC-7 exhibits a higher concentration of DEGs in the later stages of infection (40d) (Figure 4A).This suggests that LC-2 initiates a quicker response to V. dahliae infection, establishing the groundwork for its heightened resistance.Additionally, LC-2 displays a greater enrichment of genes involved in tissue development compared to LC-7 (Figure 4B), indicating a lesser impact of V. dahliae infection on the growth and development processes of LC-2.Moreover, genes associated with cell wall biogenesis and organization, serving as a dynamic barrier against pathogen invasion (Underwood, 2012), are more abundant in LC-2 than in LC-7.These findings indirectly support the superior verticillium wilt resistance observed in LC-2.
Numerous studies have indicated that co-regulated genes are often functionally related, and simultaneously recruited in coexpression network modules participating in various biological processes (van Noort et al., 2003;Romero-Campero et al., 2013;Li et al., 2022a).Here, we observed that the majority of disease resistance pathway genes are concentrated within 5 gene coexpression network modules, which are interconnected in terms of network topology and exhibit transitional functional patterns (Figure 6).Functional analysis revealed that these modules are involved in plant hormone signal transduction, activation of defense response genes and stress-responsive genes, cell wall biogenesis and thickening, and the synthesis of various compounds.This suggests that the response of eggplant to V. dahliae is the result of the integrated action of multiple processes, consistent with previous transcriptome analyses of cotton response to V. dahliae infection (Zhang et al., 2020;Xiong et al., 2021;Feng et al., 2023).Additionally, correlation analysis of the metabolome and transcriptome revealed associations between metabolites and co-expression network modules.We detected 261 metabolites that are simultaneously associated with module eigengenes and genes within key modules, involving multiple disease resistance-related pathways.Some metabolites have been shown to play roles in plant disease resistance, such as flavonoids, isoquinoline alkaloids, and cutin (Serrano et al., 2014;Zhang et al., 2022;Zhuang et al., 2023).Significant associations were observed between module genes and metabolites from the same pathway.These results indicate the reliability of the disease resistance-related modules obtained through gene co-expression network analysis.
The analysis of gene expression in plant-pathogen interaction pathways revealed significant expression differences in many disease resistance-related genes between LC-2 and LC-7.These differences reflect variations in the direct responses of the two strains to V. dahliae infection.For example, pathogenesis-related protein 1 (PR1) is known to be closely associated with plant pathogen resistance (Han et al., 2023).Among the five transcripts annotated as PR1, only one exhibits a similar differential expression pattern between LC-2 and LC-7 (Figure 5).Additionally, several members of the WRKY transcription factor family (including WRKY29, WRKY33, etc.) exhibit differential expression patterns between varieties.Numerous studies have demonstrated the crucial roles of WRKY family transcription factors in plant immunity and defense signaling processes (Eulgem and Somssich, 2007;Pandey and Somssich, 2009;Wani et al., 2021).The differential expression of these disease resistance-related genes between LC-2 and LC-7 likely contributes to the observed disparity in verticillium wilt resistance between LC-2 and LC-7.
Furthermore, the co-expression network analysis also revealed 13 key TFs that play potentially important roles in resistance (Table 2).These TFs belong to various families including AP2/ ERF, bHLH, GRAS, NAC, NF-YA, RWP-RK, HB-HD-ZIP, GARP-G2-like, and others.Serving as hub genes within the network modules, they exhibit co-expression with genes involved in plant hormone signaling pathways, defense compound biosynthesis, and defense responses (Figure 7), implying potential associations of these TFs with various biological processes underlying disease resistance.Notably, the roles of some TFs in verticillium wilt resistance have been validated.For instance, an AP2/ERF gene, GhTINY2, is strongly induced by V. dahliae in cotton, directly activating WRKY51 expression and promoting salicylic acid accumulation and signal transduction, thereby enhancing resistance to V. dahliae (Xiao et al., 2021).In this study, the AP2/ ERF gene ERF72 (SMECHR0902114) serves as the hub gene of the turquoise module, and its homolog has been demonstrated in Arabidopsis to act as a positive regulator mediating resistance to the necrotrophic pathogen Botrytis cinerea (Li et al., 2022b).Similarly, GhPAS1, a bHLH TF, is highly upregulated in cotton roots in response to V. dahliae invasion, and overexpression of GhPAS1 improves resistance to V. dahliae (Zhang et al., 2023).Two bHLH TF (SMECHR0902653 and SMECHR0400739) act in jasmonic acid were also identified as potential resistance genes in this study.Furthermore, the GRAS TF SCL14 (SMECHR0603017) were proved to play a key role in cis-jasmone induced indirect defense (Matthes et al., 2010).These findings underscore the effectiveness of our data mining approach.Furthermore, the roles of some TFs, such as the two HB-HD-ZIP family TFs Smechr0101299 and Smechr0802554, in disease resistance processes remain unclear.These transcription factors exhibit significant expression differences between LC-2 and LC-7 in response to V. dahliae.For instance, theSmechr0101299 is induced and upregulated in LC-2 upon V. dahliae infection at 20 days, whereas it remains downregulated in LC-7 (Figure 7B).Such differential expression patterns may contribute to the observed differences in disease resistance.Further investigation through methods such as transgenic approaches is needed to elucidate their roles in eggplant verticillium wilt resistance.

Conclusion
In short, our study has characterized the transcriptome and metabolome of V. dahliae in eggplant varieties with different degrees of resistance.V. dahliae infection affected the enzyme activities associated with disease resistance and altered the gene expression patterns and metabolite contents.The verticillium wiltresistant variety LC-2 exhibits higher activities of resistance-related enzymes compared to the sensitive variety LC-7.The identified DEGs and DAMs were primarily enriched in defense response, antibiotic metabolic response processes, sterol biosynthetic process, cell wall organization, plant-pathogen interaction pathways, and the disease resistance-related pathway in plant hormone signal transduction.We further identified 13 key transcription factors using a co-expression analysis.These results provide valuable insights into the molecular mechanism of eggplant resistance to verticillium wilt.
FIGURE 3 Summary of the metabolomes of LC-2 and LC-7.(A) Principal component analysis (PCA) of LC-2 and LC-7 samples.(B) The numbers of differentially accumulated metabolites (DAMs) in LC-2 and LC-7.(C, D) The statistics of KEGG enrichment of DAMs of LC-2 vs. LC-7 at 20 days (C) and 40 days (D).
(A) The numbers of differentially expressed genes (DEGs) in LC-2 and LC-7.(B) The GO enrichment results and corresponding numbers of genes for DEGs in LC-2 and LC-7.(C) Comparison and distribution of DEGs in 20-day and 40-day samples of LC-2 and LC-7, using 0-day samples as a control.(D-G) Expression profiles of 4017 LC-2specific DEGs (D), 4953 LC-7-specific DEGs (E), 351 DEGs with similar expression patterns between LC-2 and LC-7 (F), and 2861 DEGs with different expression patterns between LC-2 and LC-7 (G).The color of the heatmap represents the log2 fold change of the comparison between pairs.Red and blue respectively represent positive and negative fold change values, with higher saturation indicating higher absolute values.
FIGURE 5DEGs in KEGG pathways related to the pathogen defense response.(A) Plant-pathogen interaction (ko04626).(B) Disease resistance-related pathways in plant hormone signal transduction (ko04075).The blue rectangles represent genes, and the heatmap beside each gene shows the log2 fold change of 20-and 40-day samples for LC-2 and LC-7, respectively.The 0-day samples of LC-2 and LC-7 were used as controls to calculate the fold change for each variety.
FIGURE 6Topological structure and functional relevance of gene co-expression modules.(A) Gene co-expression sub-network constructed by five key modules (blue, green, tan, turquoise and yellow).Modules are represented by node genes with different colors.(B-F) Biological processes of blue (B), turquoise (C), yellow (D), tan (E), and green (F) modules.Biological processes are represented by nodes with module colors, and biological processes with similar functions are categorized in red circles and summarized by tags.

FIGURE 7
FIGURE 7Characterization of thirteen key hub transcription factors.(A) In the sub-network consisting of key modules, 13 key transcription factors were co-expressed with genes involved in pathways including plant hormone signaling, defense response to biotic stress, and synthesis of compounds related to defense response.They were also indirectly co-expressed with cell wall reinforcement pathway genes through network neighbors.(B) Differential expression patterns (log2 fold change) of the 13 key transcription factors The first four columns represent comparisons within the strain, using the 0-day samples of LC-2 and LC-7 as controls, while the last three columns represent comparisons between strains, using LC-7 as the control.The color saturation indicates the degree of change in fold change values.

FIGURE 9 V
FIGURE 9 V. dahliae-response patterns of the 13 transcription factors in LC-2 and LC-7.The unstressed expression level (Control: uninoculated) was assigned a value of 1.The error bar on each column represents the standard deviation (SD) of the three biological replicates (Student's t-test).20 d: 20 d after inoculation, 40 d: 40 d after inoculation.

TABLE 1
Changes in enzyme activities associated with disease resistance at different times (20 d and 40 d) after inoculation of V. dahliae.

TABLE 2
Transcription factors that act as hub genes in key modules.