Mycorrhizal symbiosis increases plant phylogenetic diversity and regulate community assembly

The intricate mechanisms shaping plant diversity and community composition are the cornerstone of ecological understanding.  Yet, the role of mycorrhizal symbiosis, the fundamental partnership between fungi and plant roots, in influencing community composition has often been underestimated.  Here, we use extensive species survey data from 1,315 terrestrial ecosystem sites to elucidate the influence of mycorrhizal symbiosis on plant phylogenetic diversity and its implications for community assembly processes.  Our findings demonstrate that increasing mycorrhizal symbiotic potential leads to greater phylogenetic dispersion within plant communities. Furthermore, we unveil a distinct dichotomy in the assembly processes governed by mycorrhizal status. Mycorrhizal species predominantly influence deterministic processes, suggesting a role in niche-based community assembly.  Conversely, non-mycorrhizal species exert a stronger influence on stochastic processes, highlighting the importance of random events in shaping community structure.  These results underscore the crucial but often hidden role of mycorrhizal symbiosis in driving plant community diversity and assembly. This study provides valuable insights into the complex mechanisms shaping ecological communities and the way for more informed conservation and management practices that acknowledge the complex interplay between symbiosis and ecological community dynamics.

Methods
Filed survey and data collection

This study utilized field survey data collected from 1315 sites across various grassland ecosystems in China.  Vegetation surveys were conducted during the peak plant growth season, specifically from mid-July to August in middle and high latitude regions, including the Qinghai-Tibetan Plateau, and from August to September in desert, subtropical, and tropical regions.  At each site, plant community data was collected using ten 1m × 1m quadrats (reduced to 0.5m × 0.5m for meadow and alpine meadow) located randomly within a 100m × 100m area.   For shrubland ecosystems, five 5m × 5m quadrats were randomly placed within the same 100m × 100m area.  Five of the ten quadrats (or all five quadrats in shrublands) were randomly selected for detailed vegetation analysis. Within these selected quadrats, plant species richness (SR) and species relative abundance were recorded.  Soil samples were collected from each quadrat using a soil core method at a depth of 20 cm.  Soil organic carbon content (SOC) were measured for each sample.  Total phosphorus (TP), total nitrogen (TN), and soil pH were interpolated from the Basic soil property dataset of high-resolution China Soil Information Grids (2010-2018).

To facilitate analysis and interpretation, the 1315 field sites were classified into four distinct ecosystem types: meadow (further categorized into lowland meadows (LM), mountain meadows (MM), and alpine meadows (AM)), steppe (divided into temperate steppe (TS) and alpine steppe (AS)), shrubland (classified as warm shrubland (WG) and tropical shrubland (TG)), and desert (DS).  This classification was based on plant community composition, climate, and prevailing environmental conditions.

Mean annual temperature (MAT, °C), mean annual precipitation (MAP, mm), and mean diurnal range (mean of monthly maximum temperature - minimum temperature, °C) for each site were obtained from the WorldClim data layers (specifically, bio_1 and bio_12) at a spatial resolution of 30 seconds × 30 seconds (approximately 1 km × 1 km at the equator) (http://www.worldclim.org/).

Remote sensing data, including Fraction of Photosynthetically Active Radiation (Fpar) were obtained from MOD15A2H Version 6 data product, slope and elevation data was extracted from the STRM 90m dataset 171 in China, based on the SRTM V4.1 database (https://www.resdc.cn/data.aspx?DATAID=123).  Fpar is defined as the fraction of incident photosynthetically active radiation, 400-700 nanometers (nm), absorbed by the green elements of a vegetation canopy.

Plant mycorrhizal status and community mycorrhizal index

To quantify the mycorrhizal status of each plant community, we calculated a mycorrhizal index representing the degree of potential mycorrhizal colonization within the community.  The mycorrhizal status of each plant species was determined using an established database of mycorrhizal associations.  To minimize potential errors during the matching process, species were matched based on their genus level (Brundrett & Tedersoo 2019).

The dominant mycorrhizal status of each community was determined based on the mycorrhizal status exhibiting the highest abundance within that community.  This approach allowed us to differentiate communities based on the predominant mycorrhizal association of their constituent plant species.

Construction of phylogenetic relationships and calculation of phylogenetic distances

Phylogenetic relationships and distances between plant species were determined using the V. PhyloMaker package and picante package in R. V. PhyloMaker generated phylogenetic hypotheses for the 1235 plant species in our study by linking them to the 'GBOTB.extended' megatree.  This megatree encompasses 74,531 species, representing all families of extant vascular plants, and serves as the largest dated phylogeny for vascular plants.  Phylogenetic distances within each community were calculated using the picante package.  We calculated two metrics: mean pairwise distance (MPD) and mean nearest taxon distance (MNTD).  To account for the relative abundance of each species within the community, species abundance was incorporated as a weighting factor in the phylogenetic distance calculations.

To mitigate the influence of species richness on community phylogenetic distances, standardized effect size metrics for MPD (SESMPD) and MNTD (SESMNTD) were calculated.  This standardization involved generating a null distribution by randomly shuffling the distance matrix labels across all taxa 999 times.  The mean of the null distribution was then used to calculate the standardized effect size.  SES values less than 0 indicate phylogenetic clustering, where species within the community are more closely related than expected by chance.  Conversely, SES values greater than 0 indicate phylogenetic overdispersion, where species are more distantly related than expected.

Statistical analyses

Comparison of plant community composition.  To compare plant community composition across different dominant mycorrhizal status, we first examined the relative abundance of different plant families within each of the community types.  We then performed non-metric multidimensional scaling (NMDS) using the metaMDS function in the vegan package, based on a Bray-Curtis distance matrix.  The stress function value was assessed to ensure model reliability.  Significant differences between groups were determined using Adonis, with significance levels denoted as follows: n.s. (not significant), *p < 0.05, **p < 0.01, and ***p < 0.001.

Assessing the influence of mycorrhizae on species richness and phylogenetic diversity.  To evaluate the role of mycorrhizae in shaping community species richness and phylogenetic diversity, we constructed a set of predictor variables encompassing climate (MAT, MAP, MDR, Fpar), geographic factors (longitude, latitude, elevation, slope), mycorrhizal status (MI), and soil properties (TN, TP, SOC, pH).

Multiple regression models were developed using the MuMIn package in R to assess the effects of these predictors on species richness and phylogenetic diversity.  A full set of models incorporating all predictors was generated and ranked based on the Akaike information criterion (AIC) using maximum likelihood estimation.  Models with ΔAIC < 2 were retained, and model averaging was employed to calculate parameter estimates and p-values.  The relative effect of each predictor was determined by calculating the ratio of its parameter estimate to the sum of all parameter estimates.  To further emphasize the importance of mycorrhizae in predicting species richness and phylogenetic diversity, we constructed comparative models excluding the mycorrhizal index while retaining all other predictors.  These reduced models were compared to their corresponding full models using AICc values, with lower AICc values indicating superior model performance.

Quantifying the relative importance of stochastic and deterministic processes in community assembly.  To assess the influence of mycorrhizae on community assembly processes, we utilized normalized stochasticity ratio (NST) analysis.  NST, an extension of the Beta-diversity metric, quantifies the relative contribution of stochastic and deterministic processes in community assembly.  NST values range from 0 to 1, with 0.5 representing an equal contribution of both processes.  NST values predominantly above 0.5 indicate dominance of stochastic processes, while values below 0.5 suggest a greater influence of deterministic processes.

We further investigated the impact of different mycorrhizal plant types on community assembly using phylogenetic bin-based null model analysis (iCAMP) (Ning et al. 2020).  This method, implemented using the iCAMP package in R, divides plant species into phylogenetic bins with significant phylogenetic signals and quantifies the contribution of each bin to deterministic (homogeneous selection (HoS), heterogeneous selection (HeS)) and stochastic (dispersal limitation (DL), homogenizing dispersal (HD), and drift (DR)) processes.  We selected the 300 most abundant plant species (representing 91.6% of total individuals), assigned them to phylogenetic bins, and determined the dominant mycorrhizal status within each bin.  The relative contribution of each bin to different assembly processes was then assessed to quantify the influence of different mycorrhizal plant types on ecological processes.