Among the 1102 + 734 upregulated DEGs in the G1 and G2 gene sets, the KEGG classification map revealed that 11 genes associated with cellular processes were annotated to the peroxisome pathway, and the MAPK signaling pathway was the most enriched in the environmental information processing signaling pathway and plant hormone signal transduction. In total, 21 genes involved in ubiquitin-mediated proteolysis were annotated to genetic information processing, and 36 genes were involved in phenylpropanoid biosynthesis. Moreover, 82 genes were associated with the plant-pathogen interaction pathway ( Figure 4A ). Carbon metabolism and plant-pathogen interactions were the most enriched pathways in the KEGG analysis. Phenylpropanoid biosynthesis, α-linolenic acid metabolism, flavonoid biosynthesis, and sphingolipid metabolism pathways were also enriched ( Figure 4B ).

In the GO classification diagram, DEGs were the most abundant metabolic process, and 120 genes related to response to stimulus, 37 genes related to signaling and detoxification, and four genes related to immune system processes were annotated. In terms of cell composition (CC), a single cell enriched the most gene entries. The binding and catalytic activities were annotated to the molecular function (MF) with the largest number of DEG entries, and 15 genes involved in antioxidant activity were annotated ( Supplementary Figure S2 ).

Owing to the large number of plant disease resistance genes involved in BP and MF, GO enrichment analysis and GO enrichment hierarchy analysis were conducted. The GO enrichment bubble diagram revealed that the tricarboxylic acid cycle had the largest number of DEG entries annotated in BP, and 10 genes involved in response to chitin were significantly enriched, which were significantly enriched with the host immune response; additionally, defense response genes were annotated ( Figure 5A ). The GO enrichment chord diagram revealed the gene expression of chitin in response to pathogens, in which DN10131 was significantly upregulated ( Figures 5B, C ), triggered the immune response of plants, and participated in the process of response to oxidative stress, which indicated that it plays an important role in plant immune signaling. In the MF bubble map, 69 genes were annotated to protein kinase activity, and the highest enrichment was protein serine/threonine kinase activity, which may act as PRRs on plant cell membranes. In total, 13 peroxidase activity genes and four phenylalanine ammonia-lyase activity genes were also enriched in small amounts ( Figure 6A ). DN1051 (GO:0003824), DN10922, and DN10613 were involved in the redox process, and DN10131 and DN10068 were upregulated in the G1 and G2 groups ( Figures 6B, C ).

Transcription factors (TFs) were predicted for the 734 co-expressed upregulated genes using a Venn diagram ( Supplementary Figure S3A ), and heatmaps were generated for the top three types of TFs with the most gene entries ( Supplementary Figure S3B ), among which the first three gene families were AP2ERF-ERF, WRKY, and C2H2. Numerous studies have identified the important roles of TFs in plant defense. One of the most studied MYB TFs is Botrytis susceptible 1 (BOS1), which restricts necrosis triggered by B. cinera (Mengiste et al., 2003). In addition, WRKY33 is phosphorylated by MPK3/MPK6 to regulate the biosynthesis of phytoalexin in response to pathogen infection (Ren et al., 2008). Arabidopsis CCCH protein C3H14 contributes to basal defense against B. cinerea mainly through the WRKY33-dependent pathway (Wang et al., 2020). WRKY transcription factors act as key regulators in plant disease resistance signaling pathways (e.g. strawberry, tomato) (Shu et al., 2021), responding to biotic stress by regulating genes involved in innate immunity, hormone signaling pathways, and phytoalexin synthesis (Lee et al., 2023).

To assess genes with the same expression pattern, a co-expression trend analysis of three sets of DEGs was conducted. The expression dynamics can be clustered into 11 expression patterns ( Supplementary Figure S4 ). According to the co-expression trend diagram, three patterns with consistent upregulated expression trends were selected for subsequent analysis ( Supplementary Figure S4 A, H, and K). Among the 189 genes whose expression was consistent, those related to plant-pathogen interactions, plant hormone signal transduction, and the MAPK signaling pathway were significantly enriched ( Figure 7A ), which was the result of interactions between P. cyrtonema and B. deweyae.

Metabolites or metabolite intermediates associated with antibacterial activity, such as phenylpropanoid biosynthesis, flavonoid biosynthesis, ubiquinone, and other terpenoid-quinone biosynthesis, and α-linolenic acid metabolism, were enriched in the other two modes. Additionally, sphingolipid metabolism was also enriched in the other two modes ( Figures 7B, C ).

When attacked by B. deweyae, the receptor protein kinase BAK1/FLS senses pathogenic effector factors and triggers the MAPK signaling immune response ( Figure 1 ), including ROS mechanisms and the activation of transcription processes to resist pathogens. This indicated that brassinosteroid-related receptor kinase 1 (BAK1) and LRR receptor-like serine/threonine protein kinase (FLS) play key roles in activating MAPK signaling as two LRRs, while WRKY TFs activate downstream genes for B. deweyae resistance, with WRKY 33 being the most annotated transcription factor, which activates the camalexin synthesis process downstream. WRKY 22 TFs mediate cell death by activating the transcription of the senescence receptor kinase FRK1 ( Supplementary Figure S5 ). The MAPK signaling pathway is an essential component of plant immunity, mediating hypersensitive response (HR) and associated cell death (Adachi et al., 2016). MPK3/MPK6 phosphorylates WRKY transcription factors including WRKY22, WRKY23, WRKY29, WRKY46 and WRKY53, which mediate the pathogen induced plant defense response (Xie et al., 2021). In our previous analysis of TFs, we propose that this conserved MAPK-WRKY regulatory module likely functions in P. cyrtonema to coordinate defense against fungal pathogens.

Based on the KEGG and GO classification information, the protein interaction network diagram nodes were edited, and the functions and interactions of each gene were displayed, among which the ubiquitin-mediated protein hydrolysis gene DN1274_c1_g1 (ko04120) showed binding relationships with multiple disease resistance genes, including the MAPK signaling pathway, α-linolenic acid metabolism pathway, and plant hormone signaling pathway. This indicated that protein ubiquitination is crucial for activating and regulating gene expression in the immune pathways of P. cyrtonema. Four genes (DN3411_c0_g1/DN21425_c0_g1/DN53217_c0_g1/DN880_c2_g1) (ko00600) involved in sphingolipid metabolism most strongly interacted with DN1274_c1_g1, and these four genes were found to be closely related. The two genes involved in MAPK signaling, including DN1158_c0_g1 (ko04626) and DN65_ c0_g1 (ko04016), also presented a high degree of correlation in the network diagram, and multiple pathways, such as chitin and sphingoesteric acid metabolism, promoted MAPK signaling ( Figure 8 ).

Real-time PCR verification of the screened DEGs was conducted, and the correlation between the RNA-seq data and relative quantitative data was evaluated. The results of the qPCR analysis of the 12 DEGs are shown in Figure 9 ; the genes confirmed included the receptor protein kinase gene in MAPK signaling pathway DN65_c0_g1 (ko04016), the DN14157_c0_g1 (ko00940) gene of phenylpropanoid biosynthesis or catalytic activity (GO:0003824), the DN1401_c0_g1 (ko00600) gene of sphingomate metabolism, and the phenylpropanoid biosynthesis DN13372_c0_g1 (ko00940) or antioxidant activity (GO:0016209).

We used the LC-QTOF platform for qualitative and quantitative metabolomics analysis of 16 samples and detected 24,245 peaks under default mode, of which 4,975 metabolites were annotated. By conducting PCA of correlations between samples to evaluate biological reproducibility within the group ( Supplementary Figure S6A ), the Spearman rank correlation coefficient (r) was used as an indicator of biological reproducibility ( Supplementary Figure S6B ). Based on the results of the PCA variable analysis and correlation assessment, the differential metabolite (DAMs) data obtained from this metabolomics study were found to be reliable and can be used for downstream DAMs screening. Based on the results of OPLS-DA ( Supplementary Figure S7 ). Permutation testing confirmed the robustness of the OPLS-DA model. All permuted R²Y and Q²Y values (blue dots) fell below the original model’s metrics (red star), indicating the model’s validity (p < 0.001). The intercepts of the regression lines (Q² < 0.05) further support non-overfitting ( Supplementary Figure S8 ). The models built in the three comparison groups were all valid and could be used for subsequent analysis.

The DAMs in each differential group were annotated using the KEGG database, and the substances with antifungal activity in the three groups were annotated to the biosynthesis of other secondary metabolites. Most of the annotated metabolites were antibiotic metabolites, such as neomycin, kanamycin, and gentamicin. Flavonoid biosynthesis, indole alkaloid, quinoline alkaloid, and alkaloid-related metabolites were also annotated to the three groups. Arachidonic acid metabolism, α-linolenic acid metabolism, and cutin, suberin, and wax biosynthesis were the most annotated DAMs in lipid metabolism. Other DAMs, such as diterpenoid biosynthesis and carotenoid biosynthesis, were annotated to the metabolism of terpenoids and polyketides ( Figure 10 ). Eight DAMs were annotated to the phenylpropanoid biosynthesis pathway in the M2. Although KEGG analysis revealed apparent enrichment of neomycin, kanamycin, and gentamicin, these compounds are not synthesized in plants. This observation likely reflects conserved primary metabolic components shared between plants and bacteria, such as enzymes and intermediates involved in amino sugar and nucleotide sugar metabolism (e.g., UDP-glucose). To improve pathway resolution and avoid such artifacts, propose future studies will employ complementary tools like MapMan or Reactome for functional enrichment analysis, which may better distinguish plant-specific metabolic processes.

The abundance score (DA score) reflects the overall changes in all metabolites in the metabolic pathway, and a positive score indicates that the expression of all annotated metabolites in the pathway is upregulated, and vice versa. The DA scores for the three comparison groups are shown in Supplementary Figure S9 . The five pathways with the greatest enrichment in each difference group were assessed by association analysis and combined with the results of enrichment network diagram analysis. In the M1 comparison group, the sulfur- relay system had the highest score, and the biosynthesis of alkaloids, plant hormone signaling, and the biosynthesis of cutin, serine, and wax were also upregulated. The enrichment network diagram showed that indole alkaloids and other biosynthesized alkaloids were the most annotated differentially expressed metabolites (with 11 metabolites in total) ( Supplementary Figure S9A ). In the M2 comparison group, β-lactam resistance had the highest score, and the metabolites associated with phenylpropanoid biosynthesis increased in score ( Supplementary Figure S9B ). In the M3 comparison group, the highest score was for α-linolenic acid metabolism, accounting for 11 metabolites, followed by microterpene biosynthesis ( Supplementary Figure S9C ). Furthermore, the bubble plot of KEGG enrichment factor for DAMs confirmed significant enrichment of cutin, serine, and wax biosynthesis ( Supplementary Figure S10A, B ) and alpha-linolenic acid metabolism ( Supplementary Figure S10C ) in the metabolome of P cyrtonema. The results of the KEGG enrichment analysis of the transcriptome revealed that α-linolenic acid metabolism was the most significant metabolite involved in disease resistance.

Subsequently, PCA was performed separately for the three differential groups, and the degree of association between the DAMs and the DEGs was evaluated by reducing the dimensionality of the DAMs to a few independent variables. The results of the analysis are shown in Supplementary Figure S12 . By comparing the pathways associated with genes in the transcriptome and the pathways associated with metabolites in the metabolome, the number of common pathways involved was obtained, as shown in Supplementary Figure S12 , with 79, 78, and 83 common metabolic pathways in the three differential groups, respectively.

The DEGs and metabolites in the three comparison groups were uniformly annotated on the KEGG pathway, and the top 30 pathways with significantly enriched differential genes/metabolites were plotted as a bar chart. Bubble plots were made using the KEGG pathways coenriched by transcriptomics and metabolomics. In the MG1 grouping bar chart, metabolic pathways related to keratin, serine, and wax biosynthesis were significantly different. As shown in the KEGG pathway bubble plot on the right ( Supplementary Figure S13A ), sphingolipid metabolism, flavonoid biosynthesis, phenylpropanoid biosynthesis, α-linolenic acid metabolism, and plant hormone signaling were significantly enriched. These pathways were also enriched in the MG2 group ( Supplementary Figure S13B ). In the MG3 group, genes and metabolites related to α-linolenic acid metabolism presented the greatest significance, whereas the DEGs and metabolites related to plant hormone signaling and keratin, serine, and wax biosynthesis presented increased significance ( Supplementary Figure S13C ).

In this study, the networks involved in Jasmonic acid (JA) synthesis and sphingolipid metabolism were annotated. Sphingolipids are a class of structurally complex lipid molecules containing long-chain sphingosine groups, with galactosyl sphingosine and sphingosine being key intermediate products in the sphingolipid metabolism pathway of P. cyrtonema ( Figure 11A ). Both metabolite levels increased after B. deweyae infected P. cyrtonema, and their gene expression also increased, which aligned with the transcriptome network diagram. Another annotated pathway was the synthesis of the plant hormone JA; we found ( Figure 11B ) that JA synthesis is influenced by abscisic acid (ABA), brassinolide, abscisate, and zeatin and is closely related to genes associated with zeaxanthin and JA. JA is an endogenous hormone that plays a very important role in plant disease resistance, and its gene expression is regulated by WRKY TFs, which can trigger HR-mediated disease resistance in plants (Shu et al., 2021). Zeatin is a novel plant hormone that can regulate various physiological mechanisms, such as growth, development, phototropism, and stress tolerance. Although abiotic stresses such as drought in wheat (Triticum aestivum) and cadmium toxicity in Desmodesmus armatus have been extensively studied (Piotrowska-Niczyporuk et al., 2024; Wang et al., 2023), research on their responses to biotic stress—particularly pathogen infection—remains relatively underexplored. In this experimental, the metabolic processes of zeatin are not significant, the network diagram showed that most genes involved in its synthesis are the same as those involved in JA synthesis, indicating that a synergistic relationship may exist between JA and zeatin, which can coregulate P. cyrtonema in the B. deweyae antidisease process, a hypothesis requiring functional validation through future interrogation of their crosstalk mechanisms.

According to the phenylpropion synthesis pathway diagram ( Figure 12 ), during B. deweyae infection, the activities of three metabolites, phenylpyruvate, p-coumaric acid, and leucopelargonidin, significantly increase in P. cyrtonema. Phenylpyruvate is the precursor for phenylalanine (Pen) synthesis, which is primarily synthesized into Pen by transferases such as aspartate aminotransferase (AST) and glutamate dehydrogenase (GDH). The phenylpropanoid biosynthesis pathway may indirectly influence the salicylate-mediated disease resistance signaling pathway. Studies have shown that p-coumaric acid has a significant inhibitory effect on the walnut anthracnose pathogen Colletotrichum gloeosporioides. In this study, Pen activated the downstream flavonoid biosynthesis pathway, where the chalcone synthase (CHS) gene plays a crucial role. In another pathway involved in flavonoid synthesis, leucopelargonidin promotes the formation of dihydrokaemplerol under the action of anthocyanin synthase (ANS), with rutin participating as a precursor in the downstream flavonoid synthesis pathway. Pen acts as a key substance in the phenylpropanoid metabolic pathway (Li et al., 2017). The phenylpropanoid metabolism-based defense responses towards pathogen attacks have been widely characterized in plants (Cui et al., 2024; Xiao et al., 2022). It can directly or indirectly activate plants to produce various secondary metabolites with antibacterial effects, such as flavonoids and phenols, B. cinerea quercetin dioxy genase (BcQdo) (catalyzing flavonoids degradation gene) involves in B. cinerea virulence towards the Panax ginseng, △BcQdo mutants showed increased flavonoids accumulation and reduced disease development (Chen et al., 2022).

Linolenic acid (LnA) and linoleic acid (LA) are the main unsaturated fatty acids in plants (Lim et al., 2017). They are natural immune inducers in plants and can act as signaling molecules to induce the natural, systematic, and lasting resistance of plants to diseases (Dedyukhina et al., 2014). In this study, the metabolite α-linolenic acid was detected in the metabolome, and its activity significantly increased with increasing inoculation time ( Figure 13 ), indirectly correlated with the plant hormone JA signaling pathway, with the metabolic levels of both showing a positive correlation. In the hormone signaling pathway, the JAZ (jasmonate ZIM-domain protein) protein complex, which is composed of the negative regulator JAZ and inositol pentophosphate molecules, has COI1 as a core member of the JAZ receptor complex (Wang, 2015). JA receptor COI 1 can specifically bind to JAZ, causing JAZ to be ubiquitinated and degraded by the proteasome, relieving the inhibition of JAZ in the transcriptional regulation of the JA pathway and thus participating in the transmission of downstream disease resistance genes (Srivastava et al., 2018). In this study, an increase in the expression of multiple JAZ genes regulated transcription and played a dominant role in the JA signal-mediated disease resistance process. Moreover, MYC2 acts as the core transcription factor of the JA signaling pathway, not only exerting positive regulatory effects on defense-related JA responses but also regulating the expression of the plant defense gene PDF1.2 (plant defensin 1.2) (Li et al., 2021), the genes COI1 and MYC2 were confirmed to participate in the JA signaling pathway, with MYC2 showing a significant association with B. squamosa resistance in onion (Allium cepa) (Lee et al., 2022). Additionally, allene oxide synthase (AOS) enzymes mediate the biosynthesis of JA from α-linolenic acid. The AOS-overexpressing line exhibited enhanced JA accumulation, which was associated with increased resistance to B. squamosa (Kim et al., 2025), suggesting that α-linolenic acid metabolism plays a crucial role in plant defense mechanisms. However, the regulatory role of α-linolenic acid in JA biosynthesis is currently inferred primarily from omics data, and its precise molecular mechanisms require further experimental dissection.

Samples were collected at 0 h, 24 h, 48 h, and 96 h after inoculation, with healthy tissues from the interface between diseased and healthy tissues (3 × 5 mm) collected in centrifuge tubes and stored in liquid nitrogen. Metabolomic samples were preserved using the same method, while transcriptomic analysis was performed on three biological replicates per group (n=3), and metabolomic analysis was performed on biological replicates per group (n=4). The repetitive groups and grouping of each sample are listed in Supplementary Table S2 .

Total RNA was extracted from the samples using an RNAprep Pure Plant Kit (Tiangen, Beijing, China). The RNA concentration and purity were measured using a NanoDrop 2000 system (Thermo Fisher Scientific, Wilmington, DE). The integrity of the RNA was assessed using an RNA Nano 6000 detection kit from Agilent Bioanalyzer 2100 Systems (Agilent Technologies, CA, USA). After confirming that the RNA was sufficient for detection, one part was used for sequencing, and the other part was used for qPCR validation.

Paired-end 150 bp (PE 150) mode sequencing was performed using the Illumina NovaSeq 6000 sequencing platform. Gene functions were annotated by comparing with the KOG/COG (Clusters of Orthologous Groups of proteins), KO (KEGG Ortholog database), and GO databases. Data were analyzed using the bioinformatics analysis workflow provided by the Biomarker Cloud Platform BMKCloud (www.biocloud.net).

Genes whose expression levels significantly differ among different samples are called DEGs, and the set of genes obtained from differential expression analysis is known as the set of DEGs, referred to as “A vs. B” (A is the control and B is the experimental treatment). The grouping of DEGs identified in this study is shown in Supplementary Table S3 .

Differential expression gene screening was performed using differential analysis software based on the count value of genes in each sample, and DESeq 2 was used (Love et al., 2014). While detecting DEGs, the screening criteria were set as a differential expression ratio (fold change) ≥ 2 and a false discovery rate (FDR) < 0.01. The greater the absolute value of log₂FC, the more significant the difference between the two groups of samples of genes with smaller FDR values.

The KEGG database resources (Kanehisa et al., 2011) (http://www.kegg.jp), and the KOBAS database (Xie et al., 2011) were used to conduct DEGs enrichment analysis. The clusterProfiler software was used to analyze the enrichment of DEGs in KEGG pathways, and the GO enrichment analysis of DEGs was performed using the clusterProfiler package based on Wallenius non-central hypergeometric distribution (Wu et al., 2021), with P < 0.05 as the threshold. Combined with the differential expression analysis results, the differential gene sequences were aligned to the STRING database (http://stringdb.org/) to obtain the protein-protein interaction (PPI) relationships of the DEGs in B. deweyae, which were then visualized using Cytoscape 3.10.1 (Shannon et al., 2003).

The DEGs used for qPCR validation were selected from the original material that passed detection for reverse transcription. The reverse transcription reaction was performed with the MightyScript First Chain cDNA Synthesis Master Mix (B639251) kit, following instructions provided by Shanghai Biotech Co., Ltd. (MightyScript). After the reaction ended, the samples were stored at –20°C until use. Primers were designed using Primer Premier 6 software with target gene band sets between 150–250 bp, and the gene sequences are shown in Supplementary Table S4 . The internal reference gene of ubiquitin-conjugating enzyme-E2-10 (UBQ-E2-10) and elongation factor 1-alpha isoform (EF-1α2) were selected to correct the relative expression levels of the genes. The Q-PCR mixture was prepared as follows: 20 μL, including 10 μL of Master premix (SGE xcel FastSYBR, B 532955, Sangon Biotech), 0.3 μL each of the upstream and downstream primers, and 2 μL of cDNA, using RNase-Free ddH₂O, was added to 20 μL, and three repetitions were set for each sample. Real-time fluorescence qPCR was performed on a Gentier 96E (Xi’an Tianlong Technology Co., Ltd.).

The one-step real-time fluorescence qPCR protocol was as follows: predenaturation at 95°C for 15 s, reaction at 95°C for 5 s, and reaction at 58°C for 30 s for 38 cycles; the melting curve program was as follows: 95°C for 15 s, followed by 60°C for 1 min. If the Ct value was between 15 and 25, the melting curve was a single peak, and if the Tm value was between 80 and 90, the target gene was considered to be successfully amplified. After obtaining the data using Gentier 96E, the software provided by Gentier 96E was used to preliminarily collate and analyze the experimental data. The 2^-ΔΔCt method was used to obtain the relative expression value, and the calculation formula was as follows: The graphs were drawn with GraphPad Prism 8.

ΔΔ Ct = (Ct of the target gene in the experimental group-Ct of the internal reference gene in the experimental group) - (Ct of the target gene in the control group-Ct of the internal reference gene in the control group).