This study was conducted in a 50-ha subtropical forest plot situated in the Heishiding nature reserve of southern China, which is located on the tropic of Cancer at 23°25′ ~ 23°29′N, 111°49′ ~ 111°55′E. Mean annual temperature and precipitation are 19.7°C and 1750 mm, respectively. The total area of the nature reserve is 4,200 ha, divided into core and experimental areas. The 50-ha forest plot was established within the 1,660 ha experimental area in 2012. All trees with a diameter at breast height (DBH) ≥ 1 cm were identified to species. In total, the plot includes about 269,000 stems of 213 woody plant species.

We analyzed an existing database of the plant–EM fungal network for the forest plot. Details regarding sampling design and taxonomic identification of root fungi (based on ITS) have been previously described (Wang et al., 2019), so we provide only a summary here. Within the 50-ha stem-mapped plot, EM fungi were quantified from 512 root samples of 43 plant species, which were collected randomly in the plot (Supplementary Figure S1). The 43 species were selected based on the phylogeny and their abundance. The 43 species included 3 species from each of Litsea and Lithocarpus (two of the most abundant genera in the plot) and 37 species from a larger set of genera where 6 species were from the same families as Litsea (Lithocarpus) and the others were from different families (Wang et al., 2019). 5–15 individuals for each plant species were randomly chosen for fine root sampling with at least three replicates taken from each individual by root tracing along different directions and pooled to generate a single sample. Root identity validation based on rbcLa fragment showed that 97% of fine-root samples of host plants were correctly traced. Root-associated fungi were identified by the internal transcribed spacer (ITS) region of fungal rDNA (details in Wang et al., 2019). The operational taxonomic units (OTUs) of root fungi were discriminated using a threshold of 97% sequence identity. EM fungi were identified based on a database of EM fungal taxa and lineages (Tedersoo and Smith, 2013). Our final database was a single quantitative interaction matrix of 862 EM fungal OTUs and 43 plant species. To account for the sampling inequality, each cell in each network matrix was filled with the mean of abundance (sequenced reads) of each OTU (species) of EM fungi on each host tree, and the number was rounded to the nearest integer.

We measured six metrics frequently used in the analysis of ecological networks: connectance, interaction evenness, interaction asymmetry (for EM fungi), nestedness, specialization and modularity. Connectance is the proportion of realized interspecific links, calculated as C = L / I J , where L is the number of nonzero cells in the binary interaction network and I and J are the numbers of plant and fungus species in the network. Interaction evenness was defined as Shannon’s index (Tylianakis et al., 2007), H = p i j log 2 p i j / log 2 F , where F is the total number of interactions in the whole network and p i j is the proportion of interactions linking plant species i and fungal OTU j. Interaction asymmetry for a given species i was calculated as A i = ∑ j d i j / k i , where k_i is the number of species interacting with i and d_ij is a measure of the symmetry of the pairwise interaction strength between i and j (Vázquez et al., 2007). Network specialization was calculated following Shannon entropy, reflecting degree of complementarity specialization of a whole plant-EM fungal network (Blüthgen et al., 2006).

We estimated the modularity of the plant–EM fungal association network using the DIRTLPAwb+ algorithm set to steps = 10⁶ searches (Beckett, 2016) in R-package bipartite (Dormann et al., 2009). Network modularity (M) ranges from 0 (low modularity) to 1 (high modularity). Nestedness was assessed using the weighted nested overlap and decreasing fill (WNODF; Almeida-Neto et al., 2008) in the R-package bipartite (Dormann et al., 2009). Network nestedness ranges from 0 (low nestedness) to 100 (perfect nestedness). To determine the significance of modularity and nestedness in observed network, we compared the observed network metric to those calculated from 1,000 null networks generated with the “swap” method (Artzy-randrup and Stone, 2005; Dormann et al., 2009). This non-sequential null model algorithm was selected because it maintains both marginal totals and connectance (number of links) in the null model networks consistent with the observed networks.

The abundance of each plant species was quantified as total number of individuals in the 50-ha forest plot. The abundance of each mycorrhizal fungal species (OTU) was calculated as the sum of the abundance of each fungal species (OTU) on each host tree (Wang et al., 2019), which was evaluated after subsampling each sample to 3,000 sequence reads to eliminate the effects of sample size (Wang et al., 2019). To quantify spatial overlap of plant and EM fungal species at the local scale, we compiled matrices of occurrence of each. We recorded the presence or absence of each species in each cell of a 20 × 20 m grid (hereafter quadrat) in the 50 ha plot. We thus obtained one spatial occurrence matrix for host plants and one for EM fungi, with species in rows and quadrat in the columns, and cells populated with ones for presences and zeros for absences.

We calculated interaction probability matrices to test whether network metrics were associated with neutral processes (e.g., relative species abundance) and niche/dispersal related processes (as implied by spatial overlap: presence of both species in a 20 × 20 m quadrat, which could be due to dispersal limitation or autocorrelated environmental conditions) or via the combined effect of both (mixed). To test the hypothesis that neutral processes dominated interactions we calculated a probability matrix (Ab) quantifying the relative interaction potential between a plant and a fungal species as the product of their relative abundances (Vázquez et al., 2009). To obtain interaction probability matrix expected from spatial processes (niche or dispersal based), we normalized the spatial overlap matrices Sp so that their elements sum to one (Vázquez et al., 2009). Thus, the probability of interaction increased with increasing spatial overlap of two species pairs, while those with no spatial overlap had zero probability of interaction. To model the combined effect of neutral and niche/dispersal related processes we calculated combined probabilities (AbSp) as the element-wise multiplication of matrices Ab and Sp, again normalizing the resulting matrices (Vázquez et al., 2009). These combined matrices represent the expected probability under the joint influence of species abundance and spatial overlap (analogous to an interaction term in regression modeling). Thus, we had three probability matrices quantifying relative abundance, spatial overlap and their combined effect: Ab, Sp and AbSp. As a null model, we also defined an equiprobable probability matrix (Null). In this matrix, all pairwise interactions had the same probability of occurrence, where each element was calculated as the inverse of the product of the numbers of plant and EM fungal species in the network (Vázquez et al., 2009). To evaluate whether these interaction probabilities inferred from species abundance, spatial overlap and null model matched the frequencies of interaction in observed networks, we also constructed an observed probability matrix (Obs), where pairwise interaction probabilities were calculated as the ratio of each element and the sum of all elements in the realized network matrix.

To evaluate the ability of interaction probability matrices to predict the network metrics (connectance, evenness, specialization, modularity, nestedness and EM fungal interaction asymmetry) we used a randomization algorithm (Vázquez et al., 2009). The algorithm assigned the total number of interactions originally observed in the network according to the three probability matrices, with the only constraint being that at least one interaction was assigned for each species (Vázquez et al., 2009).

To test the relative support for each hypothesis in explaining the pairwise interaction frequencies in the observed plant–EM fungal network (Vázquez et al., 2009), we used a maximum likelihood approach. The likelihood that a probability matrix explained the observed matrix was calculated assuming that the pairwise probability of interaction between a plant and an EM fungus followed a multinomial distribution (Vázquez et al., 2009). This likelihood was calculated using function dmultinom in the R-package stats. The comparisons between models were conducted using differential Akaike’s information criterion ( Δ AIC; Burnham and Anderson, 2004). Δ AIC was the difference in AIC for a probability matrix (Ab, Sp, AbSp and Null) with that for the best fitting model. The model with the lowest Δ AIC indicated the model providing the best fit to the observed data (Burnham and Anderson, 2004). In addition, to provide an effective way to scale and interpret the Δ AIC, we calculated the Akaike weights (w_i; Burnham and Anderson, 2004). The w_i of each model was calculated as

where R was the number of all models. The value of w_i ranges between 0 and 1, and the sum of all model weights is 1. The greater the weight of a model, the greater the support for that model as the best candidate model in the model set.

We evaluated the influence of evolutionary history on network patterns using a phylogenetic tree of plant species in a local forest plot reconstructed by four general plant DNA barcodes (rbcLa, matK, trnL and ITS2; Supplementary Figure S2) and a phylogenetic tree of EM fungal species (Supplementary Figure S3) reconstructed using ITS sequences (details for the phylogenetic reconstruction are given in Supplementary Method S1). We evaluated the impact of phylogenetic history on the assemblage of interacting partners of the species quantified using the Bray–Curtis dissimilarity index and using Mantel tests to compare the correlation between the phylogenetic distance matrix of plant/EM fungi and the community similarity matrix of their interacting partners.

We tested the phylogenetic signal on species degree (number of links for each species) and species strength (the importance of a species to its partner’s set taking into account the relative abundance of the species on each partner; Barrat et al., 2004). Two of the most widely used statistics, Pagel’s λ (Pagel, 1999) and Blomberg’s K (Blomberg et al., 2003), were conducted to test the phylogenetic signal in the two network indices. The two statistical tests were implemented within the phylosig function in the phytools R-package (Revell, 2012). Both indices assumed the classic Brownian motion (BM) evolutionary model, and their values vary from 0 to 1 for λ and from 0 to > > 1 for K. In both cases, low values (close to 0) indicate the absence of phylogenetic signal (i.e., the trait has evolved independently of phylogeny); high values (close to 1) indicate trait evolution according to BM; and, in the case of K, values >1 suggest that close relatives are more similar than expected under BM. The statistical significance of K was assessed based on a comparison of the observed phylogenetically independent contrasts and the expected contrast under 999 randomizations (Blomberg et al., 2003), whereas the statistical significance of λ was assessed based on a comparison of the likelihood of a model accounting for the observed λ with the likelihood of a model that assumes complete phylogenetic independence (Pagel, 1999).

The observed fungal network registered 4,360 links between the 862 EM fungal species (23 genera) and 43 woody plant species. We found that Cenococcum spp. and Russula spp. were the most important EM species in terms of their co-occurrences with roots of a high number of plant species (Figure 1). Cenococcum spp. and Russula spp. were each symbiont on 43 host plant species. Among plants, Ternstroemia gymnanthera, Itea chinensis and Manglietia fordiana were the species with the most interactions, associating with the roots from 23, 18 and 18 EM genera (Figure 1, 193, 147 and 92 species included), respectively. The observed plant–EM fungal network had significant strong nestedness [NODF = 10.256, CI for null model = (8.240, 8.561), p < 0.01]. For instance, most EM species habited on the root systems of a subset of the host plants of Cenococcum spp. and Russula spp. (Supplementary Figure S4), while most plant species established co-occurrence associations with a subset of the EM fungal symbionts of Ternstroemia gymnanthera, Itea chinensis and Manglietia fordiana (Supplementary Figure S4). The observed plant-EM fungal network exhibited significantly strong modular structure [M = 0.482, CI for null model = (0.058, 0.061), p < 0.01], with 14 modules identified (Supplementary Figure S5). The network was characterized by an intermediate level of specialization (H2′ = 0.42), low connectance (0.12) and high interaction evenness (0.74), indicating structure was not dominated by a small number of interactions. The average strength of EM fungi with their host plants was close to symmetry (mean asymmetry for EM = 0.17).

The top-ranked model for predicting species interaction frequencies was clearly that including the only spatial overlap of species (lowest Δ AIC and highest wAIC for the Sp model, Table 1) with all other models (Ab, AbSp and null) exceeding this by more than 100 AIC (Table 1). However, model performance differed among the network structural metrics analyzed (Figure 2). The AbSp model including both spatial overlap and species abundance performed best overall because the 95% confidence intervals of four out of six metrics (connectance, nestedness, evenness and asymmetry of EM) predicted by the Absp were more close to their observed values (Figures 2A–C,F), particularly for nestedness (Figure 2F). The AbSp model also provided a reasonable prediction for network connectance and was marginally better than the pure abundance model (Figure 2B). Although the AbSp model provided the closest prediction for interaction asymmetry for EM fungi (Figure 2A) and network evenness (Figure 2C), it was poor in absolute terms. None of the models for network specificity and modularity was an improvement on the null model (Figures 2D, E).

Phylogenetic distances of plant species were significantly negatively associated with ecological similarity of the fungal community (one-tailed Mantel test r = −0.086, p = 0.025) and phylogenetic distances of fungal species were also significantly negatively related to ecological similarity of the plant community (one-tailed Mantel test r = −0.023, p = 0.001). We found no evidence of a phylogenetic signal on species degree or species strength for plants (Table 2), while for EM fungi, results were metric-specific; Pagel’s λ suggested a significant phylogenetic signal in both degree and species strength (both p < 0.001) but this was not evident using Blomberg’s K (both p > 0.4, Table 2).

The local plant–EM fungal network in this subtropical forest was characterized by higher species number and more interactions than comparable studies (Jacquemyn et al., 2011; Montesinos-Navarro et al., 2012; Toju et al., 2014). Consistent with other mutualistic plant-fungal networks (Jacquemyn et al., 2011; Montesinos-Navarro et al., 2012; Toju et al., 2014), our network had low connectivity (Montesinos-Navarro et al., 2012). Compared with low levels of specialization (Roy-Bolduc et al., 2016), more intermediate levels of specialization as observed here (H2′= 0.42) in the plant-EM fungal network may be advantageous for host trees by increasing their chance to form a symbiosis with suitable EM fungal partners. Nonrandom associations of plants with symbiotic partners are common in biotrophic fungi, which may enhance modular structure in our plant-EM fungal network (Supplementary Figure S5). This is an important mechanism that leads to niche partitioning (reflecting ecological specialization; Jacquemyn et al., 2015). This implies that modularity and specialization are possible mechanisms allowing many EM fungal species to coexist locally on the same set of host plants.

We found that nested structures in the plant-EM fungi network where host plants of specialized EM fungi were subsets of the host plants of generalized Cenococcum spp. and Russula spp. fungi (Supplementary Figure S4). In contrast, anti-nested structures, in which co-occurring plant species share fewer fungal symbionts than expected by chance, have been reported in many smaller plant-EM fungal networks (Bahram et al., 2014; Roy-Bolduc et al., 2016; Arraiano-Castilho et al., 2020). This suggests that increasing numbers of species might favor greater network nestedness, as suggested by Bascompte et al. (2003). As the sampling effort invested in our study over a 50-ha forest extent exceeds many other studies in plant–root fungal association networks (Jacquemyn et al., 2011; Montesinos-Navarro et al., 2012; Toju et al., 2014), we believe the sampling effort is likely to have weak if any, influence on our estimates of network metrics (see Supplementary Figure S6).

Spatial co-occurrence patterns of plant and fungal species determined species interaction frequencies while species abundance or the interaction of the two had no support in our study (Table 1). Some evidence suggests that tropical mycorrhizal populations show significant spatial heterogeneity and nonrandom associations with different hosts (Husband et al., 2002). This clarifies how species interactions may be mainly determined by plant partner choice for EM fungi, which is also controlled by soil environment heterogeneity or spatial proximity (Cavender-Bares et al., 2009; Wang et al., 2019). Moreover, locally unmeasured microsite effects may confound the influence of the spatial distribution on species interaction patterns. In summary, our findings show the role of niche/dispersal-related processes (as implied by spatial overlap) in organizing symbiotic patterns of plant-EM fungal communities. This in turn points to the importance of co-occurrence as regulated by the effects of belowground environmental filtering in generating the observed links between plants and EM fungi.

Though species abundance makes only a weak contribution to the species interaction frequency of plants and EM fungi, we found that species abundance matching and spatial constraints made the highest contribution to most network metrics (nestedness included, Figures 2A–C,F). The mechanisms underlying the effect of species abundance on network structure (e.g., nestedness) could be explained by the regulation of right-skewed frequency distributions of relative species abundance on co-occurrence patterns (Canard et al., 2014). Moreover, spatial autocorrelation over short distances (in a range of 25 m, see Supplementary Figure S7) and limited dispersal capacity (Tedersoo et al., 2011; Bahram et al., 2013) may lead to aggregated distributions of EM fungi. In addition, the nonrandom distribution of host plants is also dependent to some extent on the occurrence of suitable mycorrhizal fungi (Jacquemyn et al., 2013). This may explain how spatial overlap of plants and fungi can predict network assembly (Jacquemyn et al., 2015).

However, co-occurrences of plant-EM fungi in a microsite were only a prerequisite to establishing symbiotic relationships. For example, the most dominant EM species that co-occur at a given study site are associated with distinct plants with scattered distributions (Tedersoo et al., 2010). This suggests that other biological constraints such as trait-matching of plant-fungal networks might contribute to the observed network structure (i.e., network specificity and modularity, Figures 2D,E). For example, some modules were dominated by Scleroderma spp. fungi with long-distance exploration while other modules were dominated by Russula spp. fungi with contact hyphal exploration (Supplementary Figure S5 and Zhu et al., unpublished 2022). Plant trait-based host selection behavior in EM fungi or sampling effects associated with species abundance and richness of the plant taxa may contribute to observed network specificity and modularity. Besides this, other deterministic factors such as phenology (Donatti et al., 2011; Morente-López et al., 2018) and a combination of local adaptation and competition (Valverde et al., 2020) can affect network modularity. Observed network specialization is poorly predicted by abundance, spatial overlap and plant traits (Sáyago et al., 2013), consistent with our study (Figure 2D).

Our results suggest that species abundance and spatial distribution can better predict the network structures that are regulated by both neutral processes (e.g., relative species abundance) and niche/dispersal related processes (as implied by spatial overlap), while failing in other network structural properties, including modularity, that may be mostly driven by trait-based niche processes (Donatti et al., 2011; Morente-López et al., 2018).

It has been shown that species interaction patterns of plants and root fungi are mainly constrained by the phylogeny of the host plant, including from the study plot (Toju et al., 2016; Wang et al., 2019) and weakly constrained by fungal phylogeny (Munoz et al., 2012). Similarly, we found that plant-fungal interactions were strongly associated with plant phylogeny (one-tailed Mantel test r = −0.086, p = 0.025) but weakly associated with fungal phylogeny (one-tailed Mantel test r = −0.023, p = 0.001). A lack of phylogenetic signal in species degree and interaction strength for fungal hosts (Table 2) may be related to the high vulnerability of fungal hosts and only distant relatedness between the multiple EM fungal partners supported. In contrast, phylogenetic constraints are evident for EM fungi on network characteristics, such as network nestedness (Chen et al., 2017) and compartmentalization (Vacher et al., 2008). However, the phylogenetic conservatism for EM fungi in species degree and interaction strength was only supported by λ statistics but not K statistics (Table 2). This suggests K statistics might underestimate the phylogenetic signal because it is sensitive to the number of missing plant taxa (Jardim et al., 2021).

Overall, a significant but weak phylogenetic signal on plant-EM fungal symbiotic associations suggests that the contribution of conservative traits on the network structure thus may be limited, while traits unconstrained by phylogenetic history may play important roles in regulating plant-fungal interaction patterns. Additionally, evolutionary history may influence network structure at other scales and analyzing a phylogenetically more diverse root fungal community or a larger network might lead to the consistent detection of a phylogenetic signal.