The study area was the Wangpicheon Ecosystem and Landscape Conservation Area (ELCA), located in central South Korea. This region was designated by the Ministry of Environment in 2005 as a protected area for habitat conservation for endangered species and rare wild animals and plants. The P. densiflora growing in this region are the largest in South Korea and grow with characteristically straight trunks. The total area of the conservation area is 102.84 km2, and according to forest type maps (1:5,000) from the Korea Forest Service, P. densiflora occupies the widest area, at approximately 35.37 km2 (34.40%; Supplementary Table S1). When we analyzed meteorological data for 1990–2023 using data from the Uljin Meteorological Office, the mean annual temperature was 12.9 °C and the mean annual rainfall was 1,186.64 mm, with most rainfall concentrated in the period from July to September (Figure 1b). To identify sites of multiple-tree mortality of P. densiflora within the conservation area, we conducted remote sensing using drones over 3 years from 2021 to 2023. Using the collected orthophotos (resolution 5–10 cm), we identified dead P. densiflora based on information such as color and appearance (Kim J. et al., 2017) and then extracted multiple-tree P. densiflora mortality sites based on the definition of these sites (Supplementary Figure S1a). Later, we conducted field surveys to verify the multiple-tree mortality of P. densiflora at these locations. Using these methods, we were ultimately able to identify 15 areas where multiple-tree mortality had occurred (Figure 1a). Due to the strict management system in this region, pine wilt disease and other pests are absent. When comparing climate conditions during the study period (1990–2023; Figure 1b) with historical climate data prior to 1990(Supplementary Figure S2), we attributed the observed P. densiflora mortality primarily to physiological stress induced by recent climatic changes, including more frequent and prolonged droughts, elevated temperatures, and altered precipitation patterns (e.g., decreased precipitation in winter and spring). Here, we defined climate change-related mortality as multiple-tree mortality occurring in the absence of pest or disease infestation and temporally coinciding with periods characterized by these documented climatic stressors. Indeed, no pest damage was observed during our field surveys.

Next, using the 15 secured multiple-tree mortality sites, we analyzed existing aerial photographs to identify the time of tree death and to estimate the time since tree death. For this purpose, we used aerial photographs from the South Korean National Geographic Information Institute, Kakao Map, and Google Earth Pro. From the images, we estimated the time of tree death to be the time when the tree branches first spread in a star shape, and the color began to change to white or red. By analyzing previous aerial images, we broadly categorized multiple-tree mortality of P. densiflora in the survey area into three stages (Supplementary Figures S4, S5): tree deaths before 2008 or from 2008 to 2011 were classified as Group 1 (n = 6). These were the oldest cases in this study, and as of 2024, over 15 years have passed since the multiple-tree mortality event. Tree deaths occurring from 2014 to 2017 were classified as Group 2 (n = 7). These were defined as intermediate multiple-tree mortality events that occurred 7–10 years ago, as of 2024. Finally, tree deaths occurring from 2020 to 2021 were classified as Group 3 (n = 2). These were the most recent multiple-tree mortality events, occurring approximately 3 years before the study period in 2024.

We conducted field surveys of the 15 P. densiflora multiple-tree mortality sites during the active tree-growth period from June to September 2024. We also selected and surveyed control plots comprising healthy P. densiflora forests with similar local conditions located at least 30 m from the multiple-tree mortality sites and, owing to the intervening topography and slope, positioned so that dead or declining trees were not visible from the control plots. Therefore, we conducted field surveys at 30 sites in total. Considering the local stand characteristics and number of dead trees, we used survey plot sizes of 20 m × 10 m (12 sites), 15 m × 15 m (1 site), 20 m × 20 m (1 site), and 20 m × 15 m (1 site) for the tree mortality sites. For all control sites, we conducted the surveys within 20 m × 10 m quadrats. For the vegetation survey, we investigated dominance, combining the coverage and abundance of species in the plots by stratum (tree, sub-tree, shrub layer, and herb layers) in accordance with the phytosociological survey method of Braun-Blanquet (1964). The dominance index of each species was determined using the Braun-Blanquet scale as follows: r, rare species; +, coverage <1%; 1, coverage 1–5%; 2, coverage >5–25%; 3, coverage >25–50%; 4, coverage >50–75%; 5, coverage >75–100% (Canullo et al., 2012). Diameter at breast height (DBH) was surveyed for all individuals in the plot with a DBH of at least 6 cm. For environmental factors, we used a handheld GPS 64S (Garmin, Olathe, KS, USA) and clinometer (Suunto, Vantaa, Finland) to record the elevation, slope, slope direction, and geographic coordinates, and we also recorded rock exposure and topographic position in the field. For P. densiflora with a DBH of <6 cm, we surveyed three different tree height categories (h ≤ 10 cm, 10 < h < 50 cm, h ≥ 50 cm) using the folding rule. Next, we evaluated the decay class of coarse woody debris (CWD) from P. densiflora in the plots, surveying snag (standing dead trees), downed wood (i.e., logs), and stumps separately. To evaluate the decay class of snag and log, we referred to the stand-level biodiversity monitoring protocol steps for field data collection from British Columbia (Province of British Columbia, 2009; Supplementary Figure S3), and for stumps we referred to data from Rouvinen and Kouki (2002; Supplementary Table S2). Finally, we used a fisheye lens at the center of each multiple-tree mortality survey plot to obtain crown projection data and analyze crown openness. Crown openness was analyzed using Gap Light Analyzer 2.0 (Frazer et al., 1999). The multiple-tree mortality sites have been summarized in Supplementary Table S3.

For soil physicochemical analysis, three bulk soil samples were taken from a depth of 0–10 cm within a 1-m radius of dead trees and mixed after removing the organic horizon at multiple-tree mortality site. In the control plots, soil samples were collected within a 1-m radius of living P. densiflora trees and assessed using the same method as for the multiple-tree mortality sites. We collected 15 samples each from multiple-tree mortality sites and control plots and transferred the samples to a laboratory where they were naturally dried, passed through a 2-mm sieve, and sealed. The samples were then sent to the South Korean Analysis Technology Center, a professional soil sample analysis company. The methods of analysis for each sample are described in Supplementary Table S4. In summary, we analyzed 11 chemical properties, including total nitrogen (N), exchangeable cations (K⁺, Ca²⁺, Mg²⁺, and Na⁺), available phosphorus (measured as Olsen-P), organic carbon (C), organic matter (OM), cation exchange capacity (CEC), soil pH, and electrical conductivity (EC).

For the rhizosphere soil survey, we combined analysis of past aerial photographs from multiple-tree mortality sites identified using prior remote sensing with the results of decay class assessment conducted in the field. Based on the results, we collected samples according to the time since tree death. We aimed to minimize errors due to local conditions by collecting samples from nearby survey sites in locations with similar conditions; the selected sites were located in close proximity to each other and shared similar ridge-top soil properties, microclimatic conditions (e.g., strong winds and ample sunlight), and initial stand structure. Therefore, among the total 15 multiple-tree mortality sites, we selected 10 representative plots that best met these strict homogeneity criteria for the final analysis in Section 2.4. To minimize micro-scale soil heterogeneity, three subsamples were collected from different directions around the rhizosphere of each individual tree and pooled into a single composite sample. Specifically, we collected soil tightly adhering to fine roots, which distinguishes these samples from the bulk soil collected for physicochemical analysis in Section 2.2. For Group 1 we selected three snags (decay class 6) each from three sites and collected soil samples from depths of 0–10 cm after removing the organic horizon. For Group 2 we selected three snags (decay class 4) each from five sites for soil collection. For Group 3 we selected three snags (decay class 3) each from two sites for sample collection. Of these, one site contained decay class 3 dead trees (twigs and pinecones remaining, showing the same decay condition as a multiple-tree mortality site from 3 years prior) located at the site of a previous multiple-tree mortality site from 2008. To compare with sites of tree mortality, we collected control rhizosphere soil samples from healthy P. densiflora forests in similar locations at least 50 m from the soil samples collected from multiple-tree mortality sites. In this process, we selected three healthy P. densiflora at each of the four sites and collected samples from locations where dead trees were not visible due to the terrain and slope, to ensure that the sampling was not influenced by the extending root systems or mycelia of nearby dead trees. We additionally collected four rhizosphere soil samples from one declining P. densiflora forest site (defined as the ‘decline’ group). Here, the health condition of the declining trees was evaluated based on defoliation in accordance with the method of Huo et al. (2019). The trees were found to be in the 40–80% defoliation range (Supplementary Figure S6). All the rhizosphere soil samples were placed in 50-mL conical tubes and transported to the laboratory in an icebox to maintain a low temperature. At the laboratory, the samples were stored at -20 °C. The local conditions for all collected samples showed the characteristics of a ridge environment, with strong winds and ample sunlight. The local conditions and soil environment at the collection sites are described in Supplementary Table S5.

Subsequently, after transporting the soil samples to the laboratory at a low temperature, a DNeasy^® PowerSoil^® Pro Kit (Qiagen, Hilden, Germany) was used to extract DNA from the stored samples, and a Nanodrop (ND-LITE-PR; Thermo Fisher Scientific, Waltham, MA, USA) was used for quality control. DNA sequencing was performed at the Kyungpook National University NGS Center (Daegu, Republic of Korea). A Miseq system (Illumina, San Diego, CA, USA) was used for bacterial 16S V4 amplicon sequencing, and a Hiseq3000 (Illumina, San Diego, CA, USA) was used for fungal ITS2 amplicon sequencing. The amplified DNA set was cleaned using Beckman Coulter AMPure XP (Thermo Fisher Scientific, Waltham, MA, USA) and amplified again to tag the different barcodes of each sample. The completed libraries were sequenced using each Illumina platforms using universal primer pairs (515F/805R and ITS86F/ITS4R) (Supplementary Table S6).

Finally, we used QIIME2 and DADA2 to remove chimera reads, compile phylogenetic trees, and perform classification and normalization from the DNA raw data (Callahan et al., 2016; Bolyen et al., 2019). The databases used for classification included SILVA (v 138) for bacteria and UNITE (v 8.3) for fungi (Abarenkov et al., 2010; Bokulich et al., 2018). To compile phylogenetic trees, we used SATé-enabled phylogenetic placement (SEPP) and the SILVA database for bacteria and Mafft 7 for fungi (Mirarab et al., 2012; Katoh and Standley, 2013). Furthermore, all samples were rarefied to an equal sequencing depth, and the rarefied feature tables were used for downstream statistical analysis and visualization (Supplementary Table S7).

P. densiflora CWD showed a gradual progression of decay from snags to logs and stumps (Figure 2a). Three years after tree death (Group3), the snags were mostly decay class 3, but after 7–10 years (Group2), they changed to decay class 4, and beyond 15 years (Group1), the snag forms reduced. P. densiflora saplings also showed almost no regeneration 3 years after tree death but regenerated as more time passed (Figure 2b). In addition, the extent of decay of some stumps and logs in Group 1 survey sites was so severe that it was difficult to evaluate the decay class; these were excluded from the analysis. Meanwhile, some snags observed in 2020–2021 had a relatively high decay class of 4–5, which was due to natural disasters (we observed burn marks consistent with lightning strikes).

We performed CCA between multiple-tree mortality sites, differentiating by vegetation characteristics and local conditions. The analysis revealed no statistically significant environmental variables that distinguished the different multiple-tree mortality site groups from each other (Figures 3a, b; Table 1). In contrast, the analysis comparing all multiple-tree mortality sites against the control sites (Figure 3c) showed a significant separation between the two groups (Overall model p < 0.001). The key environmental variables driving this difference were pH, Elevation, and Rock exposure (Table 1). To summarize, compared with the other groups, Group 1 was strongly associated with P. densiflora, Quercus. spp, and Rhododendron mucronulatum in the understory vegetation, and strong associations with total nitrogen and available phosphorus, which indicate soil fertility in environments with open canopies and high rock exposure, were observed (Figure 3). Particularly, as shown in Figure 3b, canopy openness tended to be inversely related to P. densiflora BA and Quercus spp., suggesting that the local environment for Group 1 was characterized by open forest gaps following the death of overstory vegetation. To investigate this relationship in more detail, we performed univariate analyses between canopy openness and the basal area (BA) of P. densiflora and Quercus spp. The simple linear regression showed no significant association between canopy openness and the BA of either P. densiflora (Estimate = 0.001, p = 0.998) or Quercus spp. (Estimate = −0.794, p = 0.169). Similarly, the Spearman partial correlation analysis, revealed weak negative but non-significant correlations for both species (Estimate = −0.21, p = 0.59 for P. densiflora; Estimate = −0.19, p = 0.63 for Quercus spp.). When we compared characteristics between multiple-tree mortality and control sites, the control sites showed similar vegetation characteristics to Group 3 and some of Group 2. Specifically, compared with Group 1, these sites showed relatively high SR, strong associations with Lindera obtusiloba and Vaccinium hirtum var. koreanum, and were closely related to soil electrical conductivity and soil pH. The IVs analysis showed some differences in dominant species among the groups (Table 2). In Group 1 (the earliest mortality sites), P. densiflora was dominant. In contrast, Fraxinus sieboldiana was dominant in the intermediate (Group 2) and most recent (Group 3) mortality sites. Notably, Group 2 showed a relatively higher IV for P. densiflora compared to Group 3, and P. densiflora was not among the top 10 species in the Group 3 understory.

For fungi, in terms of diversity, healthy P. densiflora showed lower diversity than the other groups, and diversity was highest 15 or more years after tree death (Figure 4a). A statistically significant difference was observed only in the Chao1 index between the Alive and Decline groups, and no significant differences were found among the other groups (Figure 4a). In NMDS analysis, R was 0.172, indicating slight differences between groups, but mostly overlapping trends were observed (Figure 4b). In particular, the healthy P. densiflora group showed the most similar characteristics to Group 1, which had the longest duration since tree death, while Group 3 and the Decline group showed similar trends. Group 2, meanwhile, showed mixed characteristics. For bacteria, in terms of diversity, unlike fungi communities, diversity considerably increased in the Decline group and Group 3, followed by a gradual decrease (Figure 4c). Similar to the fungi community, a statistically significant difference was found only in the Chao1 index, and only between the Alive and Decline groups. In NMDS analysis, similar to fungi, we observed a low R-value (Figure 4d). However, the healthy P. densiflora group showed characteristics similar to those of Group 1, which had the longest duration since tree death. Some of Group 2 also showed similar characteristics to the healthy P. densiflora group, but overall showed a mixed distribution.

We performed analysis using an RF model to identify the main indicator genera driving differences between the groups, and the 15 taxa with the highest importance, alongside their relative abundance ratio, are shown in Table 3. The results showed that, for fungi, several genera in the Basidiomycota and Ascomycota phyla were the most important indicator taxa, whereas for bacteria, no taxa had notably high importance. Unclassified groups among Ascomycota (fungi) and Proteobacteria (bacteria) showed relatively high abundance in healthy P. densiflora, suggesting that these might play important roles in tree health. Meanwhile, Oidiodendron spp. and Venturia spp. in the phylum Ascomycota and Umbelopsis spp. in the phylum Mucoromycota showed high relative abundance in healthy P. densiflora and lower abundance in declining P. densiflora. This suggests that a decrease in these microbes could be related to the gradual decline of healthy P. densiflora. Fungi unclassified at the genus level, as well as Russula spp. and Talaromyces spp. in the phyla Basidiomycota and Ascomycota, and Umbelopsis spp. in the phylum Mucoromycota showed relatively high abundance in Group 2.

We conducted network analyses separately for fungal and bacterial communities across five groups (Alive, Decline, and Mortality Groups 1–3). We primarily described the structural properties at the whole-network level, while using the Largest Connected Component (LCC) level for complementary analysis. In the co-occurrence network based on fungal groups, the main connectors at the phylum level were Mucoromycota, Basidiomycota, and Ascomycota, with Basidiomycota and Ascomycota making consistently high contributions to the network (Figure 5). We observed high modularity in healthy P. densiflora (whole network: 0.825, LCC: 0.563), indicating that specific microbial communities exhibited discernible network structures rather than random associations (Table 4). Meanwhile, in the Decline group, modularity was relatively high at the whole-network level (0.772) but dropped to –0.002 at the LCC level. Notably, this network contained a relatively large number of nodes and edges, and also exhibited a considerable proportion of negative correlations compared with the other groups. Group 3 showed a relative decrease in the number of nodes and edges compared to the Decline group, but exhibited modularity comparable to that of the Decline group (whole network: 0.849, LCC: –0.012). In Group 2, although the modularity of the whole-network level (0.407) was relatively low, the network at the LCC level (0.393) exhibited a modular organization where the main interactions were concentrated. In addition, the numbers of nodes and edges were relatively high at both the whole-network and LCC levels, forming a compact network structure. Meanwhile, in Group 1, modularity increased again at the whole-network level (0.802), whereas the modularity at the LCC level (–0.031) markedly decreased. In addition, the numbers of nodes and edges were similar to those in the Alive group.

In terms of connectors at the phylum level, Proteobacteria, Verrucomicrobiota, Firmicutes, Chloroflexi, and Actinobacteriota showed consistently high contributions to the network across temporal changes in bacterial network structure (Figure 6). In healthy P. densiflora stands, modularity was relatively high compared with the other groups (whole network: 0.387; LCC: 0.348). At the whole-network level, the network comprised 209 nodes and 722 edges (Table 4). Meanwhile, in the Decline group, the number of edges (1,457) and nodes (225) markedly increased at the whole-network level, whereas they decreased at the LCC level (edges: 406, nodes: 29). Furthermore, modularity showed a marked increase at the whole-network level (0.868), while it dropped sharply at the LCC level (–0.002), and, similar to fungi, this group exhibited a relatively high proportion of negative correlations compared with the other groups. In Group 3, similar to the Decline group, modularity was high at the whole-network level (0.918) but low at the LCC level (–0.080). Meanwhile, the numbers of nodes and edges were the lowest among all groups. Consistent with the fungal network results, Group 2 showed lower modularity at the whole-network level (0.664) compared with the Decline and Group 3 networks, but markedly higher modularity at the LCC level (0.582). In addition, at the LCC level, Group 2 contained more nodes and edges than that of Group 3. In Group 1, modularity (whole network: 0.647, LCC: 0.229) was lower than in Group 2 but higher than in the Decline and Group 3 groups, and the whole network comprised 177 nodes and 147 edges.

We compared the time of multiple-tree mortality and CWD decay classes in P. densiflora forests, and found that after the death of P. densiflora, the CWD generally decayed from snag forms into logs and stumps, following a similar trajectory depending on the elapsed time since tree death. Although these findings align with predictable natural decay processes, it is significant that we quantitatively demonstrated this progression at specific time intervals, tracking the conversion from snags into logs and stumps through changes in decay classes. Indeed, our observations primarily reflect this common decay sequence. However, it is important to acknowledge that this observed decay progression does not comprehensively represent all forest dynamics, as live trees can also directly fall to logs due to windthrow or other disturbances, bypassing an intermediate snag stage. Consequently, our results should be interpreted with caution, recognizing inherent variability and complexity in forest stand dynamics. Meanwhile, this dataset is particularly valuable because it serves as a reference for estimating the timing of tree death across diverse regions based on dead-tree decay classes.

Historically, various forest gap sizes have been reported to alter light availability within stands, inducing natural disturbances in the area (Canham, 1988; Bolton and D’Amato, 2011; Lee et al., 2024). According to our findings, healthy P. densiflora forests and Group 3 and some Group 2 sites, which are considered to be at an early stage after tree death, showed similar species composition and did not show major differences. The species that appeared usually grow in the understory vegetation of pine forests on ridges. Meanwhile, Group 1, at least 15 years after tree death, showed relatively favorable levels of elements that affect soil richness (Supplementary Table S5). In addition, we observed the entry and regeneration of P. densiflora and Quercus spp. in the understory vegetation. One of the shared characteristics of these sites is that as dead P. densiflora decayed, forest gaps were opened. In other words, even though the light environment improved, the entry of woody plants other than P. densiflora was slow, and this is because P. densiflora habitats are usually on ridges (or upper slopes), where strong winds, ample sunlight, and a low topographic wetness index create low soil richness and harsh conditions. These environmental factors are important in allowing P. densiflora to be dominant as an edaphic climax, suggesting the potential for continual regeneration of this species. In particular, the fact that Quercus spp., which has relatively strong shade tolerance, was established as a succession species alongside P. densiflora is similar to findings in several previous studies in temperate regions (Higo et al., 1995; Beon and Bartsch, 2003; Choung et al., 2020). This supports the notion that this is a typical transition from P. densiflora forests to P. densiflora-broadleaved mixed forests. Furthermore, many previous studies have reported that the progressive increase in CWD decay in the stand over time gradually becomes an important factor affecting the soil environment (Harmon et al., 1986; Noh et al., 2017). Indeed, considering that Group 1 showed overall favorable levels of elements affecting soil richness compared to other groups (Supplementary Table S5), we cannot exclude the possibility that these changes affected the processes of plant entry and regeneration. However, there are limitations to discussing this clearly based on the current study, given that various ecological dynamic processes act in concert, including climate and local environmental characteristics (Sturtevant et al., 1997; Laiho and Prescott, 2004). Conducting further studies to better understand the processes of CWD decay, changes in the soil environment, and the relationships with species entry would contribute to even more detailed knowledge of the processes of change in aboveground stand behaviors in P. densiflora forests. To support practical forest management decisions, future research should specifically investigate the interactive effects of CWD decay stages and soil nutrient dynamics on species establishment success. Furthermore, in sub-alpine coniferous forests in South Korea, the entry of diverse species into forest gaps that appear as dominant tree species die reportedly causes rapid changes in forest structure (Lee et al., 2024). In lowland P. densiflora forests, even though multiple-tree mortality of P. densiflora leads to further expansion of forest gaps over time, we observed regeneration of P. densiflora and slow entry of other species. This shows differences in the pattern of species entry after the death of dominant tree species in high- and low-altitude coniferous forests. Furthermore, canopy openness was not identified as a significant variable distinguishing the multiple-tree mortality sites. This might be because the statistical difference did not emerge as the forest gap between Group1 and Group2 did not show a significant difference. However, as shown in Figure 2, it is suggested that the forest gap is gradually expanding as snags decay and fall over time. Furthermore, just as Group1 and some of Group2 showed a high association with rock exposure, it is possible that the formation of gaps and substrate characteristics (e.g., soil, moss cover, rock exposure) interacted to promote the regeneration of P. densiflora saplings. In fact, a study on Picea abies by Kathke and Bruelheide (2010) also found that gap age itself did not significantly affect regeneration density; however, a clear difference was confirmed in analyses that included microsites. Therefore, although canopy openness alone could not explain the differences between sites in this study, the observed variance in P. densiflora regeneration may be the result of an interaction effect between microsites (e.g., environments with high rock exposure) and gap size. To closely verify these effects, continuous and precise monitoring is required, and it is necessary to strengthen the interpretation and verification of dynamic changes in declining stands through comparative studies across various regions.

We compared changes in the interactions and network structure between microbes and used an RF model to track differences in the composition of microbial communities between groups, thereby providing a theoretical basis for detecting the early signs of tree decline. In addition, we analyzed how microbial community composition is affected by changes in the degree of decay with time after tree death. This analysis addresses the limitations of aboveground observations and bridges the gap in existing knowledge regarding plant–microbe interactions and the resilience of forest ecosystems.

According to our results, healthy P. densiflora exhibited lower or similar bacterial and fungal diversity (Chao1 and H’) compared to declining or dead trees. This finding contrasts with studies on sub-alpine Abies koreana (Kang et al., 2024) and previous reports associating high diversity with enhanced pathogen suppression (Hu et al., 2016; Wagg et al., 2021). These discrepancies suggest that the relationship between microbial diversity and tree health is not universally linear but likely varies depending on host species and local environmental conditions. This observation supports the view that the abundance of specific functional groups is a more critical determinant of stress resistance than overall diversity (Hol et al., 2015; Kong et al., 2022). Crucially, our RF model results substantiate this perspective by identifying Oidiodendron spp. (Ascomycota) and Umbelopsis spp. (Mucoromycota) as key taxa significantly more abundant in healthy P. densiflora. Oidiodendron spp. are known ericoid mycorrhizal fungi that form symbiotic associations with plant roots and promote root growth (Perotto et al., 2012; Baba et al., 2021), while Umbelopsis spp. are reportedly effective at inhibiting phytopathogens near the roots and have positive effects on plant growth (Kolp et al., 2018; Wang et al., 2023). Therefore, forest managers should compare microbial communities between healthy and declining trees to identify taxa consistently associated with tree health. Monitoring these functionally important groups, rather than relying solely on overall diversity metrics, may provide clearer diagnostics of forest condition and more actionable guidance for management. Taken together, these rhizosphere microbial communities hold significant potential as early indicators for evaluating ecosystem resilience and informing conservation strategies (Mendes et al., 2013; Egidi et al., 2019).

In healthy P. densiflora, the rhizosphere microbial networks exhibited a relatively simple structure, characterized by fewer nodes and edges (particularly in fungi) but high modularity at both the whole-network and LCC levels. These results differ from a previous study on Pinus sylvestris growing in the southeast of Spain. Lasa et al. (2022) reported that the microbial networks of healthy P. sylvestris were more complex than those of declining individuals. In addition, several recent studies have proposed that compact networks with strong interactions are more favorable for maintaining plant health and suppressing pathogens than relatively simple networks (Tao et al., 2018; Yang et al., 2017; Fernández-González et al., 2020). Therefore, the simple network structure for healthy P. densiflora observed in our study contrasts previous studies and could be due to ecological differences, including tree species, local environment, or climatic conditions. Mendes et al. (2013) proposed that the microbial network structure of healthy plants is closely related to the concepts of core and minimal microbiomes. The core microbiome is hypothesized to maintain a minimal microbial community through functional duplication while effectively performing specific functions, such as maintaining plant health and suppressing pathogens (Mendes et al., 2011, 2013). Taken together, these findings suggest that, independent of the structural properties or forms of the network, P. densiflora regulates symbiotic relationships by selectively maintaining functionally important microbial groups (i.e., microbial communities associated with tree health), thus protecting itself from soil-mediated disease (Cook et al., 1995; Mendes et al., 2011, 2013). Moreover, after P. densiflora death, with the gradual opening of forest gaps, changes in stand productivity (Brumme, 1995; Prescott, 2002) and in the distribution of microclimatic conditions and substrate due to fallen trees (i.e., logs; Perreault et al., 2021) cannot be completely excluded. However, we did not observe statistically significant differences in the soil environment between the control and multiple-tree mortality sites (Supplementary Table S10). Even though large forest gaps were formed in the stands a long time after tree death (i.e., Group 1, over 15 years), microbial network complexity (particularly in fungal communities) resembled that of healthy P. densiflora. Thus, we did not find evidence to suggest that forest gaps and CWD had important effects on the rhizosphere soil environment.

Interestingly, compared with the Healthy group, networks in the Decline group were larger and more densely interconnected (increased nodes and edges) and showed a higher proportion of negative associations, while modularity decreased—particularly at the LCC level—indicating that the largest connected component became more intricately connected yet less compartmentalized. Similar to these characteristics, De Vries et al. (2018) showed that experimental drought made soil bacterial co-occurrence networks larger and more highly connected while reducing their modularity, and interpreted these changes as destabilising for the community. Theoretical work on microbial networks further suggests that community stability tends to be higher when interactions among community members are generally weak and the network is organized into modular, partially compartmentalized subcommunities (Coyte et al., 2015). In this context, the complex, densely connected yet weakly modular structure of the Decline network may be more closely linked to a transient, stress-related reorganization driven by competitive processes than to a stable equilibrium state. Hernandez et al. (2021) linked stress-induced reductions in network modularity and shifts in the balance between negative and positive associations to lower community stability, although in their case the relative contribution of negative associations decreased rather than increased. These similarities and contrasts suggest that such stress-related network responses may differ among host species and site conditions, highlighting the need for comparative studies across diverse tree species and regions.

Meanwhile, the microbial network in Group 2 exhibited a marked increase in both nodes and edges, characterized by high modularity. The changes in network structure could be due to exudates from mycorrhizal fungi increasing interactions between specific bacteria and fungi and altered community structure, as reported in a previous study (Hassani et al., 2018). As trees decline and die, organic matter is released through the degradation of wood by root mycorrhizal fungi (Ekblad et al., 2013). This released organic matter could promote various microbial interactions, resulting in a rapid increase in population sizes and the formation of complex networks (Ehleringer et al., 2000; Clemmensen et al., 2015). While our RF analysis indicated that certain fungal genera were relatively enriched during this stage, these findings should be viewed as exploratory associations rather than causal evidence, given the limitations of genus-level identification. Importantly, this elevated complexity was transient: by approximately 15 years after mortality, the networks reverted to a simple structure similar to that of healthy P. densiflora, with reduced numbers of nodes and edges, and whole-network level modularity tended to return to similar levels. These patterns appear to reflect a short-term reorganization of belowground interactions following mortality, possibly driven by shifts in resource availability and the transient expansion of specific microbial groups, with a subsequent trend toward community re-balancing over time. We therefore view changes in network connectivity as community-scale indicators of belowground dynamics rather than evidence for taxon-specific mechanisms. From a management perspective, tracking these temporal shifts in network complexity could offer a practical, complementary signal for assessing forest condition under climate-driven disturbance. However, to establish this approach as a diagnostic tool, future research is essential. Integrative, functionally resolved longitudinal studies are needed to establish causal links between network structure and ecosystem processes and to clarify the roles of taxa active during forest dynamics.

In our study, we found that, after the multiple-tree mortality of P. densiflora, both aboveground species composition and rhizosphere microbial network structure—particularly fungal communities—gradually tended to return toward conditions similar to those in healthy P. densiflora forests. Despite improved light availability as dead trees decayed and forest gaps expanded, the establishment of other tree species remained limited, and regeneration was consistently dominated by P. densiflora. In this chapter, we hypothesize that this regeneration process may have been driven by intimate interactions between the above- and belowground environments. Specifically, we propose that microbial communities, including mycorrhizal fungi, could have contributed to the re-establishment of P. densiflora, with resource sharing via Common Mycorrhizal Networks (CMNs) serving as one possible mechanism. In our chronosequence, RF analysis revealed that the mid-decay period (~7–10 years after death) was characterized by an increased representation of genera such as Talaromyces and Russula. Talaromyces spp. have been reported to include antifungal and plant-growth–promoting strains, whereas many Russula spp. form ectomycorrhizae that have been associated with improved host nutrient and water acquisition and increased stress tolerance (Bergemann et al., 2006; Dickie and Moyersoen, 2008; Wang et al., 2015; Filizola et al., 2019). Snags can act as reservoirs for these fungi and, together with the rhizosphere fungi of living trees, may form shared mycorrhizal networks that facilitate resource sharing between dead and living P. densiflora (Toju et al., 2013; Figueiredo et al., 2021). By ≥15 years after death (Group 1), fungal co-occurrence networks tended to resemble those in healthy P. densiflora stands at the whole-network level, exhibiting comparable complexity and modularity. We speculate that this convergence might stem from non-exclusive processes, such as a potential decrease in exudate-responsive taxa or the development of functional CMN links between snags and living trees. These patterns suggest that interactions between above- and belowground components may be linked to tree regeneration dynamics and might evolve as the decay of P. densiflora progresses. However, direct, site-specific evidence demonstrating CMN connections is required. Because our identifications are limited to the genus level, and some assignments are incomplete, we cannot assign functional roles at our sites. The current results alone are clearly limited in supporting this hypothesis and will need to be augmented by empirical research that can demonstrate the necessary conditions for the belowground connectivity between plants via the CMN (Figueriedo et al., 2021).

In addition, specific evidence is needed to confirm whether resources are transferred between plants directly via the CMN (Bever et al., 2010; Courty et al., 2010). Moreover, CMN formation does not guarantee that the effects of these interactions will always benefit plants, as the outcomes may include benefits and disadvantages depending on the circumstances, requiring further investigations (Bücking et al., 2016; Figueriedo et al., 2021). If we can develop a clearer understanding of these interactions through future research, it could help us to understand ecological resilience and processes of regeneration in coniferous forests that are declining due to environmental stress.