Ecological drivers of fine-scale distribution of arbuscular mycorrhizal fungi in a semiarid Mediterranean scrubland

Abstract Background and Aims Arbuscular mycorrhizal (AM) fungi enhance the uptake of water and minerals by the plant hosts, alleviating plant stress. Therefore, AM fungal–plant interactions are particularly important in drylands and other stressful ecosystems. We aimed to determine the combined and independent effects of above- and below-ground plant community attributes (i.e. diversity and composition), soil heterogeneity and spatial covariates on the spatial structure of the AM fungal communities in a semiarid Mediterranean scrubland. Furthermore, we evaluated how the phylogenetic relatedness of both plants and AM fungi shapes these symbiotic relationships. Methods We characterized the composition and diversity of AM fungal and plant communities in a dry Mediterranean scrubland taxonomically and phylogenetically, using DNA metabarcoding and a spatially explicit sampling design at the plant neighbourhood scale. Key Results The above- and below-ground plant community attributes, soil physicochemical properties and spatial variables explained unique fractions of AM fungal diversity and composition. Mainly, variations in plant composition affected the AM fungal composition and diversity. Our results also showed that particular AM fungal taxa tended to be associated with closely related plant species, suggesting the existence of a phylogenetic signal. Although soil texture, fertility and pH affected AM fungal community assembly, spatial factors had a greater influence on AM fungal community composition and diversity than soil physicochemical properties. Conclusions Our results highlight that the more easily accessible above-ground vegetation is a reliable indicator of the linkages between plant roots and AM fungi. We also emphasize the importance of soil physicochemical properties in addition to below-ground plant information, while accounting for the phylogenetic relationships of both plants and fungi, because these factors improve our ability to predict the relationships between AM fungal and plant communities.


INTRODUCTION
Arbuscular mycorrhizal (AM) fungi are members of the subphylum Glomeromycotina (Spatafora et al., 2016), which form mutualistic associations with >70 % of plant species (Brundrett and Tedersoo, 2018). AM fungi enhance the uptake of water and mineral nutrients by the hosts (Lambers et al., 2006;Smith and Read, 2010). In return, plant hosts provide carbohydrates to their symbiotic partners (Kiers et al., 2011). Thus, AM interactions play key roles in ecosystem functioning and the provision of services, such as plant growth and carbon sequestration (Van Der Heijden et al., 1998;Zhu and Miller, 2003;Wilson et al., 2009).
Although it is generally assumed that AM fungi-plant associations are non-specific (Brundrett, 2002;Helgason et al., 2007;Öpik et al., 2009), some studies have reported the existence of host plant preferences (Vandenkoornhuyse et al., 2002;Torrecillas et al., 2012). Such preferences would lead to non-random relationships between the diversity and composition of AM fungal communities and the distribution of plant roots in the soil (García de León et al., 2016;Xu et al., 2016;Neuenkamp et al., 2018;Šmilauer et al., 2020). Root sampling at the individual plant level has evidenced the complexity and variability of AM fungi-plant associations (Davison et al., 2015;Šmilauer et al., 2020), revealing a highly heterogeneous distribution of AM fungal communities. The heterogeneous spatial structure of AM fungal communities might be even more obvious in patchy environments (Davison et al., 2012), such as drylands, where a discontinuous distribution of the vegetation is a common feature (Aguiar and Sala, 1999;Maestre and Cortina, 2005). A better understanding of how plant community attributes (including diversity and composition) might affect the spatial structure of the AM fungal communities is particularly important in dry ecosystems because water accessibility and nutrient absorption by plants are often limited, and the AM fungal-plant associations might play a key role alleviating plant stress (Sánchez-Blanco et al., 2004;Barea et al., 2005;Khalvati et al., 2005;Allen, 2007).
Given that vegetation is more easily accessible above-ground than below-ground, most previous studies assessing the AM fungal-plant relationships at the community level have focused on the distribution of vegetation above-ground, although direct interactions take place below-ground. However, a decoupling between above-and below-ground plant community attributes has been observed in different environments (e.g. Jones et al., 2011;Träger et al., 2019;Illuminati et al., 2021), related to different plant community dynamics acting in the two compartments, such as the ability of many perennial plants to persist for long periods below-ground (Klimešová and Klimeš, 2007;Reintal et al., 2010). In this context, an obvious asymmetry between compartments is expected especially in water-limited ecosystems, where a greater root allocation is the norm (Schenk and Jackson, 2002;Mokany et al., 2006). The ability to detect AM fungal-plant associations might therefore be constrained if only plant above-ground data are used, because they might not accurately depict the below-ground plant community. Recent advances in the use of DNA metabarcoding to determine root diversity (Matesanz et al., 2019;Cabal et al., 2021) offer a means to improve the assessment of the interactions between plant roots and AM fungal communities.
Empirical studies have suggested that the relationships between AM fungi and plants above-ground might differ from those established with roots (e.g. Weigelt et al., 2021;Xia et al., 2021). Plant species vary in their ability to aerate the soil, in the physical changes induced by root growth or in the modifications of temperature and moisture under their canopy (Breshears et al., 1998). For example, AM fungal composition in forest roots can be affected by canopy-mediated light availability (Koorem et al., 2017). Moreover, plants and AM fungi can respond in a similar manner to soil physicochemical properties (Jamiołkowska et al., 2018). This can lead to misleading inferences about the importance of plant community attributes in generating patterns of AM fungal diversity. Conversely, opposite responses of plants and fungi to soil physicochemical properties could mask realized associations between plants and AM fungi. Thus, simultaneous assessment of the combined and independent effects of aboveand below-ground plant communities, together with their shared variation with the soil physicochemical properties, on AM fungal diversity is needed to provide a more complete picture of the complex nature of AM associations.
A useful approach that might help to formulate hypotheses related to the abiotic and biotic factors shaping heterogeneous AM fungal distributions is the identification of phylogenetic patterns in AM fungal community composition. The phylogenetic similarity of co-occurring AM fungal taxa has been used previously not only to infer community assembly processes (Liu et al., 2015;Vályi et al., 2016) or to identify key plant determinants of AM fungal diversity (Montesinos-Navarro et al., 2015;Valverde-Barrantes et al., 2016;Chen et al., 2017;López-García et al., 2017), but also to disentangle the existence of co-evolutionary mechanisms, such as exudate release (Brundrett, 2002;Brundrett and Tedersoo, 2018). This phylogenetic approach is based on the idea that closely related species tend to respond to their biotic and abiotic environments in a more similar manner than species selected at random (i.e. they exhibit phylogenetic signal in their ecological interactions). Therefore, species with shared evolutionary history are expected to be ecologically and functionally more similar than more distant relatives (Harvey and Pagel, 1991;Losos, 2008). For example, AM fungi exhibit a strong phylogenetic signal in hyphal growth, root colonization and spore size (Powell et al., 2009;Koch et al., 2017;Aguilar-Trigueros et al., 2019). A recent study showed that the phylogenetic structure of AM fungi can be explained by plant functional traits (López-García et al., 2017), which might also be phylogenetically structured (Valverde-Barrantes et al., 2016;Xia et al., 2021). Furthermore, variations in the abiotic environment, such as soil fertility or pH, can also drive changes in the phylogenetic structure of AM fungal communities (Liu et al., 2015;Davison et al., 2021;Li et al., 2021). Thus, a detailed analysis of the phylogenetic structure of plant and AM fungal communities might shed light on key aspects of these symbiotic relationships interacting with the abiotic environment.
Here, we assess the role of above-ground and below-ground plant community attributes, mainly composition and diversity, and soil physicochemical properties in predicting the spatial variation in the community structure of AM fungi considering the scale at which plant-plant interactions occur in a semiarid Mediterranean scrubland. To achieve this goal, we combine spatially explicit sampling of standing vegetation with nextgeneration sequencing to unravel the distributions of both roots and AM fungal communities, while accounting for variations in soil physiochemical properties. In water-limited ecosystems, the decoupling between above-and below-ground diversities is exacerbated at fine scales (Schenk and Jackson, 2002;Hiiesalu et al., 2012;Illuminati et al., 2021). Thus, Mediterranean scrublands provide an ideal model in which to test the relative importance of above-and below-ground plant community attributes as drivers of AM fungal communities. In a previous study in the same system, López-Angulo et al. (2020) found that the soil generalist fungal diversity and composition were better explained by above-ground plant community attributes compared with below-ground plant community attributes, but information on AM fungal communities and the influence of phylogeny on community assembly is still lacking. Given that above-and below-ground plant communities might drive different processes structuring AM fungi-plant associations, we hypothesize that above-and below-ground plant community attributes will jointly, but also independently, explain patterns of diversity and composition of AM fungal communities. More particularly, we ask whether assessing the more easily accessible above-ground vegetation is a reliable proxy of the linkages between plant roots and AM fungi. In addition, we test the hypothesis that phylogenetically closely related plant species might harbour more closely related AM fungal communities.

Study area and soil sampling
The study was performed in a dry Mediterranean scrubland located 50 km south-east of Madrid in central Spain (40°16ʹ08.5″N, 3°08ʹ11.1″W; 781 m a.s.l.). The mean temperature is 12.8 °C and mean annual precipitation 452 mm, with almost no precipitation during summer. The soil is calcareous and classified as Xeric Calcigypsids (Soil Survey Staff, 2014). The area is densely vegetated by dwarf scrubs dominated by Thymus vulgaris, Bupleurum fruticescens and Helianthemum cinereum interspersed with perennial grasses, such as Stipa pennata, Avenula bromoides and Koeleria vallesiana, and occasionally, with Quercus coccifera and Quercus ilex subsp. ballota.
We mapped all the above-ground perennial individuals within an 8 m × 8 m plot established in a representative area of the scrubland ( Fig. 1) in May 2016. The plot size guaranteed the inclusion of a high number of individuals (8551) and species (45 perennials) belonging to a wide range of plant families (18) (Supplementary Data Table S1). We located each individual plant by recording its rooting point, except for rosette plants, perennial herbs and tussocks, for which rooting points were assigned to the centroid of the plant, using a real-time GPS (Viva GS15; Leica, Wetzlar, Germany; absolute precision 1 cm). For each individual plant, we measured the longest diameter of its crown and the perpendicular diameter to estimate the projection area of the crown. After mapping, we established an 8 × 8 regular grid within the plot, resulting in 64 quadrats (1 m 2 ) and set 64 soil sampling points in the centre of each quadrat (Fig. 1). We set 20 additional soil sampling points close to the corners of the plot to increase the spatial resolution at finer scales (Supplementary Data Fig. S1), resulting in 0.70 m spacing between sampling points in these areas (84 soil sampling points in total). Centred around each soil sampling point, we delimited an above-ground sampling circle of 20 cm radius (84 above-ground sampling circles). For each species found within each above-ground sampling circle, we used the sum of all the intersection areas (in centimetres squared) between the projection area of the crowns of the individuals and the sampling circle as an estimation of above-ground species cover. We selected 20 cm because this was the radius that maximized the similarity in diversity between below-and above-ground plant assemblages in the same plant community (Illuminati et al., 2021).
At each soil sampling point, two contiguous soil samples were collected using steel cores (10 cm in depth by 5 cm in diameter). One soil sample was used to evaluate AM fungal diversity and soil physicochemical variables (hereafter, AMF-soil sample), and the other one to assess the below-ground plant community (hereafter, root sample). We sieved and homogenized the soil of the AMF-soil sample through a 2 mm mesh. We stored a 0.1 g subsample of each homogenized sample for molecular analyses of AM fungal diversity (see below in Arbuscular mycorrhizal fungal characterization section) and a 50 g subsample of airdried soil for soil physicochemical analyses. Specifically, we determined four variables related to nutrient stocks (organic carbon, total nitrogen, available phosphorus and potassium), two dynamic variables related to the soil microbial activity (acid phosphatase and β-glucosidase enzymatic activities) and several physical variables (pH, electrical conductivity, sand, silt and clay contents; for summary statistics, see Supplementary Data Table S2). In parallel, all plant roots from the root samples were filtered through a 1 mm mesh, thoroughly washed within 48 h after field collection and then centrifuged at 3000g for 30 s to remove excess water. We weighed the fresh root biomass and homogenized it by cutting roots into small pieces. We stored 0.1 g of root biomass per sample at −80 °C for subsequent DNA metabarcoding.

Arbuscular mycorrhizal fungal characterization
DNA isolation DNA from the AMF-soil samples was isolated using the DNeasy PowerSoil isolation kit (Qiagen, CA, USA) following the manufacturer's protocol. An extraction blank was included in each of the seven DNA extraction rounds to check for cross-contamination.
Taxonomic assignment We performed taxonomic assignment of the AM fungal OTU by querying the clustered centroids against the MaarjAM reference database (Öpik et al., 2010; accessed 3 October 2018), using the script 'assign_taxonomy.py' implemented in Qiime and the BLAST algorithm with a maximum E-value of 1 × 10 −50 and a minimum percentage identity of 90 % in order to assign the Glomeromycota sequences to an OUT (Morgan and Egerton-Warburton, 2017). An OTU table with the number of sequences of each OTU in each sample was created. OTU tables were subjected to quality filtering to remove OTUs with low frequency in the whole dataset (0.005 %; Bokulich et al., 2013) and low abundance (0.1 % threshold; Esling et al., 2015). Ten samples with no sequences and one sample with four sequences after all filtering steps were removed from further analyses.
Phylogenetic tree We built the AM fungal phylogeny from the representative sequences of each OTU. We calculated a maximum likelihood (ML) phylogeny using the GTR+I+G nucleotide substitution model [selected on the basis of corrected Akaike information criterion (AICc); ModelTest; Schliep, 2011] and 100 fast bootstrap replicates.

Root characterization
Root sequencing was performed following the pipeline described above, with exceptions specific to the group (for more details, see Matesanz et al., 2019;Illuminati et al., 2021). Specifically, a fragment of the rbcL chloroplast gene sequence (550 bp) was amplified using the primers rbcLa-F (5ʹ -ATGTCACCACAAACAGAGACTAAAGC3-3ʹ) and rbcLa-R (5ʹ-GTAAAATCAAGTCCACCRCG-3ʹ) (Levin et al., 2003;Kress et al., 2009). After sequencing on the Illumina MiSeq PE300 v.3 run, we used VSEARCH to confirm species assignments (99 % match to reference) using an in-house reference database containing the rbcL sequences of all plant species sequenced individually (see sequences in the Supplementary Data Sequences File S1). Four very closely related species pairs (Stipa pennata and Stipa tenacissima; Teucrium capitatum and Teucrium gnaphalodes; Thymus vulgaris and Thymus lacaitae; and Quercus coccifera and Q. ilex) were grouped at the genus level because their rbcL sequences were identical. Belowground species abundance was estimated as the number of reads for each species after rarefaction of each plant community (1768 sequences per sample) to equal sequence numbers. We reconstructed the rbcL phylogenetic tree from each representative sequence. We used the in-house reference database to build an ML phylogenetic tree using the GTR+I+G nucleotide substitution model and 100 fast bootstrap replicates.

Estimation of AM fungal and plant community attributes
Arbuscular mycorrhizal fungal taxonomic and phylogenetic composition and diversity Arbuscular mycorrhizal fungal taxonomic composition was assessed from the AM fungal abundance matrix. AM fungal phylogenetic composition was computed as the AM fungal phylogenetically weighted species composition matrix (Pillar and Duarte, 2010), by weighting the AM fungal abundances by their phylogenetic relationships estimated from the phylogenetic tree (Supplementary Data Methods S1; Debastiani & Pillar, 2012). AM fungal taxonomic diversity was estimated as the taxon richness in each sample (i.e. total number of OTUs). Before estimating AM fungal taxonomic diversity, AM fungal sequence reads were rarefied to the minimum number of sequences in a soil sample (1451 sequences per sample; Supplementary Data Fig. S2; Methods S1), to account for the unequal number of sequences among samples. AM fungal taxonomic diversity based on rarefied data was highly correlated with unrarefied AM fungal taxonomic diversity (r = 0.998). AM fungal phylogenetic diversity was estimated as Rao's quadratic entropy index (Rao, 1982), and one outlier was removed from subsequent statistical analysis.
Above-ground plant taxonomic composition and diversity Two measures of above-ground plant taxonomic composition were calculated as the scores of each above-ground plant sample on the first two axes of a non-metric multidimensional scaling ordination (nMDS; Supplementary Data Methods S1). nMDS is a robust method to reduce the dimensionality of the AM fungal communities in two or three variation axes . We performed the nMDS using the Bray-Curtis distance, based on the Hellinger-transformed plant cover data from each sampling circle (Barberán et al., 2015). The stress value (i.e. the parameter that shows the goodness of fit of the ordination) was 0.19. The first axis (above-ground taxonomic composition.  Fig. S3A). Above-ground plant taxonomic diversity was estimated as the plant species richness (i.e. total number of species) in each above-ground plant circle with radius 20 cm.
Above-ground plant phylogenetic composition and diversity To estimate above-ground plant phylogenetic composition, we first computed a matrix of the above-ground plant phylogenetically weighted species composition (Pillar and Duarte, 2010), using the above-ground species data matrix and the plant phylogenetic tree. Then, for each sample, two measures of plant phylogenetic composition were calculated as the scores of each sample on the first two axes of an nMDS using the Bray-Curtis distance on the Hellinger-transformed matrix. The stress value for the nMDS was 0.12. Unlike the first axis of the nMDS conducted on above-ground plant taxonomic composition, in which species belonging to the same family could be found at both ends of the taxonomic compositional gradient (e.g. Teucrium gnaphalodes and Phlomis lychnitis; Supplementary Data Fig. S3A), members of the same family are found at the same end of a phylogenetic compositional axis. The negative values of the first phylogenetic composition axis (above-ground phylogenetic composition.1) were related to grasses (Poaceae), whereas all other families (mainly forbs) were grouped towards the positive values (Supplementary Data Fig. S3B). The negative values of the second phylogenetic composition axis (above-ground phylogenetic composition.2) were related to Lamiaceae species, whereas the positive values were related to Cistaceae and Fabaceae species (Supplementary Data  Fig. S3B). The above-ground phylogenetic diversity was estimated using the Rao index.
Below-ground plant taxonomic and phylogenetic composition and diversity Taxonomic and phylogenetic composition and diversity were computed as for the above-ground community but using the number of sequences per root sample instead of plant cover.
The two below-ground plant taxonomic composition axes were negatively related to assemblages dominated by roots of Ononis tridentata, Staehelina dubia and Lithodora fruticosa (below-ground taxonomic composition.1) and of Phlomis lychnitis, Hippocrepis commutata and Lithodora fruticosa (below-ground taxonomic composition.2). Furthermore, they were positively related to assemblages dominated by roots of Linum narbonense, Euphorbia nicaeensis and Fumana thymifolia (below-ground taxonomic composition.1) and of Euphorbia nicaeensis, Linum narbonense and Salvia lavandulifolia (below-ground taxonomic composition.2) (Supplementary Data Fig. S3C). The first below-ground phylogenetic composition axis was negatively related to Cistaceae and Lamiaceae species and positively related to Poaceae species. The second below-ground phylogenetic composition axis was negatively related to Poaceae species and positively related to Fabaceae species (Supplementary Data Fig. S3D). The nMDS for the taxonomic and phylogenetic below-ground composition had stress values of 0.21 and 0.07, respectively. Below-ground plant taxonomic and phylogenetic diversities were estimated as plant species richness (i.e. total number of taxa) and the Rao index, respectively, in each sample.
Above-ground plant cover and root biomass Above-ground plant cover was estimated as the sum of the cover of all species within the circle with a radius of 20 cm, and root biomass was assessed as the total root biomass per root sample.
Soil and spatial covariates To reduce the number of soil physicochemical variables, we performed a principal components analysis (PCA) with varimax rotation. Before performing the PCA, we estimated the best number of components to retain in our analyses using the smoothing approximation of the cross-validation criterion implemented by Josse and Husson (2012). We retained the first four principal components. The first axis, explaining 22 % of the variance, was highly related to soil texture, varying from soil with high sand content at one end of the axis to soil with high clay and slit content at the other. The second axis (explaining 19 %) was mainly related to variations in soil organic carbon. The third axis (explaining 18 %) was positively correlated with nitrogen and phosphorus content (fertility). Finally, the fourth axis (explaining 16 %) was associated with variations in pH and conductivity (salinity) (for more details, see table S2 in the paper by López-Angulo et al., 2020).
To account for unmeasured spatially structured variables (i.e. spatial covariates), we generated a set of Moran's eigenvectors from the coordinates of each sampling point using distance-based Moran's eigenvector maps (dbMEM; . Initially, we fitted a trend surface (i.e. a linear regression on the x-and y-coordinates) to remove spatial trends in the four response variables (i.e. AM fungal taxonomic and phylogenetic diversity, taxonomic and phylogenetic composition). The residuals of these regressions were used to compute dbMEM eigenvectors (Borcard and Legendre, 2002). To select the dbMEM eigenvectors that significantly contributed to explain AM fungal diversity and composition, we performed a forward selection procedure with double-stopping criterion (α = 0.05, 9999 permutations; Blanchet et al., 2008) of the dbMEM eigenvectors. We repeated this procedure independently for each dependent variable (i.e. the AM fungal taxonomic and phylogenetic diversity and composition), such that a different number of dbMEM variables was selected in each case (AM fungal taxonomic diversity = 5; phylogenetic diversity = 3; taxonomic composition = 2; phylogenetic composition = 3). The x-and y-coordinates were also included as predictors in the final models.

Statistical analyses
To determine whether above-ground plant community attributes were correlated significantly with the plant community attributes estimated below-ground, Pearson correlations were conducted. We included all plant community attributes as predictors because correlation coefficients varied between 0.5 and −0.5 (Supplementary Data Figs S4 and S5). All these statistical analyses, and the subsequent ones, were performed in R (v.4.0.2; R Core Team), and a detailed description of packages used can be found in Supplementary Data Methods S1.

Arbuscular mycorrhizal fungal composition
We used partial redundancy analysis (pRDA;  to study the above-and below-ground effects of the plant community on the variation of AM fungal composition. As the response matrix, we used either the taxonomic or the phylogenetically weighted AM fungal composition matrices. As predictors, we used the above-and below-ground plant diversity and composition, root biomass and above-ground plant cover. The effect of the below-ground plant diversity and composition was tested after accounting for the effects of the above-ground plant diversity and composition, and vice versa, and root biomass and above-ground plant cover. To account for soil effects and residual spatial variation, we constrained the ordinations by soil physicochemical properties (PCAs) and spatial covariates (dbMEM eigenvectors). To avoid multicollinearity, we performed independent pRDAs for the taxonomic and phylogenetic plant community predictors. We conducted the pRDAs with forward selection using the double-stopping criteria (P < 0.05 and adjusted R 2 < global R 2 ; Blanchet et al., 2008). We tested the significance of predictors using the Monte Carlo test based on 999 permutations. Before statistical analyses, we applied a Hellinger transformation on the AM fungal and plant community matrices (Legendre and Gallagher, 2001). We applied a square root transformation on root biomass and the salinity gradient (second PCA axis of the soil PCA).

Arbuscular mycorrhizal fungal diversity
We fitted generalized linear models to both AM fungal taxonomic and phylogenetic diversity to assess the effects of the plant community attributes (above-and below-ground taxonomic and phylogenetic composition and diversity, root biomass and above-ground plant cover), controlling the effects of soil physicochemical properties (PCAs) and spatial covariates (dbMEM eigenvectors). The response of the AM fungal taxonomic diversity to predictors was evaluated considering a Poisson error distribution and logarithmic link function, whereas the response of AM fungal phylogenetic diversity was analysed using a Gaussian error distribution and identity link function. We used a model selection procedure based on minimizing the AICc to select the best predictors of AM fungal taxonomic and phylogenetic diversities. The models were ranked according to the AICc. We calculated model-averaged parameter estimates over the set of models with ∆AICc < 2 (Burnham and Anderson, 2002). We estimated 95 % confidence intervals (CIs) around model-averaged parameter estimates, considering a parameter to be significant if the 95 % CI excluded zero (Burnham and Anderson, 2002). We checked the absence of multicollinearity in all models using the variance inflation factor. In all cases, variance inflation factors were smaller than four, indicating absence of collinearity (Zuur et al., 2010). To avoid overfitting, we allowed a maximum of seven predictors in the candidate models because it is recommended that the sample size (i.e. number of soil cores) should be ten times greater than the maximum number of predictors (Harrell et al., 1984;Hair et al., 2014; Supplementary Data Methods S1).

Taxonomic description of the AM fungal and plant communities
A total of 359 507 18S rRNA gene sequence reads were assigned to 1267 AM fungal OTUs, ranging from 1451 to 13 931 reads per sample. AM fungal communities were dominated by taxa within the Glomeraceae family, representing 91 % of the sequences (Supplementary Data Fig. S6). The remaining OTUs belonged to the Acaulosporaceae, Ambisporaceae, Archaeosporaceae, Claroideoglomeraceae, Diversisporaceae, Gigasporaceae and Paraglomeraceae (Supplementary Data  Table S3). AM fungal taxonomic diversity ranged from 8 to 130 OTUs per sample (62.9 ± 22.4 OTUs; mean ± SD).
A total of 30 plant taxa, 26 identified at the species level and 4 at the genus level, were found below-ground across all soil samples (identified through DNA metabarcoding), whereas 39 plant species were detected above-ground in the area corresponding to all circles of 20 cm radius. The below-and above-ground taxonomic diversity ranged from 3 to 12 plant species (7 ± 2; mean ± SD) and from 3 to 16 (8 ± 2.7; mean ± SD), respectively. The above-ground and below-ground plant composition were significantly correlated, both taxonomically (Supplementary Data Fig. S4) and phylogenetically (Supplementary Data Fig.  S5). The above-ground plant taxonomic diversity was weakly and significantly correlated with below-ground plant taxonomic diversity (r = 0.36, P < 0.01). In contrast, the above-ground plant phylogenetic diversity was not correlated with below-ground plant phylogenetic diversity (Supplementary Data Fig. S5).

Arbuscular mycorrhizal fungal community composition
Partial pRDAs showed that the forward-selected aboveand below-ground attributes of plant communities together explained 1.2 and 5.6 % of the total variance of taxonomic and phylogenetic AM fungal composition, respectively (for forward-selected predictors, see Table 1; Supplementary Data Fig. S7), and the soil physicochemical properties (soil fertility) explained 0.6 % of the taxonomic AM fungal composition (Table 1; Supplementary Data Fig. S7). We found that belowground plant taxonomic composition and above-ground plant phylogenetic composition explained unique fractions of variation in the AM fungal taxonomic and phylogenetic composition (Table 1; Supplementary Data Fig. S8). Specifically, Glomeraceae taxa were associated with forbs, whereas the rest of the families tended to be associated with grasses ( Supplementary Data Fig. S9). Above-ground taxonomic plant diversity, but not below-ground plant diversity, was also associated with variations of AM fungal taxonomic composition (Table 1). This indicates that more species-rich plant assemblages above-ground led to significant changes in the identity of the taxa co-occurring in the AM fungal communities.
The AM fungal composition was also associated with variations in soil physicochemical properties (Table 1). Specifically, pRDA biplots showed that most AM fungal taxa belonging to Diversisporaceae and Claroideoglomeraceae families were associated with poorer soils, and the taxa belonging to Glomeraceae were associated with alkaline soils (pH ≈ 8; Supplementary Data Fig. S10). pRDA ordinations also showed that the AM fungal communities were spatially structured (Supplementary Data Table S4).

Arbuscular mycorrhizal fungal diversity
Arbuscular mycorrhizal fungal taxonomic diversity was positively associated with above-ground plant phylogenetic diversity (Fig. 2), and overall, this was consistent for all fungal families (Supplementary Data Fig. S11). AM fungal taxonomic diversity was also associated with variations in above-and below-ground plant taxonomic composition (Fig. 2), indicating that different combinations of plant species can lead to richer (higher number of OTUs) AM fungal communities. Furthermore, the two nMDS axes of above-ground phylogenetic plant composition affected the AM fungal taxonomic diversity, and, in the opposite direction, the AM fungal phylogenetic diversity (Fig. 2). This result indicated that areas with above-ground cover of Lamiaceae species were related to more taxa-rich but less phylogenetically diverse AM fungal communities, matching the increase in Glomeraceae taxa (Supplementary Data Fig. S11). The abundance of Poaceae, Fabaceae and Cistaceae species led to poor (in terms of the number of OTUs) but more phylogenetically diverse AM fungal communities (Fig. 2). This relationship was probably attributable to the finding that these plant groups decreased Glomeraceae taxa and increased Claroideoglomeraceae and Paraglomeraceae ( Supplementary Data Fig. S11).
The AM fungal taxonomic diversity decreased with soil fertility (1.1 % of explained variance), whereas AM fungal phylogenetic diversity decreased with soil alkaline pH (7.2 % of explained variance; Supplementary Data Fig. S12). The spatial eigenvectors representing unmeasured spatially structured factors exerted positive effects on AM fungal taxonomic diversity and negative effects on AM fungal phylogenetic diversity (Supplementary Data Fig. S12). Table 1

. Significance level (P-values) and adjusted coefficient of determination (adjusted R 2 ) of the effects of above-and below-ground plant attributes (taxonomic and phylogenetic diversity and composition) and soil physicochemical properties (pH and fertility) on the arbuscular mycorrhizal fungal taxonomic and phylogenetic composition and diversity.
Arbuscular mycorrhizal fungal composition Arbuscular mycorrhizal fungal diversity Taxonomic  Phylogenetic  Taxonomic  Phylogenetic Plant predictor set Adjusted R 2 P-value Adjusted R 2 P-value Adjusted R 2 P-value Adjusted R 2 P-value The significance of predictors on the AM fungal composition was tested against 999 Monte Carlo permutations in partial redundancy analyses. The significance of predictors on the AM fungal diversity was tested using a z-statistic after estimating 95 % confidence intervals around model-averaged parameter estimates using generalized linear models. The reported adjusted R 2 represents the unique variance explained by each variable after variance partitioning (see shared variance in Supplementary Data Figs S7 and S8). Only variables significant in at least one model are shown in the table. Composition.1 and composition.2 are variables calculated as the scores of each sample on the two first axes of non-metric multidimensional scaling ordinations (nMDS).

DISCUSSION
Our findings provide evidence of a heterogeneous distribution of AM fungal communities at the spatial scale in which plants interact in a Mediterranean scrubland. We found that the plant community attributes, soil physicochemical properties and spatial variables jointly explained a significant fraction of the local distribution of AM fungal communities. Thus, we demonstrated that using sets of different predictors related to aboveand below-ground vegetation, soil and space provided a more complete picture of the complex nature of AM associations. Several studies have also shown a high spatial variation of the AM fungal communities at fine scales (Horn et al., 2017;Avio et al., 2020). However, unlike these studies, we found that plant community attributes predicted unique fractions of AM fungal community structure and composition. Importantly, these results point out that above-and below-ground plant community attributes that predicted unique fractions of AM fungal community structure were decoupled, as previously shown in the same plant community (Illuminati et al., 2021) and in other environments (Schenk and Jackson, 2002;Reintal et al., 2010;Hiiesalu et al., 2012). This above-ground plant influence, together with a positive effect of root biomass and negative effect of soil fertility on AM fungal taxonomic diversity, has also been documented recently for the taxonomic diversity of key soil fungi in the same study system (López-Angulo et al., 2020), suggesting that fungi belonging to subphylum Glomeromycotina respond in a similar manner to general ecological drivers of soil fungal diversity (i.e. nutrient gradient and habitat availability) (Legay et al., 2014). Nevertheless, AM fungal diversity and composition were also more directly associated with the below-ground compartment of the plant community (particularly with the below-ground plant composition) than soil fungal diversity and composition (López-Angulo et al., 2020), confirming the expected direct relationship between AM fungi and plant roots.
Many studies to date have used above-ground information to establish relationships between plant and AM fungal communities, obviating the root fraction and thus, assuming a certain symmetry between both community compartments (Horn et al., 2017;Van Geel et al., 2018;Bittebiere et al., 2020). The recent development of high-throughput DNA metabarcoding has offered a means to determine root diversity (Matesanz et al., 2019;Cabal et al., 2021), providing opportunities for a better and more realistic understanding of the interactions between plants and AM fungi. Despite these methodological advances to identify plant species below-ground, our results show that above-ground information such as plant cover, which is easy to estimate visually, might also be an indicator of the linkages between plants and AM fungi. The rationale behind the strong association between above-ground plant diversity and these obligate root symbionts might be that the larger sampling size above-ground provides an extensive characterization of the entire plant community, including that below-ground. On the contrary, soil core sampling might fail to capture all species recorded in the above-ground sampling rings. Although the radius for the above-ground sampling circles maximizes the similarity between both plant compartments, significant discrepancies remain between them in this community (Illuminati et al., 2021). For example, plant shoot and root systems can exhibit large differences in lateral spread and temporal turnover dynamics (Schenk and Jackson, 2002;Sun et al., 2016). In this sense, our above-ground sampling might have provided a longer-lasting picture of the plant community, because above-ground effects might be more stable temporally than those below-ground, given that root distributions change faster over time than the above-ground parts (Peek et al., 2005;Sun et al., 2016). These different temporal dynamics could lead to AM fungal communities being coupled differently with both compartments of plant communities. The AM fungal diversity and composition in our Mediterranean scrubland were explained mainly by variations in the above-ground plant composition, showing that certain plant species can stimulate the performance of particular AM fungal taxa. These findings, which denote some degree of preferential association between plants and AM fungi, concur with results of microcosms and observational studies conducted in other systems (Johnson et al., 2004;Hausmann and Hawkes, 2009;Alguacil et al., 2019;Šmilauer et al., 2020). Furthermore, the fact that plant composition, but not plant diversity, affected AM fungal taxonomic diversity suggests that the identity of plant species influences the assembly of AM fungal communities more directly than the variety or the number of plant species (Trinder et al., 2009;Neuenkamp et al., 2021). Another explanation for the lack of relationship between below-ground plant and AM fungal taxonomic diversity would be the occurrence of common mycorrhizal networks (Van Der Heijden and Horton, 2009;van der Heijden et al., 2015). In these networks, few AM fungal taxa colonize different plant roots, interconnecting different plant species and thus leading to a decoupling between AM fungal and plant diversity, especially at this fine scale.
When the phylogenetic relationships were considered in the analyses, above-ground plant phylogenetic composition (first axis of nMDS), which separated grasses (Poaceae) from forbs (Supplementary Data Fig. S3B), was the only predictor explaining a significant and unique fraction of variation in the AM fungal phylogenetic composition and diversity (Supplementary Data Fig. S8). This result is in line with other studies showing the existence of a phylogenetic signal for AM associations (Yang et al., 2012;Montesinos-Navarro et al., 2015;Chen et al., 2017;López-García et al., 2017). Our findings indicated that AM fungal taxa belonging to family Glomeraceae tended to be associated with forbs, whereas other AM fungal families (especially Paraglomeraceae, Ambisporaceae, Claroideoglomeraceae and Archaeosporaceae) were more related to grasses. This was probably the reason why AM fungal communities related to grasses were more phylogenetically diverse although they exhibited a lower number of OTUs. This agrees with Šmilauer et al. (2020), who found that plant species (mainly Poaceae) negatively affecting the richness of the AM fungal communities also positively affected the AM fungal phylogenetic diversity (mean nearest taxon phylogenetic distance).
In addition, our linear models showed that Lamiaceae species were positively associated with a larger number of OTUs dominated by Glomeraceae taxa, whereas Fabaceae and Cistaceae species supported a larger phylogenetic diversity of AM fungal communities. These results could be explained by the relative fitness differences and the niche differences among competing AM fungal taxa, as postulated by the coexistence theory (Chesson, 2000). Plant assemblages dominated by Lamiaceae, the most diverse and abundant plant family in our scrubland, led to more phylogenetically clustered AM fungal communities but with a larger number of OTUs, dominated by Glomeraceae taxa. This pattern might be the result of asymmetric competition (Mayfield and Levine, 2010) through the exclusion of poor competitors belonging to distantly related AM fungal lineages produced by Glomeraceae taxa. Members of the Glomeraceae establish a faster (Hart and Reader, 2002), more extensive (Maherali and Klironomos, 2007) and more permanent (Hart and Reader, 2002;López-García et al., 2014) mycelium in roots compared with other AM fungal families, and they might also exhibit higher competitive ability for carbon uptake (Yang et al., 2014). In contrast, in soils dominated by grasses and legumes, niche differentiation might favour AM fungal communities composed of more distantly related taxa, in turn enhancing coexistence among AM fungi. This is supported by the fact that grasses are expected to be less dependent on mycorrhizal symbiosis owing to the high efficiency of their fibrous roots to absorb nutrients (Javaid, 2008). Accordingly, legumes might be associated with phylogenetically overdispersed AM fungal communities because the functional traits needed to associate with N-fixing plants are not phylogenetically conserved (López-García et al., 2017). These opposite patterns of phylogenetic and taxonomic diversity suggesting competitive exclusion were also found in a study of an alpine meadow, in which an increase in soil fertility (N + P) caused a loss of taxonomic diversity (particularly of Glomus species), but an increase in genus taxonomic diversity (Liu et al., 2015). Our findings highlight that changes in environmental conditions, including those mediated by specific groups of plants, can shape the outcome of the direct biotic interactions among AM fungi (Thonar et al., 2014;Knegt et al., 2016).
Finally, in our study, the variance of AM fungal phylogenetic diversity and composition explained by the spatial eigenvectors was higher than the variance explained by the plant attributes and soil physicochemical properties. This suggests that the AM fungal assemblage exhibited strong phylogenetic and spatial structuring. Other spatially autocorrelated abiotic factors, together with biotic interactions among soil microorganisms, could be mechanisms that drive the spatial patterns in the phylogenetic structure of AM fungal communities (Maherali and Klironomos, 2007;Knegt et al., 2016). Furthermore, the large unexplained variation of AM fungal diversity and composition observed in our models probably reflects unmeasured deterministic or neutral processes that are not spatially structured at the spatial scales considered in our sampling. Such unexplained variation might include the effect of topographic or microclimatic environmental variations, such as soil moisture and temperature heterogeneity (Kivlin et al., 2011;Teste et al., 2020). Therefore, the inclusion of other abiotic variables that might influence AM fungi and plants in parallel could increase the shared variation between plants and the soil environment.

Conclusions
We provide evidence that the plant community attributes, soil physicochemical properties and spatial variables jointly affect the fine-scale distribution of AM fungal communities in a Mediterranean scrubland. Importantly, our results highlight that, despite the decoupling found between the above-ground and below-ground plant compartments, the more easily accessible above-ground vegetation might serve as an indicator of the active interactions of AM fungi with plant roots. However, agreeing with our hypothesis, both below-and above-ground plant community attributes also explained unique fractions of the variation in AM fungal composition and diversity at the fine scale in which plants interact. Thus, we emphasize the importance of considering the above-ground plant distribution and soil environmental factors, together with data from the below-ground plant compartment, because they clearly improve our ability to predict the relationships between AM fungal and plant communities. As hypothesized, we show the importance of incorporating the phylogenetic relatedness of mutualistic partners into metrics of diversity and composition to improve our understanding of the ecological and evolutionary mechanisms underlying the assembly of AM fungal communities.

SUPPLEMENTARY DATA
Supplementary data are available online at https://academic. oup.com/aob and consist of the following.
Sequences File S1: species in-house reference database containing the rbcL sequences of all plant species individually sequenced. Methods S1: R packages used for specific applications. Figure S1: core sampling design. Figure S2: rarefaction curve of arbuscular mycorrhizal fungal communities for each soil sample. Figure S3: non-metric multidimensional scaling ordinations showing patterns of variation in plant above-ground and below-ground taxonomic and phylogenetic composition. Figure S4: pairwise scatterplots of the plant community taxonomic attributes. Figure S5: pairwise scatterplots of the plant community phylogenetic attributes. Figure S6: proportion of sequence reads of arbuscular mycorrhizal fungi at the family level. Figure S7: Venn diagrams showing variance partitioning results of AM fungal taxonomic and phylogenetic composition, and taxonomic and phylogenetic diversity explained by plant community attributes, soil properties and spatial covariates. Figure S8: Venn diagrams showing variance partitioning results of AM fungal taxonomic and phylogenetic composition, and taxonomic and phylogenetic diversity explained by the below-ground and above-ground taxonomic and phylogenetic plant community attributes. Figure S9: redundancy analysis biplots showing the relationship of AM fungal taxonomic composition and phylogenetic composition to the below-ground and above-ground taxonomic and phylogenetic plant community attributes. Figure S10: redundancy analysis biplots showing the relationship of AM fungal taxonomic composition and phylogenetic composition to the soil properties and spatial covariates. Figure S11: linear relationships between the AM fungal taxonomic and the taxonomic and phylogenetic plant community attributes. Figure S12: effect of soil physicochemical properties and the spatial covariates. Table S1: list of plant species found in our study plot. Table S2: summary statistics for soil physicochemical properties. Table S3: list of arbuscular mycorrhizal fungi virtual taxa indicating taxonomic ranks, included known species and unidentified cultures. Table S4: ANOVA-like results based on redundancy analysis testing the effect of the forward-selected soil physicochemical properties and spatial covariates on the arbuscular mycorrhizal fungi taxonomic and phylogenetic composition. Table S5: results of the AICc-based model selection based on linear models testing the response of arbuscular mycorrhizal fungi diversity to plant community attributes controlling for the effects of above-ground plant cover, root biomass and soil physicochemical properties.