Transcriptomics of Three Larix Species in Response to Geographically Distinct Bursaphelenchus xylophilus Strains in China

TL;DR

This study explores how three larch species respond at the gene level to different strains of a deadly pine wood nematode in China.

Contribution

The study reveals species-specific and strain-specific transcriptional responses to geographically distinct PWN isolates in larch trees.

Findings

Genes related to oxidative stress and secondary metabolite production were differentially expressed in larch species infected with PWN.

The northern Fushun strain elicited a stronger and more distinct defense response compared to the southern Changde strain.

Transcriptomic analysis identified potential genetic targets for improving resistance in larch trees.

Abstract

Pine wood nematode (PWN, Bursaphelenchus xylophilus) is fatal to the pine trees around the world. Its northward and westward expansion in China endangers Larix spp., yet its molecular response remains understudied. We conducted transcriptomic analysis (RNA-seq) on three economically important larch species (Larix principis-rupprechtii, L. olgensis, and L. kaempferi) infected by geographically distinct PWN isolates (northern Fushun and southern Changde strains) at 1 and 3 days post-inoculation. Comparative RNA-seq analysis of 36 samples revealed that genes such as oxidative stress, and secondary metabolite production were differentially expressed in Larix spp. upon infection by the PWNs. Furthermore, compared to the Changde strain, infection with the Fushun PWN strain can elicit a consistently stronger and more distinct transcriptional defense response across all tested larch species. These results provide insights into the molecular mechanisms of plant defense against PWNs, offering genetic target for resistance breeding and informing the development of targeted control measures against this pathogen.

Figures

Full text

1. Introduction

The pine wood nematode (Bursaphelenchus xylophilus, PWN) is a globally invasive pest that severely threatens coniferous forest ecosystems by causing pine wilt disease (PWD). The disease has caused significant ecological and economic losses globally [1]. Worldwide, PWD is posing a threat to more than 900 million acres of coniferous forest, with a continued spread that could precipitate irreversible ecosystem collapse [2]. According to recent monitoring data, PWD is spreading northward, transcending the traditionally recognized constraints [3,4]. Northern PWN populations exhibit increased cold tolerance and reproductive output compared to southern ones, traits that likely facilitate their northward spread [5,6]. Therefore, elucidating the pathogenic mechanisms of PWD and developing effective management strategies have become imperative.

China, particularly in the northeast, is rich in coniferous forest resources, including Pinus. bungeana, P. tabuliformis, P. koraiensis, P. sylvestris, alongside Picea asperata and Larix spp. [7]. A total of 17 species within the genera Pinus and Larix in China have been identified as natural hosts of PWN [8,9], which renders these forests highly vulnerable to PWD, with profound implications for potential ecosystem devastation and substantial economic costs. Artificial inoculation experiments have further validated this host status, demonstrating that PWN can infect and induce wilt symptoms in seedlings of L. olgensis [10], L. principis-rupprechtii [11], and L. kaempferi [12]. Extensive research, primarily focused on pine species, has begun to elucidate the molecular defense responses triggered by PWN infection. These responses encompass the modulation of phytohormone signaling pathways (e.g., jasmonic acid, salicylic acid, and ethylene) [13], the biosynthesis of defensive secondary metabolites such as terpenoids and phenylpropanoids [14], cell wall reinforcement via lignification [15], and the activation of oxidative stress responses [16]. However, systematic and comparative studies on the defense mechanisms across different host species—particularly among various Larix species, which are economically and ecologically vital in threatened regions—remain notably limited. Our preliminary observations have revealed a clear phenotypic contrast: under controlled inoculation conditions, L. kaempferi exhibits delayed symptom onset compared to the more susceptible L. principis-rupprechtii. Furthermore, PWN populations from distinct geographical origins (e.g., the northern Fushun strain versus the southern Changde strain) appear to differ in virulence.

Therefore, this study employed comparative transcriptomics to analyze and contrast the early defense responses of three key larch species (L. principis-rupprechtii, L. olgensis, and L. kaempferi) to infection by PWN isolates of distinct geographical origins (the northern Fushun strain and the southern Changde strain). By integrating gene annotation, we aimed to characterize species-specific expression patterns of key differentially expressed genes, offering genetic target for resistance breeding and identification the molecular targets against the PWN.

2. Materials and Methods

2.1. Materials

The northern strains (Fushun strain, FS) were collected from infected Pinus koraiensis in Fushun, Liaoning Province, and the southern strains (Changde strain, CD) were isolated from infected P. massoniana in Changde, Hunan Province. Both isolates were routinely cultured and maintained on Botrytis cinerea grown on potato dextrose agar(PDA) medium at 25 °C. Briefly, B. cinerea was first cultured on PDA plates until the mycelia fully covered the surface. Subsequently, surface-sterilized PWNs were inoculated onto the fungal lawn and subcultured under laboratory conditions in a 25 °C incubator to obtain sufficient nematodes for experimentation. All nematode isolates are stored in the Beijing Key Laboratory for Forest Pest Control of Beijing Forestry University, China.

Three-year-old seedlings of L. principis-rupprechtii (Lp), L. olgensis (Lo), and L. kaempferi (Lk) with uniform growth conditions and genetic backgrounds were selected and grown in a greenhouse of Beijing Forestry University, which was set at 25 °C with 80% relative humidity. All seedlings were planted in uniform pots 14.5 cm × 19 cm × 18 cm and domesticated for a few months in order to eliminate the influence from abiotic stress factors.

2.2. Inoculation of Host Plants

Three-year-old seedlings of L. principis-rupprechtii, L. kaempferi, and L. olgensis with uniform growth vigor were selected for the study. The seedlings were cultivated in a soil mixture consisting of 40% (v/v) nutrient substrate and 60% (v/v) native soil.

The bark-grafting inoculation method was employed following established protocols [17]. For each tree species, 40 seedlings were inoculated with one of the two PWN strains, and 10 seedlings were mock-inoculated with sterile water as controls. The inoculation site was selected 1–2 cm above the soil surface. A sterile scalpel was used to create a “T”-shaped wound (approximately 0.3 cm wide and 0.5 cm long) on the stem ensuring the incision penetrated the bark and reached the xylem. A 200 μL pipette tip was then inserted into the wound to contact the exposed xylem surface.

PWNs were collected from 7-day-old B. cinerea cultures using a modified Baermann funnel technique. Nematodes were washed three times with sterile distilled water and concentrated to a final density of 10,000 individuals per 50 μL, as determined by counting under a stereomicroscope. A 50 μL aliquot containing approximately 10,000 nematodes (for treatment groups) or 50 μL of sterile distilled water (for controls groups) was dispensed into the pipette tip. Finally, the tip and the wound area were securely sealed with Parafilm^®^ to prevent leakage and desiccation. Stem samples were collected at 1 and 3 days post-inoculation for transcriptomic analysis.

2.3. Sampling

Stem tissues were collected at 1 and 3 days post-inoculation (dpi) for transcriptomic analysis. For each type of Larch species and treatment combination (FS strain, CD strain, or sterile water control), three young seedings with uniform growth conditions were selected as biological replicates. Sample groups were named according to the convention: species_treatment_time (e.g., Lp_1_FS for L. principis-rupprechtii sampled 1 day after FS strain inoculation). At each time point, a 2 cm stem segment encompassing the inoculation site was excised using a sterile blade, immediately frozen in liquid nitrogen and stored at −80 °C prior to RNA extraction. Three independent biological replicates were performed for each condition, with each replicate consisting of pooled tissue from three individual seedlings to minimize individual variation.

2.4. RNA Extraction

Total RNA was isolated from the collected stem tissues (see Section 2.3) with TRIzol reagent (Invitrogen, Carlsbad, CA, USA). Subsequent quality control involved three measures as follows: RNA concentration was quantified using a NanoDrop 8000 spectrophotometer (NanoDrop, Wilmington, DE, USA); RNA integrity was assessed on an Agilent 2100 Bioanalyzer (Agilent, Waldbronn, Germany); and the absence of degradation or contamination was confirmed by 1% agarose gel electrophoresis. Only samples passing all quality thresholds were used for library construction.

2.5. Preparation of Transcriptome Sequencing Libraries

The total RNA (3 μg) was used for constructing transcriptome sequencing libraries using the Illumina kit (San Diego, CA, USA). Briefly, this was followed by an A-tailing step and the ligation of Illumina sequencing adapters. The libraries were then size-selected with AMPure XP beads to enrich for fragments of the desired size and subsequently amplified by PCR to amplify the adapter-ligated products.

2.6. Transcriptome Analysis and Sequencing Assembly

Following quality confirmation, the qualified RNA samples were sent for sequencing at Shanghai Meiji Biomedical Technology Co., Ltd. (Shanghai, China), in accordance with their standard protocols for library construction and sequencing. Raw paired-end sequencing reads were generated in FASTQ format. De novo transcriptome assembly was performed using Trinity (v2.13.2). The initial assemblies were subsequently refined and filtered using TransRate to remove low-quality and poorly supported transcripts. Assembly completeness was assessed using BUSCO (v5) against the embryophyta_odb10 dataset.

Gene expression levels were quantified using RSEM, which estimates transcript abundances based on the number of reads uniquely mapped to each assembled transcript via Bowtie2. Functional annotation of all assembled unigenes was performed by sequence alignment against six major public databases: NCBI non-redundant protein (NR), Swiss-Prot, Pfam, EggNOG, Gene Ontology (GO), and Kyoto Encyclopedia of Genes and Genomes (KEGG), using BLASTX with an e-value threshold of 1e-5. Statistical summaries of the annotation results from each database were compiled. DEGs were identified using the DESeq2 package with a threshold of false discovery rate FDR < 0.05 & |log2FC| ≥ 1.

2.7. Differential Gene GO, KEGG Analysis

Functionally annotated DEGs were subsequently analyzed for enrichment of specific biological functions and pathways. Gene Ontology (GO) enrichment analysis was performed using the Goatools software, applying Fisher’s exact test to evaluate statistical significance. To correct for multiple hypothesis testing, raw p-values were adjusted using the Benjamini–Hochberg method to control the FDR. GO terms with an adjusted p-value (q-value) < 0.05 were considered significantly enriched. Similarly, Kyoto Encyclopedia of Genes and Genomes (KEGGs) pathway enrichment analysis was conducted to identify metabolic or signaling pathways significantly overrepresented among the DEGs compared to the entire transcriptomic background, using Fisher’s exact test with FDR correction (q-value < 0.05) for statistical assessment.

2.8. qRT-PCR

To validate the reliability of the RNA-seq data, qRT-PCR was performed on selected genes. Gene-specific primers were designed using Primer Premier 3.0 (PREMIER Biosoft) (Supplementary Table S5). cDNA was synthesized from total RNA using the Hifair^®^ II 1st Strand cDNA Synthesis SuperMix (Yeasen), which includes a gDNA removal step. Quantitative PCR was carried out in technical triplicates using 2× HQ SYBR Mix (Zoman) on a Bio-Rad CFX Connect Real-Time PCR system, controlled by CFX Manager 3.0 software. Relative gene expression levels were calculated using the 2^−ΔΔCt^ method [18].

2.9. Statistical Analysis and Drawing

Statistical analyses were performed using SPSS Statistics (version 17.0; SPSS Inc., Chicago, IL, USA) and Origin (version 8.0; OriginLab Corporation, Northampton, MA, USA). Prior to analysis, the assumptions of normality (assessed using the Shapiro–Wilk test) and homogeneity of variances (assessed using Levene’s test) were verified for all datasets. For comparisons of multiple group means (e.g., among different treatment groups within a single species and time point), a one-way analysis of variance (ANOVA) was conducted. When the ANOVA indicated a statistically significant overall effect (p < 0.05), post hoc pairwise comparisons were performed using Duncan’s multiple range test to identify specific group differences. Significance levels are denoted as * p < 0.05, ** p < 0.01, *** p < 0.001, and **** p < 0.0001. All graphical representations of the data were generated using Origin.

3. Results

3.1. RNA Sequencing Quality

A total of 36 cDNA libraries were constructed and sequenced. High-quality clean data were obtained for all samples, with an average output of 6.63 Gb per library. The minimum Q20 and Q30 scores were 98.47% and 95.28%, respectively, confirming that the sequencing quality met the required standards (Table S1). All clean reads were pooled and subjected to de novo assembly using Trinity, yielding 228,721 unigenes and 404,237 transcripts. The assembly demonstrated good continuity, with an N50 length of 1622 bp for unigenes. To assess the completeness of the transcriptome assembly, BUSCO analysis was performed. The results showed a complete BUSCO (C) score of 78.0%, consisting of 70.6% single-copy and 7.4% duplicated genes, indicating a relatively reliable assembly suitable for downstream analysis (Table S2).

Subsequently, clean reads from each sample were mapped back to the assembled transcriptome for gene expression quantification. The mapping rates across all 36 samples ranged from 74.13% to 83.47%, demonstrating effective utilization of the sequencing data (Table S3). For functional annotation, the assembled unigenes were searched against major public databases. In total, 148,921 unigenes (65.11% of all unigenes) were annotated in at least one database. Specifically, 114,307 (49.98%), 112,159 (49.04%), and 109,902 (48.05%) unigenes were annotated in the NR, Swiss-Prot, and Pfam databases, respectively. Furthermore, 83,947 (36.70%) and 70,590 (30.86%) unigenes were assigned GO terms and KEGG pathways, respectively, providing a substantial basis for functional interpretation (Table S4).

3.2. RNA Sequencing (RNA-Seq) Profiles and Identification of Differentially Expressed Genes (DEGs)

Analysis of DEGs among different groups can help reveal the molecular mechanism by which Larix resists PWD. To elucidate the inherent resistance mechanisms of different Larix spp., all samples were divided into the following two major categories: a control category and a treatment category. The control category comprised healthy samples (since its injected substance was sterile water, it was assumed to be in a healthy state in the subsequent samples), including the comparisons Lk_1_CK_vs_Lp_1_CK, Lk_3_CK_vs_Lp_3_CK, Lo_1_CK_vs_Lp_1_CK and Lo_3_CK_vs_Lp_3_CK (Figure 1 and Figure 2). The treatment category consisted of nematode-inoculated samples, containing the comparisons Lk_1_FS_vs_Lk_1_CD, Lk_3_FS_vs_Lk_3_CD, Lo_1_FS_vs_Lo_1_CD, Lo_3_FS_vs_Lo_3_CD, Lp_1_FS_vs_Lp_1_CD and Lp_3_FS_vs_Lp_3_CD. A total of 55,292 DEGs were detected in the Lk_1_CK_vs_Lp_1_CK comparison, accounting for approximately 24.2% of all expressed unigenes, with 29,871 upregulated and 25,421 downregulated. This comparison yielded the highest number of DEGs, indicating substantial basal transcriptional divergence between L. kaempferi and L. principis-rupprechtii. The overall distribution of DEGs across all comparisons is presented in Figure 1A.

To investigate the similarities and differences in gene expression across transcriptomic, Venn diagrams were used to visualize the distribution of DEGs among the different groups. The healthy control and different treatment categories contained 12 and 4 uniquely upregulated DEGs, and 10 and 5 uniquely downregulated DEGs, respectively (Figure 3). Functional annotation based on GO and Swiss-Prot indicated that core responsive DEGs across comparisons were significantly associated with oxidative stress, signal transduction, biosynthesis of secondary metabolites, and energy metabolism. Notably, among the DEGs associated with secondary metabolism, several were identified as encoding key enzymes such as terpene synthases and cytochrome P450s. The coordinated regulation of these pathways underscores a conserved transcriptional reprogramming involved in defense responses in Larix species following infection with the PWN.

3.3. DEGs Significantly Enriched in Healthy Category

To elucidate the inherent transcriptional differences among the three Larix species, we analyzed DEGs between healthy control samples. Distinct expression profiles, characterized by substantial fold changes, revealed clear species-specific patterns (Table 1). In the Lk_1_CK_vs_Lp_1_CK-up group, five oxidative stress-related genes were identified among the ten genes. These included photosynthesis-related components such as chlorophyll a-b binding protein and multiple peroxidases (e.g., Peroxidase 9). Conversely, genes for secondary metabolism and defense were strongly suppressed, such as gamma-humulene synthase and allene oxide synthase. Furthermore, in the Lk_3_CK_vs_Lp_3_CK-up set encoding cysteine proteinases (FC up to 7285.05) and an EF-hand domain protein (FC = 16,053.27) became prominent.

Within the Lo_1_CK_vs_Lp_1_CK-up and Lo_1_CK_vs_Lp_1_CK-down sets, genes involved in redox and transport were highly induced, such as probable inositol oxygenase (FC = 7411.48) and a ZEB2-regulated ABC transporter (FC = 4469.36). This was accompanied by a marked suppression of many pathogenesis-related genes, exemplified by endochitinase 4 (FC as low as 5.02 × 10^−6^). The transcriptional divergence intensified (Table 1).

Collectively, these comparisons highlight substantial basal transcriptional differences among the species. L. kaempferi exhibited higher expression of genes related to photosynthesis and specific peroxidases. L. olgensis showed elevated expression of genes involved in redox metabolism and energy production, some with exceptionally high fold changes.

3.4. DEGs Significantly Enriched in Different Treatment Category

Building on the comparative analysis of the healthy control groups aimed at revealing host resistance mechanisms, we next analyzed the treatment groups to investigate how PWN strains from different geographic origins affect different Larix species. In L. kaempferi, early response (Lk_1_FS_vs_Lk_1_CD) was characterized by the specific upregulation of genes such as glutamate decarboxylase 4 (FC = 3772.44), leucoanthocyanidin reductase (FC = 3307.12), and glucokinase (FC = 2650.21). Conversely, several genes associated with oxidative stress response, including catalase (FC = 7.80 × 10^−4^) and peroxidase 54 (FC = 1.38 × 10^−3^), were suppressed in response to FS inoculation compared to the CD strain. At a later stage (Lk_3_FS_vs_Lk_3_CD), genes like peroxiredoxin Pen c 3 (FC = 6009.96) were markedly induced, whereas multiple genes encoding ribosomal proteins (e.g., 50S ribosomal protein L4, 40S ribosomal protein S8) and tubulin subunits were downregulated (Table 2). Collectively, these results indicate that in L. kaempferi, the FS strain modulates specific metabolic and antioxidant pathways while concurrently repressing fundamental cellular biosynthesis processes, compared to the CD strain.

In the Lo_1_FS_vs_Lo_1_CD-up set, defense-related genes were enriched, such as multiple peroxidases (e.g., Peroxidase 53, Peroxidase N1), an LRR receptor-like serine/threonine-protein kinase (EFR), an endochitinase (Endochitinase 4), and phenylalanine ammonia-lyase. Conversely, the Lo_1_FS_vs_Lo_1_CD-down set showed suppression of genes involved in primary metabolism, including a putative alpha,alpha-trehalose-phosphate synthase (FC = 5.01 × 10^−4^), catechol O-methyltransferase A (FC = 8.55 × 10^−4^), and peroxiredoxin PRX1 (FC = 8.74 × 10^−4^). A distinct pattern was observed at a later time point. The Lo_3_FS_vs_Lo_3_CD-up set was dominated by the marked upregulation of transcription initiation factor IIB (FC = 7213.24). Simultaneously, the Lo_3_FS_vs_Lo_3_CD-down set featured pronounced downregulation of secondary metabolism-associated genes, such as xanthohumol 4-O-methyltransferase (FC = 2.01 × 10^−4^), catechol O-methyltransferase (FC = 7.96 × 10^−4^), and alpha terpineol synthase (FC = 1.58 × 10^−3^) (Table 3). Collectively, these patterns indicate that the response of L. olgensis to the FS strain involves an early (Day 1) induction of defense-related genes (e.g., peroxidases, hydrolases, and phenylpropanoid pathway components), followed by a later (Day 3) shift characterized by the suppression of specific secondary metabolic pathways and the upregulation of key transcriptional initiation factors.

The Lp_1_FS_vs_Lp_1_CD-up gene set was characterized by the dominance of genes associated with oxidative stress response, most notably catalase-peroxidase (FC = 69,220.11) and superoxide dismutase [Cu-Zn] (FC = 40,867.03). Lp_1_FS_vs_Lp_1_CD-down set comprised genes involved in secondary metabolism (e.g., alpha terpineol synthase, FC = 2.38 × 10^−3^) and defense (e.g., endochitinase 4, FC = 3.19 × 10^−3^). At a later time point (Lp_3_FS_vs_Lp_3_CD), both up- and down-regulated gene sets contained genes related to oxidative stress and hydrolysis. For instance, catalase-1 (FC = 640.35) was upregulated, while several peroxidases (e.g., Peroxidase 4, Peroxidase 53) were downregulated (FC ranging from 1.27 × 10^−3^ to 1.52 × 10^−3^). Notably, key oxidative stress-related genes remained differentially expressed on day 3, albeit with lower fold changes compared to day 1 (Table 4). These expression dynamics in L. principis-rupprechtii indicate an early and potent induction of oxidative stress-responsive genes upon inoculation with the FS strain, which occurs alongside the modulation of other metabolic and defense-related pathways.

3.5. GO Analysis

GO enrichment analysis provided a broad functional overview of the biological processes affected by PWN inoculation. The top 20 most significantly enriched GO terms across all comparisons are summarized in Figure 3.

Analyses under healthy control comparisons revealed pronounced species-specificity. The transcriptional differences between L. kaempferi and L. principis-rupprechtii at both time points were characterized by the enrichment of broad metabolic categories, primarily oxidoreductase activity (Molecular Function, MF) and biosynthetic process (biological process, BP) (Figure 4A,B). In contrast, differences between L. olgensis and L. principis-rupprechtii at day 1 were strongly associated with translation and ribosomal structure/function (Figure 4C), a profile that shifted by day 3 to more closely resemble that of the Lk vs. Lp comparison (Figuure 4D).

The three Larix species exhibited distinct response patterns to PWN inoculation from different geographical sources. In both L. kaempferi and L. principis-rupprechtii, the most enriched high-level GO terms were similar, most notably oxidoreductase activity (MF) and organic acid metabolic process (BP) (Figure 4E,F,I,J). Conversely, L. olgensis exhibited a clear temporal progression: the early response (day 1) was characterized by terms related to membrane intrinsic components (Cellular Component, CC), while the later response (day 3) shifted towards processes such as organic substance metabolic process (BP) (Figure 4G,H).

3.6. KEGG Analysis

KEGG pathway enrichment analysis revealed that inoculation with the FS strain triggered extensive metabolic reprogramming across the three Larix species. At one day 1, pathways including oxidative phosphorylation, fatty acid degradation, and pyruvate metabolism were significantly enriched in multiple comparisons (Figure 5). The consistent enrichment of these core energy and catabolic pathways implies a rapid and coordinated shift in host energy metabolism following pathogen challenge. Furthermore, variation in the enrichment profiles among the three Larix species was evident, pointing to distinct early transcriptional responses to the same PWN strain.

3.7. Verification of Gene Expression by qRT-PCR

To independently verify the transcriptome findings, candidate DEGs associated with oxidative phosphorylation, pyruvate metabolism, and fatty acid degradation were subjected to qRT-PCR. A correlation analysis between the qRT-PCR and theRNA-seq results showed that nine genes displayed consistent expression patterns, thereby validating the accuracy of our sequencing data (Figure 6 and Figure 7). Compared to Lp_1_CK, most genes in three Larix spp. exhibited a marked upregulation in expression on the first day following nematode infection, particularly in L. kaempferi, such as cytochrome c (TRINITY_DN118629_c0_g1), ATP synthase (TRINITY_DN5077_c0_g2), inorganic pyrophosphatase (TRINITY_DN40521_c0_g1), phosphoenolpyruvate carboxykinase (TRINITY_DN181650_c0_g1), alcohol dehydrogenase (TRINITY_DN134703_c0_g1, TRINITY_DN4004_c0_g3, and TRINITY_DN8551_c3_g1), phosphoenolpyruvate carboxylase (TRINITY_DN4239_c0_g1), and thiolase (TRINITY_DN48265_c0_g1). However, the expression of these genes declined pronouncedly by the third day. Overall, most genes in L. kaempferi showed higher expression levels upon infection with the FS nematode strain compared to the CD strain, such as oxidative stress, secondary metabolites synthesis, and detoxification, which indicated its higher resistance.

4. Discussion

A comparative transcriptomic analysis was performed to investigate the species-specific defense mechanisms of three Larix spp. against PWN infection as well as the pathogenicity divergence between geographically distinct PWN populations in China. The findings revealed that following the PWN challenge, all three Larix species exhibited marked differential expression in genes broadly associated with oxidative stress and metabolism. Notably, the transcriptional changes induced by the FS strains were more pronounced and distinct compared to those induced by the CD strains across all hosts.

Upon PWN infection, nematodes migrate through and feed on plant resin canals, disrupting water conductivity in stems and subsequently reducing transpiration and photosynthesis in needles, ultimately leading to systemic oxidative damage. The induction of oxidative stress represents a conserved plant defense mechanism, in which ROS function as signaling molecules to activate downstream defense responses. Consistent with studies in Pinus species challenged by PWN [16,19,20,21], genes encoding key antioxidant enzymes such as peroxidases and superoxide dismutases were strongly differentially expressed in infected Larix seedlings. These antioxidant enzymes work in concert to eliminate excess ROS and mitigate oxidative damage. However, our transcriptomic analysis revealed marked interspecific and temporal differences in the oxidative stress responses among the three larch species. In L. principis-rupprechtii, inoculation with the FS strain triggered an extreme upregulation of oxidative stress-related genes at 1 dpi, most notably catalase-peroxidase (FC = 69,220.11) and superoxide dismutase [Cu-Zn] (FC = 40,867.03) (Table 4). By 3 dpi, however, the expression levels of these genes had decreased substantially. In contrast, L. kaempferi exhibited a more modulated early response: some typical antioxidant genes, including catalase and peroxidase 54, showed lower expression by FS strains compared to CD strains at 1 dpi (Table 2). Meanwhile, L. olgensis displayed a delayed but sustained response, with strong upregulation of peroxidase 53 (FC = 3931.93) and other oxidative stress-related genes peaking at 3 dpi (Table 3). These findings demonstrate that while all three larch species mount oxidative stress responses to PWN infection, the timing and magnitude of these responses differ substantially, suggesting a species-specific regulatory pattern.

Secondary metabolites of coniferous trees are important components of their defense responses which may participate in host recognition of the nematode and impede its spread and reproduction. Our transcriptomic analysis revealed species-specific patterns in the regulation of secondary metabolism-related genes following PWN infection.

In L. principis-rupprechtii, genes involved in the synthesis of potential defense compounds were significantly downregulated at 1 dpi, including members of the terpene synthase family (e.g., alpha-terpineol synthase, FC = 2.38 × 10^−3^) and chitinase class I genes (e.g., endochitinase 4, FC = 3.19 × 10^−3^) (Table 4). In contrast, L. kaempferi exhibited strong induction of genes involved in specific metabolic pathways at 1 dpi, such as glutamate decarboxylase 4 (FC = 3772.44) and leucoanthocyanidin reductase (FC = 3307.12) (Table 2). Meanwhile, L. olgensis showed early upregulation of phenylalanine ammonia-lyase (PAL, FC = 528.66) at 1 dpi (Table 3), accompanied by GO enrichment for membrane-associated components. By 3 dpi, the response in L. olgensis shifted toward broader metabolic reprogramming, characterized by strong upregulation of genes such as transcription initiation factor IIB and enrichment for diverse metabolic processes (Table 3, Figure 4G,H). The crucial role of terpenoid and phenylpropanoid pathways in the defense of pine trees against PWN has been well established. Upregulation of genes involved in the terpene backbone biosynthesis pathway following PWN inoculation has been reported in several pine species, including P. massoniana [16], P. pinaster [13,21,22], P. pinea [21], and P. yunnanensis [22]. Moreover, multiple terpene synthase genes were more highly expressed in resistant genotypes of P. massoniana compared to susceptible ones after PWN infection [13]. For instance, Trindade et al. [23] reported the differential upregulation of α-pinene synthase in resistant P. pinaster and of limonene cyclase in susceptible P. densiflora, demonstrating that species- or genotype-specific expression of terpenoid synthesis genes is associated with host resistance. Furthermore, in resistant P. massoniana genotypes, genes such as α-pinene synthase (PmTPS4) and longifolene synthase (PmTPS21) encode enzymes capable of synthesizing a repertoire of mono- and sesquiterpenes [14], many of which have been shown to possess repellent or nematicidal activity against PWN, thereby contributing to pine resistance.

The adaptation of PWN pathogenicity to different pine species is one of the key factors affecting the spread of the disease. In this study, in addition, the FS strains can trigger a more significant response mechanism in different Larix species. A large number of domestic and foreign studies have shown that there are differences in the pathogenicity of PWN from different geographic sources [24,25], and PWN isolates from different host sources exhibit different levels of virulence to different tree species. For example, a comparison study involving eight southern strains and eight northern strains found that the pathogenicity of northern PWN on P. tabuliformis was significantly higher than that of southern strains, with further characterized by larger average body size and higher reproductive capacity [26]. In addition, comparative analyses identified Liaoning strains as significantly more virulent than those from Nanjing or Chongqing, based on assessments of seedling kill, migration, and reproduction, based on assessments of seedling kill, migration, and reproduction [5]. It is important to note that this study focused on the host transcriptome; therefore, the specific nematode-derived effectors or mechanisms that drive these differential host responses remain to be elucidated and represent a key direction for future research.

5. Conclusions

Our comparative transcriptome analysis revealed distinct species-specific defense strategies among larch species in response to infection by the PWN, particularly involving genes related to secondary metabolite synthesis. Furthermore, infection with the northern-origin (Fushun, China) PWN strain elicited a consistently stronger and more distinct transcriptional defense response across all tested larch species compared to the southern-origin (Changde, China) strain. Collectively, our study provides molecular evidence for the divergent defense mechanisms that likely contribute to the differential susceptibility observed among these larch species. The finding that northern PWN isolates trigger a heightened host transcriptomic response also offers a novel perspective for further investigating the interaction dynamics between geographically distinct nematode populations and their conifer hosts.