Jin Zhou
Ziyan Liao
Zhidong Liu1,3
Youhua Chen- 1China-Croatia “Belt and Road” Joint Laboratory on Biodiversity and Ecosystem Services, Chengdu Institute of Biology, Chinese Academy of Sciences, Chengdu, China
- 2Key Laboratory of Bio-Resources and Eco-Environment of Ministry of Education, College of Life Sciences, Sichuan University, Chengdu, Sichuan, China
- 3University of Chinese Academy of Sciences, Beijing, China
Urbanization not only profoundly alters landscape profiles, ecosystems and vertebrate faunal diversity but also disturbs microbial communities by increasing stochasticity, vulnerability, biotic homogenization, etc. However, because of the buffering effect of host species, microbial communities are expected to be influenced by both host species and urbanization stresses. Therefore, the impacts of urbanization on animals’ microbial symbionts could be more complex and uncertain. In this study, we quantified the urbanization degree of sampling sites and surveyed the gut and skin microbes of three amphibian host species in different sites in urban parks and nearby villages of Chengdu, Southwest China. Furthermore, a co-occurrence network analysis, the phylogenetic normalized stochasticity ratio and Sloan neutral community models were applied to infer the impact of urbanization on symbiotic microbial communities. For the three host species, urbanization increased the diversity of symbiotic microbes and the number of keystone microbial taxa. However, the negative effects of such increased diversification were evident, as the community stochasticity and co-occurrence network structure vulnerability also increased, while the network structure complexity and stability were reduced. Finally, the community stochasticity had positive associations with the network vulnerability, implying that the existence of many transient symbiotic rare microbial taxa in urban parks makes the symbiotic microbial community structure more fragile. Conclusively, urbanization increased the symbiotic microbial diversity at the cost of community stability; the results provide a new perspective for better understanding the complex triangulated environment–host–microbe relationship.
1. Introduction
Urbanization is one of the most direct manifestations of human activities (Rothwell et al., 2015; Zhai et al., 2020) and has been a major cause of the biodiversity crisis, which has remarkably reduced the suitable habitats available for wild animals to survive and breed (Grimm et al., 2008; McDonald et al., 2008; Fenoglio et al., 2020; Li et al., 2022). Various studies have reported the effects of urbanization on macroecological biodiversity patterns (Kondratyeva et al., 2020; Li et al., 2022; Pan et al., 2022). However, few studies (Lokmer et al., 2020; Barnes et al., 2021) have reported the impacts of urbanization on the biological diversity and composition of microorganisms associated with the organs of animal hosts.
Cities are expected to expand due to increase in the population through immigration from neighboring wild areas and the aggregation of supporting infrastructure, particularly in developing countries (Henderson, 2002; Xu et al., 2014; Zhai et al., 2020). Therefore, urbanization is inevitable and will be intensified. Over the past decade, Chengdu has been the fastest-growing city in Southwest China, with a total population of over 20 million. Consequently, suitable habitats for wild animals to inhabit within city boundaries have been greatly reduced (Kong et al., 2010). Furthermore, as a result of increased urbanization, these ‘green spaces’ are continuously fragmented and isolated, leaving wild animals vulnerable to human disturbance (Gaston, 2010) (e.g., pollution, roadkill, and illegal injury) and making it less likely to avoid population decline or local extinction through migration (Marzluff and Rodewald, 2012; Bendik et al., 2014). This equally applies to microorganisms of host species. Because of the difficulty of finding suitable habitats for animal host species in cities, microbial communities of the hosts are also greatly selected and modulated by urban environments in comparison to those in rural and wild habitats (Lokmer et al., 2020; Barnes et al., 2021).
Over one-third of amphibian species are threatened by urbanization (Hamer and McDonnell, 2008), which may be because amphibians are the most sensitive taxa in response to climate and urban environments compared with others such as birds and mammals (Hamer and McDonnell, 2008; Magle et al., 2012; Callaghan et al., 2021). Small patches and fragmented habitats caused by urbanization can remarkably affect the home range of animals and, as a result, their reproduction and survival, particularly for amphibians (Hinam and Clair, 2008; Decout et al., 2012; Amburgey et al., 2021). However, as amphibians are small-bodied animals, they can survive in small patches of forested and humid areas of ‘green spaces’, provided that these areas are not heavily polluted or disturbed (Hamer and McDonnell, 2008). The impact of urbanization on the biodiversity (Callaghan et al., 2021; Yang, L. et al., 2022) and genetic structure (Wei et al., 2021) of amphibian populations has been extensively studied, but microbial data records associated with amphibian hosts in urban settings are lacking.
Symbiotic microorganisms are critical for amphibian survival (Becker et al., 2015) and adaptation to rapid environmental changes (Xu et al., 2020b; Barnes et al., 2021). Consequently, it is crucial to elucidate the response of the amphibian symbiotic microbial community to rapid urbanization, which is valuable for wildlife conservation in urban settings (Wei et al., 2021). However, urban ecosystems present distinctive microhabitat and microclimate conditions in comparison to wild regions (Zhang et al., 2022). Therefore, it is expected that the microorganism community structure and diversity of amphibian hosts must also alter and adapt to the urban environment. However, to the best of our knowledge, few studies have reported on the potential distinctive features (e.g., stochasticity and stability) of the microbial community structure of animal hosts in urban settings (see concept Figure 1A), especially for amphibians, which are indicator species for assessing urbanization’s ecological impact (Guzy et al., 2012).
Figure 1. Community composition of microbes of amphibians in the urban parks and wild area habitats. (A) Schematic illustration of the effects of urbanization on the community network stability and assembly mechanisms of amphibian symbiotic microorganisms. (B) Diagram of grouping. (C) Map of sampling sites in city parks and wild areas of Chengdu city in China based on the land use types of the dataset. The green points represent urban land, including the 51 (urban land), 52 (rural residential areas), and 53 (other construction land) land use types. The purple and blue points represent the U (urban park samples) and W (wild area samples) groups, respectively. (D) Comparison of the taxonomic composition of the top 10 phyla and family between the U and W groups. BG, Bufo gargarizans; FM, Fejervarya multistriata; PN, Pelophylax nigromaculatus. G, gut; SK, skin; W_G, gut microbiota of amphibians in wild areas; W_SK, skin microbiota of amphibians in wild areas; U_G, gut microbiota of amphibians in urban parks; U_SK, skin microbiota of amphibians in urban parks.
The process of community assembly is generally used to uncover the formation mechanisms of microbial diversity and composition differences, which include deterministic processes based on niche theory and stochastic processes based on neutral theory (Stegen et al., 2013, 2015; Ning et al., 2020). Theories based on niches contend that species distribution is primarily controlled by the deterministic processes involve with abiotic factors (e.g., environmental factors) and biotic factors (e.g., biological interactions) (Stegen et al., 2013; Ning et al., 2020). In contrast, neutral theories view all species as ecologically equivalent, with stochastic processes largely determining species dynamics and patterns, such as birth/death, speciation/extinction, and migration (Stegen et al., 2013; Ning et al., 2020). Furthermore, ecological stability is vital for microbial communities to adapt to rapid environmental changes (Yuan et al., 2021; Liu et al., 2022). Generally, microbial network analysis is used to evaluate the interspecific interactions and ecological stability or complexity of microbial communities (Yuan et al., 2021; Liu et al., 2022).
In summary, we aimed to answer the following questions: (i) Does urbanization affect the diversity and community structure of symbiotic microbes? (ii) What is the actual influence of urbanization on the stochasticity and stability of amphibian symbiotic microbial communities? (iii) What is the relationship between the stochasticity and stability of symbiotic microbial communities under urbanization gradient? To this end, we chose three amphibian host species (Pelophylax nigromaculatus, abbr: PN; Fejervarya multistriata, abbr: FM; Bufo gargarizans, abbr: BG) that are widely distributed in the city parks and surrounding wild areas of Chengdu, southwest China, to elucidate the community assembly and diversity patterns of microbial taxa in the skins and guts of the host species. We compared the microbial differences between amphibians in urban and wild environments, between host species and between different organs of the host species. This work explores the urbanization-host animal-microbial relationship with the aim of providing new perspectives on bioenvironmental interactions in the context of rapid global land use change.
2. Materials and methods
2.1. Sample collection and sequencing
From April to August 2021, we collected 194 microbial samples from three frog species in downtown Chengdu city and surrounding wild areas, which included four city park sampling sites and eight wild sampling sites (that is, exurban regions) (Figure 1B; Supplementary Table S1). For each sampling site, we searched for frog species alongside ditches, brooks or ponds. We captured adult individuals by wearing sterile gloves and took the microbial samples from the skin by using cotton swabs after rinsing three times with sterile water to remove the potential transient bacteria (Lauer et al., 2007). A 2-ml aseptic centrifuge tube was used to collect gut microbial content after the frog was euthanized and dissected. After returning to the laboratory, all gut and skin samples were immediately stored at –80°C. Experimental approval was obtained from the Animal Ethical and Welfare Committee of Chengdu Institute of Biology, Chinese Academy of Sciences (permit no. CIBDWLL2022008).
Microbial DNA was extracted using the MN NucleoSpin 96 Soil kit (MACHEREY-NAGEL) according to the manufacturer’s instructions. The V3-V4 regions of the 16S rRNA gene were amplified using the universal primer pair 338F/341F and 806R/805R (Yao et al., 2019; Gutierrez-Villagomez et al., 2021). Briefly, PCR thermocycling was carried out under the following conditions: 95°C for 5 min, followed by 30 cycles at 95°C for 30 s, 50°C for 30 s, 72°C for 60 s and a final extension at 72°C for 5 min. The PCR products were further purified with e.Z.N.A. TM Cycle-Pure Kit (Omega).
The raw sequence data were processed using the QIIME2 software package (version QIIME2-2021.2) (Bolyen et al., 2019). The paired reads (2 × 250 bp paired-end mode) generated by the Illumina platforms of MiSeq and Nova6000 were demultiplexed and filtered by vsearch and quality filtered with the following quality-control criteria: sequences of length < 200 bp and average quality score < 25 were removed; no ambiguous bases were allowed. High-quality sequence reads with >97% identity were defined as operational taxonomic units (OTUs) with the vsearch cluster-features-de-novo (Rognes et al., 2016). Then, the samples were rarefied to the same sequence depth (33,187 bacterial sequences per sample). To minimize inclusion of sequencing errors, the clustered feature tables were further filtered using QIIME2 feature-table filter-features (Revise: p-min-samples 2 --p-min-frequency 10) (Zhang et al., 2020). Taxonomy assignment to the OTUs was performed using the Silva v138 database and the Naive Bayes classifier (Quast et al., 2012). Subsequently, the sequences identified as chloroplasts or mitochondria were removed. The final dataset retained 19,042 OTUs for the downstream analyses. In our study, we classified the collection of microbial samples from urban park samples and wild area samples into U and W groups, respectively (Figure 1B). We also separated the microbial samples into two groups for three host species in the analyses. Group 1 contained all the microbes from all the organs (gut and skin) of the three host species, including W_BG, U_BG, W_FM, U_FM, W_PN, and U_PN (Figure 1C). Group 2 distinguished the microbial communities from the skin (abbr: SK) and gut (abbr: G) ecosystems, including W_BG_G, W_BG_SK, U_BG_G, U_BG_G _SK, W_FM_G, W_FM_SK, U_FM_G, U_FM_SK, W_PN_G, W_PN_SK, U_PN_G, and U_PN_SK (Figure 1B). In this study, Group 1 was the primary research target. The analyses were conducted and compared at the OTU, family and phylum levels.
2.2. Measures of diversity and community structure
All result analysis and visualization were performed in R v4.1.3.1 For the alpha diversity, we computed the richness, Shannon index and phylogenetic diversity for each community and used the Wilcoxon rank-sum test and the false discovery rate (FDR) corrected method to identify significant differences between the U and W groups.
For the beta diversity, we compared the difference in beta diversity built on the Bray–Curtis distance among the six groups by using Kruskal–Wallis tests with multiple testing correction Behrens–Fisher using the nmpc function in the ‘nmpc’ R package (Munzel and Hothorn, 2001). Additionally, PERMANOVA and ANOSIM tests (number of permutations: 999) based on the Bray–Curtis distance were used to analyze the differences in gut and skin microbes between the Group U and Group W amphibians. Nonmetric multidimensional scaling (NMDS) was used to visualize the results using the ‘vegan’ R package (Oksanen et al., 2002).
We assessed whether neutrality is the principal driver of the microbial community structure in urban and wild environments. To determine the strength of the neutral stochasticity, we utilized the Sloan neutral community model (NCM) (Sloan et al., 2006; Chen et al., 2019). This method evaluates the relationship between occupancy and species abundance, in which R2 indicates the overall fit to the neutral model (Sloan et al., 2006; Chen et al., 2019). In this model, the parameters Nm and m represent the metacommunity size times migration and migration rate, respectively, where higher m values indicate less dispersal limitation (Burns et al., 2016; Chen et al., 2019). To further assess the relative importance of stochastic versus deterministic processes in the community assembly, the phylogenetic normalized stochasticity ratio (pNST, ranging from 0 to 1) index was calculated using the pNST function in the ‘NST’ package in R (Ning et al., 2019, 2020). The pNST was used with 50% as the boundary point between stochastic assembly (>50%) and deterministic assembly (<50%). The total differences in the pNST between the U and W groups were tested by Kruskal–Wallis tests with multiple testing correction (Behrens–Fisher) using the nmpc function in the ‘nmpc’ R package (Munzel and Hothorn, 2001). In addition, the relationships between the pNST and urbanization, robustness and vulnerability were estimated with Spearman’s correlation analyses.
To evaluate the levels of OTU species co-occurrence patterns, the co-occurrence networks were constructed based on robust correlations with Spearman’s correlation coefficients (r) of >0.6 or <−0.6 and false discovery rate-corrected p-values of <0.01 using the ‘WGCNA’ R package (Langfelder and Horvath, 2008). To minimize the impacts of rare OTUs in the dataset, only OTUs with relative abundance >0.01% and the bacterial taxa that were detected in at least 10% of all samples per group were included in the analysis (Jiao et al., 2021). Subsequently, visualization of networks and the calculation of network properties were performed using the ‘igraph’ R package (Csardi and Nepusz, 2006). Additionally, we extracted the subnetwork topological characteristics of each sample by using the subgraph functions in the ‘igraph’ R package (Csardi and Nepusz, 2006; Ma et al., 2016). In this study, the network characterizations of the node numbers, edge numbers, average degree, degree centralization, graph density, average path length, clustering coefficient, graph modularity, and betweenness centralization were calculated to estimate the network structure differences in different groups. In addition, the robustness based on the method of 50% of nodes was randomly removed with 100 repetitions of the simulation, and the vulnerability was also calculated to assess the network stability (or ecological stability) for the symbiotic microbiota of amphibians from urban park and wild area habitats (Yuan et al., 2021; Liu et al., 2022). The total differences in the network characterization in different groups were tested by using Kruskal–Wallis tests with multiple testing correction (Behrens–Fisher) using the nmpc function in the ‘nmpc’ R package (Munzel and Hothorn, 2001). Furthermore, the relationships between the urbanization and robustness and vulnerability were assessed with Spearman’s correlation analysis.
To decipher changes in keystone taxa between Group 1 levels, we identified the putative role of each node based on two network topological characteristics of the within-module connectivity (Zi) and the among-module connectivity (Pi) following the previous classification: peripherals (Zi < 2.5 and Pi < 0.62), connectors (Zi < 2.5 and Pi > 0.62), module hubs (Zi > 2.5 and Pi < 0.62) and network hubs (Zi > 2.5 and Pi > 0.62) (Olesen et al., 2007). The connectors, module hubs and network hubs are generally defined as keystone taxa, which represent key species that may play crucial roles in maintaining the stability of the microbial community structure (Röttjers and Faust, 2019; Yuan et al., 2021; Liu et al., 2022).
2.3. Impacts of abiotic factors
We tested the impacts of abiotic variables on the microbial community structure and diversity. To do so, we utilized a set of variables, including 19 bioclimatic variables such as the latitude, longitude, urban_degree (degree of urbanization), elevation, PH, AT (air temperature) and AH (air humidity) (Supplementary Table S2). Bioclimatic variables were extracted from the WorldClim database.2 PH was measured using a portable HANNA HI 9913001/1 instrument, and AT and AH were detected every 10 min by using an automatic temperature and humidity recorder (Elitech GSP-8A) at each sampling site. In addition, we quantified the degree of urbanization by calculating the total proportion of three land use types, i.e., urban land (code 51), rural residential areas (52) and other construction land (53) to all land use classes within a 10-km radius buffer around each sampling site. Land use class data from 2020 (approximate to our sampling time) were obtained from the Resource and Environment Science and Data Center.3
Because the utilization of all variables could lead to multicollinearity, to reduce the risk of overfitting, we first utilized random forest (RF) to identify potentially highly correlated and important variables, which can indicate the importance of variables by the percentage increase in the MSE (mean squared error) (Jiao et al., 2021). In the RF models, the environmental variables are referred to as predictors for symbiotic microbial richness. Then, we selected the top 10 most important environmental variables for collinearity analysis using the check_collinearity function in the ‘performance’ R package (Lüdecke et al., 2021). Finally, we retained five environmental variables, the PH, AT, urban_degree, bio_3 (isothermality) and bio_18 (precipitation of warmest quarter) (Supplementary Table S2), for subsequent analysis. In this work, the relative abundance of the top 10 phyla and families and network characteristic correlations with environmental variables were calculated using the corrplot function in the ‘corrplot’ package (Wei et al., 2017). Moreover, redundancy analysis (RDA) was applied to disentangle the relationship between the symbiotic microbes of amphibians and environmental variables with the rda function of the ‘vegan’ package (Oksanen et al., 2002).
3. Results
3.1. Overview of the amphibian symbiotic bacterial composition and diversity
In terms of the community composition of host microbes, at the phylum level, for the skin and gut microbiota of amphibians, the bacterial composition was primarily dominated by Proteobacteria, Firmicutes, Actinobacteriota, and Bacteroidota. Specifically, for the skin microbiota of the three amphibians, the relative abundance of Proteobacteria, Firmicutes, and Actinobacteriota in wild areas (W_SK) was significantly higher than that in urban parks (U_SK) (Figure 1D). Conversely, for the gut microbiota, the relative abundance of Proteobacteria in wild areas (W_G) was significantly lower than that in urban parks (U_G), while the relative abundance of Firmicutes in the W_G group was significantly higher than that in the U_G group (Figure 1D). Furthermore, in contrast to the skin microbiota, Desulfobacterota and Verrucomicrobiota were also major phyla in the gut microbiota. Analogously, the relative abundance of Desulfobacterota in the W_G group was significantly higher than that in the U_G group, in addition to Bufo gargarizans (Figure 1D).
Interestingly, at the family level, the relative abundance of families in different groups presented remarkable microbial community composition differences (Figure 1D). For skin microbiota, Enterobacteriaceae, Comamonadaceae, and Burkholderiaceae were major dominant families. Additionally, in addition to Comamonadaceae and Burkholderiaceae, the relative abundance of Caulobacteraceae and Sphingomonadaceae of the three amphibians in the wild areas (W_SK) was also significantly higher than that in the urban parks (U_SK), whereas Enterobacteriacea had a higher relative abundance in the U_SK group (Figure 1D). For gut microbiota, Enterobacteriaceae was the dominant family in the urban parks (U_G), while Lachnospiraceae, Erysipelotrichaceae, Ruminococcaceae, and Desulfovibrionaceae were the major dominant families in the wild areas (W_G; Figure 1D).
In regard to the biological diversity of the microbes, at the alpha diversity level, the overall microbial diversity was consistently high in the U groups for the species-level richness (Figure 2A), Shannon index (Figure 2B), and phylogenetic diversity (Figure 2C). However, when comparing the alpha diversity difference for the skin and gut of host species, in general, skin microbes presented higher diversity (Supplementary Figure S2). This might be attributed to the fact that skin microbes are jointly influenced by host species and external abiotic environments. However, at the beta diversity level, the beta diversity built on the Bray-Curtis distance of the three amphibians was significantly higher in the W group than in the U group (Kruskal–Wallis tests, p < 0.001; Figure 2D; Supplementary Figure S3). Non-metric multidimensional scaling (NMDS) analysis was used to further disentangle this dissimilarity. The results show that there was clear segregation between the U and W groups, which presented significant differences (PERMANOVA, p < 0.001; ANOSIM, R = 0.12, p < 0.001; Figure 2E).
Figure 2. Microbial diversity of amphibians in the urban park and wild area habitats. Comparison of the alpha diversity between Group U (urban park samples) and Group W (wild area samples) amphibian symbiotic microorganisms, including the richness (A), Shannon index (B), and PD phylogenetic diversity index (C). The difference in the alpha diversity was tested using the Wilcoxon rank-sum test and false discovery rate (FDR) at the corrected level of p < 0.05. (D) Differences in the beta diversity between the U and W groups based on the Bray–Curtis distances. The different superscripts indicate significant differences between groups (p < 0.05, Kruskal–Wallis tests and Behrens–Fisher post hoc tests). (E) Nonmetric multidimensional scaling (NMDS) analysis based on the Bray–Curtis distance of the relative abundance of operational taxonomic units (OTUs). Differences are denoted as follows: ∗p < 0.05; ∗∗p < 0.01; ∗∗∗p < 0.001. W_BG, W_FM, and W_PN represent the symbiotic microbiota of Bufo gargarizans, Fejervarya multistriata, and Pelophylax nigromaculatus in the wild areas, respectively. U_BG, U_FM, and U_PN indicate the symbiotic microbiota of Bufo gargarizans, Fejervarya multistriata, and Pelophylax nigromaculatus in the urban parks, respectively.
3.2. Urbanization influences the symbiotic bacterial community assembly of amphibians
Here, the NCM model was used to effectively evaluate the relationship between the occurrence frequency and relative abundance of OTUs. The results show that the explained proportion of variation (R2) of the symbiotic microbial community in the urban parks (U group: U_BG, R2 = 0.626; U_FM, R2 = 0.602; U_PN, R2 = 0.679) versus the wild areas (W group: W_BG, R2 = −0.027; W_FM, R2 = 0.357; W_PN, R2 = 0.214) present a remarkable difference. The greater R2 values of the U group indicate that stochastic processes played a more important role in the U group than in the W group (Figure 3A). This further suggests that urbanization enhanced the ecological stochasticity of the microbial community assembly. Furthermore, the Nm and m values were higher for microbial taxa in the U group (U_BG, Nm = 2,477, m = 0.078; U_FM, Nm = 1931, m = 0.061; U_PN, Nm = 1,583, m = 0.05) than in the W group (W_BG, Nm = 702, m = 0.023; W_FM, Nm = 919, m = 0.03; W_PN, Nm = 693, m = 0.023) (Figure 3A). These results illustrate that species dispersal of symbiotic microbes was higher in the U group than in the W group, indicating that the microbial taxa of the W group had a higher chance of dispersal limitation. Additionally, similar results were found in the gut and skin symbiotic microbe taxa between the U and W groups (Group 2 level) (Supplementary Figure S4).
Figure 3. Microbial assembly mechanisms of amphibians in urban park and wild area habitats. (A) Neutral community model (NCM) applied to estimate the influence of random dispersal and ecological drift on the community assembly of the symbiotic microbiota of amphibians from urban parks and wild areas. Nm is the metacommunity size times immigration, R2 is the fit to the neutral model, and m indicates the migration rate. The solid yellow line indicates the best fit to the neutral model, while a dashed yellow line indicates 95% confidence intervals around the prediction. The different colors are used for OTUs that occur more or less frequently than predicted by the NCM. (B) The ecological stochasticity in the assembly of the amphibian symbiotic microbiota community evaluated by the phylogenetic normalized stochasticity ratio (pNST). The value of 0.5 is regarded as the boundary between a stochastic (>0.5) and deterministic (<0.5) assembly. The different superscripts indicate significant differences between groups (p < 0.05, Kruskal–Wallis tests and Behrens–Fisher post hoc tests). (C) Regression relationships between the pNST and urban_degree (urban degree of sampling sites). (D,E) Regression relationships between the pNST and the robustness and vulnerability, respectively. W_BG, W_FM, and W_PN represent the symbiotic microbiota of Bufo gargarizans, Fejervarya multistriata, and Pelophylax nigromaculatus in the wild areas, respectively. U_BG, U_FM and U_PN indicate the symbiotic microbiota of Bufo gargarizans, Fejervarya multistriata, and Pelophylax nigromaculatus in the urban parks, respectively.
The pNST model was further used to quantify the relative importance of stochastic and deterministic processes of the symbiotic microbe community assembly. The results show that the pNST value of the W group was <0.5, while that of the U group was >0.5, which indicates that stochastic processes and deterministic processes were the predominant processes for the W and U groups, respectively (Figure 3B). Moreover, significant differences in the pNST values were found between the U and W groups (Kruskal–Wallis tests, p < 0.001; Figure 3B). Importantly, the pNST was significantly positively correlated with urban_degree (Figure 3C), which further identifies the effect of urbanization on microbial community assembly processes. These results indicate that urbanization significantly alters microbial community assembly processes and increases stochastic processes in amphibian symbiotic microbes. These results are in line with previous NCM findings.
3.3. Urbanization influences the network structure of the amphibian symbiotic bacterial community
Using network analyses, we further identified important evidence for community simplicity, with a well-fitted power-law model (R2 = 0.995) for the degree distribution of the total network, which indicates that the network was scale-free and nonrandom (Barberan et al., 2012). The edge numbers were higher in the U groups than in the W group, whereas the node numbers were lower (Figure 4A). Importantly, the co-occurrence instances (edge/node ratio) in the U group (U_BG, edge/node = 19.32; U_FM, edge/node = 3.26; U_PN, edge/node = 5.97) were higher than those in the W group (W_BG, edge/node = 86.70; W_FM, edge/node = 27.11; W_PN, edge/node = 36.18), which indicates that the microbial taxa of the W group tend to co-occur more tightly and frequently (Figure 4A; Shu et al., 2021). The average degree and graph density, representing the network complexity (Wu et al., 2016), decrease from the W group to the U group, indicating that the W group had higher network complexity than the U group (Figure 5). The degree centralization shows similar patterns, suggesting that the network structures of the U group are more ‘even’ networks (Wu et al., 2016). Compared with the U group, the W group had a higher clustering coefficient and shorter average path length, indicating that the small-world properties of the W group were typical (Figure 5). Additionally, the betweenness centralization, which reflects the role of a node as a bridge between components of network (Banerjee et al., 2019), was higher in the U group than in the W group, excluding the FM amphibian (Figure 5).
Figure 4. Co-occurrence networks of symbiotic microbial taxa of amphibians in the urban park and wild area habitats. (A) Network diagram of the symbiotic microbiota of amphibians in urban parks and wild areas. The robustness (B) was calculated as the proportion of the remaining species in the community after random node 50% removal, and the vulnerability (C) was determined by the maximum node vulnerability in each community. (D) The regression relationships between urban_degree (urban degree of sampling sites) and robustness and vulnerability. The difference in robustness was tested using the Wilcoxon rank-sum test and false discovery rate (FDR) at the level of corrected p < 0.05. Differences are denoted as follows: ∗p < 0.05; ∗∗p < 0.01; ∗∗∗p < 0.001. W_BG, W_FM, and W_PN represent the symbiotic microbiota of Bufo gargarizans, Fejervarya multistriata, and Pelophylax nigromaculatus in the wild areas, respectively. U_BG, U_FM, and U_PN indicate the symbiotic microbiota of Bufo gargarizans, Fejervarya multistriata, and Pelophylax nigromaculatus in the urban parks, respectively.
Figure 5. Comparison of different network characterizations between the U and W groups, including the node number (nodes_num), edge number (edges_num), average degree (degree), degree centralization, graph density, average path length, clustering coefficient, graph modularity, and betweenness centralization. The different superscripts indicate significant differences between groups (p < 0.05, Kruskal–Wallis tests and Behrens–Fisher post hoc tests). W_BG, W_FM, and W_PN represent the symbiotic microbiota of Bufo gargarizans, Fejervarya multistriata, and Pelophylax nigromaculatus in the wild areas, respectively. U_BG, U_FM, and U_PN indicate the symbiotic microbiota of Bufo gargarizans, Fejervarya multistriata, and Pelophylax nigromaculatus in the urban parks, respectively.
As indicators of the network stability, the robustness and vulnerability (Yuan et al., 2021) were calculated. Our results show that the robustness was significantly higher (all p < 0.001) in the W group than in the U group (Figure 4B), while the vulnerability had the opposite trend (Figure 4C), suggesting that the U group had lower network stability than the W group. More importantly, the correlation results show that urbanization (urban_degree) was remarkably negatively correlated with the robustness (R = −0.82, p < 0.01), whereas vulnerability presented the opposite trend (R = 0.94, p = 0.017) (Figure 4D), which further suggests that urbanization significantly altered the robustness and vulnerability of the network. The above findings strongly indicate that urbanization profoundly affects microbial associations, significantly reduces complexity and stability, and enhances the vulnerability of the symbiotic microbial community.
The keystone taxa were identified to further estimate the effect of urbanization on the network structure of symbiotic microbes by using Zi–Pi relationships, with the connectors, module hubs and network hubs regarded as keystone taxa. The results reveal that the periphery had the most nodes in all groups, and the numbers of keystone taxa were higher in the U group (U_BG, n = 713; U_FM, n = 729; U_PN, n = 1,088) than in the W group (W_BG, n = 99; W_FM, n = 178; W_PN, n = 124) (Figures 6A–C), which demonstrates that urbanization notably increased the numbers of keystone taxa. Additionally, we found that the Acidobacteriota, Bacteroidota, Firmicutes, Actinobacteriota, Proteobacteria, and Planctomycetota play critical roles in maintaining the ecological stability of commensal microbes (Figures 6A,B). Dechloromonas at the genus level plays a dominant role in symbiotic microbes under urban settings (Figure 6A). Importantly, the degree of keystone taxa in group U was significantly lower than in group W, which may imply that many rare taxa play a more important role in regulating the ecological stability of symbiotic microbes in the urban setting (Figures 6A,B).
Figure 6. The keystone taxa distribution of amphibian symbiotic microbiota in urban parks and wild area habitats. The phylogenetic tree of keystone species of connectors (A) and module hubs (B). Colours of the tip point indicate the phyla distribution. Blue, black, and orange dots indicate the keystone taxa that occurred in both urban and wild settings, only in urban setting (U group), and only in wild areas (W group), respectively. The outermost bar plot represents the degree of keystone taxa, and the annotations correspond to genus-level species. (C) Identifying keystone taxa across habitats based on their topological roles in networks. The pie chart shows the proportion of each keystone taxa. n denotes the number of keystone taxa, including the module hubs, network hubs, and connectors. W_BG, W_FM, and W_PN represent the symbiotic microbiota of Bufo gargarizans, Fejervarya multistriata, and Pelophylax nigromaculatus in the wild areas, respectively. U_BG, U_FM, and U_PN indicate the symbiotic microbiota of Bufo gargarizans, Fejervarya multistriata, and Pelophylax nigromaculatus in the urban parks, respectively.
3.4. Impacts of urbanization and environmental variables on the symbiotic bacterial community of amphibians
The RF analysis show that bio_3, AT, urban_degree and PH were the four most important variables for symbiotic microbial richness (Supplementary Figure S1A). Specifically, correlations between the symbiotic bacterial composition and urbanization (urban_degree) and environmental factors were detected. At the phylum level (Figure 7A), the urban_degree, AT, bio_18, bio_3 and pH were significantly correlated (all p < 0.05) with the relative abundance of Proteobacteria, Firmicutes, and Desulfobacterota, respectively. At the family level (Figure 7B), the urban_degree and all variables presented significant correlations with the relative abundance of the top 10 families, excluding Enterobacter australis, for which urban_degree presented a remarkable positive correlation.
Figure 7. The relationship between the symbiotic microbes of amphibians and environmental variables. Correlations between environmental variables and the top 10 phyla (A) and families (B). (C) Correlations between network characterizations and environmental variables. The size and color of the circles show the strength and sign of the correlation, respectively. Differences are denoted as follows: ∗p < 0.05; ∗∗p < 0.01; ∗∗∗p < 0.001. (D) Redundancy analysis (RDA) of the community structure of symbiotic microbiota amphibians and environmental variables, including bio_18 (precipitation of warmest quarter), bio_3 (isothermality), PH, AT (air temperature), and urban_degree (degree of urbanization). W_BG, W_FM, and W_PN represent the symbiotic microbiota of Bufo gargarizans, Fejervarya multistriata, and Pelophylax nigromaculatus in the wild areas, respectively. U_BG, U_FM, and U_PN indicate the symbiotic microbiota of Bufo gargarizans, Fejervarya multistriata, and Pelophylax nigromaculatus in the urban parks, respectively. The full names of the environmental factors are given in Supplementary Table S2.
Additionally, we examined the correlations between the network topological features and environmental variables (Figure 7C). The results show that AT and bio_18 were significantly positively correlated with nodes_num, the average path length and graph modularity, and the correlations with the degree, graph density and degree centralization presented the opposite trend. Moreover, PH was negatively correlated with nodes_num, the average path length and graph modularity, while it was positively correlated with the degree, graph density and degree centralization. Importantly, significant positive correlations between the degree and graph density and urban_degree were observed, indicating that urbanization remarkably affects the complexity of microbial communities. Furthermore, correlations between urban_degree and edges_num and the degree centralization presented similar patterns, while urban_degree was significantly negatively correlated with nodes_num, the average path length, graph modularity and betweenness centralization.
The RDA further revealed the effect of urbanization and environmental variables on the symbiotic bacterial community (Figure 7D). The results show that five variables were correlated with the microbial community, especially PH, AT and urban_degree, with the RDA axis accounting for the total variance and RDA axis 2 accounting for 20.64%. The PH, urban_degree and bio_3 were positively correlated with the microbial community in the U group, while bio_18 and AT were negatively correlated. Moreover, the W group showed opposite correlations.
Finally, the relationship between the pNST and robustness was further identified. The results show that the pNST was remarkably positively correlated with the network vulnerability (R = 0.89, p = 0.033; Figure 3D), yet the robustness was slightly negatively correlated with the pNST (R = −0.66, p = 0.18; Figure 3E), which could indicate a certain degree of correlation between the stochasticity and network stability.
4. Discussion
The environmental disturbances, high spatial heterogeneity, and habitat fragmentation caused by urbanization (Gaston, 2010; Kondratyeva et al., 2020) can remarkably affect the composition and diversity of an animal’s symbiotic microbiome and adaption to complex urban settings (Lokmer et al., 2020; Barnes et al., 2021). In this study, we found that urbanization increased the richness and homogenization of microbes in amphibian hosts. First, in the alpha diversity analyses, we found that the microbial richness, abundance-based diversity (Shannon) and phylogenetic diversity tended to be higher in urban green spaces than in wild habitats (Figure 2A). This suggests that urbanization significantly increases the diversity of amphibian commensal microbes, while altering their kinship and phylogenetic diversity. Importantly, these results align with previous ecological theory that more diverse communities are resistant to invasion and further protect the host (Costello et al., 2012) as urban wildlife has more enhanced immune function, pathogen prevalence and stress than wildlife in wild areas (Murray et al., 2019). However, in our beta analyses, we found that the beta diversity of microbial taxa on the skin, in the gut or both organs generally presented a marked decrease in urban host specimens. Furthermore, ordination analyses (NMDS) showed that the microbial composition from urban samples tended to group together, but the microbial taxa from wild samples tended to be disperse in the ordination space. Therefore, it can be concluded that the microbial taxa of amphibian hosts presented high heterogeneity in wild areas and high homogeneity in urban areas. One explanation for this could be that the urban environments have similar local species pools and high dispersal rates, which drive the homogenization of symbiotic microorganisms (Figure 3A). This may suggest that there is remarkable habitat variation between sample sites in the wild, which leads to significant differentiation of symbiotic microorganisms. In contrast, urban habitats tend to have more similar habitats, which results in similar environmental pressures and symbiotic microbial convergence. More importantly, our results showed that a high alpha diversity does not necessarily imply high community compositional differences or high biotic heterogeneity for host-specific microorganisms.
Furthermore, we revealed the sensitivity and importance of amphibian skin as a first line of defense (Cramp et al., 2014; Xu et al., 2020a) in terms of community structural differences and community assembly. For instance, our results showed that the relative abundance of the remaining families (Figure 1D) and the three alpha indices (Supplementary Figure S2) of the skin microbes differed more between urban settings and wild areas than the gut microbes. At the community assembly (Supplementary Figure S4), we found that urbanization notably shifts differential patterns of stochasticity (R2) and dispersal rates (m) on the skin and gut symbiotic microbes. These results suggest that the environmental stressors imposed by urbanization markedly alter the conditions and niche spaces of the skin and gut, with the skin microbes becoming more sensitive (Xu et al., 2020a) and having a higher rate of dispersal and stochasticity versus the gut (Xu et al., 2022). Importantly, these findings will enhance our understanding of the mechanisms underlying the effects of urbanization on amphibians and contribute to the conservation of sensitive taxa.
Previous studies have revealed that stochastic processes dominate the assembly of skin microbial communities of amphibians along land use gradients, especially drift (Barnes et al., 2021). The stochasticity includes three assembly processes: homogeneous dispersal (as high dispersal rates can homogenize communities and cause high taxonomic turnover), dispersal limitation (as low dispersal rates could enhance community taxonomic alters) and ecological drift (Stegen et al., 2015; Ning et al., 2020). In this study, we found that urbanization increased the microbial stochasticity of amphibian hosts by transforming the dominant deterministic processes in wild areas into dominant stochastic processes in urban settings (Figure 3B). This change could be attributed to higher heterogeneity of conditions and the environmental filtering ability of different sampling sites in wild areas versus urban parks, while the sampling sites in urban parks have similar environmental stress and more numbers of rare or keystone taxa (Figure 6C). Interestingly, our results also verified that the amphibians in urban settings had lower abundance of existence of transient unique or rare symbiotic microbial taxa, such as more network nodes (Figures 4A, 5) or keystone taxa (Figures 6A,B) with no or low degree and a higher proportion of the remaining other families with lower abundance (Figure 1D).
The urban habitat is a separate local species pool, containing many new or unique microbial taxa, which provides a greater diversity of symbiotic taxa for the animals that live there (Barnes et al., 2021; Dillard et al., 2022). Previous studies show that rare taxa are susceptible to local exclusion by drift or competition with other dominant taxa (Melbourne and Hastings, 2008). However, in this study, the passive dispersal could play crucial roles in community assembly (Ai et al., 2013; Jia et al., 2018), as the diverse and distinctive urban local species spool may bring more low-abundance rare or unique taxa to the amphibians, thus causing high turnover and increasing the stochasticity (Stegen et al., 2015). These results align with the previous studies in which the assembly of rare species is mediated by stochastic processes (Jia et al., 2018). More importantly, this is further corroborated by our results that the symbiotic microbiome in urban parks has higher dispersal rates, indicating that homogeneous dispersal can play a profound role in the community assembly in urban settings. This may be due to frequent disturbances (such as human activities) and highly connected river network systems (Fan et al., 2022), which increase the exposure of amphibians to environmental microbes (especially water microbes) and could further increase dispersal rates. However, this may also introduce more new microbial taxa, thereby increasing the speciation, diversification, priority effects and colonization of new species, which may increase stochasticity (such as drift or homogeneous dispersal) (Vellend, 2010; Zhou and Ning, 2017).
In this study, stochastic features were also indicated by the random combination of occasional unique or rare microbial taxa at the OTU level in different locations and host individuals in urban samples, which is indicated by the phylum-and family-level compositional analyses of OTUs and the network analyses (Figures 4A, 6A,B). These microbial taxa increased the node number in the networks in urban areas, but they had few edges connecting other taxa in the networks (Figures 6A,B). Furthermore, although the average degree and graph density representing the network complexity over the urbanization gradient were markedly reduced, the network modularity, representing functional units (Luo et al., 2006) or niches (Eiler et al., 2012), and keystone taxa numbers exhibited an increasing trend (Figure 6C). This suggests that although urbanization reduces the complexity of the network, it increases the functional units to cope with various risks or pressures brought by the urban environment.
Urbanization can significantly change the ecological stability and complexity of environmental microorganisms (Yuan et al., 2021), but the impact on the ecological stability of symbiotic microorganisms has not been studied, which may be because it is a more complex triangular relationship involving the host, symbiotic microorganisms and environment. Here, we examine the network stability and complexity of a symbiotic community by calculating the robustness, vulnerability and topological properties. We found that urbanization reduces the network complexity and ecological stability, which may be one of the reasons for the increased randomness, since the network structure of microorganisms generally represents a nonrandom association (Barberan et al., 2012).
The key is that urban ecosystems are always disturbed and manipulated by human activities and pollution (Zhai et al., 2020; Yang, Y. et al., 2022). For the host species, they are completely different from wild habitats. For example, water bodies in city rivers, ditches or ponds are usually polluted by sewage emissions. Moreover, they are typically managed by local administrative units for the purpose of improving water quality. As such, microhabitat conditions are easily disturbed or altered because of human interference. Two typical pieces of evidence are the remarkable abiotic factor differences in water bodies and air features of urban settings: PH and AT showed significant correlations with symbiotic microbial composition and network topological properties. Additionally, we found that bio_18 and bio_3 also presented remarkable correlations. In summary, these results suggest the important roles of current and historical environmental factors in driving the dissimilarity of symbiotic communities between urban and wild settings.
More importantly, we quantified the degree of urbanization and revealed that it is significantly positively correlated with the stochasticity and vulnerability and significantly negatively correlated with the complexity and robustness. At the same time, we verified that the stochasticity has a significant negative correlation with the ecological stability, which is an interesting and important result. Consequently, our findings provide valuable and new insight into understanding the complex triangular relationship involving the host, environment, and symbiotic microbes.
It must be noted that there are some uncertainties (discussed in Supplementary Text S1 and Figure S5) in the analysis. For future research, it would be interesting to explore the impacts of urbanization on the host species community structure and how the impacts affect symbiotic microbiota. We argue that such analyses would clarify the triangular relationship between hosts, microbes and abiotic environments. In addition, it is important to test the transmission processes of microbial communities hosted within species at various trophic levels between food webs, which may have potential excavation value for understanding interspecific relationships in the context of global change. Furthermore, rare taxa play a pivotal role in regulating the community assembly and ecological stability of environmental microorganisms (Liang et al., 2020; Shu et al., 2021), but their specific regulatory mechanisms in amphibian symbiotic microbes along an urbanization gradient have not been quantitatively and deeply studied, especially the functional mechanisms and assembly processes of each taxon.
Data availability statement
The data presented in the study are deposited in the National Center of Biotechnology Information (NCBI: https://www.ncbi.nlm.nih.gov/) repository, accession number PRJNA889361.
Ethics statement
The animal study was reviewed and approved by the Animal Ethical and Welfare Committee of Chengdu Institute of Biology, Chinese Academy of Sciences.
Author contributions
JZ analyzed the data, wrote the manuscript, performed the research, investigation, and methodology. ZiL wrote the manuscript. ZhL, XG, and WZ: investigation. YC: wrote the manuscript, methodology, and funding acquisition. All authors contributed to the article and approved the submitted version.
Funding
This work was supported by the Strategic Priority Research Program of the Chinese Academy of Sciences (XDB31000000) and the Second Tibetan Plateau Scientific Expedition and Research Program, Grant Number 2019QZKK0303 the National Key Research and development Program of China (2020YFE0203200).
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2022.1108662/full#supplementary-material
Footnotes
References
Ai, D., Chu, C., Ellwood, M. D. F., Hou, R., and Wang, G. (2013). Migration and niche partitioning simultaneously increase species richness and rarity. Ecol. Model. 258, 33–39. doi: 10.1016/j.ecolmodel.2013.03.001
Amburgey, S. M., Miller, D. A., Rochester, C. J., Delaney, K. S., Riley, S. P., Brehme, C. S., et al. (2021). The influence of species life history and distribution characteristics on species responses to habitat fragmentation in an urban landscape. J. Anim. Ecol. 90, 685–697. doi: 10.1111/1365-2656.13403
Banerjee, S., Walder, F., Büchi, L., Meyer, M., Held, A. Y., Gattinger, A., et al. (2019). Agricultural intensification reduces microbial network complexity and the abundance of keystone taxa in roots. ISME J. 13, 1722–1736. doi: 10.1038/s41396-019-0383-2
Barberan, A., Bates, S. T., Casamayor, E. O., and Fierer, N. (2012). Using network analysis to explore co-occurrence patterns in soil microbial communities. ISME J. 6, 343–351. doi: 10.1038/ismej.2011.119
Barnes, E. M., Kutos, S., Naghshineh, N., Mesko, M., You, Q., and Lewis, J. D. (2021). Assembly of the amphibian microbiome is influenced by the effects of land-use change on environmental reservoirs. Environ. Microbiol. 23, 4595–4611. doi: 10.1111/1462-2920.15653
Becker, M. H., Walke, J. B., Murrill, L., Woodhams, D. C., Reinert, L. K., Rollins-Smith, L. A., et al. (2015). Phylogenetic distribution of symbiotic bacteria from Panamanian amphibians that inhibit growth of the lethal fungal pathogen Batrachochytrium dendrobatidis. Mol. Ecol. 24, 1628–1641. doi: 10.1111/mec.13135
Bendik, N. F., Sissel, B. N., Fields, J. R., O’Donnell, L. J., and Sanders, M. S. (2014). Effect of urbanization on abundance of Jollyville plateau salamanders (Eurycea tonkawae). Herpetol. Conserv. Biol. 9, 206–222.
Bolyen, E., Rideout, J. R., Dillon, M. R., Bokulich, N. A., Abnet, C. C., Al-Ghalith, G. A., et al. (2019). Reproducible, interactive, scalable and extensible microbiome data science using QIIME 2. Nat. Biotechnol. 37, 852–857. doi: 10.1038/s41587-019-0209-9
Burns, A. R., Stephens, W. Z., Stagaman, K., Wong, S., Rawls, J. F., Guillemin, K., et al. (2016). Contribution of neutral processes to the assembly of gut microbial communities in the zebrafish over host development. ISME J. 10, 655–664. doi: 10.1038/ismej.2015.142
Callaghan, C. T., Liu, G., Mitchell, B. A., Poore, A. G. B., and Rowley, J. J. (2021). Urbanization negatively impacts frog diversity at continental, regional, and local scales. Basic Appl. Ecol. 54, 64–74. doi: 10.1016/j.baae.2021.04.003
Chen, W., Ren, K., Isabwe, A., Chen, H., Liu, M., and Yang, J. (2019). Stochastic processes shape microeukaryotic community assembly in a subtropical river across wet and dry seasons. Microbiome 7:138. doi: 10.1186/s40168-019-0749-8
Costello, E. K., Stagaman, K., Dethlefsen, L., Bohannan, B. J., and Relman, D. A. (2012). The application of ecological theory toward an understanding of the human microbiome. Science 336, 1255–1262. doi: 10.1126/science.1224203
Cramp, R. L., McPhee, R. K., Meyer, E. A., Ohmer, M. E., and Franklin, C. E. (2014). First line of defence: the role of sloughing in the regulation of cutaneous microbes in frogs. Conserv. Physiol. 2:cou012. doi: 10.1093/conphys/cou012
Csardi, G., and Nepusz, T. (2006). The igraph software package for complex network research. Int. J. Complex Syst. 1695, 1–9.
Decout, S., Manel, S., Miaud, C., and Luque, S. (2012). Integrative approach for landscape-based graph connectivity analysis: a case study with the common frog (Rana temporaria) in human-dominated landscapes. Landsc. Ecol. 27, 267–279. doi: 10.1007/s10980-011-9694-z
Dillard, B. A., Chung, A. K., Gunderson, A. R., Campbell-Staton, S. C., and Moeller, A. H. (2022). Humanization of wildlife gut microbiota in urban environments. eLife 11:e76381. doi: 10.7554/eLife.76381
Eiler, A., Heinrich, F., and Bertilsson, S. (2012). Coherent dynamics and association networks among lake bacterioplankton taxa. ISME J. 6, 330–342. doi: 10.1038/ismej.2011.113
Fan, Y., Zheng, J., Deng, L., Rao, W., Zhang, Q., Liu, T., et al. (2022). Spatiotemporal dynamics of microplastics in an urban river network area. Water Res. 212:118116. doi: 10.1016/j.watres.2022.118116
Fenoglio, M. S., Rossetti, M. R., Videla, M., and Baselga, A. (2020). Negative effects of urbanization on terrestrial arthropod communities: a meta-analysis. Glob. Ecol. Biogeogr. 29, 1412–1429. doi: 10.1111/geb.13107
Grimm, N. B., Faeth, S. H., Golubiewski, N. E., Redman, C. L., Wu, J., Bai, X., et al. (2008). Global change and the ecology of cities. Science 319, 756–760. doi: 10.1126/science.1150195
Gutierrez-Villagomez, J. M., Patey, G., To, T. A., Lefebvre-Raine, M., Lara-Jacobo, L. R., Comte, J., et al. (2021). Frogs respond to commercial formulations of the biopesticide Bacillus thuringiensis var. israelensis, especially their intestine microbiota. Environ. Sci. Technol. 55, 12504–12516. doi: 10.1021/acs.est.1c02322
Guzy, J. C., McCoy, E. D., Deyle, A. C., Gonzalez, S. M., Halstead, N., and Mushinsky, H. R. (2012). Urbanization interferes with the use of amphibians as indicators of ecological integrity of wetlands. J. Appl. Ecol. 49, 941–952. doi: 10.1111/j.1365-2664.2012.02172.x
Hamer, A. J., and McDonnell, M. J. (2008). Amphibian ecology and conservation in the urbanising world: a review. Biol. Conserv. 141, 2432–2449. doi: 10.1016/j.biocon.2008.07.020
Henderson, V. (2002). Urbanization in developing countries. World Bank Res. Obs. 17, 89–112. doi: 10.1093/wbro/17.1.89
Hinam, H. L., and Clair, C. S. (2008). High levels of habitat loss and fragmentation limit reproductive success by reducing home range size and provisioning rates of northern saw-whet owls. Biol. Conserv. 141, 524–535. doi: 10.1016/j.biocon.2007.11.011
Jia, X., Dini-Andreote, F., and Falcao Salles, J. (2018). Community assembly processes of the microbial rare biosphere. Trends Microbiol. 26, 738–747. doi: 10.1016/j.tim.2018.02.011
Jiao, S., Zhang, B., Zhang, G., Chen, W., and Wei, G. (2021). Stochastic community assembly decreases soil fungal richness in arid ecosystems. Mol. Ecol. 30, 4338–4348. doi: 10.1111/mec.16047
Kondratyeva, A., Knapp, S., Durka, W., Kühn, I., Vallet, J., Machon, N., et al. (2020). Urbanization effects on biodiversity revealed by a two-scale analysis of species functional uniqueness vs redundancy. Front. Ecol. Evol. 8:73. doi: 10.3389/fevo.2020.00073
Kong, F., Yin, H., Nakagoshi, N., and Zong, Y. (2010). Urban green space network development for biodiversity conservation: identification based on graph theory and gravity modeling. Landsc. Urban Plan. 95, 16–27. doi: 10.1016/j.landurbplan.2009.11.001
Langfelder, P., and Horvath, S. (2008). WGCNA: an R package for weighted correlation network analysis. BMC Bioinform. 9, 1–13. doi: 10.1186/1471-2105-9-559
Lauer, A., Simon, M. A., Banning, J. L., André, E., Duncan, K., and Harris, R. N. (2007). Common cutaneous bacteria from the eastern red-backed salamander can inhibit pathogenic fungi. Copeia 2007, 630–640. doi: 10.1643/0045-8511(2007)2007[630:CCBFTE]2.0.CO;2
Li, G., Fang, C., Li, Y., Wang, Z., Sun, S., He, S., et al. (2022). Global impacts of future urban expansion on terrestrial vertebrate diversity. Nat. Commun. 13:1628. doi: 10.1038/s41467-022-29324-2
Liang, Y., Xiao, X., Nuccio, E. E., Yuan, M., Zhang, N., Xue, K., et al. (2020). Differentiation strategies of soil rare and abundant microbial taxa in response to changing climatic regimes. Environ. Microbiol. 22, 1327–1340. doi: 10.1111/1462-2920.14945
Liu, S., Yu, H., Yu, Y., Huang, J., Zhou, Z., Zeng, J., et al. (2022). Ecological stability of microbial communities in Lake Donghu regulated by keystone taxa. Ecol. Indic. 136:108695. doi: 10.1016/j.ecolind.2022.108695
Lokmer, A., Aflalo, S., Amougou, N., Lafosse, S., Froment, A., Tabe, F. E., et al. (2020). Response of the human gut and saliva microbiome to urbanization in Cameroon. Sci. Rep. 10:2856. doi: 10.1038/s41598-020-59849-9
Lüdecke, D., Ben-Shachar, M., Patil, I., Waggoner, P., and Makowski, D. (2021). Performance: an R package for assessment, comparison and testing of statistical models. J. Open Source Softw. 6, 3139. doi: 10.21105/joss.03139
Luo, F., Zhong, J., Yang, Y., and Zhou, J. (2006). Application of random matrix theory to microarray data for discovering functional gene modules. Phys. Rev. 73:31924. doi: 10.1103/PhysRevE.73.031924
Ma, B., Wang, H., Dsouza, M., Lou, J., He, Y., Dai, Z., et al. (2016). Geographic patterns of co-occurrence network topological features for soil microbiota at continental scale in eastern China. ISME J. 10, 1891–1901. doi: 10.1038/ISMEJ.2015.261
Magle, S. B., Hunt, V. M., Vernon, M., and Crooks, K. R. (2012). Urban wildlife research: past, present, and future. Biol. Conserv. 155, 23–32. doi: 10.1016/j.biocon.2012.06.018
Marzluff, J., and Rodewald, A. (2012). Conserving biodiversity in urbanizing areas: nontraditional views from a bird’s perspective. Cities Environ. 1:360. doi: 10.1002/9781118279519.ch20
McDonald, R. I., Kareiva, P., and Forman, R. T. (2008). The implications of current and future urbanization for global protected areas and biodiversity conservation. Biol. Conserv. 141, 1695–1703. doi: 10.1016/j.biocon.2008.04.025
Melbourne, B. A., and Hastings, A. (2008). Extinction risk depends strongly on factors contributing to stochasticity. Nature 454, 100–103. doi: 10.1038/nature06922
Munzel, U., and Hothorn, L. A. (2001). A unified approach to simultaneous rank test procedures in the unbalanced one-way layout. Biom. Z. 43, 553–569. doi: 10.1002/1521-4036(200109)43:5<553::AID-BIMJ553>3.0.CO;2-N
Murray, M. H., Sánchez, C. A., Becker, D. J., Byers, K. A., Worsley-Tonks, K. E. L., and Craft, M. E. (2019). City sicker? A meta-analysis of wildlife health and urbanization. Front. Ecol. Environ. 17, 575–583. doi: 10.1002/fee.2126
Ning, D., Deng, Y., Tiedje, J. M., and Zhou, J. (2019). A general framework for quantitatively assessing ecological stochasticity. Proc. Natl. Acad. Sci. U. S. A. 116, 16892–16898. doi: 10.1073/pnas.1904623116
Ning, D., Yuan, M., Wu, L., Zhang, Y., Guo, X., Zhou, X., et al. (2020). A quantitative framework reveals ecological drivers of grassland microbial community assembly in response to warming. Nat. Commun. 11:4717. doi: 10.1038/s41467-020-18560-z
Oksanen, J., Blanchet, F. G., Kindt, R., Legendre, P., Minchin, P. R., O’Hara, R., et al. (2002). Vegan: community ecology package. Anim. Behav. 63, 1037–1046. doi: 10.1006/anbe.2002.3022
Olesen, J. M., Bascompte, J., Dupont, Y. L., and Jordano, P. (2007). The modularity of pollination networks. Proc. Natl. Acad. Sci. 104, 19891–19896. doi: 10.1073/pnas.0706375104
Pan, Y., Hersperger, A. M., Kienast, F., Liao, Z., Ge, G., and Nobis, M. P. (2022). Spatial and temporal scales of landscape structure affect the biodiversity-landscape relationship across ecologically distinct species groups. Landsc. Ecol. 37, 2311–2325. doi: 10.1007/s10980-022-01477-x
Quast, C., Pruesse, E., Yilmaz, P., Gerken, J., Schweer, T., Yarza, P., et al. (2012). The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res. 41, D590–D596. doi: 10.1093/nar/gks1219
Rognes, T., Flouri, T., Nichols, B., Quince, C., and Mahé, F. (2016). VSEARCH: a versatile open source tool for metagenomics. PeerJ 4:e2584. doi: 10.7717/peerj.2584
Rothwell, A., Ridoutt, B., Page, G., and Bellotti, W. (2015). Feeding and housing the urban population: environmental impacts at the peri-urban interface under different land-use scenarios. Land Use Policy 48, 377–388. doi: 10.1016/j.landusepol.2015.06.017
Röttjers, L., and Faust, K. (2019). Can we predict keystones? Nat. Rev. Microbiol. 17:193. doi: 10.1038/s41579-018-0132-y
Shu, D., Guo, Y., Zhang, B., Zhang, C., Van Nostrand, J. D., Lin, Y., et al. (2021). Rare prokaryotic sub-communities dominate the complexity of ecological networks and soil multinutrient cycling during long-term secondary succession in China's loess plateau. Sci. Total Environ. 774:145737. doi: 10.1016/j.scitotenv.2021.145737
Sloan, W. T., Lunn, M., Woodcock, S., Head, I. M., Nee, S., and Curtis, T. P. (2006). Quantifying the roles of immigration and chance in shaping prokaryote community structure. Environ. Microbiol. 8, 732–740. doi: 10.1111/j.1462-2920.2005.00956.x
Stegen, J. C., Lin, X., Fredrickson, J. K., Chen, X., Kennedy, D. W., Murray, C. J., et al. (2013). Quantifying community assembly processes and identifying features that impose them. ISME J. 7, 2069–2079. doi: 10.1038/ismej.2013.93
Stegen, J. C., Lin, X., Fredrickson, J. K., and Konopka, A. E. (2015). Estimating and mapping ecological processes influencing microbial community assembly. Front. Microbiol. 6:370. doi: 10.3389/fmicb.2015.00370
Vellend, M. (2010). Conceptual synthesis in community ecology. Q. Rev. Biol. 85, 183–206. doi: 10.1086/652373
Wei, X., Huang, M., Yue, Q., Ma, S., Li, B., Mu, Z., et al. (2021). Long-term urbanization impacts the eastern golden frog (Pelophylax plancyi) in Shanghai City: demographic history, genetic structure, and implications for amphibian conservation in intensively urbanizing environments. Evol. Appl. 14, 117–135. doi: 10.1111/eva.13156
Wei, T., Simko, V., Levy, M., Xie, Y., Jin, Y., and Zemla, J. (2017). Package ‘corrplot’. Statistician 56:e24.
Wu, L., Yang, Y., Chen, S., Zhao, M., Zhu, Z., Yang, S., et al. (2016). Long-term successional dynamics of microbial association networks in anaerobic digestion processes. Water Res. 104, 1–10. doi: 10.1016/j.watres.2016.07.072
Xu, L., Chen, H., Zhang, M., Zhu, W., Chang, Q., Lu, G., et al. (2020b). Changes in the community structure of the symbiotic microbes of wild amphibians from the eastern edge of the Tibetan plateau. Microbiology 9:e1004. doi: 10.1002/mbo3.1004
Xu, H. J., Li, S., Su, J. Q., Nie, S., Gibson, V., Li, H., et al. (2014). Does urbanization shape bacterial community composition in urban park soils? A case study in 16 representative Chinese cities based on the pyrosequencing method. FEMS Microbiol. Ecol. 87, 182–192. doi: 10.1111/1574-6941.12215
Xu, L., Xiang, M., Zhu, W., Zhang, M., Chen, H., Huang, J., et al. (2020a). The behavior of amphibians shapes their symbiotic microbiomes. mSystems 5:e00626-20. doi: 10.1128/mSystems.00626-20
Xu, L., Zhou, J., Zheng, P., Wan, X., Zhu, L., Huang, J., et al. (2022). Seasonal variation significantly influenced the stochasticity of community assembly of amphibian symbiotic bacteria. Environ. Microbiol. 24, 5734–5748. doi: 10.1111/1462-2920.16224
Yang, Y., Chen, H., Al, M. A., Ndayishimiye, J. C., Yang, J. R., Isabwe, A., et al. (2022). Urbanization reduces resource use efficiency of phytoplankton community by altering the environment and decreasing biodiversity. J. Environ. Sci. 112, 140–151. doi: 10.1016/j.jes.2021.05.001
Yang, L., Zhao, S., and Liu, S. (2022). A global analysis of urbanization effects on amphibian richness: patterns and drivers. Glob. Environ. Chang. 73:102476. doi: 10.1016/j.gloenvcha.2022.102476
Yao, Q., Yang, H., Wang, X., and Wang, H. (2019). Effects of hexavalent chromium on intestinal histology and microbiota in Bufo gargarizans tadpoles. Chemosphere 216, 313–323. doi: 10.1016/j.chemosphere.2018.10.147
Yuan, M. M., Guo, X., Wu, L., Zhang, Y., Xiao, N., Ning, D., et al. (2021). Climate warming enhances microbial network complexity and stability. Nat. Clim. Chang. 11, 343–348. doi: 10.1038/s41558-021-00989-9
Zhai, T., Wang, J., Fang, Y., Qin, Y., Huang, L., and Chen, Y. (2020). Assessing ecological risks caused by human activities in rapid urbanization coastal areas: towards an integrated approach to determining key areas of terrestrial-oceanic ecosystems preservation and restoration. Sci. Total Environ. 708:135153. doi: 10.1016/j.scitotenv.2019.135153
Zhang, Q. Q., Liu, Z. H., Liu, L. L., Hu, G., Lei, G. L., Wang, Y., et al. (2020). Prebiotic maltose gel can promote the vaginal microbiota from BV-related bacteria dominant to lactobacillus in rhesus macaque. Front. Microbiol. 11:594065. doi: 10.3389/fmicb.2020.594065
Zhang, L., Yang, L., Zohner, C. M., Crowther, T. W., Li, M., Shen, F., et al. (2022). Direct and indirect impacts of urbanization on vegetation growth across the world’s cities. Sci. Adv. 8:eabo0095. doi: 10.1126/sciadv.abo0095
Keywords: biotic homogenization, disturbance, abiotic gradient, community patterns, ecological niches
Citation: Zhou J, Liao Z, Liu Z, Guo X, Zhang W and Chen Y (2023) Urbanization increases stochasticity and reduces the ecological stability of microbial communities in amphibian hosts. Front. Microbiol. 13:1108662. doi: 10.3389/fmicb.2022.1108662
Edited by:
Lifeng Zhu, Nanjing University of Chinese Medicine, ChinaCopyright © 2023 Zhou, Liao, Liu, Guo, Zhang and Chen. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Youhua Chen, 
Y2hlbnloQGNpYi5hYy5jbg==
†These authors have contributed equally to this work