The Composition and Diversity of Soil Bacterial and Fungal Communities Along an Urban-To-Rural Gradient in South China

Soil microbes are of great significance to driving the biogeochemical cycles and are affected by multiple factors, including urbanization. However, the response of soil microbes to urbanization remains unclear. Therefore, we designed an urban-to-rural gradient experiment to investigate the response of soil microbial composition and diversity to urbanization. Here, we used a high-throughput sequencing method to analyze the biotic and abiotic effects on soil microbial composition and diversity along the urban-to-rural gradient. Our results showed that soil bacterial diversity was the highest in urban areas, followed by suburban areas, and was the lowest in exurbs; however, fungal diversity did not vary significantly among the three areas. Plant traits, i.e., tree richness, shrub richness, the number of tree stems, diameter at breast height of trees, and soil properties, i.e., pH, soil organic carbon, soil exchangeable calcium and magnesium, and soil water content, were only significantly influenced bacterial diversity, but not fungal diversity. The effect of trees and shrubs was higher than that of herbs on microbial composition. Soil organic carbon, pH, soil available nitrogen, soil exchangeable calcium, and magnesium were the major soil factors influencing the soil bacterial and fungal composition. Soil properties had a greater influence on bacterial than on fungal composition at genus level, while plant traits contributed more to fungal than to bacterial composition at genus level. Our study suggests that the urban-to-rural gradient affect the composition and diversity of bacterial community as well as the fungal composition, but not the fungal diversity.


Introduction
Acceleration of urbanization globally has led to an explosive increase in global urban population. In 2018, the global urban population of 4.2 billion accounted for 55% of the total population worldwide [1]. Despite rapid economic development, urbanization has resulted in a series of ecological challenges, such as urban heat islands [2], increased atmospheric CO 2 concentration [3], nitrogen deposition [4], and decreased biodiversity [5]. However, urban forests play a positive role under the enormous pressure of urbanization. Urban forests are associated with multiple ecological benefits such as absorption of carbon and release of oxygen, regulation of urban microclimate, and diminishing heat island effect, thus alleviating the impacts of urbanization [6].
As an essential component of urban forest ecosystem, soil microbes participate in many ecological processes and play an important role in urban forest ecosystem. Soil microbes drive biogeochemical cycles [7], via litter decomposition [8,9], catalyzing the turnover of soil carbon and nutrients [10] and alleviation of changes induced by urbanization. Therefore, subtle changes in soil microbes may reflect We selected the following study sites: Qifeng Park (113°45′50″-113°46′26″ E, 23°00′50″-23°01′04″ N) in the urban areas, the Tongsha Ecological park (113°45′54″-113°46′45″ E, 22°57′45″-23°58′25″ N) in the suburban areas about 6 km from the city center, and the Dalingshan Forest Park (113°45′56″-113°46′46″ E, 22°52′28″-22°53′09″ N) in the exurban areas about 16 km from the city center. The population densities of the three sites decreased from urban to exurban areas [33]. Qifeng Park is located in the center of Dongguan with a high degree of disturbance. Plantation is the main vegetation type. Tongsha Ecological Park is located in the suburb of Dongguan. The vegetation types include mainly plantations and natural secondary forests. The existing subtropical evergreen broad-leaved forests are mainly artificial, and the natural secondary forests are mainly composed of evergreen shrubs and small trees. Dalingshan Forest Park is located in the exurb of Dongguan. The vegetation types mainly include secondary forests and plantations. The secondary forest is well preserved due to effective long-term protection (Table 1).

Experimental Design, Plant Census, and Soil Sampling
The experiment was carried out in Dongguan in January 2018. Nine independent plots with identical size (20 × 20 m) were randomly designated within each site for vegetation census and soil sampling. The space between each two plots was more than 20 m. In each plot, we recorded all trees The population densities of the three sites decreased from urban to exurban areas [33]. Qifeng Park is located in the center of Dongguan with a high degree of disturbance. Plantation is the main vegetation type. Tongsha Ecological Park is located in the suburb of Dongguan. The vegetation types include mainly plantations and natural secondary forests. The existing subtropical evergreen broad-leaved forests are mainly artificial, and the natural secondary forests are mainly composed of evergreen shrubs and small trees. Dalingshan Forest Park is located in the exurb of Dongguan. The vegetation types mainly include secondary forests and plantations. The secondary forest is well preserved due to effective long-term protection (Table 1).

Experimental Design, Plant Census, and Soil Sampling
The experiment was carried out in Dongguan in January 2018. Nine independent plots with identical size (20 × 20 m) were randomly designated within each site for vegetation census and soil sampling. The space between each two plots was more than 20 m. In each plot, we recorded all trees (diameter at breast height (DBH) ≥ 3 cm) including the species name, DBH, and the number of stems. Four 25 m 2 (5 × 5 m) subplots and five 1 m 2 (1 × 1 m) subplots were selected in each plot for investigation of shrubs and herbs, respectively. In the shrub quadrats, the individual species name and their number were recorded. In the herb quadrats, the species name and cover were recorded. Plant traits used for analyses included diameter at breast height of trees (DBH), the number of tree stems (TS), tree richness (TR), shrub richness (SR), and herb richness (HR) ( Table S1).
For soil sampling, after removing the litter layer, six samples were collected in an "S" shape (0-10 cm) and thoroughly mixed. Then, we removed impurities such as plant roots and stones to provide a composite soil sample in each plot. The composite soil sample was divided into two subsamples: one (>6 g) was immediately transferred into the portable ice box designed for microbial identification and another (>500 g) carried in a sealed bag for the determination of soil chemical properties. We also obtained soil samples from each plot using a cut ring (50 mm diameter) to measure the soil water content. All soil samples were collected within a single day.

Soil Physicochemical Properties
Soil water content (SWC) was calculated gravimetrically via oven drying to constant mass at 105 • C. Soil pH was measured via glass electrode method in a 1:2.5 (v/v) soil:water suspension. We determined soil organic carbon (SOC) using the K 2 Cr 2 O 7 oxidation method and soil available nitrogen (AvN) via alkaline hydrolysis diffusion method. Soil available potassium (AvK), soil exchangeable calcium (Ca 2+ ) and soil exchangeable magnesium (Mg 2+ ) were determined via flame atomic absorption spectrophotometry after extraction with 1 mol·L −1 ammonium acetate. Soil available phosphorus (AvP) was determined using the molybdenum-antimony colorimetric method after extraction with 0.03 mol·L −1 ammonium fluoride and 1 mol·L −1 hydrochloric acid [34] (Table S1).

Soil Microbial Communities
Soil microbial DNA was extracted from fresh soil using a MOBIO PowerSoil ® DNA Isolation Kit (MOBIO laboratories, Carlsbad, CA, USA) for the corresponding sample according to the manufacturer's instructions. The concentration and purity of the extracted DNA were quantified using the NanoDrop One (Thermo Fisher Scientific, MA, USA). The 16S rRNA (for bacteria) and ITS3 (for fungi) genes of distinct regions were amplified using a specific primer with 12 bp barcode. Primers were synthesized by Invitrogen (Invitrogen, Carlsbad, CA, USA). PCR reactions, containing 25 µL 2x Premix Taq (Takara Biotechnology, Dalian Co. Ltd., China), 1 µL of each primer (10 mM) and 3 µL DNA (20 ng/µL) template in a volume of 50 µl, were amplified by thermocycling: 5 min at 94 • C for initialization; 30 cycles of 30 s each for denaturation at 94 • C, 30 s annealing at 52 • C, and 30 s extension at 72 • C; followed by a 10 min final elongation step at 72 • C.
PCR products were analyzed and purified by 1% agarose gel electrophoresis for further experiments. The selected PCR products were mixed in equidensity ratios according to the GeneTools Analysis Software (Version4.03.05.0, a division of Synoptics, Cambridge, England). The PCR product mixtures were purified with EZNA Gel Extraction Kit (Omega, USA). Sequencing libraries were finally generated using NEBNext ® Ultra™ DNA Library Prep Kit for Illumina ® (New England Biolabs, MA, USA), and sequenced on an IlluminaHiseq2500 platform to generate 250 bp paired-end reads (Guangdong Magigene Biotechnology Co.,Ltd. Guangzhou, China).
The paired-end raw reads were filtered to obtain the high-quality clean reads according to the Trimmomatic (V0.33, http://www.usadellab.org/cms/?page=Trimmomatic, Aachen, Germany). The paired-end clean reads were merged according to the degree of overlap between the paired-end reads by FLASH (V1.2.11, https://ccb.jhu.edu/software/FLASH/, MD, USA), and the spliced sequences designated as Raw Tags. Using Mothur software (V1.35.1, http://www.mothur.org, MA, USA), the sequences were assigned to each sample based on their unique barcode and primer, after which the barcodes and primers were removed to obtain the effective Clean Tags. Usearch software (V10, http://www.drive5.com/usearch/, CA, USA) was used to select operational taxonomic units (OTU) by combining the reads of clustered OTUs with 97% similarity. Finally, the normalized (subsampled) OTU table was obtained according to the sample with the least sequences. The OTU table with annotations of fungal taxonomy was then used to analyze ecological groups of fungi with FUNGuild software (MN, USA) [35]. The OTU table with annotations of bacterial taxonomy was then used with FAPROTAX (Vancouver, Cananda) software to find the bacterial functional group for each OUT [36].

Statistical Analyses
The Chao1 index and Shannon index of bacteria (Chao1 B and Shannon B ) and fungi (Chao1 F and Shannon F ) were calculated with QIIME (V1.9.1, http://qiime.org, CA, USA) based on bacterial and fungal OTUs. The Chao1 index indicated the microbial species richness, and the Shannon index represented both the microbial species richness and evenness.
Kruskal-Wallis nonparametric test was conducted to analyze the differences in soil properties, plant traits, microbial diversity, and the relative abundance of the main bacterial and fungal phyla and genera along the urban-to-rural gradient. The differences of both bacterial and fungal composition (OTUs) along the urban-to-rural gradient and the differences of both bacterial and fungal functional groups along the urban-to-rural gradient were detected by permutational multivariate analysis of variance (PerMANOVA) with 999 permutations. Redundancy analysis (RDA) was used to determine the correlation between the composition of bacterial and fungal communities, and soil properties and plant traits. Also, we used the Spearman correlation analyses to determine the relationship of bacterial and fungal diversities, composition, and bacterial and fungal functional groups with soil properties and plant traits. Furthermore, the contribution of soil properties and plant traits to the composition of bacterial and fungal community was identified via variation partitioning analysis (VPA) [37]. All analyses were conducted in R (3.4.4, Vienna, Austria) software with 'agricolae', 'psych', and 'vegan' packages. p value < 0.05 was considered statistically significant.

Changes in Plant and Soil Variables
Plant and soil properties greatly varied along the urban-to-rural gradient displaying different trends in these variations ( Table 2). Among plant traits, TR, SR, and TS significantly increased along the urban-to-rural gradient (p < 0.01), whereas HR and DBH showed the opposite trend. All the soil properties varied significantly along the urban-to-rural gradient (p < 0.05) except AvK. The soil was acidic and the pH decreased along the urban-to-rural gradient. SOC, AvN, and SWC were the highest in exurban areas. However, Ca 2+ and Mg 2+ had the highest values in urban area soils. AvP in the suburban areas was significantly higher than in urban and exurban areas (p < 0.01).

Communities of Soil Bacteria and Fungi
A total of 82,683 bacteria OTUs and 27813 fungi OTUs were identified in soil samples. The main bacterial phyla (relative abundance > 1%) across all samples were Proteobacteria, Acidobacteria, Actinobacteria, Verrucomicrobia, WPS-2, Planctomycetes, Chloroflexi, Bacteroidetes, Gemmatimonadetes, and Latescibacteria. Proteobacteria and Acidobacteria were the most abundant bacterial phyla constituting 40.94% and 25.73%, respectively ( Figure S1A, Table S2). Except for the Acidobacteria, Chloroflexi, Planctomycetes, and Verrucomicrobia, the relative abundance of other main bacterial phyla differed significantly along the urban-to-rural gradient (p < 0.05, Table S2). Additionally, the relative abundance of Actinobacteria, Proteobacteria, and WPS-2 increased significantly along the urban-to-rural gradient, whereas those of Bacteroidetes, Gemmatimonadetes, and Latescibacteria decreased along the gradient (p < 0.05, Table S2). At the genus level, there were 20 main bacterial genera (relative abundance > 1%, Figure 2A). The relative abundance of all the main bacterial genera, except for Occallatibacter, exhibited a significant difference along the urban-to-rural gradient (p < 0.05, Table S2). Nine genera with the most abundance showed different trends. Bryobacter, from Acidobacteria, exhibited the lowest abundance in the exurban areas; however, another genus belonging to Acidobacteria, Candidatus_Solibacter, showed the highest abundance in the suburban areas (p < 0.05, Table S2). Within the Proteobacteria, the abundance of three genera (Acidibacter, Bradyrhizobium, and Roseiarcus) significantly increased from urban to exurban areas. The abundance of other two genera (Burkholderia-Caballeronia-Paraburkholderia and Rhodoplanes) in the urban areas was significantly lower than in suburban and exurban areas (p < 0.05, Table S2). Acidothermus in Actinobacteria also increased in the abundance along the gradient, and ADurb.Bin063-1 in Verrucomicrobia was predominant in suburban areas (p < 0.05, Table S2). The PerMANOVA analyses of bacterial OTUs indicated that the composition of bacterial community significantly differed along the urban-to-rural gradient (p < 0.001, Table S4). exurban areas. The abundance of other two genera (Burkholderia-Caballeronia-Paraburkholderia and Rhodoplanes) in the urban areas was significantly lower than in suburban and exurban areas (p < 0.05, Table S2). Acidothermus in Actinobacteria also increased in the abundance along the gradient, and ADurb.Bin063-1 in Verrucomicrobia was predominant in suburban areas (p < 0.05, Table S2). The PerMANOVA analyses of bacterial OTUs indicated that the composition of bacterial community significantly differed along the urban-to-rural gradient (p < 0.001, Table S4). Among the fungal communities, the main phyla (relative abundance > 1%) were Ascomycota (38.77%), Basidiomycota (50.71%), and Zygomycota (4.12%, Figure S1B, Table S3). Only the relative abundance of Zygomycota showed a significant difference along the urban-to-rural gradient (p < 0.05, Table S3). At the genus level, the fungal genera varied greatly along the urban-to-rural gradient ( Figure 2B, Table S3). There were 19 main fungal genera (relative abundance > 1%), and the relative abundance of almost half of which significantly differed along the urban-to-rural gradient (p < 0.05, Table S3). The relative abundance of Aspergillus, Staphylotrichum, and Trichocladium belonging to Ascomycota, and Tomentella and Trichosporon included in Basidiomycota was significantly higher in the urban than in exurban areas. Howerver, Oidiodendron, from Ascomycota, and Abundisporus, Entoloma, and Veluticeps from Basidiomycota had the highest abundance in exurban areas (p < 0.05, Table S3). The PerMANOVA analyses of fungal OTUs also indicated that the composition of fungal community significantly differed along the urban-to-rural gradient (p < 0.001, Table S4).
Both indices (Chao1 and Shannon) characterizing the diversity of bacterial communities decreased significantly along the urban-to-rural gradient ( Figure 3A,B). At the same time, diversity of fungal communities did not vary significantly along the gradient ( Figure 3C,D). analyses of fungal OTUs also indicated that the composition of fungal community significantly differed along the urban-to-rural gradient (p < 0.001, Table S4).
Both indices (Chao1 and Shannon) characterizing the diversity of bacterial communities decreased significantly along the urban-to-rural gradient ( Figure 3A and 3B). At the same time, diversity of fungal communities did not vary significantly along the gradient ( Figure 3C and 3D).

Bacterial and Fungal Functional Groups
FAPROTAX assigned 108,854, 128,508, and 153,970 squences to 62, 61, and 62 bacterial functional groups in urban, suburban, and exurban areas, respectively. The main six abundant groups were influenced by urban-to-rural gradient ( Figure 4A). The relative abundance of chemoheterotrophy, aerobic_chemoheterotrophy, cellulolysis, and iron_respiration was higher in exurban areas than in suburban and urban areas. However, the relative abundance of predatory_or_exoparasitic increased in urban areas ( Figure 4A). The PerMANOVA analyses indicated that the composition of bacterial functional groups significantly differed along the urban-to-rural gradient (p < 0.01, Table S5).

Bacterial and Fungal Functional Groups
FAPROTAX assigned 108,854, 128,508, and 153,970 squences to 62, 61, and 62 bacterial functional groups in urban, suburban, and exurban areas, respectively. The main six abundant groups were influenced by urban-to-rural gradient ( Figure 4A). The relative abundance of chemoheterotrophy, aerobic_chemoheterotrophy, cellulolysis, and iron_respiration was higher in exurban areas than in suburban and urban areas. However, the relative abundance of predatory_or_exoparasitic increased in urban areas ( Figure 4A). The PerMANOVA analyses indicated that the composition of bacterial functional groups significantly differed along the urban-to-rural gradient (p < 0.01, Table S5).
FUNGuild assigned 32,878, 51,000, and 39,390 sequences to 59, 54, and 58 fungal functional guilds in urban, suburban, and exurban areas, respectively, with the "confidence ranking" of "highly probable" or "probable". There were 15 main guilds ( Figure 4B). For example, the relative abundance of wood saprotroph and ericoid mycorrhizal fungi was the highest in exurban areas than in urban and suburban areas. The guilds of undefined saprotroph, plant pathogen, animal pathogen, and animal pathogen-endophyte-plant saprotroph-soil saprotroph were most abundant in urban areas ( Figure 4B). The PerMANOVA analyses also indicated that the composition of fungal functional guilds were significantly influenced by the urban-to-rural gradient (P < 0.01, Table S5).

Correlation of Characteristics of Soil Bacterial and Fungal Communities with Plant and Soil Variables
According to Spearman correlation analyses, DBH, TS, TR, SR, HR, pH, SOC, Ca 2+ , Mg 2+ , and SWC significantly influenced six bacterial functional groups (p < 0.05, Table S6). For fungal functional groups, the abundance of wood saprotroph was negatively correlated with DBH and pH,  FUNGuild assigned 32,878, 51,000, and 39,390 sequences to 59, 54, and 58 fungal functional guilds in urban, suburban, and exurban areas, respectively, with the "confidence ranking" of "highly probable" or "probable". There were 15 main guilds ( Figure 4B). For example, the relative abundance of wood saprotroph and ericoid mycorrhizal fungi was the highest in exurban areas than in urban and suburban areas. The guilds of undefined saprotroph, plant pathogen, animal pathogen, and animal pathogen-endophyte-plant saprotroph-soil saprotroph were most abundant in urban areas ( Figure 4B). The PerMANOVA analyses also indicated that the composition of fungal functional guilds were significantly influenced by the urban-to-rural gradient (p < 0.01, Table S5).

Correlation of Characteristics of Soil Bacterial and Fungal Communities with Plant and Soil Variables
According to Spearman correlation analyses, DBH, TS, TR, SR, HR, pH, SOC, Ca 2+ , Mg 2+ , and SWC significantly influenced six bacterial functional groups (p < 0.05, Table S6). For fungal functional groups, the abundance of wood saprotroph was negatively correlated with DBH and pH, while positively with TS, TR, SR, SOC, and AvN. The abundance of ericoid mycorrhizal fungi was correlated with plant traits and soil properties except for AvK. The abundance of animal pathogen and animal pathogen-endophyte-plant saprotroph-soil saprotroph was positively correlated with Ca 2+ and Mg 2+ . The abundance of plant pathogen was positively and significantly associated with pH, but negatively with TS, TR, SOC, and AvN (p < 0.05, Table S6).
According to Spearman correlation analyses, the diversity of bacterial communities (Chao1 B , Shannon B ) was strongly correlated with almost all plant and soil variables (Table 3). However, the fungal diversity was not influenced by plant and soil variables (Table 3). The Redundancy analyses (RDA) and Spearman correlation analyses indicated that plant traits including TS, DBH, TR, SR, and HR significantly affected the bacterial composition. SR was the most important factor, while HR had the least effect ( Figure 5B and Figure S2B, Tables S7, S8, and S11). Among soil properties, pH, SOC, Ca 2+ , and Mg 2+ were identified as the most important predictors of bacterial composition ( Figure 5A and Figure S2A, Tables S7, S8, and S11). AvP did not influence the bacterial community significantly (Table S11). The changes in the abundance of Proteobacteria, Bacteroidetes, and WPS-2 were related to plant traits (SR, TR, TS, and DBH) and soil properties (pH, Ca 2+ , Mg 2+ , SOC, and SWC) ( Figure S2A,B, Table S7). Moreover, the abundance of Acidibacter and Roseiarcus in Proteobacteria was positively associated with SWC, AvN, SOC, SR, TR, and TS, but negatively correlated with pH, Ca 2+ , Mg 2+ , and DBH; however, the effect of the observed factors on the other three genera (Ellin6067, Haliangium, and MND1) belonging to Proteobacteria was opposite ( Figure 5A,B, Table S8).  The plant and soil properties had no effect on fungal community at the phylum level (p > 0.05, Figures S2C and S2D). However, plant traits including TS, DBH, TR, SR, and HR significantly affected the fungal composition at the genus level. SR and TS influenced fungal composition the most at the genus level ( Figure 5D, Tables S10 and S11). Among soil properties, SOC, pH, AvN, Ca 2+ , and Mg 2+ most significantly affected the fungal composition at the genus level ( Figure 5C, Tables S10 and S11). Similarly, AvP did not influence the fungal community significantly (Table S11). In addition, only the abundance of Zygomycota was positively and significantly correlated with SOC, AvN, SWC, and SR, but negatively with pH and DBH (Table S9). SOC, AvN, SWC, TS, TR, and SR positively influenced the aboundance of Oidiodendron, from Ascomycota, whereas pH, Ca 2+ , Mg 2+ , and DBH had negative influence. However, SOC, AvN, SWC, TS, TR, and SR negatively influenced The plant and soil properties had no effect on fungal community at the phylum level (p > 0.05, Figure S2C,D). However, plant traits including TS, DBH, TR, SR, and HR significantly affected the fungal composition at the genus level. SR and TS influenced fungal composition the most at the genus level ( Figure 5D, Tables S10 and S11). Among soil properties, SOC, pH, AvN, Ca 2+ , and Mg 2+ most significantly affected the fungal composition at the genus level ( Figure 5C, Tables S10 and S11). Similarly, AvP did not influence the fungal community significantly (Table S11). In addition, only the abundance of Zygomycota was positively and significantly correlated with SOC, AvN, SWC, and SR, but negatively with pH and DBH (Table S9). SOC, AvN, SWC, TS, TR, and SR positively influenced the aboundance of Oidiodendron, from Ascomycota, whereas pH, Ca 2+ , Mg 2+ , and DBH had negative influence. However, SOC, AvN, SWC, TS, TR, and SR negatively influenced the aboundance of Staphylotrichum in Ascomycota, pH, Ca 2+ , Mg 2+ , and DBH had positive influence ( Figure 5C,D, Table S10).
Variation partitioning analysis (VPA) indicated that 76% of variation in bacterial community at the phylum level ( Figure S3A), and 77% at the genus level were attributed to plant and soil variables ( Figure 6A). Soil properties contributed to higher variation in bacterial community compared with plant traits (Figure 6 and Figure S3A). However, a large proportion of variation in fungal community could not be explained by plant and soil variables ( Figure 6B and Figure S3B). More importantly, plant and soil variables did not significantly affect fungal community at the phylum level ( Figure S3B), but significantly affected fungal community at the genus level ( Figure 6B). The result showed that soil properties had a greater effect on bacterial than on fungal composition at genus level. By contrast, plant traits contributed more to fungal than to bacterial composition at genus level. the aboundance of Staphylotrichum in Ascomycota, pH, Ca 2+ , Mg 2+ , and DBH had positive influence ( Figures 5C and 5D, Table S10). Variation partitioning analysis (VPA) indicated that 76% of variation in bacterial community at the phylum level ( Figure S3A), and 77% at the genus level were attributed to plant and soil variables ( Figure 6A). Soil properties contributed to higher variation in bacterial community compared with plant traits (Figures 6A and S3A). However, a large proportion of variation in fungal community could not be explained by plant and soil variables (Figures 6B and S3B). More importantly, plant and soil variables did not significantly affect fungal community at the phylum level ( Figure S3B), but significantly affected fungal community at the genus level ( Figure 6B). The result showed that soil properties had a greater effect on bacterial than on fungal composition at genus level. By contrast, plant traits contributed more to fungal than to bacterial composition at genus level.

The Composition and Diversity of Soil Microbial Community Varied Along the Urban-To-Rural Gradient
The composition and diversity of soil microbial communities are closely related to urbanization [38]. Our results confirmed that the urban-to-rural gradient affected the composition of soil microbial community, with varying degrees of influence on bacterial and fungal communities (Figures 2, 3, and S1, Tables S2 and S3). Along the urban-to-rural gradient, the relative abundance of the main bacterial phyla and genera differed significantly (Figures 2 and S1, Table S2). The relative abundance of main fungal genera differed significantly along the urban-to-rural gradient (Figures 2, Table S3). In addition, the PerMANOVA analyses showed that the composition of both bacterial and fungal OTUs had significant difference along the urban-to-rural gradient (Table S4). This result indicates that bacteria respond to urbanization via changes in relative abundance at the phylum and genus level, while fungi respond to urbanization via changes in relative abundance at the genus level. Further, we found that bacterial diversity decreased significantly along the urban-to-rural gradient, whereas the fungal diversity did not change significantly (Figure 3). Our results are similar to the studies of Barrico et al. [29]. With the intensity of urbanization, a few main bacteria OTUs are suppressed, while the survival of other OTUs were promoted in specific habitats. In addition, disturbances can increase environmental heterogeneity, creating more diverse niches for coexistence of other OTUs. Therefore, the diversity of bacteria is higher in urban areas with large disturbance. Moreover, due to the high functional redundancy of

The Composition and Diversity of Soil Microbial Community Varied Along the Urban-To-Rural Gradient
The composition and diversity of soil microbial communities are closely related to urbanization [38]. Our results confirmed that the urban-to-rural gradient affected the composition of soil microbial community, with varying degrees of influence on bacterial and fungal communities (Figure 2, Figure 3, and Figure S1, Tables S2 and S3). Along the urban-to-rural gradient, the relative abundance of the main bacterial phyla and genera differed significantly (Figure 2 and Figure S1, Tables S2). The relative abundance of main fungal genera differed significantly along the urban-to-rural gradient ( Figure 2, Table S3). In addition, the PerMANOVA analyses showed that the composition of both bacterial and fungal OTUs had significant difference along the urban-to-rural gradient (Table S4). This result indicates that bacteria respond to urbanization via changes in relative abundance at the phylum and genus level, while fungi respond to urbanization via changes in relative abundance at the genus level. Further, we found that bacterial diversity decreased significantly along the urban-to-rural gradient, whereas the fungal diversity did not change significantly (Figure 3). Our results are similar to the studies of Barrico et al. [29]. With the intensity of urbanization, a few main bacteria OTUs are suppressed, while the survival of other OTUs were promoted in specific habitats. In addition, disturbances can increase environmental heterogeneity, creating more diverse niches for coexistence of other OTUs. Therefore, the diversity of bacteria is higher in urban areas with large disturbance.
Moreover, due to the high functional redundancy of bacteria community [39], higher bacterial diversity contributes to maintaining the ecosystems stability.
According to the bacterial and fungal functional groups analyses, our results indicated that the urban-to-rural gradient significantly impacted the composition of bacterial and fungal functional groups. The relative abundance of chemoheterotrophy, aerobic_chemoheterotrophy, cellulolysis, wood saprotroph, and ericoid mycorrhizal fungi was the highest in exurban areas. Chemoheterotrophic bacteria and aerobic_chemoheterotrophic bacteria utilize the soil organic carbon which is higher in exurban areas. Because the secondary forest is well preserved in exurban areas, there are more litter and coarse woody debris containing lignin, cellulose, and hemicellulose in exurban areas. Wood saprotroph can decompose lignin, cellulose, and hemicellulose; therefore, the relative abundance of wood saprotroph fungi was higher in exurban areas. In addition, the spearman correlation analyses showed wood saprotroph fungi was highly correlated with SR and TR, wich showed the highest value in exurban areas (Table S6). Ericoid mycorrhizal fungi, limited of Ericaceae [40], have a strong ability to take up phosphorus directly from organic matters [41]. Our study showed that ericoid mycorrhizal fungi were negatively correlated with available phosphorus (Table S6), and the available phosphorus was the lowest in exurbs (Table 1). Moreover, the Ericaceae plants are abundant in exurban areas, especially Rhododendron moulmainense. Therefore, the relative abundance of ericoid mycorrhizal fungi was the highest in exurban areas. Moreover, the higher relative abundance of plant and animal pathogens in urban areas indicates that urbanization can increase the potential risk of pathogen infection.
Our results indicated that the urban-to-rural gradient significantly affected the composition of soil bacterial communities at the level of phyla, genera, and OTUs and fungal communities at the level of genera and OTUs. The bacterial diversity was significantly influenced by the urban-to-rural gradient, while the fungal diversity was not.

Plant and Soil Variables Significantly Affect Soil Bacterial Diversity, but Not Fungal Diversity
Along the urban-to-rural gradient, differences in plant and soil variables explain changes in soil microbial communities [29]. Our results showed that soil bacterial diversity was the highest in urban areas and the lowest in exurban areas. However, fungal diversity did not vary significantly along the urban-to-rural gradient ( Figure 3) because urbanization-induced changes in plant and soil variables (TR, SR, TS, DBH, pH, SOC, AvN, Ca 2+ , Mg 2+ , SWC) significantly affect bacterial diversity but not fungal diversity. Generally, soil microbes interact with plants, with synergistic as well as positive or negative feedback effects [42]. Therefore, we expect a significant correlation between microbial diversity and plant diversity. However, we found that only bacterial diversity was significantly related to TR, SR, TS, DBH, while fungal diversity was not correlated with plant traits. Hu et al. [20] also found that fungal diversity was weakly related to plant diversity in the forests of five climate regions in China; similarly, plant diversity did not determine fungal diversity at a global scale [43]. However, Hiiesalu et al. [44] found that plant diversity determined fungal diversity at a regional scale in temperate pine forests. This finding indicates that the relationship between fungal diversity and plant communities is very complex and depends on spatial scales.

Soil Properties Affect Bacterial Composition More Than Fungal Composition, While Plant Traits Affect Fungal Composition More Than Bacterial Composition
The changes in soil microbial communities along the urban-to-rural gradient and the relationship between soil microbial diversity and plant and soil variables are described above. However, the link between soil microbial composition and plant and soil variables remains unclear. In this study, we found that plant and soil variables explained large variation in the bacterial community composition, with 76% and 77% at the phylum and genus level, respectively ( Figure 6A and Figure S3A). However, plant and soil variables did not significantly affect fungal community at the phylum level and explained 30% of variation in fungal community composition at the genus level ( Figure 6B and Figure S3B).
Both RDA and VPA showed that explanatory variables contributed to significant variation in fungal composition at the genus level ( Figure 5C,D and Figure 6B). Studies have shown that fungal community composition is mainly affected by climatic factors (rainfall and temperature) [18]. The heat island effect caused by urbanization raises the urban temperature significantly higher than that of suburban and exurban areas [2]. Therefore, temperature may be an important factor affecting the composition of fungal community in urban forest soils.
We found that TR and SR had a greater impact on soil bacterial and fungal composition, whereas HR had less impact (Tables S7-S10). Therefore, woody plants had a stronger effect than herbs on soil microbial community. The study of Hu et al. [20] reported similar results. We also found that the effect of plant diversity on soil fungal composition was stronger than on bacterial composition at genus level ( Figure 6). Fungi are the major organisms that can decompose woody litter containing lignin, which accounts for 60-75% of woody litter affected by plant diversity [45]. Moreover, mycorrhizal fungi develop a symbiosis with plants [46]. Therefore, plant diversity explained more variations in fungal composition (13% of variation at genus level, Figure 6B) than in bacterial composition (5% of variation at genus level, Figure 6A). However, only a small proportion of variation in fungal composition was explained by plant diversity. Perhaps plant composition greatly influences fungal composition [17,20].
Moreover, soil pH and SOC largely affect soil bacterial community composition. Soil pH, SOC, and AvN were the important factors influencing soil fungal community composition (Table S11), which was similar to the results of other studies [21,23]. Notably, we found strong effects of Ca 2+ and Mg 2+ on soil bacterial and fungal composition (Table S11). Barrico et al. [29] also showed that Mg 2+ was an important factor affecting the fungi community. Our research sites were located in urban parks and subjected to human disturbance, such as artificial fertilization. Since calcium and magnesium are essential trace elements in plants, MgSO 4 is often used by gardeners for the management of garden plants, releasing large amounts of magnesium within a short time [47]. Treatment of ornamental plants in urban parks with artificial fertilizers raises the soil content of magnesium and calcium. Soil microbes facilitate the absorption of trace elements by plants. Therefore, in urban forests, due to human disturbance (fertilizers, etc.), the content of soil Ca 2+ and Mg 2+ is changed, resulting in a strong correlation between soil microbial community and levels of Ca 2+ and Mg 2+ . Overall, these results support the hypothesis that plant and soil variables significantly affected the composition of soil bacterial and fungal composition. The influence of trees and shrubs on the composition of soil microbial community was larger than that of herbs. Soil properties exert larger effect on bacterial than on fungal composition at genus level, while plant variables had larger effect on fungal than on bacterial composition at genus level.

Conclusions
The results demonstrated that urbanization affected the composition and diversity of soil bacteria and fungi. The composition of soil bacterial communities at the level of phyla, genera, and OTUs and fungal communities at the level of genera and OTUs significantly varied along the urban-to-rural gradient. The bacterial diversity decreased significantly along the urban-to-rural gradient. However, fungal diversity did not significantly differ along the urban-to-rural gradient. Moreover, the composition of both soil bacterial and fungal communities varied in response to plant and soil variables. Soil properties had a greater effect on bacterial than on fungal composition at genus level. By contrast, plant traits contributed more to fungal than to bacterial composition at genus level. Trees and shrubs had a higher impact than herbs on the microbial community composition. Among soil properties, in addition to pH, SOC, and AvN, the changes in soil Ca 2+ and Mg 2+ caused by urbanization had a significant impact on soil bacterial and fungal communities.
Supplementary Materials: The following are available online at http://www.mdpi.com/1999-4907/10/9/797/s1, Figure S1: Relative abundance of bacterial phyla (A) and fungal phyla (B) along the urban-to-rural gradient at the Dongguan city, southern China, Figure S2: RDA plot showing influence of soil properties and plant traits on the distribution of main phyla in soil bacterial (A and B) and fungal (C and D) communities, Figure S3: The proportions of variation at phylum level in composition of soil bacterial community (A) and soil fungal community (B) composition explained by soil and plant variables, Table S1: Abbreviations and corresponding full names of plant traits and soil properties in this paper, Table S2: Relative abundance of bacterial composition along the urban-to-rural gradient at phylum and genus level, Table S3: Relative abundance of fungal composition along the urban-to-rural gradient at phylum and genus level, Table S4. The PerMANOVA analyses of bacterial and fungal OTUs along the urban-to-rural gradient, Table S5: The PerMANOVA analyses of bacterial and fungal functional groups along the urban-to-rural gradient; Table S6. Spearman correlation of bacterial and fungal functional groups with plant and soil variables, Table S7: Spearman correlation of bacterial composition with plant and soil variables at phylum level, Table S8: Spearman correlation of bacterial composition with plant and soil variables at genus level, Table S9: Spearman correlation of fungal composition with plant and soil properties at phylum level, Table S10: Spearman correlation of fungal composition with plant traits and soil properties at genus level, Table S11: The results of RDA monte carlo tests of soil bacterial and fungal community compositions.
Author Contributions: L.Z. and X.T. conceived and designed the experiments. X.T. and L.K. performed field surveys and collected the data. X.T. analyzed the data, prepared the figures and tables, and wrote the first draft. Z.S., X.L. and L.Z. reviewed and edited the manuscript.