The study was conducted at Dahuofang Forest Station in Fushun City, Liaoning Province, Northeast China (124°16′E, 41°58′N). The dominant tree species are P. koraiensis, Pinus tabuliformis (Chinese pine), and Larix olgensis (Changbai larch) in the forest station (Yu et al., 2015). This region was found natural occurrence of PWD in 2017 and recorded as the northernmost boundary of the PWD epidemic area in China (Yu and Wu, 2018). It has a typical continental monsoon climate within the mid-temperate zone. The mean annual temperature in this region is 5.5 °C, and the mean annual precipitation is 661.22 mm as well as mean annual evaporation is 1,411.78 mm according to CGIAR-CSI Global Aridity and PET Database (Zomer et al., 2022). The pine forests of P. koraiensis were established in the 1960s as pure stands, and all the pine individuals in the forests were similar in age (ca. 60 years old) and size (Hou et al., 2025).

The study was conducted in early October 2021 when the symptoms of newly infected pines with PWD could be visually determined (Ma et al., 2020; Hou et al., 2025). In the forest station, we selected four forest stands of P. koraiensis that follow a natural chronosequence of PWD, and a healthy L. olgensis stand. The selection was made based on unmanned aerial vehicle remote sensing data of discolored standing trees in the experimental site (Figure 1). In the stand of P. koraiensis where all the needles of each individual tree were reddish-brown, the trees were classified as fully killed by PWD (abbreviated as Killed Pk). In the neighboring stands where the pines had the majority needles of mixed colors of red and green, and these pines were considered diseased due to PWD (abbreviated as Diseased Pk). The healthy stand of P. koraiensis was composed of pine trees that had entirely green needles (abbreviated as Healthy Pk). A piece of triangle-shaped tissue was taken from the trunk of sampled individuals in each stand at breast height using a chain saw. The tissues were immediately transported to the laboratory and extracted using the “Baerman funnel method”. The microscopic examination verified the absence of B. xylophilus in the tissue of healthy pines and presence of B. xylophilus in the tissue of diseased and killed pines. In order to compare the soil communities of resistant vs. susceptible conifer species, we also selected a stand of larch species L. olgensis (abbreviated as Healthy Lo) in the vicinity of the killed stands of P. koraiensis (Fang et al., 2022; Wang et al., 2024). In addition, we also selected an open field nearby that had been created by clear-cutting the killed P. koraiensis (abbreviated as Clear-cut Pk) in the previous year. The inclusion of this field in the design was used to understand the consequences of PWD for soil microbial and nematode communities. It is noted that the ideal design was to select multiple stands that experienced similar levels of disease to avoid the pseudoreplication of sampling in the study. However, since the infection of P. koraiensis by PWD only recently occurred and the epidemic area was very small, such design was very difficult in practice. Therefore, we selected these forest stands that were spaced a minimum of 5 km in distance to minimize the effects of pseudoreplication on sampling errors. This selection was also to prevent the natural infection between the stands by the PWD at the time of sampling, given that the maximal spread distance of the vectored beetles was 5 km (Hou et al., 2025).

Within each forest stand, eight individuals of P. koraiensis or L. olgensis were randomly selected and the individuals were at least 30 m in distance between each other) and the surface soil under their canopy was collected (see below for details). Therefore, a total of 40 soil samples were collected (5 sampling fields × 8 samples per field) in the current study. All experimental sites were located on a hilly terrain orienting south to north, with slopes ranging from 40° to 50°. The sampled soils are of dark brown forest soils that have averaged pH 6.59, contents of 5.21 g·kg−1 total carbon, 0.38 g·kg−1 total nitrogen, and 0.78 g·kg−1 total phosphorus (Hou et al., 2025).

Before sampling, the understory vegetation and litter around each tagged tree were removed using a sterilized shovel. Topsoil at a 0–15 cm depth was collected from each of the four directions ca. 45 cm from the trunk using a soil sampler (38 mm in diameter). The four soil cores from each individual tree were thoroughly mixed to form a composite soil sample and placed in a new plastic bag. All the collected samples were stored in an ice box and immediately transported to the laboratory at the Shenyang Institute of Technology, Shenyang, China. The shovel and soil corer were sterilized using 75% alcohol between sampling events. In the laboratory, each soil sample was divided into three subsamples for microbial DNA extraction, nematode extraction, and physicochemical analyses.

Soil nematodes were extracted from 100 g fresh soil sample using a modified Baermann funnel method (Viglierchio and Schmitt, 1983). The extracted soil nematodes were preserved in formalin for later identification. The first 100 encountered nematodes in each sample were identified using a microscope (OLYMPUS CKX53, Tokyo, Japan). If the sample contained fewer than 100 nematodes, all the nematodes were identified. Nematodes were determined at the genus level according to identification guides of Bongers (1988) and Zhang et al. (2013). The nematode density was calculated to be total number of individuals per 100 g dry soil. All nematodes were further categorized into plant feeders (PF), bacterial feeders (BF), fungal feeders (FF), carnivores (CA), and omnivores (OM) according to their feeding habits (Yeates et al., 1993).

All samples were analyzed for pH, soil water content (SWC), total carbon (TC), total nitrogen (TN), and total phosphorus (TP). Soil pH was determined at a soil-to-water ratio of 1:2.5 (w/v) by the pH meter S210–K (Mettler Toledo, Germany). SWC was estimated by oven-drying the samples at 105 °C to constant weight. TC and TN were analyzed using a CN elemental analyzer (Elementar Vario EL III, Elementar, Germany). TP was measured by colorimetric analysis after digestion with sulfuric acid and perchloric acid.

Soil microbial genomic DNA was extracted using the PowerSoil DNA Isolation Kit (Omega Bio-tek, Norcross, GA, USA) following the manufacturer's protocol. The primers 338F (5′-ACTCCTACGGGAGGCAGCAG-3′)/806R (5′-GGACTACHVGGGTWTCTAAT-3′) were used to amplify V3-V4 regions of the bacterial 16S rRNA gene and the primers ITS1F (5′-CTTGGTCATTTAGAGGAAGTAA-3′)/ITS2R (5′-GCTGCGTTCTTCATCGATGC-3′) were used to amplify the fungal ITS1 region (Crowther et al., 2014). Genomic DNA from each sample was independently amplified and the PCR products were pooled and purified using the QIAquick Gel Extraction Kit (Qiagen). Subsequently, the concentration of purified amplicons was quantified using the Picogreen assay (Invitrogen) on a TBS-380 microplate reader, and equimolar concentrations were mixed. Finally, the pooled samples were subjected to paired-end sequencing (2 × 250 bp) on the Illumina Novaseq 6000 platform at Majorbio Bio-Pharm Technology Co., Ltd. (Shanghai, China). One bacterial sample from the clear-cut soil stand was excluded due to failed amplification.

The initial sequence processing was conducted using QIIME2 (version 2021.2) (Bokulich et al., 2018). Samples were first demultiplexed using the Demux plugin, followed by adapter and Barcode sequence removal using the Cutadapt plugin (Martin, 2011). Subsequently, the DADA2 plugin was employed for filtering, denoising, and quality control, resulting in the generation of the Amplicon Sequence Variants (ASV) table (Callahan et al., 2016). For bacterial sequences, taxonomic classification was performed using the SILVA (V138.1) database (Quast et al., 2013). For fungal sequences, the UNITE (V8.2) database (Kõljalg et al., 2013) was used for species annotation of ASV sequences. Finally, data were rarefied across samples for subsequent biodiversity analysis and community structure comparisons.

FAPROTAX was used to predict bacterial ecological functions by analyzing gene sequences and inferring their functional potential, such as metabolic functions and aromatic compound degradation (Louca et al., 2016; Sansupa et al., 2021). For the fungal data, we used FUNGuild to assign fungal sequencing data into three trophic groups (symbiotroph, saprotroph, and pathotroph), by comparing them against a database of known fungal genes associated with specific nutritional modes (Nguyen et al., 2016). Sequences that were assigned more than one trophic strategy in FUNGuild were excluded from our analysis (Deng et al., 2023). For each soil sample, the relative abundances were calculated for each functional or trophic group.

To explore the co-occurrence patterns of bacterial, fungal, and nematode communities in different forest types, correlation-based network analyses were conducted for each community as well as inter-kingdom network analysis. Only bacterial ASVs, fungal ASVs, and nematode genera present in at least 25% of the samples were included to ensure correlation reliability. To select highly significant associations, overall networks were constructed using correlations that were significant (p < 0.05) with strong positive (r > 0.6) or strong negative (r < −0.6) coefficients. Additionally, to simplify networks for better visualization and analysis, we used ASVs with relative abundance higher than 0.01% to construct the microbial networks for different soil origins (Khatri-Chhetri et al., 2024). Relationships among microbial and nematode communities were visualized in Gephi (version 10.1) (Bastian et al., 2009) using the Fruchterman Reingold layout (Fruchterman and Reingold, 1991). Various network structural attributes including node level topological properties such as degree and density were calculated in package “igraph”. Nodes (ASVs of bacteria and fungi, and genera of nematodes) are the essential units of a network, while edges represent links between nodes.

The α-diversity indices including species richness and Shannon index (H'), as well as β-diversity were calculated using the “vegan” package (Oksanen et al., 2013). We employed one-way ANOVA to examine the effects of PWD chronosequence on diversities of soil microbial and nematode communities, in which α-diversity indices were included as response variables and “soil origin” as a fixed factor. Tukey's HSD (Honest Significant Difference) post hoc tests were used to compare α-diversity among soil origins. Principal coordinate analysis (PCoA) was used to visualize the separations of the communities. Effects of PWD on soil microbial and nematode community composition were assessed by permutation multivariate analysis of variance (PERMANOVA) with Bray-Curtis dissimilarity (999 permutations) in “vegan” package (Oksanen et al., 2013). Canonical Correspondence Analysis (CCA) was performed to examine the relationships between soil biota and soil physiochemical properties. Prior to the CCA, the Detrended Canonical Analysis (DCA, 999 permutations) was performed to assess the explained variance of each environmental factor, using the “envfit” function in “vegan” package (Dong et al., 2024). Differential relative abundances of taxa and functions were analyzed using the Kruskal-Wallis test with Benjamini-Hochberg correction in Statistical analysis of taxonomic and functional profiles (STAMP) software (Parks et al., 2014). Kruskal-Wallis test was alternatively used in the “agricolae” package when the assumptions were not met even after data transformation. The p-values of the Kruskal-Wallis test were adjusted using the Benjamini-Hochberg method to control the false discovery rate (FDR). All statistical analyses and data visualizations were performed using R version 4.4.1 (R Core Team, 2024).

The species richness of bacterial community in soil of healthy P. koraiensis was 20.8% lower than in clear-cut soil and 17.3% lower than in soil of diseased pines, respectively (Figure 2A). Similarly, the species richness of fungal community in soil of healthy P. koraiensis was lower than in clear-cut soil (by 24.1%) and in soil of killed pines (by 24.5%) (Figure 2B). The Shannon index of both bacterial and fungal communities was also significantly lower in the soil of healthy pines than in soils of other origins including in soil of larch, but the differences were relatively small (all <5.0%) (Figures 2A, B). However, this difference in bacterial diversity disappeared when P. koraiensis was killed by PWD (Figure 2A). On the contrary, the diversity of soil nematode community did not differ across different soil origins (Figure 2C).

The first two axes of PCoA explained 35.1% and 10.4% of the variance in soil bacterial community, respectively. The composition of bacterial community in the soil of healthy P. koraiensis differed from that in soil of larch, as well as from the compositions in soil of diseased or killed P. koraiensis trees, although the latter two soils did not differ (PERMANOVA: R² = 0.47, p < 0.001; Figure 3A and Supplementary Table 1). A similar pattern was observed for both soil fungal community (8.6% and 6.8% for PC1 and PC2, respectively; PERMANOVA: R² = 0.21, p < 0.001; Figure 3B and Supplementary Table 1) and nematode community (17.4% and 12.8% for PC1 and PC2, respectively; PERMANOVA: R² = 0.27, p < 0.001; Figure 3C and Supplementary Table 1). The results of CCA revealed that edaphic factors including pH, SWC, TC, TN, and TP all significantly contributed to the difference in the compositions of microbial and nematode communities of these soil origins, except TC that was not related to the composition of nematode communities (Figure 3 and Supplementary Table 2). However, only the pH and SWC among these soil factors of diseased and killed P. koraiensis were higher than healthy P. koraiensis (Supplementary Figure 2).

The relative abundances of Proteobacteria and Gemmatimonadota was lower while the abundance of Methylomirabilota was higher in soil of diseased or killed P. koraiensis than healthy pines (Figure 4A and Supplementary Figure 3A). At genus level, the relative abundance of Mycobacterium was higher, while the relative abundance of Microlunatus and Vicinamibacteraceae were lower in soil of healthy P. koraiensis than diseased or dead P. koraiensis (Supplementary Figure 3D). The relative abundance of phylum Basidiomycota was lower but the abundance of Mortierellomycota was higher in soil of diseased or killed P. koraiensis than in soil of healthy pines (Figure 4B and Supplementary Figure 3B). At genus level, the relative abundance of Solicoccozyma was higher but that of Sebacina and Cryptococcus was lower in soil of diseased or killed P. koraiensis than in soil of healthy pines (Supplementary Figure 3E). We also found a lower relative abundance of genus Aphelenchoides in soil of killed P. koraiensis than in soil of healthy pines.

In addition, the relative abundances of Methylomirabilota, Chloroflexi, and Firmicutes in bacterial community were also lower in soil of healthy P. koraiensis than in soil of larch (Figure 4A and Supplementary Figure 3A). We also observed lower relative abundances of fungal genera Mortierella, Trichocladium, and Penicillium in the soil of healthy P. koraiensis than in soil of larch (Supplementary Figure 3E). In nematode community, the relative abundance of Pararotylenchus was lower but abundance of Aphelenchoides was higher in soil of healthy P. koraiensis than in soil of larch (Figure 4C and Supplementary Figure 3C).

Bacterial ASVs related to functions in chemoheterotrophy and aerobic chemoheterotrophy had the highest relative abundance across all soil samples (Figure 5A). We found a lower abundance of bacterial ASVs of aromatic compound degradation but higher abundances of cellulolytic, predatory or exoparasitic bacterial taxa in the soil of healthy P. koraiensis than in soil of diseased or killed P. koraiensis (Figure 5A and Supplementary Figure 4A). In addition, we observed a lower relative abundance of bacteria related to nitrate reduction and aromatic compound degradation in soil of healthy P. koraiensis than in soil of larch (Figure 5A and Supplementary Figure 4A). The relative abundance of symbiotrophs in the fungal community, e.g., ectomycorrhizal fungi (EcM) dominated by Agaricomycetes, was significantly higher in soil of healthy P. koraiensis than in soil of diseased or killed P. koraiensis (Figure 5B and Supplementary Figure 4B). On the contrary, the relative abundance of saprotrophs was lower in the soil of healthy P. koraiensis than in the soil of larch (Figure 5B and Supplementary Figure 4B). Regarding the nematode community, the relative abundances of plant feeders seemed to be lower in the soil of diseased P. koraiensis than in the soil of healthy P. koraiensis (Figure 5C).

The numbers of nodes and edges in the inter-kingdom, bacterial, and fungal community networks were overall lower in soil of forest stands compared to those in soil of clear-cut field (Figures 6A–C). The inter-kingdom network complexity decreased following the progression of PWD, with the proportion of bacterial nodes decreasing but the proportion of fungal nodes increasing as the severity of PWD intensified (Figure 6A and Supplementary Table 3). The numbers of nodes and edges in bacterial community networks decreased but those numbers in fungal community networks increased over the progression of PWD (Figures 6B, C). The number of edges in nematode community networks was lower in soil of diseased P. koraiensis than in soil of healthy P. koraiensis, but this difference disappeared when trees were killed over time (Figure 6D). The average degree and density of the nematode community network first decreased in diseased pines but surpassed the original values in healthy pines at the later stage when the trees were killed (Supplementary Table 4).