Impacts of Urbanization Undermine Nestedness of the Plant–Arbuscular Mycorrhizal Fungal Network

Cities are prone to ecological problems, yet the impacts of rapid global urbanization on the feedback between above- and belowground subsystems remain largely unknown. We sampled the roots of 8 common herbaceous plants within the Fifth Ring (urban areas) and in Jiufeng National Forest Park (rural areas) in Beijing (China) to assess the impacts of urbanization on the network of plant-arbuscular mycorrhizal (AM) fungal associations. Using Illumina MiSeq sequencing, 81 AM fungal OTUs were identified in 78 herb root samples. The Shannon, Simpson, and Pielou indices of root AM fungi in urban areas were significantly higher than those in rural areas. In this study, a significantly nested mycorrhizal association network was observed in rural areas (NODF = 64.68), whereas a non-nested pattern was observed in urban areas (NODF = 55.50). The competition index C-score (0.0769) of AM fungi in urban areas was slightly lower than that in rural areas (0.1431), and the species specialization (d’) of 8 host plants and fungal dissimilarity among 8 host plants in urban areas were significantly lower than those in rural areas. Convergent associations among hosts may be an important factor influencing this non-nested pattern of the plant-AM fungi network in urban areas. Generalists, rather than specialists, were enhanced during the establishment of mycorrhizal associations in urban areas. Our results suggest that reduced selectivity of host plants, and generalist promotion and specialist reduction of AM fungi during urbanization may contribute to the non-nested network of plant-AM fungal associations.


INTRODUCTION
Urban expansion is occurring worldwide (Angel et al., 2005) and may affect aboveand belowground system feedback. Interactions between plants and mycorrhizal fungi are tremendously complex, and the multispecies networks resulting from these associations have consequences for plant growth and productivity (Garg and Chandel, 2010;Hacquard and Schadt, 2015). Habitat fragmentation, isolation, heavy metal pressure, high nutrient availability, and plant diversity changes that occur with rapid urbanization could present great ecological risks to plant-mycorrhizal fungal associations (Pickett et al., 2001;Jacquot-Plumey et al., 2003;Hoch et al., 2019). Currently, urban mycorrhizal studies are mainly focused on understanding the essential role of AM fungi in plant growth (Newbound et al., 2010;Chagnon and Brisson, 2017) and their diversity along urbanization gradients (Bainard et al., 2011;Martinova et al., 2016). Knowledge of how AM fungi are associated with co-occurring plant species could offer novel insights into the community-scale processes of plant-AM fungus symbioses (Chagnon et al., 2012). However, whether the network consists of host plants and their associated fungi are affected by urbanization remains largely unknown.
AM fungi assist strict associations with plants to compensate for the loss of saprotrophic capabilities (Bonfante and Genre, 2010). Associations between AM fungi and plants follow nonrandom assemblages (Montesinos-Navarro et al., 2012;Sepp et al., 2019). Roots colonized with different AM fungal species exhibit different alleviation effects on Arizona cypress seedlings under fuel pollutants (Aalipour et al., 2019). The interspecific interactions in communities with a nested network structure in which specialist species interact with subsets of species that interact with generalist species have been found to be highly robust to species extinctions (Bascompte et al., 2006;Bascompte, 2010). Significantly nested network topologies are observed in plant-AM fungal associations (Chagnon et al., 2012;Montesinos-Navarro et al., 2012). Moreover, partner selection in mycorrhizal associations could be mediated by competition, resource abundance, and microhabitat conditions (Verbruggen et al., 2012;van der Heijden et al., 2015;Lekberg and Waller, 2016) and is thus context-dependent. Experiments conducted in urban green roofs confirm that nutrient conditions can affect AM fungal association establishment (Verbruggen et al., 2012;Chaudhary et al., 2019). Some experiments have revealed no significant colonization difference between urban parks and mature forests (Karliński et al., 2014). Under similar resources, optimal foraging theory predicts selection for generalization in choosing partners (Waser et al., 1996), -the latter conflicting with the plants' interest in specialized partners. The selective environmental filters in urbanized areas differ from those in rural areas (Knapp et al., 2008). To adapt to the urban environment, urban plants present a homogenization of the spatial plant composition (Ramalho et al., 2018). However, for mycorrhizal associations, whether urbanization raises positive or negative impacts on the nestedness of the association network remains poorly understood.
To address this research gap, we studied 8 widely distributed herbs in both urban and rural areas of Beijing (China) to investigate the impact of urbanization on the pattern of plant-AM fungal networks. In consideration of the findings of previous studies (Ramalho et al., 2018), such as the nested pattern of plant-AM fungal associations and the homogenization of the spatial plant composition in urban areas, we hypothesized that (1) the urban greenspace plant-AM fungal network possess lower levels of nestedness than the rural network and that (2) selection for generalists is the main factor accounting for the change in nestedness during urbanization. The results deepen our understanding of above-and belowground system feedback in the context of rapid global urbanization.

Study Site and Sampling
This study was conducted in Beijing (39 • 28 -41 • 05 N, 115 • 25 -117 • 30 E), which is located in the northern part of the North China Plain. Beijing has a typical continental monsoon climate, with a mean annual precipitation of 571.8 mm and a mean annual temperature of 10-12 • C. This area has cinnamon soil with a loamy texture. Urban area is located within the Fifth Ring of Beijing. Rural area is located at Jiufeng National Forest Park (40 • 03 54 N, 116 • 05 45 E), with a forest coverage of 96.2% and approximately 30 km from the urban area of Beijing (Figure 1).
Fieldwork was performed from July to August 2013. Based on field investigation, 10 common local herbaceous species, 10 genera, 7 families (He et al., 1992), easily colonized by AM fungi were selected in urban areas. Yet, Trigonotis peduncularis and Taraxacum mongolicum failed to be sampled in the Jiufeng National Forest Park. Thus, 8 herbs Crepidiastrum sonchifolium, Potentilla supina, Cirsium setosum, Setaria viridis, Viola philippica, Solanum nigrum, Plantago depressa, and Inula japonica in rural and urban areas were used in network construction. For each herb, 5 spatially separated 20 m * 20 m plots were established in urban areas and in rural areas (Figure 1). Within each 20 m * 20 m plot, 3 1 m * 1 m subplots were established, and 3 individuals of each herbaceous species and 5 soil cores (3.5 cm diameter * 20 cm deep) were collected in each subplot. Thus, fifteen soil cores were collected from each 20 m × 20 m plot and pooled to yield a composite soil sample. The roots of the sampled plants were rinsed with tap water and then with sterilized distilled water. To reduce the workload, for each herb species, the roots of 9 plant individuals collected from a 20 m × 20 m plot were pooled to yield a composite root sample. In total, 80 root samples were collected in rural and in urban areas. All root samples were frozen and stored at −80 • C for subsequent DNA extraction.

DNA Extraction and Sequencing
The root samples were pulverized and mixed with liquid nitrogen using a mortar that was cleaned with HNO 3 and sterilized distilled water and sterilized in an autoclavation. Total DNA was extracted from 0.05 g of each sample using a MOBIO PowerPlant DNA extraction kit (MO Bio Laboratories, United States). The final DNA quality and contents was measured by NanoDrop (Thermo Fisher Scientific, United States). The total DNA was diluted 10 times, 20 times, and 50 times as template for the primary PCR amplifications.
The SSU regions of the nuclear ribosomal RNA genes of Glomeromycota were amplified from DNA extracts using a semi-nested PCR protocol with primer pairs AML1/AML2 (Lee et al., 2008) and NS31/AML2 that combined with adapter sequences and barcode sequences (Simon et al., 1992). Primary PCR carried out in 25 µL reaction volumes, each containing 1 µL of DNA template (10 times, 20 times, or 50 times DNA dilutions), 2.5 µL of buffer, 2 U of Taq DNA polymerase, 4 mM dNTPs, and each primer at 4 µÌ. The following cycling parameters were used: 94 • C for 5 min; 35 cycles at 94 • C for 30 s, 58 • C for 45 s, and 72 • C for 1 min; and a final step at 72 • C for 10 min. PCR products of 3 replications were mixed and diluted 10 times with double-distilled H 2 O as the DNA template for the second PCR amplifications. The second PCR carried out in 50 µL reaction volumes, each containing 5 µL of DNA template, 5 µL of 10 × buffer, 2 U of Taq DNA polymerase, 10 mM dNTPs, each primer at 8 µÌ, and double-distilled H 2 O to a final reaction volume of 50 µL. The   second PCR was performed under the conditions as the first one except that 30 cycles was used. The PCR products were examined on ethidium bromide-stained 2% agarose gels by electrophoresis and visualized under UV light. To minimize heterogeneity, three replicates of the second PCR products from each sample were pooled together. This DNA mix was purified with a AxyPrepDNA Gel Extraction Kit (AXYGEN, United States) and equal concentrated amplicons were sequenced

Bioinformatics Analyses
The sequence reads were processed with a QIIME toolkit (Caporaso et al., 2010). First, sequence reads were quality-filtered and demultiplexed according to the following criteria: minimum length ≥ 200 bp (excluding barcode and primer sequences); ambiguous bases ≤ 0; homopolymer length ≤ 10 bp; maximum number of primer or barcode mismatches ≤ 0; and minimum mean quality score ≥ 30 in a window of 50 nt. The 80 root samples yielded a total of 1 757 618 sequences after quality control, with an average sequence length of 291 bp. Using the de novo approach, we detected 9 232 chimeric sequences, and the remaining non-chimeric sequences yielded 164 initial OTUs based on 97% sequence similarity using the USEARCH algorithm. The OTUs sequences were then subjected to a BLASTN search of the SILVA database (Pruesse et al., 2007). Through comparison using the SILVA database, 1 323 639 sequences were identified as Glomeromycota. However, the pot-U1 and pot-R4 samples only yielded 3 and 7 Glomeromycota sequences and were excluded from the subsequent processing (Supplementary Table S1). Clusters representing < 0.0001% of the total reads were removed in subsequent analysis to avoid α diversity overestimation due to sequencing errors (Bokulich and Mills, 2013). For the comparative analyses, the OTU table was rarefied without replacement to 4 227 reads (the lowest number of sequences among all the samples) per sample. Multiple sequence alignment and construction of the maximum likehood (ML) phylogenetic tree were performed using MEGA6.0 (Tamura et al., 2013). The robustness of this phylogenetic tree was evaluated using 1 000 bootstrap replications. The raw reads have been submitted to the NCBI (National Center for Biotechnology Information) Sequence Read Archive (SRA) database under accession numbers PRJNA64452.

Statistical Analyses
Using the vegan package (Oksanen et al., 2019) of R software (https://cran.r-project.org), we calculated Shannon's diversity . We used the picante package (Kembel et al., 2010) to calculate Faith's PD phylogenetic diversity (PD), where S is the number of OTUs in each sample and p i representes a single OTU sequence abundance ratio of each sample's total abundance. A linear mixed effects model (LMM) in the nlme package (Pinheiro et al., 2018) was used to measure differences in plant heavy metal contents, plant individual biomass, relative abundance of fungal OTUs, and fungal alpha diversity between the urban and rural samples, with the land use type (urban/rural) set as a fixed effect and plant species set as a random effect. Difference of soil properties in rural areas and those in urban areas were tested by one-way ANOVA.
The impacts of urbanization on AM fungal community structure were measured by Adonis (permutational multivariate analysis of variance), ANOSIM (analysis of similarity), and MRPP (multiresponse permutation procedure) tests via the vegan package, with 999 random permutations. The homogeneity of the variances were checked using the function betadisper in vegan package. Fungal dissimilarity among 8 host species were measure by Bray-Curtis distance (BC). BC = n k = 1 x ik −x jk n k = 1 x ik +x jk where i and j represented different host species, x was the sequence number of the k-th OTU, and n was the sum of the total number of OTUs in all the samples. Difference between the fungal dissimilariy (BC) among host plants in rural areas and that in urban areas was tested by the paired sample t test.
The original observed data consisted of a matrix with herbaceous plants in rows and AM fungi being aggregated within each plant species in columns. Therefore, the network in this study represents plant-symbiotic fungal associations (Chagnon et al., 2012;Toju et al., 2016). The NODF (nestedness metric based on paired overlap and decreasing fill) was used to measure the nestedness of mycorrhizal interactions   where marginal total of row j (column l) less than marginal total of row i (column k) DF paired = 0; marginal total of row j (column l) no less than marginal total of row i (column k) DF paired = 100; PO is percentage of co-present 1's in column (row) pairs to those in a column k (row i). In this study, all AM fungi in rural or urban areas were considered during nesting calculation, for no significant alpha and beta diversity difference being observed between soil AM fugnal pools (Lin et al., 2020). NODF calculations and significances tests using 50 replicates of column-random row (EE) nulls and 50 replicates of fixed column-fixed row (FF) nulls were conduceted in ANINHADO software (Guimarães and Guimarães, 2006). EE null: presences are randomly assigned to any cell within the matrix. FF null: the probability of a cell a ij show a presence is P i C P j R /2, in which P i is the number of presences in the row i; P j is the number of presences in the column j; C is the number of columns; R is the number of rows. The NODF_c is propsed to compare nestedness level betweeInclude the following n netoworks by Song et al. (2017). NODF_c = NODF/max (NODF)/ C * log(S) , where Frontiers in Microbiology | www.frontiersin.org max(NODF) is the maximum nestedness in a network with given rows, columns and links; C is the network connectance; S is the geometric mean number of species in the network. NODFc of observed network and significances tests using 999 replicates of EE nulls were measured via maxnodf package (Christoph and Benno, 2020) and vegan package.
Using bipartite package (Dormann et al., 2018), we calculate the degree, d' , H 2 ' , and C.score metrics to describe the network topology based on the mycorrhizal interaction matrices. Degree is the number of links per species, and values range from 1 (specialist) to the number of hosts (generalist) (Gonzalez et al., 2010). The significances of plant-AM fungal links in networks and d' of AM fungi were tested by using 999 replicates of swap.web null in bipartite package. Standardized effect size (z-score),

obs−mean(nulls) sd(nulls)
, where obs is the observed links (or d' value) of a fungus; mean(nulls) is the mean value of the links (or d' value) from 999 swap.web randomizations; sd(nulls) is the standard deviation of the links (or d' value) from randomizations. Zscore > 1.96 (<-1.96) indicate that the observed links or d' value of a fungus is 0.05 significantly higher (lower) than that of randomizations.
Species specialization (d') & overall degree of specialization (H 2 ') are quantitative indices used to describe the degree of interaction specialization,: d' characterizes specialization at the species level, while H 2 ' characterizes specialization at the network level (Blüthgen et al., 2006). The d' index emphasizes the intensity of a species deviating from a random sampling, calculated as the Kullback-Leibler distance. Values range from 0 (generalist) to 1 (specialist). D' is calculated first by finding d i : d i = c j = 1 p ij ln p ij q j , where c is the range of resources; p' ij is the interactions of species i/ the sum of performances of species i; q j is the sum of interactions of resource j/ the total number of interactions in the matrix.
where d max is given analytically, and d min is found heuristically. H 2 ' is the extension of d' for the entire association network. In this study, H 2 ' calculates how strongly the observed interactions deviating from random expectation with the given species' marginal totals. C.score calculates the (normalised) mean number of checkerboard combinations. Values range from 0 (aggregation, no repelling forces between species) to 1 (evidence for disaggregation, e.g. through competition).

Characterization of Soil Properties
Soil properties were significantly different between rural and urban sites ( Table 1). The pH, total nitrogen (TN), magnesium (Mg), nickel (Ni), copper (Cu), and lead (Pb) contents in urban soil were significantly higher than those in rural soil (P < 0.05). No significant differences between soil organic carbon (SOC), available phosphorus (AP), chromium (Cr), cadmium (Cd), zinc (Zn), and manganese (Mn) contents were observed between urban sites and rural sites. Herbaceous richness in urban areas was significantly lower than that in rural areas (P < 0.05).
The magnesium (Mg) and lead (Pb) contents of urban plants were significantly higher than those of rural plants (P < 0.05). The individual biomass, chromium (Cr), and cadmium (Cd) levels of urban plants were not significantly different from those of rural plants (Figure 2).
The Shannon, Simpson, and Pielou indices of the root AM fungi in the urban areas were significantly higher than those in the rural areas ( Table 2). The communities of urban samples were significantly different from those of rural samples ( Table 3). Dissimilarities in the fungal community and Glomus Group II among 8 host plants in urban areas showed significantly lower values than those in rural areas (Figure 4).

Impact of Urbanization on Network Structure
The plant-AM fungal network in the urban area showed a higher connectance (0.5895) level and lower H2' (0.3972) and C.score (0.1431) values than those in the rural area (connectance = 0.5093, H2' = 0.2619, C.score = 0.0769) ( Table 4). The association network in urban areas showed a significantly lower NODF value (55.50) than that in rural areas (64.68) (Figure 5). The NODF of the rural network was significantly higher than those of the EE nulls and FF nulls (P < 0.05), while the urban network exhibited a non-nested pattern (P > 0.05) ( Table 4). The NODFc value of the urban network (0.3284) was lower than that of the rural network (0.4446) and did not differ from that of the simulated EE nulls (P > 0.05) (Figure 5).
The LMM analysis showed that the d' values of the rural plants were significantly higher than those of the urban plants (P < 0.05) ( Table 5). The pattern of significant plant-AM fungal links differed between the rural and urban networks (Figure 6). Urbanization enhanced the percentage of generalist OTUs (degree = 8) and reduced the proportion of specialist OTUs (degree = 1) (Figure 7).

The Nestedness of Mycorrhizal Networks Differs Between Rural and Urban Areas
In the rural areas, the network of mycorrhizal associations showed a significantly nested pattern (Figures 4, 5). This asymmetric interactive pattern could contribute to the maintenance of biodiversity (Bascompte et al., 2006;Jabot and Bascompte, 2012). In a previous study, a 10-plant-47-AM fungi network in a hemiboreal forest also showed a significantly nested web structure (Chagnon et al., 2012). However, experiments carried out by Toju et al. (2014) suggest a non-nested 33-plant-10-AM fungal network in a temperate secondary forest. This non-nested pattern could be attributed to using non-AM specific primers and nonherbaceous hosts and only 14 out of 33 host plants were found to form associations with AM fungi. Thus, fewer co-shared fungal OTUs and lower network nestedness might be observed. Some meta-community   analyses suggested that soil nutrient conditions could be a driving force of nestedness (Verbruggen et al., 2012). Habitat heterogeneity could result in nonrandom species spatial distribution in pollination networks (Olesen et al., 2007). In this study, plant roots were sampled in a large geographic and coincidental area. No significant correlation was observed between the rank of nestedness of plants and beta diversity of AM fungi (Supplementary Table S2). Thus, the effect of sampling heterogeneity was not the essential factor contributing to this nested pattern of rural networks.  In contrast to the observation that network structure shows some resistance to fluctuations in species and interactions (Petanidou et al., 2008;Montesinos-Navarro et al., 2012), urbanization disrupted the significantly nested pattern of the plant-AM fungal network (Figure 5, Table 4). NODF exhibits positive correlations with network connectance (Almeida-Neto et al., 2008). Bascompte et al. (2003) suggested that an increased number of interactions within a given number of species yields a higher relative nestedness level. However, a similar pattern was only observed in the NODF EE nulls, and the rural network with a lower connectance level showed a higher NODF value than that of the urban network (Table 4). Thus, a large relative nestedness difference existed between the urban and rural networks. The connectance and network size independent metric NODFc also indicated that a less nested pattern was observed in urban areas (0.3284) than in rural areas (0.4446).
In this study, the significantly increased d' of host plants ( Table 5) and significantly lower fungal beta diversity among 8 host plants in the urban areas (Figure 4) than those in the rural areas, indicate a reduced selectivity of urban plants in choosing partners. It is hypothesized that the numbers of nondominant microbial species decline as soil incompatibility increases (Encinas-Viso et al., 2016). In this study, significantly higher soil Pb, Ni and Cu and plant Pb contents were observed in urban areas than those in rural soils and plants (Figure 2, Table 1). The toxicity of soil heavy metals may hinder the ranges of partner selection. Due to similar adaptations to the urban environment, urban plants show a homogenization of the spatial community composition (Knapp et al., 2008;Ramalho et al., 2018), and a similar convergent pattern was also observed in root AM fungal associations.
The frequency of generalist fungal OTUs but not that of specialist fungal OTUs was slightly enhanced in the urban areas (Figure 7). AM fungi do not strictly host specific taxa (Sanders, 2003), and habitat fragmentation during urbanization could obstruct the diffusion and colonization of specialists (Telford et al., 2006). Carbon resource or root niche heterogeneity catalyzes expansion in resource use (niche expansion) and diversification of specialist AM fungi (Dickie, 2007;Yu et al., 2018). Filters during urbanization, such as heavy metal stresses, particle deposition, and habitat fragmentation, contribute to the reduction of plant richness ( Table 1) and homogenization of plant communities (Knapp et al., 2008;Ramalho et al., 2018), thus hindering the plants from associating with specialist AM fungi. Nonetheless, if specialist fungi constitute a small fraction of the overall fungal richness or are rare, nestedness is difficult to detect (Bahram et al., 2014). All these processes may explain the reduction in nested interaction patterns at the network level ( Figure 5, Table 3).

Changes in the Strength of Mycorrhizal Associations Under Urbanization
Significantly higher Shannon's diversity and Simpson's diversity of root AM fungi were observed in urban areas than in rural areas ( Table 2). Many studies have reported the significant role of AM fungi in improving host plant resistance to heavy metal toxicity and the disease and insect resistance of host plants (Hildebrandt et al., 2007;Khade and Adholeya, 2008;Joy, 2013;Schneider et al., 2016). In this study, significantly higher soil nutrient availability and soil Mg, Ni, Cu, and Pb contents were observed in urban areas ( Table 1). Plants are expected to be more dependent on mycorrhizal symbioses in urban areas than in rural areas due to the frequent exposure of plants to biotic and abiotic stresses in urban areas (Pourrut et al., 2011;Alzetta et al., 2012). Thus, the trade-off strategy (Lipson et al., 2009) supports that plants with similar biomass in urban areas yield a higher diversity of AM fungal communities under environmental pressures than those in rural areas ( Table 1, Table 2).
In this study, lower demands of nutrients and higher demands of environmental stress tolerance in urbanized areas suggest a significantly different fungal composition ( Table 3). The relative abundance of Glomus Group II was significantly enhanced during urbanization (Figure 3); this group often appears in degraded ecosystems and has a strong adaptability to environmental changes ( Supplementary Table S3). Furthermore, the beta diversity of Glomus Group II and AM fungi among the 8 host species in urban areas was significantly lower than that in rural areas (Figure 4). Therefore, selective filters of environmental stress resistance or tolerance may play an important role in the establishment of root mycorrhizal associations (Knapp et al., 2008). The patterns of significant plant-AM fungal links in the networks were quite different between the rural and urban areas (Figure 6). Competition is a driving force of niche differentiation (Wright, 2002). In this study, urbanization reduced the level of competition among AM fungi, with a lower C-score value (0.0769) obtained for the urban network than for the rural network (0.1431) ( Table 4). Fewer specialist OTUs and more generalist OTUs were observed in the urban areas (Figure 7). Thus, plants in urban areas may tend to cooperate with many kinds of available, robust AM fungi to reduce the risk of species loss (Table 1 and Figures 5, 7), other than diffusion and selectivity among host species.

CONCLUSION
The rural mycorrhizal association network shows a significantly nested network structure and is highly robust to species extinctions. Heavy metal stress, habitat fragmentation, and host diversity changes during urbanization may hinder fungal specialists and enhance fungal generalists in the network of mycorrhizal associations. Urban plants may tend to cooperate with many kinds of generalist AM fungi to compensate the risks of diffusion isolation and species extinctions. All these reductions in selectivity contribute to the non-nested plant-AM fungal network in urban areas.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm. nih.gov/, PRJNA64452.

AUTHOR CONTRIBUTIONS
LL performed methodology, formal analysis, writing original draft, and writing review and editing. YC performed methodology, investigation, and formal analysis. GX performed methodology. YZ performed methodology. SZ performed methodology. KM performed methodology and writing review and editing. All authors contributed to the article and approved the submitted version. Supplementary Table 2 | Linear regression of the plant rank of nestedness to sampling spatial distance, environmental distance, and fungal variation in urban and rural areas.
Supplementary Table 3 | Description of sequences corresponding to the closest matches from the reference database of the representative sequence of each OTU in Glomus Group II.