Ectomycorrhizal fungal network complexity determines soil multi-enzymatic activity

Soil functioning is intrinsically related to the structure of associated biological communities. This link is barely understood in the multispecies context of soil microbial communities, which often requires complex analytical approaches to discern structural and functional roles of microbial taxa inhabiting the soil. To investigate these ecological properties, we characterized the assembly and soil functioning contribution of ectomycorrhizal (ECM) fungal communities through co-occurrence network analysis. Co-occurrence networks were inferred from ECM root tips of Cistus albidus, Quercus faginea and Q. ilex on a regional scale, in Mediterranean mixed forests. Soil enzymatic activities related to carbon and nutrient cycling were also measured, and soil functionality outcomes related to ECM fungal network structure were evaluated on the community to taxon levels. Network complexity relied on habitat characteristics and seasonality, and it was linked to different dominant ECM fungal lineages across habitats. Soil enzymatic activities were habitat-dependent, driven by host plant identity and fungi with reduced structuring roles in the co-occurrence network (mainly within Thelephorales, Sebacinales and Pezizales). ECM fungal co-occurrence network structure and functioning were highly context-dependent, pointing to divergent regional fungal species pools according to their niche preferences. As increased network complexity was not related to greater soil functionality, functional redundancy might be operating in Mediterranean forest soils. The revealed differentiation between structural and functional roles of ECM fungi adds new insights into the understanding of soil fungal community assembly and its functionality in ecosystems.

1 Introduction

The structure of biological communities is a key element for predicting the potential role of biodiversity in ecosystem functioning (Wagg et al., 2021). The ways through which species are assembled, e.g., via environmental filters, disturbances or even demographic stochasticity, can lead their functions within an ecosystem (Mazel et al., 2018; Wang et al., 2022). Several studies have attempted to link the structure of biological communities to their functionality, such as those relating taxonomic or phylogenetic diversity to functional traits of species or to multifunctionality indices (Delgado-Baquerizo et al., 2016; Le Bagousse-Pinguet et al., 2019; Liu et al., 2022). Some of the most widespread tools to associate community metrics with ecosystem functions are based on correlations across variables or regression models that allow one to infer ecosystem function responses through community structure predictors, e.g., species richness or phylogenetic structure (Pérez-Izquierdo et al., 2017; Bastida et al., 2019; Krapu and Borsuk, 2020). However, most community structure metrics do not consider how the occurrences of species, and even their abundance, are affected by the presence or absence of other species and to what extent this species co-occurrence might explain ecosystem functions (Wagg et al., 2019). Within this context, network analysis has emerged as a promising analytical approach to studying multispecies co-occurrence and, hence, to evaluating assembly rules, niche differentiation and expected functionality of integrated species assemblages (Mandakovic et al., 2018; Zhan et al., 2021; Goberna and Verdú, 2022).

Co-occurrence network analysis is particularly useful on data obtained from next-generation sequencing at the level of operational taxonomic units (OTUs), species hypothesis (SH) or amplicon sequence variants (ASVs) as species equivalents (Barberán et al., 2012; Barroso-Bergadà et al., 2020). However, this approach is affected by biases, e.g., due to the sampling differences or intrinsic biotic and abiotic agents triggering indirect interactions that may result in spurious associations among species (Gloor et al., 2017; Barner et al., 2018; Barroso-Bergadà et al., 2020). To accurately infer species associations and overcome these constraints, several tools have been developed by applying distribution modeling or integrating correlation-based algorithms, such as the SParse InversE Covariance Estimation for Ecological Association Inference (SPIEC-EASI) (Kurtz et al., 2015; Morueta-Holme et al., 2016; Watts et al., 2019). These analytical frameworks have proven useful for characterizing environmental microbial communities, usually involving hyperdiverse communities with key implications for ecosystems (Pauvert et al., 2019; Wagg et al., 2019; Barroso-Bergadà et al., 2020; Luo et al., 2022).

Among soil microbial communities, ectomycorrhizal (ECM) fungi exhibit high taxonomic and functional diversity in forest ecosystems (Nguyen et al., 2020). ECM fungi mainly form symbiotic interactions in the roots of tree and shrub lineages and are known to promote positive plant–soil feedbacks (Kadowaki et al., 2018; Bahram et al., 2020; Liang et al., 2020; Segnitz et al., 2020). In the ECM symbiosis, the plant delivers photosynthetic carbon to the fungus in exchange for nutrients that the fungus obtains through the external mycelium spreading in the soil (Smith and Read, 2008). ECM fungi are capable of scavenging nutrients from soil organic matter by producing exoenzymes (Pritsch and Garbaye, 2011) which vary across ECM fungal lineages, promoting contrasting nutrient cycling and soil functional outcomes (Miyauchi et al., 2020; Lebreton et al., 2021; Prieto-Rubio et al., 2023). In addition, the C-enriched environment generated by the ECM symbiosis, i.e., mycorrhizosphere or hyphosphere (sensu van der Heijden et al., 2015), may harbor microorganisms (e.g., bacteria and saprotrophic fungi) with enzymatic capabilities to decompose organic matter, including ECM symbiotic structures and mycelia (Azcón-Aguilar and Barea, 2015; Cheeke et al., 2017; Lebreton et al., 2021; Nguyen, 2023). Hence, as the ECM environment seems key for carbon and nutrient cycling, understanding the structure of their co-occurrence networks might bring new insights into how ECM fungal species are assembled within the communities and how they contribute to explaining forest soil functioning (Barroso-Bergadà et al., 2020; Wagg et al., 2021, Martin and van der Heijden, 2024).

On the basis that the structure of ECM fungal communities is driven by habitat conditions and host plant community composition (Glassman et al., 2017; Pérez-Izquierdo et al., 2018; Wu et al., 2018), we firstly hypothesized that ECM fungal co-occurrence network structure would be dependent on host plant identity and spatial–temporal factors (habitat and seasonality). The key role of ECM fungal communities in soil functioning and soil biodiversity maintenance (Courty et al., 2010; Miyauchi et al., 2020) was the origin of our second hypothesis: the richer and/or more complex the ECM fungal co-occurrence network, the greater the soil functionality (i.e., carbon and nutrient cycling). This would be based on the associated fungal taxa with contrasting profiles of enzymatic activity leading to functional overdispersion as a result of niche segregation (Götzenberger et al., 2012; Wagg et al., 2019).

To test these hypotheses, we characterized the ECM fungal co-occurrence networks of Mediterranean mixed forests of southern Spain, at individual host plant levels from three representative plant species, and measured different potential soil extracellular enzymatic activities as a proxy of carbon and nutrient mobilization. We aimed to (i) infer the ECM fungal co-occurrence networks and study their structure and complexity across contrasting scales (forest, plot, host plant individual and ECM fungal taxon), (ii) determine the impact of ECM fungal co-occurrence network complexity on soil functionality and (iii) estimate the relationship between co-occurring ECM fungal taxa and the associated enzymatic activities in the soil.

2 Materials and methods

2.1 Study sites

This study was conducted in Mediterranean mixed forests located in the natural parks of Sierras de Cazorla, Segura y Las Villas (38.29° N, −2.57° W, referred to hereafter as Segura) and Monte la Sierra (37.64°N, −3.73° W, referred to hereafter as Jaén), in southern Spain. These protected areas have been largely monitored for plant–plant and plant–microbe interactions to better understand topics concerning Mediterranean forest community dynamics (from Siles et al., 2008, to Garrido et al., 2023, and Pajares-Murgó et al., 2023). With average altitudes of 1300 and 1000 m in Segura and Jaén, respectively, the climate in both regions is continental Mediterranean, showing annual mean temperatures of 13–15 °C and rainfall of 550–600 mm, with cool wet winters and summer droughts. Soils are characterized by calcareous limestone and low nutrient retention. The canopy layer mixes Quercus ilex L. and Q. faginea Lam. with Pinus halepensis Mill. in Jaén and with Q. pyrenaica Willd., Pinus nigra subsp. salzmanii J. F. Arnold and Pinus pinaster Ait. in Segura. At both sites, the understory includes Juniperus spp. and tree and shrub species of Cistaceae, Oleaceae, Rosaceae and Lamiaceae, among others (Pulgar et al., 2017; Alcántara et al., 2018). To compare ECM fungal community network functionality outcomes across the study sites, we selected those ECM plants that were common at both study sites (i.e., Cistus albidus L., Q. faginea and Q. ilex).

2.2 Experimental design and sampling

Four 50×50 m plots were selected per forest, and samplings were carried out in fall (2016) and spring (2017) (Alcántara et al., 2018). A total of 92 adult plant individuals were sampled, 32 for Q. ilex, 31 for Q. faginea and 29 for C. albidus (47 in Jaén and 45 in Segura, and 48 in fall and 44 in spring), with at least 10 m distance between the given plots. After removing the litter layer, we collected three subsamples ($\mathrm{10}\times \mathrm{10}\times \mathrm{20}$ cm holes) with secondary roots tracked from the target plant and surrounding soil and further combined into a single composite sample for each individual (Prieto-Rubio et al., 2022). Once in the lab, roots were separated from the surrounding soil and washed with tap water over 2 and 0.5 mm sieves, and all ECM root tips were collected under a Carl Zeiss Stemi 2000 stereomicroscope (Rincón et al., 2014). Before physical–chemical analyses, the collected soil in each sample was homogenized, air-dried and 2 mm sieved. The soil pH was measured in solutions (1:5, w:v in H₂O), the gravimetric moisture (GM) was calculated as the difference in soil weight before and after drying the samples at 105 °C for 48 h, and the soil organic matter (SOM) was determined by measuring weight loss on ignition at 400 °C for 4 h, following the methodology in Walkley and Black (1934) (see the details of the soil physical–chemical analyses in Prieto-Rubio et al., 2022).

2.3 ECM fungal community and soil multi-enzymatic activity

To frame the ECM fungal community in co-occurrence networks and relate it to potential soil enzymatic activities, we used the raw data (ITS1 sequences by Illumina MiSeq) of ECM fungal communities present at root tips and the enzymatic activities of soils collected in the rhizosphere of C. albidus, Q. ilex and Q. faginea, both at Jaén and Segura. The taxonomy of fungal amplicon sequence variants was assigned using the UNITE database v7.2 (Cole et al., 2014; Abarenkov et al., 2010) and clustered at 97% similarity into OTUs. After filtering the obtained fungal OTUs using the FUNGuild database v1.1 (Nguyen et al., 2016), the final output resulted in 449 OTUs with 6 582 941 reads associated with the ECM lifestyle. Two independent ECM abundance matrices, one per study forest, Jaén or Segura, were used for subsequent network analyses. Further details on molecular and bioinformatics analyses are provided as Supplement S1.

Nine potential extracellular enzymatic activities related to C, N and P cycling in forest soils were measured: β-glucosidase, cellobiohydrolase, β-xylosidase, β-glucuronidase and laccase related to C cycling; acid and alkaline phosphatase activities promoting P cycling; and chitinase and leucine-aminopeptidase activities as proxies of N cycling. Laccase was determined through a photometric assay and the rest through fluorogenic assays using a Victor microplate reader (Perkin-Elmer Life Sciences, Massachusetts, USA) (see details in Supplement S1). A multi-enzymatic activity index (MAI) was calculated from the enzymatic activities' dataset with a principal component analysis (PCA)-based procedure described in Meyer et al. (2018). Briefly, each enzymatic activity was standardized into z scores with the decostand function in the vegan R package (Oksanen et al., 2022). Then, a PCA was run with the prcomp function over the z-score dataset in the stats R package, and PC axes were multiplied by the proportion of variation that each captured (i.e., weighted PC scores). To ensure that higher scores corresponded to higher enzymatic activities, the signs of the scores on each PC axis were inspected and the whole axis was multiplied by −1 when necessary (reoriented PC scores, Meyer et al., 2018; Wagg et al., 2019). The MAI of each sample was calculated by summing its reoriented weighted PC scores.

2.4 Co-occurrence networks in ECM fungal communities

In our study, we followed the recommendations of Goberna and Verdú (2022) concerning the use of a co-occurrence concept for soil microbial communities, instead of other studies that assume direct association and/or interaction between the organisms. To infer the co-occurrence of ECM fungi within the communities, we adjusted the methodology as described by Wagg et al. (2019), which allows associations among pairs of OTUs to be sought and framed under the conditional independence concept (Kurtz et al., 2015). The OTU abundance matrix was rescaled by the proportion of minimum sequencing depth (39 221 sequences in Jaén and 4076 in Segura; see Prieto-Rubio et al., 2022, for more details on sequencing yielding), and OTUs that appeared <1 % in relative abundance were discarded. We obtained two curated datasets with 47 (Jaén) and 66 (Segura) ECM fungal OTUs that were used to infer the co-occurrence matrix of each ECM fungal community by applying the SPIEC-EASI algorithm with the spiec.easi function in the SpiecEasi R package (Kurtz et al., 2015). This is a robust correlation-based method for ecological network inference that differs from others in discarding spurious associations due to the lack of independence, i.e., correlations among OTUs that might be indirectly affected by other OTUs occurring in the network (Kurtz et al., 2015). The networks were constructed with the Meinshausen and Bühlmann (2006) neighborhood selection method. The minimum lambda ratio (i.e., a scalar tuning parameter in high-dimensional data) was 0.01, and 50 values of lambda for every 100 cross-validation permutations were fitted to evaluate the most parsimonious networks by using the StARS model selection criterion (Liu et al., 2010).

From the global adjacency co-occurrence matrices of each forest (spiec.easi objects), we extracted the presence–absence adjacency co-occurrence matrix per sample by pruning the global matrix to only include the OTUs in each sample. At the host plant individual level, we determined the OTU richness, the number of links (i.e., the number of co-occurrences detected among ECM fungal taxa in each sample) and the ratio number of links and OTU richness that were used as proxies of network complexity in each sample (Box 1). We then built weighted networks for the global adjacency matrices with the graph.adjacency function in the igraph R package (Csardi and Nepusz, 2006), with edge weight representing the co-occurrence frequency of ECM fungal OTUs. The adjacency matrices were used to graphically represent the networks with the Gephi software v.0.9.2 (Bastian et al., 2009). Both the Fruchterman–Reingold distribution and ForceAtlas2 attraction–repulsion algorithms were used in the simulations (Fruchterman and Reingold, 1991; Jacomy et al., 2014). We obtained the network graph and calculated different network centrality metrics at the ECM fungal OTU level (degree, eigenvector, closeness and betweenness centralities) and the hub score as analogous to a keystone taxa index (see Box 1). These metrics were used to measure the contribution of each OTU in structuring the co-occurrence network at site level. Network centrality metrics were performed through several algorithms in Gephi, all of them scaled from 0–1 (Box 1). Further, we determined the modules (i.e., sub-communities integrating the ECM fungal co-occurrence network per site) to which the OTUs belonged by using the algorithm described in Blondel et al. (2008) and implemented in Gephi.

Box 1 Ectomycorrhizal (ECM) fungal co-occurrence network parameters analyzed in this study. OTU: operational taxonomic unit.

2.5 Statistical analyses

2.5.1 ECM fungal network structure and habitat characteristics (hypothesis 1)

The response of network complexity to the factors studied was evaluated for each forest (Jaén and Segura) separately. Seasonality was fitted as a fixed factor and host plant identity nested in plots as a random factor by means of linear mixed-effects models (LMMs), using the lme4 R package (Bates et al., 2015). The significance of explanatory factors was tested with the Anova and rand functions from the lmerTest and car R packages, respectively. The variance explained of both fixed and random factors was evaluated with the r.squaredGLMM function in the MuMIN R package. In addition, the relationships of soil properties (log-transformed pH, SOM and GM) and network complexity parameters were assessed through Spearman correlation tests with Bonferroni adjustment using the psych and ggcorrplot R packages. The relationships between OTU-level network metrics were also evaluated by following the same correlation and adjustment tests.

To determine the importance of OTU abundance for network structuring (Carvalhais and Dennis, 2021), its effect on network centrality metrics (Box 1) was analyzed by means of permutational ANOVA (PERMANOVA) (adonis function, vegan R package) using Bray–Curtis dissimilarity with 999 permutations. Besides, to evaluate whether some modules were more central than others within the ECM fungal networks, we tested centrality metrics of OTUs with the corresponding module ID using PERMANOVA.

2.5.2 ECM fungal networks and soil multi-enzymatic activity (hypothesis 2)

To detect whether the structure of ECM fungal co-occurrence networks explained soil functionality (i.e., carbon and nutrient cycling), the response of the MAI was evaluated with LMMs by considering OTU richness, the number of links or network complexity (i.e., metrics determined at the host plant individual level) and seasonality as fixed factors. The response of detrended network complexity, which discards the intrinsic variation derived from OTU richness, was also evaluated with the same model and factors. The plant species nested in the plot was considered a random factor in all LMMs.

The ECM fungal OTUs potentially affecting each soil enzymatic activity were detected by the elastic net (ENET) regularization method (Zou and Hastie, 2005; Lassalle et al., 2019; Wagg et al., 2019). Briefly, the ENET regularization combines LASSO (least absolute shrinkage and selection operator) and Ridge penalty-based regression modeling, allowing us to avoid overfitting and the correlation effects between predictors (i.e., the OTUs embedded in the co-occurrence networks) and to minimize the influence of those that explain the variations in the response variables to a lesser extent (i.e., each of the soil enzymatic activities) (Zou and Hastie, 2005). As a previous step before model training, the abundance matrix of ECM fungal OTUs by host plants was standardized. Significant effects were determined by the standardized effect size | SES| >1.96 (i.e., the comparison of the observed effect of the given OTU with the simulated ones – 999 iterations – divided by the standard deviation of the simulated ones). From the ENET results, we recorded all the OTUs predicting enzymatic activities (i.e., those that contributed to predicting a given enzymatic activity) and determined the OTU richness that explained enzymatic activities per host plant individual. Further, we could determine the number of co-occurring ECM fungal OTUs that supported a common enzymatic activity (i.e., positive or negative fungal links) and the functional network complexity by dividing positive or negative fungal links by the total OTU richness (Box 1). In addition, as fungal network properties can be affected by the species or OTU number, especially when this decreases (Berry and Widder, 2014), the detrended network parameters by ECM fungal OTU richness allowed us to evaluate how the co-occurrence network structure yielded multi-enzymatic activity responses by discarding the variance intrinsically derived from the OTU richness. All functional network parameters recorded in each study forest (i.e., positive or negative OTU richness, positive or negative fungal links, and positive or negative complexity) and their detrended values for the total OTU richness were used in linear regression models (LMs) to explore their effects on MAI responses by fitting a separate LM for each functional parameter with the lm function in the stats R package.

3 Results

3.1 Structure of ECM fungal co-occurrence networks

Contrasting results in network parameters were observed between the study forests (Table 1). In Jaén, seasonality significantly affected ECM fungal network complexity, i.e., lower in spring than in fall (Table 1, Fig. S1 in the Supplement), while in Segura, this effect was detected in the OTU richness and the number of fungal links, in both cases higher in spring than in fall. Host plant species identity affected the detrended network complexity at the genus level in Segura (Table 1) (Cistus and Quercus spp. respective scores of 0.16 and 0.22, $F_{\mathrm{1},\mathrm{43}}=\mathrm{3.88}$, p=0.05).

Table 1 Effects of spatial–temporal factors and host plant species on overall ectomycorrhizal fungal network properties inferred from spiec.easi.

Linear mixed models were separately performed for each study forest: Jaén and Segura. Network properties (see Box 1) were log-transformed; seasonality was fitted as a fixed factor and host plant species nested in the plot as a random factor. The effects of fixed and random factors were determined by the F and χ² statistics, respectively, and the variance explained by both factor pools (R²). Significance level: “${}^{***}$” p<0.001, “${}^{**}$” p<0.01, “^*” p<0.05, “^.” p<0.10 and “ns” – nonsignificant (significant effects in bold).

When the forest was analyzed as a fixed factor, the ECM fungal network complexity was revealed to be significantly higher in Segura than in Jaén (Fig. S2). In fact, the network complexity and soil properties of each site were differentially correlated (Fig. S3). In general, soil pH was positively related to network parameters in Jaén and negatively in Segura, while higher SOM or GM values were associated with lower network complexity in both Jaén and Segura. The number of modules of the inferred co-occurrence networks varied between forests (Fig. 1, Table S1 in the Supplement). In Jaén, ECM fungal OTUs were distributed across 33 modules, 27 % of them composed of more than one OTU, and none of them was more central than the others (PERMANOVA, $F_{\mathrm{1},\mathrm{41}}=\mathrm{0.27}$, p=0.67, R²=0.01). By contrast, ECM fungal OTUs were distributed across 27 modules in Segura, 44 % of them with two or more OTUs, and some (e.g., modules 2, 18 or 22) were more central than the rest (PERMANOVA, $F_{\mathrm{1},\mathrm{60}}=\mathrm{5.31}$, p=0.02, R²=0.08) (Fig. 1, Table S1).

Figure 1 Co-occurrence networks of ectomycorrhizal (ECM) fungal communities in each study forest: Jaén and Segura. Graphs are inferred through the Gephi software by using the Fruchterman–Reingold and ForceAtlas algorithms. In each site network, each node number represents an ECM fungal OTU (see Table S1 for the taxonomic assignment and network metrics for each forest), and the node size is related to the hub score index (i.e., a proxy of keystone taxa). Co-occurring OTUs are connected by edges whose width indicates the co-occurrence strength (co-exclusions were not detected among ECM fungal taxa across the study sites). Within each forest ECM fungal community, a shared node color indicates a separate module (legends) inferred by Gephi (Table S1). Nodes in light-grey color indicate non-co-occurring OTUs in each forest network.

Negative correlations between the average abundance of ECM fungal OTUs and network centrality metrics (degree, closeness, eigenvector and betweenness centralities) were observed at both sites (Fig. 2) (PERMANOVA $F_{\mathrm{1},\mathrm{45}}=\mathrm{4.64}$, p=0.02, R²=0.09 in Jaén; $F_{\mathrm{1},\mathrm{64}}=\mathrm{8.23}$, p<0.001, R²=0.11 in Segura). Correlations between centrality metrics were all positive and overall greater in Jaén than in Segura (Fig. 2). In addition, hub scores (i.e., a proxy of keystone taxa; Box 1) were observed in specific Gephi modules, particularly in modules 4, 24 and 28 in Jaén and in modules 2, 18 and 22 in Segura (Fig. 1, Table S1). In Jaén, a total of seven hub OTUs were detected, most belonging to Thelephoraceae (six OTUs) (e.g., Tomentella with four) and one to Russulaceae (Table S1). However, in Segura, a widespread taxonomic diversity of hub OTUs (23 in total) across different fungal families and genera, e.g., Inocybe (7), Tomentella (5), Sebacina (4), Tuber (1), Tylospora (1), Amphinema (1), Pustularia (1), Humaria (1) and Hydnum (1), was observed (Table S1).

Figure 2 Correlations between average ECM fungal OTU abundance and network parameters (Box 1). Network metrics calculated for each ECM fungal OTU are inferred through the Gephi software. Correlation tests are calculated with the Spearman method using the Bonferroni adjustment for the dataset of each study forest: Jaén and Segura. Significant correlations (p<0.05) are plotted per site, indicating positive (red) or negative (blue) correlations.

3.2 Functionality of ectomycorrhizal co-occurrence networks

When functional ECM fungal co-occurrence network parameters were tested as explanatory variables of the multi-enzymatic activity MAI, a number of significant effects – positive OTU richness, detrended positive and negative OTU richness, functional network complexity – were observed, but only in Jaén (Table 2). In Segura, only positive OTUs were recorded, although they did not significantly affect the soil MAI (Table 2); no effect of the general ECM fungal network parameters (i.e., fungal links and network complexity) on the MAI was detected (Table S2).

Table 2 Effects of functional network parameters inferred from elastic net (ENET) regularization models on the soil multi-enzymatic activity index (MAI) per study forest: Jaén and Segura.

Notes: each network metric was fitted as a predictor of the soil MAI in linear models when recorded from the given forest. Detrended functional network metrics were also fitted as predictors. Significant values are highlighted in bold: “${}^{**}$” p<0.01, “^*” p<0.05, “${}^{***}$” p<0.10 and “ns” >0.10. Nr: not recorded.

The number of OTUs that explained at least one soil enzymatic activity differed greatly between the study forests, i.e., 76.6 % and 9.1 % of their respective total OTUs in Jaén and Segura, respectively (Fig. 3, Table S3). In both forests, ECM taxa that explained soil enzymatic activities showed low centrality scores in the co-occurrence networks (Tables S1, S3). Additionally, keystone taxa (i.e., those showing hub scoring, Fig. 1) also recorded low effects on the enzymatic activities (Fig. 3). In Jaén, 36 OTUs positively or negatively predicted soil enzymatic activities, particularly those related to C cycling (Fig. 3a, Table S3a). Most of these OTUs belonged to Thelephorales (23.40 %), Sebacinales (21.28 %) and Pezizales (12.77 %) and in many cases affected more than one enzymatic activity, as was also the case with the unique OTU belonging to Boletales (Fig. 3, Table S3a). Fungi within the Pezizales order showed a predominantly positive effect, particularly on enzymatic activities related to C (β-glucosidase and laccase) and P (alkaline phosphatase) cycling (Fig. 3, Table S3a).

Figure 3 Co-occurrence ECM fungal networks and contributions of fungal OTUs to soil extracellular enzymatic activity in Jaén and Segura. Co-occurring OTUs are marked with bold blue lines in the graphs and classed by fungal orders and keystone taxa with asterisks. Heatmaps show the estimated values (positive in red and negative in blue) of ECM fungal OTUs, explaining at least one soil enzymatic activity by elastic net (ENET) regularization models. The estimated values are 0.1-scaled. We used order as the main resolving taxonomic category to visualize trends of ECM fungi in predicting enzymatic activities within a network context and incorporating phylogenetic relationships. The ECM fungal OTU taxonomic assignment and estimated values are listed in Table S3.

In Segura, only six OTUs positively affected soil enzymatic activities (half of them either one or two enzymatic activities), particularly those related to C cycling (Fig. 3, Table S3b). The greater estimated values were detected for Peziza succosella and species belonging to Thelephoraceae, which positively explained β-glucosidase and β-xylosidase activities (Fig. 3, Table S3b).

4 Discussion

ECM fungal co-occurrence networks in Mediterranean forest soils reveal a strong habitat dependency, supporting our first hypothesis. Network metrics differ between the two study forests by season and also by the host plant identity (i.e., Cistus vs. Quercus spp.) in Segura. When investigating the ECM fungal co-occurring network properties further, the abundance of occurring taxa appears to be a key parameter for network structuring since, for a given ECM fungal OTU, centrality metrics diminish with increasing OTU abundance. The co-occurrence network structure does not have a significant effect on the soil multi-enzymatic activity MAI, but we do find relationships when functional network metrics are analyzed (e.g., OTUs and their co-occurrences explaining at least one soil enzymatic activity), partially supporting our second hypothesis. Our study reveals that fungal OTUs within a co-occurrence network context may explain soil enzymatic activities (mostly those related to carbon cycling) and that major effects are directed by taxa belonging to the orders Thelephorales, Sebacinales and Pezizales.

4.1 Ectomycorrhizal fungal networks differed between Mediterranean mixed forests

ECM fungal co-occurrence networks in the study on Mediterranean mixed forests showed a strong habitat dependency. Environmental factors such as summer drought, wildfire recurrence or soil nutrient retention have been identified as important long-term habitat filters affecting Mediterranean forest community adaptations (Pérez-Valera et al., 2017; Pérez-Izquierdo et al., 2020; Hernandez et al., 2021). These filters could contribute to the regional changes in fungal OTU richness and network parameters observed in this study. In addition, soil properties were related to the outcomes associated with co-occurrence network complexity. Co-occurrence networks revealed differential responses to soil pH ranging from 6.9 to 8.3 (averaging 7.7±0.1) and from 4.9 to 8.2 (averaging 6.6±0.2) in Jaén and Segura, respectively, suggesting local adaptation to homogeneous neutral–alkaline soils in Jaén, in contrast to those more variable in Segura, which also depended on the host plant species identity, e.g., higher soil pH under C. albidus than Quercus spp. (Prieto-Rubio et al., 2023). On the basis that pH is a well-known filter for soil microbial community assembly (Glassman et al., 2017; Tripathi et al., 2018), our results suggest that neutral pH contributed to diminished network complexity, showing an increased number of potential single-OTU modules (Jaén). By contrast, neutral–acidic conditions could imply a wider range of available niches for ECM fungi, which would explain the increasing fungal network complexity observed in Segura.

The seasonality effect on the network parameters varied in each study forest, suggesting that this effect is dependent on the habitat characteristics. This habitat dependency of seasonal outcomes associated with nutrient availability and fungal communities in forest soils has been shown previously (Siles and Margesin, 2017; Rincón et al., 2015; Ji et al., 2021). In fact, the altitudinal variations between Jaén and Segura (∼300 m) may imply variations in physical–chemical and SOM properties, as shown previously in other Mediterranean forest ecosystems (Gutiérrez-Girón et al., 2015; Centenaro et al., 2023). These altitudinal variations also span the limit between meso- and supra-Mediterranean bioclimatic belts (∼1200 m a.s.l.) and lead to changes in the physiological responses of host plant communities (Rivas-Martínez, 1987) and likely also in the ECM fungal co-occurrence networks. Further studies including a broader altitudinal zonation would be needed to confirm this causal effect on Mediterranean forest soil communities.

4.2 The dominance of certain ECM fungi was key in network structuring

The richness and abundance of ECM fungal taxa were key to understanding how co-occurrence networks were structured. As argued by Wagg et al. (2019), the network complexity detrended by OTU richness discards the intrinsic variation due to OTU richness. This parameter revealed a habitat effect: that is, in the low-complexity network of Jaén, detrending network metrics helped find fine-tuned community outcomes that would be obscured by the low OTU richness, e.g., detecting seasonal variations. By contrast, in the more complex network of Segura, those seasonal variations were not detected but showed a host plant effect differing at the genus level (Cistus vs. Quercus) (as previously observed when phylogenetic turnover outcomes of the ECM fungal communities were evaluated, Prieto-Rubio et al., 2022). The network composition also varied between the studied forests but revealed common fungal orders (Pezizales, Sebacinales and Thelephorales) that dominated the network in Jaén (∼80 % of the total OTUs), whilst others (e.g., Agaricales, Boletales or Russulales) rarely occurred. In Segura, the abundance of the common fungal orders found in Jaén decreased, allowing other fungal lineages to increase their network representativeness (e.g., OTUs belonging to Agaricales and Atheliales). This contrasting pattern among fungal lineages might suggest that ECM fungal communities conserved a regional species pool with lineages likely differing in their ecological niche space, as observed in arbuscular mycorrhizal fungal communities by Davison et al. (2021).

The negative correlation found between fungal abundance and centrality metrics suggests that local dominance plays a key role in ECM fungal network structuring (Lamanna et al., 2016; Prieto-Rubio et al., 2022). This property of ECM fungi is consistent with that observed in arbuscular mycorrhizal fungal communities (Davison et al., 2022), and its extent could also be influenced by the environmental context (Wang et al., 2022). In Jaén, the edaphic environment (neutral–alkaline soils with high organic matter content and homogeneous under all host plant species, Prieto-Rubio et al., 2023) could limit niche breadth for ECM fungi, likely leaving a more important role to dominance and explaining the observed lower community richness and network complexity (Davison et al., 2022), whilst in Segura, where soil properties were more differentiated across host plant species (Prieto-Rubio et al., 2022, 2023), the increasing niche breadth for ECM fungi could mitigate the role of dominance (Davison et al., 2022). This would facilitate the coexistence among more fungal taxa, hence explaining the higher observed community richness and network complexity in this forest (Chomicki et al., 2019). Moreover, the outcomes concerning keystone taxa could also be affected by dominance, as revealing that these OTUs were not locally abundant confirmed that other mechanisms besides competitive exclusion, such as those related to demographic stochasticity (e.g., priority effects), could impact the co-occurrence network structure (Kennedy et al., 2009; Debray et al., 2022). Further studies decoupling the role of niche processes from that of stochasticity in interaction or co-occurrence networks will bring novel understandings of fungal community assembly.

4.3 ECM fungal network structure predicted enzymatic activities in forest soils

Contrary to our expectations, the overall network properties barely explained the soil MAI. However, when functional network properties inferred from ENET regularization were tested, we found that an increase in the number of ECM fungal OTUs related positively to any soil enzymatic activity, leading to positive soil MAI values. These patterns were mainly observed in Jaén, where a greater number of ECM fungal taxa explained any soil enzymatic activity compared with those in Segura. The effects of ECM fungal networks on enzymatic activities configuring the soil MAI were also visualized through the co-occurrence networks. Across the study forests, different ECM fungal lineages could be linked to the enzymatic activities, as previously observed in soil microbial communities from different ecotones (Trivedi et al., 2016; Banerjee et al., 2018a). The number of ECM fungal OTUs explaining soil enzymatic activities substantially differed between the study forests, with a greater number of OTUs explaining enzymatic activities in Jaén than in Segura. In accordance with microbial networks in grasslands (Wagg et al., 2019; Luo et al., 2022), we observed that highly complex ECM fungal networks presented fewer taxa explaining enzymatic activity rates in forest soils, i.e., as if increasing network complexity would lead to functional redundancy, likely explained by the functional complementarity among OTUs conforming to the network (Wang et al., 2021). This result supporting the ideas of functional redundancy (i.e., enzymatic activity conservation across ECM fungi, Baldrian and Kohout, 2017) fits well into the holobiont theory (Zobel et al., 2024). Indeed, from the holobiont perspective, this redundancy ensures that essential functions carried out by microbiota are maintained even if one or more species are lost or disrupted, providing resilience against environmental changes, disturbances or the loss of specific microbial species (Vandenkoornhuyse et al., 2015). Alternatively, other soil microbial groups could be contributing to P and N mobilization to a much greater extent than ECM in these forest soils (Ward et al., 2021; Xun et al., 2021).

Furthermore, we found that the effects of ECM fungal taxa on soil enzymatic activities varied across ECM fungal lineages. In Jaén, positive or negative effects on soil enzymatic activities related to C mobilization were shown in taxa belonging to Agaricales, Atheliales, Cantharellales, Pezizales and Russulales. These ECM fungal lineages are known to conserve exoenzyme pools related to C cycling (Miyauchi et al., 2020). We could also show that OTUs belonging to Pezizales mainly drove P mobilization, whilst taxa contributing to explaining N cycling were few but phylogenetically widespread across ECM fungal lineages. Our results also show that keystone ECM fungal taxa were not major determinants of soil enzymatic activities, especially in Segura, where few effects were detected despite the higher richness and network complexity. This suggests that keystone ECM fungal taxa play community assembly-related roles rather than those associated with biogeochemical functions, as initially assumed for these microbial groups in soil networks (Banerjee et al., 2018b). These findings show the need in further studies to integrate other microbial guilds (e.g., saprotrophic bacteria and fungi) into co-occurrence networks to deepen the soil community attributes that explain carbon and nutrient mobilization.

5 Conclusions

This study shows the importance of integrating co-occurrence network analyses to disentangle ectomycorrhizal fungal community structure and function. We conclude that ECM co-occurrence network outcomes are highly dependent on habitat characteristics and that increased network complexity could imply functional redundancy in forest soils. The abundance of ECM fungal taxa is key to understanding the network structure, supporting the potential role of dominance in the ECM fungal community assembly in Mediterranean forests. In this co-occurrence network context, we found that several ECM fungal lineages explained soil functioning through enzymatic activities. In addition, keystone ECM fungal taxa were not necessarily those that profoundly explain soil enzymatic activities, showing a main influence on the network structuring. Incorporating other soil microbial guilds (e.g., saprotrophs, pathogens or arbuscular mycorrhizas) into co-occurrence network analyses will help to accurately discern structural and biogeochemical roles within forest soil communities.