Chieh-Qua (B. hispida Cogn. var. Chieh-qua How) inbred line “feicui” (FC), a common cultivar in Southern China, was provided by the Guangzhou Academy of Agricultural Sciences (Guangzhou Academy of Agricultural and Rural Sciences). The wilt-resistant lines were natural mutations isolated from “feicui” ( Supplementary Figure 1 ). The seedlings were grown in 32-well plates filled with an aseptic organic substrate at 28°C/20°C day/night temperatures in a greenhouse under a 16-h light/8-h dark photoperiod for approximately 20 days until the second true leaf stage.

The F. oxysporum provided by the Plant Protection Research Institute Guangdong Academy of Agricultural Sciences was cultivated on PDA solid medium at 28°C in the dark for 4 days, then cultured in potato dextrose broth on a shaker at 180 rpm at 28°C for 3 days. The spore suspension was diluted to 1×10⁵ spores per milliliter with sterile distilled water.

Chieh-qua’s leaves from both wilt-resistant and wilt-susceptible lines isolated from FC were chosen as the subjects of this study ( Table 1 ). Seedlings were infected with F. oxysporum by irrigation of the roots with a fungal spore suspension (3 mL per seedling) at the second true leaf stage. The true leaves of FC seedlings were harvested at 48 h after inoculation for RNA sequencing (RNA-seq) and metabolomic analysis. Three biological replicates were performed for each treatment, with 10 seedlings for each replicate, and each sample weighed approximately 3 g. Three sample groups were obtained: CK (control group, wilt-susceptible lines without pathogen inoculation), GB (wilt-susceptible lines with pathogen inoculation), and KB (wilt-resistant lines with pathogen inoculation). Two differential comparative analysis groups were established: wilt-susceptible (GB) vs. control (CK), as well as wilt-resistant (KB) vs. GB.

The total RNA from wax gourd leaf samples subjected to different treatments was extracted using the CTAB (cetyltrimethylammonium bromide) method (Gasic et al., 2004), with the following steps: The leaf samples were ground into powder in liquid nitrogen and transferred to 2-mL centrifuge tubes containing 1 mL of preheated CTAB extraction buffer (the CTAB extraction buffer was preheated in a water bath at 65°C and supplemented with 2% mercaptoethanol). The mixture of liquid and powder was vortexed thoroughly and then incubated in a water bath at 65°C for 5 min before adding an equal volume of chloroform/isoamyl alcohol (volume ratio of 24:1). The mixture was centrifuged at 1,200 rpm for 30 min, and the supernatant was transferred to new 2-mL centrifuge tubes. The chloroform/isoamyl alcohol extraction and centrifugation steps were repeated once more, and the pellet was discarded. The supernatant was mixed with 4 mol/L LiAc and incubated at 4°C for 2 h before centrifugation at 1,200 rpm for 10 min. The supernatant was discarded, and the pellet was washed three times with 70% ethanol. After discarding the supernatant, the pellet was air-dried in a laminar flow hood. The RNA was dissolved in 30 μL of DEPC-treated double-distilled water after treatment with DNase I (Invitrogen) to remove genomic DNA contamination.

Equal amounts of RNA from each sample’s three biological replicates were used to construct cDNA libraries. RNA purity (OD_260/280 and OD_260/230) was measured using a NanoPhotometer spectrophotometer, and RNA concentration was accurately measured using a Qubit 2.0 fluorometer. RNA integrity was assessed using an Agilent 2100 Bioanalyzer. Subsequently, cDNA libraries were constructed, followed by sequencing using the Illumina HiSeq platform.

The transcriptome sequencing data were initially processed using fastp (v0.19.5) (Chen, 2023) to remove low-quality sequences and adapters. Subsequently, HISAT2 (v2.2.1) (Kim et al., 2019) was employed to map the filtered reads (in fastq format) to the reference genome of wax gourd (Xie et al., 2019). The resulting BAM files were sorted using SAMtools (v1.18) (Danecek et al., 2021). Transcript assembly and quantification were performed using StringTie to assemble the mapped transcripts and generate a transcript annotation file (GTF format). STAR (v2.7.10b) (Dobin et al., 2013) was utilized to map the filtered fastq files to the reference genome with the newly assembled transcript GTF file. Finally, RSEM (Grabherr et al., 2011; Dobin et al., 2013) was used to quantify the transcripts. The transcript reads count was used to represent transcript resolution. The expression levels of genes were quantified from the transcriptome data mapped to the reference genome using featureCounts (v2.0.1) (Liao et al., 2014). Differential expression analysis of genes and transcripts with differential abundance between the resistance and susceptible lines relative to the control (KB vs. CK and GB vs. CK, respectively) was conducted using the R package DESeq2 (Wang et al., 2010). Transcripts and genes with |log2FoldChange| >1 and p _adj < 0.05 were considered differentially abundant or differentially expressed. The expression levels of genes and the abundance of transcripts were normalized using FPKM and presented accordingly.

The open reading frame of LOC120087936 and LOC120075251 was amplified by PCR and inserted into the pBI121 vector. The primers used for plasmid construction are listed in Supplementary Table 1 . The recombinant constructs, as well as empty plasmids, were transformed into Agrobacterium tumefaciens strain GV3101 using the freeze–thaw method, and then they were transiently expressed in N. benthamiana leaves through infiltration, as described previously (Ma et al., 2012).

Pathogenicity test assays were performed as previously described, with some modifications (Li et al., 2021). After 36 h of infiltration, the infiltrated plants were sprayed with fresh spore suspension (1 × 10⁵ conidia/mL) of F. oxysporum and transferred in a growth chamber under long-day conditions (LD, 25°C, 16 h light/8 h dark, light intensity of 150 μmol m⁻² s⁻¹). After a 1-week incubation, the plants were photographed using a digital camera. The experiment was repeated three times, and each treatment used six seedlings of N. benthamiana.

For the fungal biomass assay, a well-washed piece of infected N. benthamiana leaf (approximately 1 cm²) was used for DNA extraction using a Fungal DNA kit (Omega, United States) according to the manufacturer’s protocol. DNA-based qPCR was performed with 2× iTaq™ Universal SYBR Green Supermix (Bio-Rad, Hercules, CA, USA). Relative fungal biomass was calculated as a ratio (FoEF1α/NtEF-1α) represented by the equation 2 ^{[CT(NtEF-1α)-CT(FoEF1α)]} as previously described (Park et al., 2012). The primer pairs for qPCR are listed in Supplementary Table 1 .

The upregulated and downregulated DEGs selected from the two comparison groups (GB vs. CK and KB vs. GB) were subjected to Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis, respectively, and the internal function enricher of the R language package cluster Profiler (Wu et al., 2021) was used to perform functional enrichment analysis (the threshold is p < 0.05, q < 0.05). Copy the output GO column and Qvalue column to Revigo (http://revigo.irb.hr/) (Supek et al., 2011) for processing and draw using CirGO software (https://github.com/IrinaVKuznetsova/CirGO) (Kuznetsova et al., 2019) GO enrichment analysis circle plot. KEGG enrichment analysis results were visualized with GraphPad.

The RGAugury pipeline (Li et al., 2016) was used to identify RGAs in the wax gourd genome, including four major categories: NBS, RLK, RLP, and TM-CC. Log2foldchanges of the differentially expressed RGA were used to draw the ridge plot using the R package ggridges and ggplot2 (Villanueva and Chen, 2019). The protein sequences of plant transcription factors (TFs) were downloaded from the plantTFDB4.0 database (Jin et al., 2017). Subsequently, the protein sequences of wax gourd were subjected to a comparative analysis against the downloaded plant TF protein sequences using the diamond blastp (Buchfink et al., 2021). Proteins from wax gourd exhibiting sequence similarity below this 1E-5 were considered as the TFs of wax gourd.

The Pearson correlation coefficient (PCC) between TFs and RGAs was calculated, and their significance was determined using the Benjamini–Hochberg (BH) method. TF–RGA gene pairs with a |PCC| > 0.8 and a q-value < 0.001 were considered co-expressed. The co-expression network was plotted using Cytoscape (Smoot et al., 2011).

The SNP data used in this study were sourced from Xie et al. (2019). To assess genetic differentiation among populations, the sliding window approach implemented in VCFtools (Danecek et al., 2011) was employed. Specifically, windows of 2,500 base pairs (bp) in size were moved along the RGA gene intervals and their 4,000-bp upstream and downstream regions with step sizes of 50 bp and 250 bp, respectively.

The samples of both control and infected groups of wilt disease-resistant varieties were subjected to vacuum freeze-drying. Subsequently, they were ground into a powder using a grinding mill (MM400, Retsch) at a frequency of 30 Hz for 1.5 min within liquid nitrogen. Following this step, 100 mg of the powdered sample was dissolved in 1 mL of extraction solution (70% methanol solution). The dissolved samples were stored in a refrigerator at 4°C overnight, during which they were vortexed three times to enhance the extraction yield. Afterward, the samples were centrifuged at 10,000g for 10 min, and the supernatant was collected. The filtered samples were passed through a microporous filter membrane with a pore size of 0.22 μm and stored in sample vials for subsequent liquid chromatography–tandem mass spectrometry (LC-MS/MS) analysis.

Based on the Metware database (Metware Biotechnology Co., Ltd, Wuhan, Hubei, China) and publicly available metabolite information databases, substance qualitative analysis was conducted using secondary mass spectrometry data. During the analysis, isotope signals were excluded, as well as duplicate signals originating from K+ ions, Na+ ions, NH4+ ions, and fragments of larger molecules themselves (Chen et al., 2013). Metabolite quantification is accomplished through the use of a triple quadrupole mass spectrometer employing a multiple reaction monitoring (MRM) mode. After obtaining mass spectrometry data for different samples, peak areas of all metabolite spectra are integrated, and peak integration is corrected for the same metabolite across different samples (Fraga et al., 2010). Based on the results of partial least squares–discriminant analysis (PLS-DA), we can initially screen for metabolites that exhibit differences between different varieties or tissues. Additionally, we can further refine our selection by incorporating p-values or fold change values from univariate analysis. The combination of fold change and the variable importance in projection (VIP) from the OPLS-DA model is used to identify differential metabolites. The selection criteria are as follows: Metabolites with a fold change ≥2 or ≤0.5 are chosen. A fold change of 2 or more or 0.5 or less indicates significant differences between the control and experimental groups. Building upon the above criteria, metabolites with a VIP value ≥1 are selected. VIP values represent the strength of the impact of the intergroup differences for the corresponding metabolites in the model’s discriminative classification of group samples. Generally, metabolites with a VIP value of ≥1 are considered to be significantly different.

The PCC between the expression level of DEGs and the content of the differential metabolites was calculated, and their significance was determined using the BH method. DEGs and differential metabolites with a |PCC| > 0.8 and a q-value < 0.0005 were considered significant correlations.

F. oxysporum is a widely prevalent fungal pathogen in plants (Fravel et al., 2003), and the cultivation industry of chieh-qua is also affected by this disease. Therefore, understanding the corresponding mechanisms of chieh-qua against F. oxysporum infection is of significant value for controlling the pest and disease encountered during the chieh-qua cultivation process. In this study, RNA-seq was performed on samples from the CK, GB, and KB. Compared to CK, there were 1,912 upregulated DEGs and 2,818 downregulated DEGs in GB ( Figure 1A ). Compared to the GB inoculated with the F. oxysporum, 3,448 genes were highly expressed in KB, with 3,072 DEGs being downregulated. In the two comparison groups, there were 4,071 (56.7%) commonly shared DEGs, with 660 (9.2%) specific DEGs in GB vs. CK and 2,449 (34.1%) specific DEGs in KB vs. GB ( Figure 1B ). KEGG enrichment analysis of DEGs revealed enrichment of secondary metabolite biosynthesis pathways involved in biological defense functions among the upregulated genes after inoculation in the GB samples. ( Figure 1C ). At the same time, the rapid synthesis of secondary metabolites is activated by upstream signal regulation. The transport of plant hormones from their synthesis sites to target organs and their subsequent binding to receptors play crucial roles in plant physiology (Dermastia, 2019). The significantly upregulated genes were enriched in GO terms such as “hormone signal transduction” in the GB vs. CK groups, indicating that the infection process may activate multiple metabolite synthesis pathways. Simultaneously, many essential basic functions related to genes associated with processes such as photosynthesis appear to be downregulated ( Figure 1D ). Compared to GB, genes that were highly expressed in KB were enriched in pathways related to energy and growth, such as photosynthesis and nitrogen metabolism. This suggests that KB samples can maintain a better growth status after infection ( Figure 1E ). However, genes that were downregulated in KB were enriched in certain metabolic synthesis-related GO terms, such as phenylpropanoid biosynthesis ( Figure 1F ). This may be because these secondary metabolites do not participate in KB’s adaptation to F. oxysporum infection.

To explore the functions of DEGs between KB (wilt-resistant lines) and GB (wilt-susceptible), LOC120087936 and LOC120075251, which exhibited high expression levels in KB, were selected for further analysis. Indeed, LOC120087936 encodes a homolog of Arabidopsis MILDEW RESISTANCE LOCUS O 13 (AtMLO13), which belongs to a large family of seven-transmembrane domain proteins that are specific to plants and are involved in conferring resistance to biotrophic powdery mildew fungus in barley. LOC120075251 encodes a homologous version of Arabidopsis MDIS1-INTERACTING RECEPTOR LIKE KINASE2 (MIK2), which is a receptor heteromer involved in responding to various environmental stresses, including cell wall integrity sensing, salt stress tolerance, and resistance to F. oxysporum (Julkowska et al., 2016; Van der Does et al., 2017; Engelsdorf et al., 2018; Chaudhary et al., 2020). The plant transient expression vector pBI121 plasmid was used to create constructs of 35S: LOC120087936 and 35S: LOC120075251. The A. tumefaciens harboring constructs were infiltrated into at least four leaves (per seedling) of N. benthamiana for transient expression. After 36 h of infiltration, two transient expressed lines were evaluated regarding their resistance to F. oxysporum. At 7 days after the wounded leaves were sprayed with F. oxysporum, the control tobacco leaves (CK, pBI121 plasmid only) were wilted and shorter. In contrast, the leaves of overexpression of LOC120087936 and LOC120075251 were growing well and appeared healthy ( Supplementary Figure 2A ). To further determine whether the expression of LOC120087936 and LOC120075251 affected the fungal growth in planta, we estimated the relative fungal biomass in the infected leaves by DNA-based quantitative PCR (q-PCR). The assays showed that the relative fungal biomass was lower in the expressed LOC120087936 and LOC120075251 plants compared with CK ( Supplementary Figure 2B ). These results showed that LOC120087936 and LOC120075251 play an essential role in the pathogenicity of F. oxysporum. The differential expression results in this study provide a resource for the functional validation of F. oxysporum-resistant genes.

GO enrichment analysis of upregulated and downregulated DEGs in GB vs. CK provided more comprehensive information for assessing gene functions related to F. oxysporum infection. Specifically, the upregulated DEGs in GB vs. CK were mainly enriched in pathways such as the auxin-activated signaling pathway, adaptation to oxidative stress, and cinnamic acid biosynthetic process ( Figure 1G ). These enriched secondary metabolite synthesis and regulatory pathways further underscore the important role of secondary metabolites in the adaptation to F. oxysporum infection. Similarly, the downregulated DEGs in GB vs. CK were mainly enriched in pathways such as sulfate assimilation, lipid metabolic process, glycerol ether metabolic process, response to light stimulus, and photosynthetic electron transport in photosystem I ( Figure 1H ), consistent with the KEGG enrichment results.

The upregulated DEGs in KB vs. CK were primarily enriched in pathways related to protein folding, transport, xenobiotic transmembrane transport, response to light stimulus, response to wounding, and flavonoid biosynthetic process ( Figure 1I ). Conversely, the downregulated DEGs in KB vs. CK were mainly enriched in pathways such as DNA replication initiation, L-phenylalanine catabolic process, multicellular organism development, auxin-activated signaling pathway, and defense response to fungus ( Figure 1J ). These findings provide further insights into the molecular mechanisms underlying the adaptation of KB to F. oxysporum infection, highlighting the complex interplay of various biological processes and pathways involved in the host–pathogen interaction.

Previous studies have demonstrated that the plant MAPK signaling pathway and plant hormone signaling transduction, including salicylic acid (SA), are key factors in regulating plant immunity. Here, we compared the gene expression changes of genes in the MAPK pathway in GB samples before and after infection ( Figure 2 ). BAK1 showed significantly upregulated expression after infection. Studies have shown that BAK1 acts as a positive regulator in the MAPK signaling pathway (Chinchilla et al., 2007). Its upregulation activates downstream MEKK1, which is the starting point of several pathways formed by downstream MAPKs (Zipfel et al., 2004). Positive feedback responses activate various downstream pathways, such as the activation of camalexin synthesis through WRKY33, which is a secondary metabolite that inhibits bacterial and fungal infection (Koprivova et al., 2019). During the plant’s adaptation to microbial infection, various defense genes are regulated by ethylene (Ecker and Davis, 1987). The results of this study found that genes COPA, MPK3, and ERF1 on this pathway showed upregulation after infection, revealing the role of the ethylene-regulated pathway in the adaptation to F. oxysporum infection. Similarly, genes ANP1, MPK7/14, and PR1 involved in the activation of cell death, H₂O₂ production, and pathogen defense pathways were upregulated after F. oxysporum infection. SA is a plant hormone, and several genes related to its signal transduction, including NPR1 (Nonexpresser of PR Genes 1), TGA (Transcription Factor GATA), and PR1 (Pathogenesis-Related Gene 1), showed a similar expression pattern, significantly upregulated in GB compared to CK. These results not only suggest that the SA signaling pathway may play a role in the adaptation to F. oxysporum infection but also provide insights for further experiments to study its effects on downstream secondary metabolite synthesis.

RGAs are an important class of disease resistance genes. In this study, 664 RGAs were identified in the wax gourd, with 82, 36, 441, and 105 genes belonging to the NBS, RLP, RLK, and TM-CC gene families, respectively. As RGAs associated with plant resistance, differential expression analysis can identify candidate genes for further screening and application. From the distribution of RGAs’ log2FoldChanges in GB vs. CK ( Figure 3A ) and KB vs. GB ( Figure 3B ), it can be observed that some members of the TM-CC gene family not only respond to infection in GB but also show differences between GB vs. KB. A total of 23 TM-CC genes are differentially expressed in both GB vs. CK and GB vs. KB, indicating that these genes not only respond to F. oxysporum infection but also exhibit a stronger adaptation in resistant varieties. The number of differentially expressed RLP, TM-CC, and NBS in KB vs. GB is significantly fewer than that of the RLKs (177), with 17, 33, and 33 genes overlapping with those responding to F. oxysporum infection in GB, respectively.

To explore potential TFs regulating these RGAs, the PCC between TFs and differentially expressed RGAs was calculated. Using a threshold of |PCC| > 0.8 and q-value < 0.001, a total of 119 TFs were significantly correlated with RGAs ( Figure 3C ). Among them, the most abundant TFs were bHLH, which were mostly associated with the regulation of RLKs. There were 6, 2, 40, and 8 NBS, RLP, RLK, and TM-CC genes, respectively, potentially regulated by TFs. These results provide a data foundation for understanding the regulatory patterns of resistance genes in chieh-qua in the adaptation to F. oxysporum infection.

The study of Xie et al. divided the wax gourd population into cultivar1, cultivar2, landrace, and wild, as well as provided resequencing information (Xie et al., 2019), which laid the foundation for our analysis of selection signals for differentially expressed RGAs in the populations. The calculation of pi values for all RGAs showed that the diversity of RGAs in the landrace population was similar to that in the wild population ( Figure 4A ). However, the diversity of the TMCC, RLP, RLK, and TM-CC families in the cultivar2 population was the lowest, indicating a reduction in the genetic diversity of resistance genes in the cultivar2 population during breeding selection. Correspondingly, the Fst values between cultivar2 and wild were 0.3339, 0.4608, 0.4259, and 0.4732, respectively, indicating the possible presence of selection signals in the RGA regions between wild and cultivar2 populations.

To further screen for RGAs that may be under selection, this study analyzed the distribution of Fst values between populations for all differentially expressed RGA genes and their upstream and downstream 4-kb regions. Eight RGA genes (LOC12006786, LOC120073192, LOC120075251, LOC120075583, LOC120079812, LOC120083870, LOC120085463, and LOC120090401) had Fst values greater than 0.5 between cultivated (cultivar1 and cultivar2) and wild populations, indicating significant genetic differences between cultivated and wild populations within the 4-kb range upstream and downstream of these genes ( Figure 4B ). These results provided more candidate genes for resistance research.

Gene alternative splicing can lead to the existence of multiple transcriptional isoforms. In this study, a total of 71,771 transcripts were identified, originating from the expression of 21,004 genes. Comparing the differential expressed transcripts with that of genes can provide new insights into the role of gene alternative splicing in the adaptation of chieh-qua to F. oxysporum infection. Firstly, compared to CK, 1,912 upregulated and 2,818 downregulated DEGs were identified in GB, along with 4,158 upregulated and 4,891 downregulated differential transcripts ( Figure 5A ). The proportions of upregulated and downregulated DEGs, as well as non-DEGs in GB vs. CK were 13.4%, 9.1%, and 77.5%, respectively ( Figure 5B ). Among the 2,086 genes without differential expression but with differential transcripts, 152 genes had both upregulated and downregulated transcripts. Among the 16,274 genes with no differential expression, transcripts from 14,188 genes also showed no differential expression between GB and CK. Additionally, among the transcripts from the 2,086 genes, the percentages of transcripts with upregulation, downregulation, and both upregulation and downregulation between GB and CK were 57%, 36%, and 7%, respectively ( Figure 5C ). This indicates that different transcripts from the same gene may exhibit different expression patterns. Among the upregulated DEGs in GB vs. CK, 97.3% of genes also had upregulated transcripts ( Figure 5D ), while among the downregulated DEGs, 97.4% of genes also had downregulated transcripts ( Figure 5E ). Genes associated with functions such as Photosystem II and chloroplast envelope showed consistent expression patterns between transcripts and genes ( Supplementary Figure 3 ). Comparatively, 3,448 upregulated and 3,072 downregulated DEGs were identified in KB compared to GB, along with 7,052 upregulated and 6,695 downregulated differential transcripts ( Figure 5F ). The proportions of upregulated DEGs, downregulated DEGs, and non-DEGs in GB vs. CK were 14.6%, 16.4%, and 69%, respectively ( Figure 5G ). The percentages of upregulated DEGs, downregulated DEGs, and genes with both upregulated and downregulated transcripts were similar to those in GB vs. CK ( Figures 5H, I ). However, downregulated transcripts enriched in pathways such as auxin polar transport and auxin homeostasis in upregulated genes deserve further attention ( Supplementary Figure 4 ). Among the genes with no differential expression between KB and GB, there were 2,576 transcripts that were significantly differentially abundant, with proportions of upregulation, downregulation, and both upregulation and downregulation being 36%, 53%, and 11%, respectively ( Figure 5J ).

LOC120069995 is a gene encoding a protein with the McbC_SagB-like_oxidoreductase functional domain. The gene structure and expression patterns of its transcripts indicate that LOC120069995 is not differentially expressed in both GB vs. CK and KB vs. GB ( Figure 6A ). The gene has five exons; its transcript MSTRG.25886.12 has five exons, with upregulated expression in GB vs. CK and downregulated expression in KB vs. GB ( Figure 6B ); MSTRG.25886.7 has four exons, with downregulated expression in GB vs. CK and upregulated expression in KB vs. GB ( Figure 6C ); MSTRG.25886.7 has two exons, with no differential expression in both GB vs. CK and KB vs. GB ( Figure 6D ). These results provided insights into the role of gene alternative splicing in the adaptation of chieh-qua to F. oxysporum infection.

In this study, metabolites from GB and CK samples were analyzed for metabolomics using ultra-performance liquid chromatography (UPLC) and MS/MS. Through annotation in the Metware database, a total of 672 known metabolites were detected ( Supplementary Tables 2 , 3 , 4 ). Among them, 417 metabolites with annotation information showed differential contents between GB and CK, including 44 upregulated and 373 downregulated metabolites ( Figure 7A ). These metabolites belong to various categories such as flavanone, terpene, alkaloids, organic acids and derivatives, and phenylpropanoids. Metabolites of the organic acid and derivative types, such as oxalic acid and isochlorogenic acid B, exhibited the highest upregulation in content after F. oxysporum infection ( Figure 7B ). These results indicated the significant role of these substances in the adaptation of chieh-qua to F. oxysporum infection. Meanwhile, metabolites of the organic acid and derivative types, such as citraconic acid, showed a decrease in content after F. oxysporum infection, suggesting a potential inhibitory effect of F. oxysporum infection on their synthesis.

To further explore the relationship between genes and these metabolites, we performed a correlation analysis between DEGs involved in plant MAPK signaling pathway and plant SA signaling transduction identified in the GB vs. KB comparison group, and differential metabolites. The results revealed 14 DEGs and 26 metabolites (belonging to 11 classes) showing a positive correlation ( Figure 7C ). In particular, genes such as LOC120071745 (Histidine-containing phosphotransfer protein 3), LOC120077686 (BRASSINOSTEROID INSENSITIVE 1-associated receptor kinase 1), and LOC120070643 (Respiratory burst oxidase homolog protein B) were associated with multiple metabolites in the positive correlation network. Lipids, flavanones, and other substances play important roles in the interaction between plants and pathogens (Aseel et al., 2019; Jeon et al., 2020). This study identified 647 genes, including LOC120067029 (ATPase family AAA domain-containing protein 1-like), LOC120068100 (zinc finger CCCH domain-containing protein 15-like), and LOC120080874 (proteinase-activated receptor 2), whose expression was significantly positively correlated with these metabolites. Furthermore, 41 DEGs and 72 metabolites (belonging to 17 classes) showed a negative correlation ( Figure 7D ), among which organic acids and derivatives had 19 metabolites connected to other genes in the negative correlation network. Organic acids and derivatives may have more connections in the plant MAPK signaling pathway and plant SA signaling transduction process.