Disruption of Energy Metabolism, Locomotion, and Stress Tolerance by BxPDHB Silencing in Bursaphelenchus xylophilus: Implications for RNAi-Based Control of Pine Wilt Disease

1. Introduction

Pine wilt disease, caused by the pine wood nematode (Bursaphelenchus xylophilus), is one of the most destructive forest diseases worldwide and is listed by the Food and Agriculture Organization of the United Nations and by national plant-protection authorities as an internationally quarantined pest of major concern [1]. Since its transboundary and intercontinental spread, the disease has killed hundreds of millions of conifers across Eurasia and parts of the Americas, posing serious threats to the health of coniferous forest ecosystems, to forest carbon-sink stocks, and to timber trade, and causing substantial losses in direct economic value and ecosystem services [2]. Two key biological characteristics underlie the continued global dispersal and epidemic establishment of pine wilt disease. On the one hand, strong stress-adaptive capacities—including resistance to desiccation, low temperature, and oxidative stress—enable the nematode to survive and establish infection sites in the nutrient-poor, toxic secondary-metabolite environment of the host xylem [3,4]. On the other hand, its locomotory and migratory capacities—active within-tree dispersal as well as long-distance transport via vector insects (cerambycid beetles, especially Monochamus spp.)—endow it with the potential to rapidly colonize new host space and achieve transregional spread [5,6]. The coordinated interplay of these two capacities is a key intrinsic factor underlying the rapid onset and strong dispersal characteristic of pine wilt disease.

Previous studies have delineated relatively comprehensive molecular maps of these two dimensions independently. With respect to stress adaptation, B. xylophilus responds to low temperature, oxidative stress, and host terpenoid pressure through transcriptional reprogramming, upregulating cytochrome P450 family genes for detoxification and cold tolerance [7,8], and UDP-glycosyltransferases (UGTs) and glutathione S-transferases (GSTs) for xenobiotic clearance and antioxidation [9,10]; sulfate permease (SULP) family members participate in sulfur assimilation and cellular homeostasis [11], and autophagy-related genes (e.g., BxATG1, BxATG8) mediate survival regulation under adverse conditions [12]. With respect to locomotion and behavior, the sodium-dependent neurotransmitter transporter (SNF) family maintains motor coordination by recycling neurotransmitters [13], and G-protein-coupled receptors (e.g., BxGPCR17454) influence behavioral decisions through signal transduction [14]. Overall, these studies have largely focused on the effects of individual genes on specific behavioral or stress phenotypes (e.g., chemotaxis, locomotion speed, survival rate) and thus exemplify a “single gene–single phenotype” paradigm. However, the upstream coupling nodes that simultaneously govern the two major ecological traits of “high locomotory capacity” and “strong stress resistance” remain unidentified at the gene level, and cross-dimensional coupling studies are conspicuously lacking.

Importantly, locomotion and stress tolerance share a common underlying platform at the level of energy metabolism: locomotory behavior depends on neurotransmission and muscle contraction, which continuously consume ATP, whereas stress responses such as antioxidation and detoxification likewise depend on energy supply and reducing equivalents. Mitochondrial energy metabolism, therefore, constitutes the most plausible upstream hub linking locomotory output and stress defense, suggesting that a regulatory factor located at the core of energy metabolism is most likely to establish coupling between locomotion and stress tolerance. On this basis, we focused on a highly specialized dispersal window in the life history of B. xylophilus—the JIII dispersal-stage juvenile. This stage is the critical period during which nematodes emerge from diseased wood and actively migrate toward, and aggregate at, the beetle pupal chambers [15,16]; it is both a stress-tolerant phase that must cope with desiccation, low temperature, and host defenses and a core link in achieving spatial dispersal and transmission. The JIII stage is therefore the only developmental stage that simultaneously requires high locomotory output and high stress resistance, and it represents a natural biological model for dissecting the “locomotion–stress resistance” coupling. Accordingly, we integrated bioinformatics approaches—weighted gene co-expression network analysis (WGCNA) and gene set enrichment analysis (GSEA)—together with random forest machine learning to systematically identify metabolic hub genes specifically and highly expressed at the JIII stage. Candidate targets were then functionally validated by double-stranded RNA (dsRNA)-mediated RNA interference (RNAi), and their regulatory effects on energy conversion, locomotory frequency, substrate homeostasis, and redox defense were assessed across multiple dimensions. The aim of this study was to identify a candidate hub gene involved in the coordinated regulation of locomotory capacity and stress tolerance, thereby providing a candidate molecular target and a theoretical basis for developing novel green control strategies centered on RNAi from the perspective of interfering with the energy-metabolism pathway.

2. Materials and Methods

2.1. Transcriptomic Data Acquisition and Bioinformatics Analyses

Gene expression data for different developmental stages of B. xylophilus (egg, J2, J3, J4, JIII, and adult) were obtained from the NCBI Gene Expression Omnibus (GEO) database under accession number GSE242182 [17]. Gene set enrichment analysis (GSEA) [19] was performed using the clusterProfiler package [18] to compare the transcriptional profiles of the JIII stage with those of the propagative J2 stage, so as to identify candidate classical metabolic pathways that were significantly enriched and specifically highly expressed at the JIII stage. A weighted gene co-expression network was constructed using the WGCNA package in R [20]. The raw dataset comprised 15,884 genes across 18 samples; after filtering out genes with low expression variability (standard deviation ≤ 0.5), 5,843 genes were retained for subsequent analysis. The pickSoftThreshold function was used to evaluate scale-free fit and mean connectivity over a soft-threshold power range of 1–30, and a power of 30 was ultimately selected to construct the adjacency matrix. Based on the selected power, hierarchical gene clustering and module identification were performed using default parameters (minModuleSize = 30, mergeCutHeight = 0.25), yielding 17 gene modules (the grey module comprised genes that could not be assigned to any module). Using each developmental stage as a trait, Pearson correlations between module eigengenes and traits, together with their p-values, were computed for module–trait association analysis. To explore the functional pathways of the module genes, KEGG pathway enrichment analysis was performed on the target gene sets using TBtools, with the genome-wide set of 6,280 genes obtained by KEGG annotation serving as the background gene set. Fisher’s exact test was used to assess enrichment significance, and an adjusted p-value (Benjamini–Hochberg correction) < 0.05 was set as the threshold for significant enrichment. Significantly enriched KEGG pathways were retained for subsequent functional analysis.

Gene sequences within the specified module were extracted, and a protein–protein interaction (PPI) network was constructed based on the STRING database (Caenorhabditis elegans, default confidence) [21]. Cytoscape software [22] was used to calculate node degree centrality and identify hub genes, and the top 15 hub genes were selected. The core genes identified from the PPI network were then intersected with the candidate core genes enriched by GSEA. Subsequently, a random forest classification model was constructed using the randomForest (v4.7-1.2) package in R (v4.3.2) [23], with developmental stages (adult, DL3, egg, L2, L3, L4) as class labels. The caret package was used for five-fold cross-validation to assess model performance, with parameters set as follows: ntree = 500, and mtry tuned within {2, 21, 40}; the optimal parameters were selected by maximizing accuracy (mtry = 40; model accuracy = 0.853, Kappa = 0.814). Based on the optimal model, variable importance measures were calculated for each gene and used to identify key genes.

2.2. B. xylophilus Propagation and Target Gene Silencing

The B. xylophilus strain and the fungal food source (Pestalotiopsis sp.) used in this study were provided by the Key Laboratory of Major Pest Control for Ecological Public-Welfare Forests, Fujian Province. Nematodes were inoculated onto potato dextrose agar plates overgrown with Pestalotiopsis sp. and cultured at 25 °C in the dark for 7 days. Nematodes were harvested using a modified Baermann funnel method [24]. The collected nematode suspension was centrifuged at 3000×g for 2 min, and the supernatant was discarded to collect the nematode pellet for subsequent experiments.

For synchronously cultured B. xylophilus, 10,000 nematodes per group were separately soaked in nuclease-free water, in a solution of 1000 ng/μL EGFP double-stranded RNA (dsRNA), and in an equal concentration of target-gene dsRNA, serving as the blank control, negative control, and treatment group, respectively. Following the RNAi soaking method established by Maeda et al. (2001) and widely applied to nematode gene-function studies [25], each sample was incubated at 25 °C in the dark for 48 h. After treatment, nematodes were washed several times with nuclease-free water to remove residual dsRNA.

Specific dsRNA was designed against the region spanning nucleotides 143–643 (501 bp; GC content 51%) of the coding sequence of the target gene BxPDHB (BXYJ_LOCUS10935). During design, target-sequence features such as secondary structure, repetitive segments, and palindromic sequences were considered, and dsRNAEngineer was used to predict potential off-target sites (29 predicted targets, of which 15 were high-scoring); the segment with lower off-target risk was selected for subsequent RNAi experiments. To exclude cross-interference with other genes in the B. xylophilus genome, the target sequence was subjected to whole-genome homologous-sequence screening by BLAST, confirming its specificity to BxPDHB with no significant homologous off-target sequences. In addition, dsRNA of the non-homologous EGFP sequence was used as a negative control to evaluate off-target effects and non-specific interference.

2.3. Determination of Gene Silencing Efficiency

Total RNA was extracted from nematodes in each group using the RC112 total RNA extraction kit (Vazyme Biotech, Nanjing, China) according to the manufacturer’s instructions. After confirming that the concentration, integrity, and purity of the extracted total RNA were acceptable, first-strand cDNA was synthesized using the Vazyme first-strand cDNA synthesis kit (Cat. No. R312). Real-time quantitative PCR (RT-qPCR) was performed on a quantitative nucleic acid amplification instrument at the Forestry Research Center, Fujian Agriculture and Forestry University, using the Vazyme Q312 fluorescent quantitative kit. The reaction program was set according to the manufacturer’s instructions: 95 °C pre-denaturation for 30 s, followed by 40 cycles of 95 °C for 10 s and 60 °C for 30 s (annealing and extension). After amplification, a melting-curve analysis was performed (95 °C for 15 s, 60 °C for 15 s, and 95 °C for 15 s). Three independent biological replicates were set, each with three technical replicates. The relative expression of the target gene was analyzed using the 2^(−ΔΔCt) relative quantification method of Livak and Schmittgen [26], with endogenous β-actin used as the reference gene. The primers used for RT-qPCR were synthesized by Sangon Biotech (Shanghai, China); primer sequences are listed in Table S1.

2.4. Survival Rate and Head-Swing Frequency Assays

The effect of target-gene silencing on the survival rate of B. xylophilus was assessed using a soaking assay. Nematodes cultured on Pestalotiopsis sp. plates were collected by the Baermann funnel method and transferred to 24-well cell culture plates. One thousand active nematodes were added to each well and separately soaked in nuclease-free water, in a 1000 ng/μL EGFP dsRNA solution, and in an equal concentration of target-gene dsRNA solution. The plates were incubated at 25 °C in a constant-temperature incubator in the dark for 48 h. After incubation, nematodes were washed five times with nuclease-free water and centrifuged at 3000×g for 2 min. Each treatment comprised three independent biological replicates, each with three technical measurements (nine values in total). Survival status was observed and survival rate calculated 48 h after interference treatment.

Head-swing behavior was defined as the worm bending into an S shape and swinging its head once to the left or right before returning to the initial position [27]. One hundred active nematodes per group were collected and transferred to a 24-well plate, and a nematode suspension was prepared with nuclease-free water at a water depth of 3 mm. For the control groups, 0.3 mL nuclease-free water and 0.3 mL EGFP dsRNA dilution were added, whereas the treatment group received an equal volume of the target-gene dsRNA solution. The plates were incubated at 25 °C in a constant-temperature incubator in the dark for 48 h, after which the nematodes were rinsed with nuclease-free water. At 48 h after treatment, under a stereomicroscope, 10 actively moving nematodes were randomly selected from each group, and the number of head swings within 1 min was recorded. Each treatment comprised three independent biological replicates, each with three technical measurements.

2.5. Measurement of Intracellular Reactive Oxygen Species (ROS)

Intracellular ROS production was measured using the cell-permeable non-fluorescent probe 2′,7′-dichlorodihydrofluorescein diacetate (DCFH-DA), following the protocol of Chávez et al. (2014) as applied to nematode ROS detection [28]. One thousand active nematodes were collected and resuspended in 500 μL nuclease-free water in a 1.5 mL centrifuge tube; the target-gene dsRNA working solution was added to a final working concentration of 1000 ng/μL. The samples were incubated at 25 °C in the dark for 48 h and then centrifugally washed five times with nuclease-free water. The treatment and control groups were incubated with the probe reaction solution at a final concentration of 50 μmol/L at 25 °C in the dark for 2 h. After the reaction, nematodes were washed five times with M9 buffer to thoroughly remove background free fluorescence. Worms were immobilized using a 60 μg/mL levamisole solution. Fluorescence signal intensity in the nematodes was observed using a Zeiss LSM 880 ultra-high-resolution laser scanning confocal microscope. Each treatment comprised three independent biological replicates, each with three technical measurements (30 worms were randomly examined per group), and the mean fluorescence intensity of the worm body was calculated using ZEN analysis software.

2.6. Measurement of Intracellular ATP Content

ATP content was quantified using the ATP detection kit (Cat. No. K2315; ApexBio, USA) according to the manufacturer’s instructions. Before measurement, B. xylophilus was collected by the Baermann funnel method and washed three times with nuclease-free water. The treatment group received the target-gene dsRNA solution at a final concentration of 1000 ng/μL, whereas the control groups received an equal volume of nuclease-free water or GFP-dsRNA solution. After incubation at 25 °C in the dark for 48 h, the samples were centrifuged at 3000×g for 2 min and washed five times with nuclease-free water. An appropriate volume of pre-cooled ATP detection buffer was added to the collected nematode pellet for ice-bath homogenization. The homogenate was then centrifuged at 13,000×g for 5 min at 4 °C, and the supernatant was transferred to a 10 kDa ultrafiltration tube and centrifuged at 13,000×g for 10 min at 4 °C for deproteinization; the lower filtrate was collected as the test sample. The kit components were equilibrated to room temperature, and a series of standard working solutions was prepared by gradient dilution of the standard. In the standard wells, different concentrations of standard solution were added and brought to a volume of 50 μL with detection buffer; in the sample wells, the test sample was added and brought to 50 μL with detection buffer; blank wells contained no sample. Subsequently, 50 μL of freshly prepared fluorescent detection reaction solution was added to each well and incubated at room temperature in the dark for 30 min. Finally, the absorbance of each well at 570 nm was measured with a multifunctional microplate reader, and the absolute ATP content of each sample was calculated from the standard curve. Each treatment comprised three biological replicates, each with three technical replicates.

2.7. Measurement of Pyruvate Content

Pyruvate content was quantified using the pyruvate quantitative detection kit (Cat. No. AKAC002M; Beijing Boxbio Biotechnology Co., Ltd., Beijing, China) according to the manufacturer’s instructions. An appropriate volume of extraction solution was added to the collected nematode pellet for homogenization; after ice-bath treatment, the homogenate was incubated at room temperature for 30 min and then centrifuged at 8000×g for 10 min at room temperature; the supernatant was collected as the test sample. Measurement, standard, and blank wells were set up in the reaction microplate: the measurement wells received 75 μL of test supernatant, the standard wells received 75 μL of serially diluted standard pyruvate solutions, and the blank wells received 75 μL distilled water. To each well, 25 μL of Reagent 1 was first added; after thorough vortex mixing, the plate was incubated at room temperature for 2 min. Then 125 μL of Reagent 2 was added, mixed thoroughly, and incubated at room temperature for 15 min. After the reaction, the absorbance at 520 nm was measured with a microplate reader, and the absorbance difference between each sample and the blank group was calculated. The measured absorbance values were substituted into the regression equation to calculate the absolute pyruvate content of each sample. Each treatment comprised three independent biological replicates, each with three technical measurements.

3. Results

3.1. Identification of JIII-Stage-Specific Co-Expression Modules by WGCNA and Topological Hub-Gene Mining

To systematically identify key upstream master genes governing the survival and locomotion of the third-stage dispersal juveniles (JIII) of B. xylophilus, we first performed weighted gene co-expression network analysis (WGCNA) on the transcriptomic data of different developmental stages. The results showed that, among the identified gene co-expression modules, the purple module exhibited the highest correlation with the JIII stage (r = 0.97, P < 0.01), and the genes in this module showed a highly significant stage-specific high expression in the third-stage dispersal juveniles, indicating that this module’s genes are deeply involved in the development and physiological regulation of the dispersal juveniles (Figure 1a). KEGG pathway enrichment analysis of the 144 genes in this module revealed that these genes were mainly enriched in three categories of pathways—ribosome translation, redox and detoxification, and mitochondrial energy metabolism (Figure 1b and Table S2). Based on the pairwise interactions among the 144 genes in the co-expression module, we constructed a protein–protein interaction (PPI) network. Topological analysis showed that the network spontaneously clustered into three physically relatively independent yet functionally tightly coupled functional clusters, corresponding to the three pathway categories above; the mitochondrial energy-metabolism cluster was located at the topological core of the network, serving as a hub that connects the ribosome-translation and redox-and-detoxification modules (Figure 1c). To quantify the criticality of each node in information transmission across the network, we calculated the betweenness centrality of each node and, in descending order, selected the top 15 node genes as core hub genes maintaining the scaffold stability of the PPI network (Figure 1c, highlighted in yellow).

3.2. Multi-Dimensional Cross-Validation with GSEA and Random Forest and Identification of the Core Effector Target

To independently validate the above results from the perspective of dynamic developmental processes and further narrow the target range, we performed gene set enrichment analysis (GSEA) to assess the coordinated changes of gene sets during the transition from J2 to JIII metamorphosis. The analysis revealed that, during the J2-to-JIII developmental transition, genes in ten KEGG pathways were coordinately upregulated; notably, these ten pathways overlapped highly with the pathways enriched by WGCNA modules and, through convergent routes, pointed to the same three core biological themes—ribosome translation, redox and detoxification, and mitochondrial energy metabolism (Figure 2a and Table S3). The core gene sets from the same pathways identified by GSEA were intersected with the 15 hub genes identified in the previous step, yielding ten candidate genes. To evaluate the differential contribution of these ten genes to JIII phenotype maintenance, a random forest machine-learning algorithm was applied, using developmental stage as the dependent variable, to rank the feature importance of the candidate genes. The results showed that three genes—BXYJ_LOCUS8868, BXYJ_LOCUS14538, and BXYJ_LOCUS10935—had feature-importance scores greater than 80 and ranked among the top candidates (Figure 2b). Given that BXYJ_LOCUS10935 encodes the pyruvate dehydrogenase β subunit, an important component of the pyruvate dehydrogenase complex and a key rate-limiting gate controlling the entry of glycolysis into the tricarboxylic acid (TCA) cycle, BxPDHB was ultimately selected as the target for subsequent studies (hereafter referred to as BxPDHB).

3.3. Silencing Efficiency of the Specific Interference Fragment

To determine the overall efficacy of the interference experiments, the relative expression of the target gene was compared among the blank control (nuclease-free water), negative control (non-specific random-sequence dsRNA), and experimental (specific-targeting dsRNA) groups by one-way analysis of variance (ANOVA). The results showed a highly significant overall difference among the three groups (F = 38.27, P < 0.001), indicating that the different treatments exerted an overall regulatory effect on the transcriptional level of the target gene and thus provided the statistical prerequisite for subsequent pairwise comparisons. We then used Tukey’s multiple-comparison test to dissect the specific differences among groups. In the pairwise comparison, there was no significant difference in target-gene relative expression between the blank control and the negative control (P > 0.05) (Figure 3a and Table S4). This negative result confirmed that the transfection procedure, buffer environment, and introduction of exogenous non-specific nucleic acid fragments did not produce unexpected background interference on the basal transcriptional activity of the target gene, effectively ruling out potential confounding by off-target effects on the subsequent phenotypic readouts.

In sharp contrast, the transcript level of the target gene in the experimental group (specific-interference-fragment treatment) showed a sharp and highly consistent downregulation relative to both control groups. Multiple-comparison analysis showed that the differences between the experimental group and both control groups reached a highly significant statistical level (P < 0.01) (Figure 3a and Table S4), indicating that the specific interference fragment could efficiently overcome the imperfect RNAi response barrier in the nematode and achieve strong suppression of target-gene transcriptional activity.

Taken together, these RT-qPCR results confirmed, within a rigorous statistical framework, that the interference fragment exerts a sequence-dependent and highly efficient silencing effect on the target gene BxPDHB. The successful establishment of this interference system provided a highly specific test material for subsequent phenotype-rescue experiments, locomotion assessment, and metabolomic validation, and is an important prerequisite for in-depth dissection of the core function of this gene in regulating the mitochondrial energy gate and the dispersal adaptability of JIII juveniles.

3.4. Effects of BxPDHB Silencing on Energy Metabolism and Locomotory Behavior

One-way ANOVA showed a highly significant overall difference among the blank control, negative control, and interference-treatment groups (F = 16.4, P < 0.0001). Tukey’s post-hoc multiple comparisons showed no significant difference in head-swing frequency between the blank control and the negative control (P > 0.05) (Figure 3b and Table S4), excluding non-specific interference of the transfection procedure and exogenous non-specific nucleic acids on nematode locomotory activity. By contrast, the head-swing frequency of the interference-treatment group was highly significantly reduced compared with both control groups (P < 0.001), with a ~50.00% reduction relative to the blank control. This physical behavioral phenotype confirmed the indispensable key-driving role of the target gene in regulating the dispersal locomotion of the nematode.

To clarify the biochemical basis of the impaired locomotion at the cellular energy-metabolism level, we determined the absolute ATP content of nematodes in each group using a bioluminescence method. One-way ANOVA showed a highly significant overall difference in ATP content among the three groups (F = 16.23, P < 0.001). Multiple comparisons showed no significant difference in ATP content between the blank control and the negative control (P > 0.05), indicating that the transfection procedure itself did not exert a non-specific impact on the cellular energy load of the nematode. However, the ATP content of the interference-treatment group showed a significant decrease, differing significantly from the negative control (P < 0.01) and highly significantly from the blank control (P < 0.001) (Figure 3c and Table S4), with a ~16.91% decrease relative to the blank control. This quantitative energy indicator confirmed that knockdown of the target gene directly limited the efficiency of ATP synthesis within B. xylophilus cells.

Concurrently, as the direct catalytic substrate of the pyruvate dehydrogenase complex, pyruvate content was measured synchronously in each group. One-way ANOVA showed a highly significant overall difference in pyruvate content among the three groups (F = 89.06, P < 0.0001). Multiple comparisons showed no significant difference in pyruvate content between the blank control and the negative control (P > 0.05), whereas the pyruvate content of the interference-treatment group was highly significantly increased (P < 0.0001) (Figure 3d and Table S4), with a ~4.61-fold increase relative to the blank control. This abnormal substrate accumulation confirmed at the molecular level that BxPDHB silencing blocked the catalytic conversion of pyruvate to acetyl-CoA, disrupting the basal carbon-metabolic homeostasis of the nematode.

3.5. Effects of BxPDHB Silencing on Oxidative-Stress Defense and Survival Rate

To investigate the effect of target-gene silencing on the gross survival capacity of the nematode, we first measured the survival rate of each treatment group. One-way ANOVA showed a highly significant overall difference in survival rate among the three groups (F = 43.53, P < 0.0001). Multiple comparisons showed that the negative control did not exhibit a decreased survival rate relative to the blank control (P > 0.05), indicating that the in vitro transfection procedure and buffer environment did not stress the physiological state of the nematode. By contrast, after specific target-gene silencing, the survival rate of the interference-treatment group showed a sharp and highly significant decrease relative to both control groups (P < 0.0001 for both comparisons) (Figure 4a and Table S4), with a ~26.02% decrease relative to the blank control. This marked deterioration of the survival phenotype preliminarily indicated that the gene plays an important defensive-barrier role in maintaining the basic survival capacity of B. xylophilus under ex vivo stress conditions. To further explain the mechanism underlying the decreased survival rate at the level of oxidative stress, we performed in situ tracing and quantitative analysis of intracellular ROS levels using a cell-permeable fluorescent probe combined with laser scanning confocal microscopy. One-way ANOVA (F = 985.20, P < 0.0001) and multiple-comparison analysis showed that the mean fluorescence intensity of nematodes in the blank control (CK) and the negative control (ds-EGFP) groups was both low, with no significant difference between the groups, indicating that the transfection procedure and introduction of exogenous non-specific sequences did not perturb the normal ROS homeostasis of the nematode. However, the mean fluorescence intensity of the specific-targeting interference group (ds-BxPDHB) showed a sharp and highly significant increase (Figure 4b and Table S4), with a ~4.26-fold increase relative to the blank control. This result confirmed, at the quantitative molecular level, that knockdown of the target gene directly induced pathological accumulation of intracellular ROS, revealing the indispensable core role of this gene in maintaining the normal oxidative-stress defense mechanism of the nematode.

Confocal imaging showed that the blank control and negative-control groups exhibited only weak, scattered background green fluorescence in the nematodes, indicating that the control groups were in a normal redox physiological balance (Figure 5a,b). By contrast, the interference-treatment group exhibited uniform, intense green fluorescence throughout the body (Figure 5c). This imaging evidence confirmed at the molecular level that knockdown of the target gene induced pathological accumulation of intracellular ROS in B. xylophilus, revealing that oxidative-stress damage resulting from gene loss is one of the important causes of the sharp decrease in survival rate.

4. Discussion

4.1. PDHB as a Key Metabolic Regulator at the JIII Dispersal-Juvenile Stage

Elucidation of the stage-specific metabolic regulatory mechanisms of the JIII dispersal juveniles (DL3) of B. xylophilus will help to clarify its energy-allocation strategy during dispersal adaptation. In this study, WGCNA revealed that the 144 genes with the highest correlation with the DL3 (JIII) stage (r = 0.97, P < 0.01) were significantly enriched in functional pathways including ribosome translation, redox and detoxification, and mitochondrial energy metabolism, suggesting that the JIII stage may involve coordinated regulation of protein synthesis, oxidative-homeostasis maintenance, and energy supply.

Similar stage-specific metabolic adaptation strategies are found in other nematodes. For example, the dauer stage of Caenorhabditis elegans maintains cellular homeostasis and enhances environmental adaptation through global metabolic reprogramming and selective gene-expression regulation [29,30], a feature shared with the JIII stage of B. xylophilus. However, the two stages may differ in their energy strategies: dauer is a survival state that mainly reduces metabolic rate for long-term energy conservation, whereas the JIII dispersal juveniles must maintain strong locomotory output to complete dispersal and may instead depend on higher mitochondrial energy supply. It is this coexistence of “high locomotory output and high stress-resistance demand” that renders the JIII stage an ideal model for dissecting the “locomotion–stress resistance” coupling.

Further analysis showed that, within the PPI network, the ribosome-translation and redox-and-detoxification modules were connected by genes related to mitochondrial energy metabolism, forming a functional association network centered on energy metabolism. Mitochondria not only undertake ATP production but also serve as an important regulatory platform for the generation of biosynthetic precursors, maintenance of redox homeostasis, and cellular stress responses [31]. The high enrichment of energy-metabolism functions in the JIII-stage gene network suggests that this stage must meet the physiological demands of locomotion, dispersal, and environmental adaptation by enhancing energy-metabolism capacity.

Through the integration of PPI-network centrality analysis, GSEA core-gene-set enrichment, and random forest feature-importance ranking, we ultimately identified PDHB as a candidate key regulatory factor. PDHB encodes the β subunit of the pyruvate dehydrogenase complex (PDHc), which catalyzes the oxidative decarboxylation of pyruvate to acetyl-CoA, a key step linking glycolysis and the TCA cycle [32]. From the perspective of metabolic-network architecture, PDHB is located at the junction of glycolysis and mitochondrial oxidative metabolism: upstream, it accepts the pyruvate generated by glycolysis, and downstream, it affects the entry of acetyl-CoA into the TCA cycle and subsequent oxidative phosphorylation. Thus, PDHB dysfunction may simultaneously affect carbon-utilization efficiency and mitochondrial energy output, thereby affecting ATP-dependent physiological processes such as locomotory behavior and oxidative-stress regulation.

4.2. Cascading Effects of PDHB Silencing on Energy Metabolism and Locomotion

RNAi experiments further validated the functional role of PDHB in energy-metabolism regulation. After PDHB silencing, pyruvate content in the nematode significantly increased while ATP levels markedly decreased, exhibiting a typical “substrate accumulation–energy-product deficiency” pattern. This result suggests that, when PDHB function is impaired, the conversion of pyruvate to acetyl-CoA is blocked, the carbon source entering the TCA cycle is reduced, and the energy output of oxidative phosphorylation is consequently decreased. Previous studies have confirmed that the PDH complex is a key node linking glycolysis and mitochondrial oxidative metabolism, and that its functional changes can significantly affect cellular carbon-flow distribution and energy homeostasis [33]. Therefore, the ATP decrease caused by PDHB silencing in this study may not be limited to a single metabolic reaction but may result in an overall reduction in energy-supply capacity.

The energy-metabolic disorder further propagated to the individual behavioral level. After PDHB silencing, the head-swing frequency of the nematodes significantly decreased, indicating that their locomotory capacity was markedly impaired. Nematode locomotion depends on the coordination of neural signal transmission, muscle contraction, and energy metabolism, among which the actin–myosin interaction driving muscle contraction continuously consumes ATP [34]. The reduced ATP supply caused by PDHB silencing may lower muscle-contraction efficiency and affect the locomotory output supported by the neuromuscular system. Previous studies have shown that mitochondrial dysfunction and impaired energy metabolism can lead to reduced locomotory capacity in nematodes, indicating that energy status is closely related to locomotory performance [35]. Therefore, the locomotory decline observed in this study may be a combined result of energy-supply deficiency and impaired neuromuscular function caused by PDHB loss.

From an ecological-adaptation perspective, the high expression of PDHB at the JIII stage may be related to the special energy demands of this stage. The JIII dispersal juveniles must maintain strong migratory capacity to complete the ecological transition of the dispersal phase, a process that depends on sustained energy supply [36]. Enhanced PDHB expression may promote pyruvate oxidative metabolism and improve mitochondrial energy-conversion efficiency to meet the higher locomotory demands of the dispersal stage. In this study, PDHB silencing disrupted this metabolic state, leading to reduced ATP supply accompanied by weakened locomotory capacity, further supporting the possibility that PDHB participates in regulating the energy-adaptation process during the dispersal stage of B. xylophilus.

4.3. Effects of PDHB Silencing on Oxidative-Stress Homeostasis and Survival

Mitochondrial oxidative phosphorylation, while producing ATP, is also an important source of intracellular reactive oxygen species (ROS). Under normal conditions, mitochondria maintain the balance between ROS production and clearance through electron-transport-chain regulation, antioxidant enzyme systems, and redox-buffering mechanisms [37]. Notably, the PDH complex itself is also considered a potential source of mitochondrial ROS generation, and abnormalities in its electron-transfer process may promote increased oxidative pressure [38]. In this study, after PDHB silencing impaired PDHc function, pyruvate oxidation was blocked and mitochondrial metabolic flow was remodeled, which may alter the redox state of the electron transport chain; at the same time, the reduced ATP supply may weaken energy-dependent antioxidant defenses such as the glutathione cycle [39]. Together, these two aspects ultimately led to pathological ROS accumulation. Therefore, the oxidative stress caused by PDHB silencing is likely the combined result of “electron-flow imbalance promoting ROS generation” and “ATP deficiency weakening ROS clearance,” although the relative contribution of each remains to be further investigated.

Confocal detection showed that, after PDHB silencing, the ROS fluorescence signal in the nematodes was markedly enhanced, and the mean fluorescence intensity was significantly higher than that of the control groups, indicating that reduced PDHB function induced elevated intracellular ROS levels. Similar phenomena have been observed in other systems; for example, PDH complex dysfunction can lead to increased intracellular H2O2 levels [40]. By contrast, the blank and negative controls exhibited only low background fluorescence, indicating that the experimental procedures themselves did not cause significant oxidative stress and that ROS accumulation was mainly related to the loss of PDHB function. This result suggests that, in addition to participating in energy-metabolism regulation, PDHB may also be involved in maintaining the redox homeostasis of the nematode.

The combined effects of enhanced oxidative stress and energy-metabolic disorder further affected the survival capacity of the nematode. In this study, the survival rate was significantly decreased after PDHB silencing. Notably, the negative-control survival rate was slightly lower than that of the blank control, suggesting that the RNAi treatment procedure may have imposed a certain degree of non-specific physiological stress; however, compared with the large decrease in survival rate caused by PDHB silencing, this effect was small, and the loss of PDHB function remained the major factor responsible for the significantly reduced survival. It is worth noting that, as a key node linking glycolysis and mitochondrial oxidative metabolism, the effect of PDH complex dysfunction on cellular energy and redox homeostasis has been confirmed in a variety of organisms, including nematodes [32,33,38], indicating a relatively high degree of species-level conservation of this metabolic hub. The phenotypic cascade of “substrate accumulation–energy depletion–oxidative stress–reduced survival” observed upon BxPDHB silencing is consistent with this mechanism. It should be noted that the direct causal relationship between BxPDHB silencing and reduced survival, as well as its generality across different life stages of B. xylophilus, still needs to be further confirmed through rescue experiments and multi-stage validation.

4.4. Implications for Forestry Pest Control

The results of this study have certain potential implications for the green control of pine wilt disease. First, PDHB silencing significantly reduced survival rate and inhibited locomotory capacity in mixed life stages (Figure 3 and Figure 4), suggesting that this gene not only serves the specialized demand of the JIII dispersal stage but may also maintain basal energy metabolism in other active life stages. If its general importance is further confirmed across more developmental stages, an RNAi strategy targeting PDHB could theoretically interfere simultaneously with population growth in propagative stages and the dispersal migration of the JIII stage, offering greater comprehensive control potential than strategies targeting stage-specific genes. Second, RNAi technology is gradually overcoming delivery bottlenecks; previous studies have confirmed that nanocarriers such as nanoclays, chitosan, and layered double hydroxides can protect dsRNA from degradation in the complex environment of pine phloem and enhance cell-penetration efficiency [42,43]. In the future, BxPDHB-dsRNA could be combined with these carriers and applied in the field via trunk injection or root application. In addition, the catalytic activity of the pyruvate dehydrogenase complex can be interfered with by commercially available pyruvate dehydrogenase kinase (PDK) inhibitors [44], providing a potential chemical-skeleton reference for developing small-molecule nematicides with both rapid action and systemic effects. Taken together, this study provides a theoretical basis and a molecular target for building an integrated control strategy from the dimension of “energy-metabolism-hub interference,” with PDHB as the molecular target and combining RNAi biopesticides with chemical synergists.

4.5. Limitations and Future Perspectives

This study has certain limitations: the transcriptomic data used were obtained from a public database with a limited sample size; the PPI network was constructed by cross-species mapping based on the C. elegans STRING database, which may introduce annotation bias. In addition, rescue experiments were not performed in the RNAi validation, and whether ROS accumulation is a direct result of mitochondrial damage or a secondary effect of energy depletion remains to be determined. In future studies, the network could be corrected by combining interaction information from the B. xylophilus genome itself, and the feasibility of PDHB as a control target could be further validated through rescue experiments, antioxidant intervention, and field RNAi delivery assessment.

5. Conclusions

In summary, by combining multi-omics analysis with functional validation, this study identified the β subunit of the pyruvate dehydrogenase complex (PDHB) as an important metabolic regulator at the JIII dispersal-juvenile stage of B. xylophilus (pseudonymously BxPDHB). The specific high expression of PDHB at the JIII stage may be related to the increased demand for energy supply and environmental adaptation at this stage. RNAi-mediated silencing of PDHB disrupted the entry of pyruvate into the mitochondrial oxidative-metabolism pathway, resulting in decreased ATP levels accompanied by reduced locomotory capacity, increased ROS accumulation, and decreased survival rate. These results indicate that PDHB may influence the locomotory performance and oxidative-stress homeostasis of B. xylophilus during its dispersal stage by regulating energy-metabolism status, revealing a potential link among energy metabolism, behavioral output, and environmental adaptation. This study provides new molecular evidence for understanding the metabolic-adaptation mechanism during the dispersal stage of B. xylophilus and provides a theoretical basis for further exploring dispersal regulation and potential control targets in this nematode.