Associations between human impacts and forest soil microbial communities

We are living in a new epoch—the Anthropocene, in which human activity is reshaping global biodiversity at an unprecedented rate. Increasing efforts are being made toward a better understanding of the associations between human activity and the geographic patterns in plant and animal communities. However, similar efforts are rarely applied to microbial communities. Here, we collected 472 forest soil samples across eastern China, and the bacterial and fungal communities in those samples were determined by high-throughput sequencing of 16S rRNA gene and internal transcribed spacer region, respectively. By compiling human impact variables as well as climate and soil variables, our goal was to elucidate the association between microbial richness and human activity when climate and soil variables are taken into account. We found that soil microbial richness was associated with human activity. Specifically, human population density was positively associated with the richness of bacteria, nitrifying bacteria and fungal plant pathogens, but it was negatively associated with the richness of cellulolytic bacteria and ectomycorrhizal fungi. Together, these results suggest that the associations between geographic variations of soil microbial richness and human activity still persist when climate and soil variables are taken into account and that these associations vary among different microbial taxonomic and functional groups.


Introduction
Ecosystems are being rapidly degraded due to increasing human population size, economic growth, urbanization, and pollution Vitousek et al., 1995). These unprecedented human activities are causing exceptional rates of modification in atmospheric composition, climate, soil fertility, and land use (Bai et al., 2017) and progressively reshaping the spatial distribution and community structure of various organisms (Gossner et al., 2016) as well as their associated ecosystem services (Mitchell et al., 2015). Therefore, understanding the associations between human impacts and the distribution of biomes has attracted considerable attention and will provide important insights for conservation planning (Cincotta et al., 2000).
There is an increasing body of literature supporting the idea that geographic patterns of species assemblages are associated with human impact variables such as human population density (Luck, 2007), per capita gross domestic product (Marques et al., 2019), and pollution (Barker and Tingey, 2012). Human population density is widely identified as a surrogate measure of human impacts because human activities such as human-induced biological invasions and land-use changes tend to be more intense in more densely populated regions (Laurance et al., 2002;McKee et al., 2004). Previous studies have shown that there is a negative relationship between human population density and species richness of birds (Koh et al., 2006). However, there is also some evidence supporting that human population density is positively related to species richness of amphibians and mammals (Luck et al., 2004). Given that geographic pattern of species richness is frequently found to be environmentally dependent (Stein et al., 2014) and human activity is widely identified as a major cause of current environmental changes (Bai et al., 2017;Fang et al., 2018;Song et al., 2018), it is expected that human activity shapes the geographic variation of species richness through its associations with environmental factors (D'agata et al., 2014;Torres-Romero and Olalla-Tárraga, 2015).
Soil microbes represent most of the biodiversity on Earth (Whitman et al., 1998) and play pivotal roles in ecosystem functioning (Falkowski et al., 2008), as well as in maintaining the aboveground plant diversity (Bardgett and Van der Putten, 2014) and productivity (Van der Heijden et al., 2008). Thanks to the development of high-throughput sequencing methods, the last decade has witnessed extensive efforts in exploring the key environmental factors structuring the geographic variations of soil microbial communities (Fierer and Jackson, 2006;Tedersoo et al., 2014). However, the effect of human demography and its concomitant associations with the environment have been largely ignored. The environmental pollution associated with increasing human activity can change soil properties and lead to degradation (Bai et al., 2017), and this might alter microbial communities. Indeed, microbial studies in human-modified landscapes have provided evidence that human activity has similar impacts on soil microbial communities (Rodrigues et al., 2013;Gossner et al., 2016;Mueller et al., 2016;Reese et al., 2016;Schmidt et al., 2017). For example, it has been shown that urbanization had no impact on the richness of soil bacteria and fungi (Reese et al., 2016) but dramatically reduced the richness of ectomycorrhizal fungi (Schmidt et al., 2017). In general, soil microbial community composition has been found to be more similar across space in human-modified than undisturbed landscapes (Rodrigues et al., 2013;Schmidt et al., 2017). These studies and others provide clues on how, as with macroorganisms (Torres-Romero and Olalla-Tárraga, 2015), the geographic variations of microbial communities are associated with human activity, and more importantly, different soil microbial taxonomic and functional groups might be linked with human impacts idiosyncratically.
Here, we collected 472 forest soil samples across eastern China (Wang et al., 2019), and the bacterial and fungal communities in those samples were determined by highthroughput sequencing of the 16S rRNA gene and internal transcribed spacer (ITS) region, respectively. We compiled human impact variables including demography, urbanization, economic growth, and pollutant emissions and variables related to temperature and precipitation to represent climate conditions as well as variables related to physicochemical properties (e.g., pH and moisture) to represent soil conditions. After statistically taking climate and soil variables into account, we aimed to explore how human activity was associated with the geographic variations of bacterial and fungal richness and how these associations varied among different microbial functional groups.

Study sites
The detailed information about soil sampling and soil property measurement, as well as the basic description of fungal communities, has been previously described in Wang et al. (2019). Briefly, we collected 472 soil samples from forests scattered across eastern China (Wang et al., 2019). These forests span a wide range of vegetation types including tropical forests, temperate mixed coniferous broad-leaved forests, temperate deciduous broad-leaved forests, subtropical forests, and cold temperate coniferous forests (Wang et al., 2019). All soil samples were kept on ice for transport to the laboratory. For each soil sample, a subsample was stored at À20 C for DNA extraction, and the remaining soils were used for soil physical and chemical analyses. High-throughput sequencing of the ITS region was used to profile the soil fungal community. In this study, we included the fungal data set from a previous study (Wang et al., 2019) and further characterized the bacterial communities.

Molecular and bioinformatic analyses of bacterial communities
Total genomic DNA was extracted using the MoBio Power-Soil DNA Isolation Kit following the manufacturer's instructions. We amplified the V4 hypervariable region of the 16S rRNA gene using the 515-F (GTGCCAGCM GCCGCGGTAA) and 806-R (GGACTACHVGGGTWTCTAAT) primer pair (Caporaso et al., 2012). The reverse primers included an error-correcting 12-bp barcode unique to each sample to permit demultiplexing. PCR products were pooled in equimolar concentrations, purified by E.Z.N.A.1 Gel Extraction Kit (Omega BioTek, Doraville), and sequenced on a 2 Â 300 bp paired-end Illumina MiSeq platform at Sun Yat-sen University, Guangzhou, China.
Raw sequences were processed with Mothur (Schloss et al., 2009) for quality filtering and merging of pair-end reads. Reads were first trimmed to remove short and low-quality sequences (trim.seqs command: minlength ¼ 300, maxhomop ¼ 12, qwindowsize ¼ 50, qwindowaverage ¼ 20, maxambig ¼ 0, pdiffs ¼ 2, bdiffs ¼ 1). These quality-filtered reads were subsequently merged into paired-end sequences (make.contigs command with default parameters). Assembled contigs without an exact match to one of the barcode sets or primers were discarded. The high-quality sequences were clustered into operational taxonomic units (OTUs) based on a 97% similarity threshold with the UPARSE algorithm (Edgar, 2013), using the "-cluster_otus" command in USEARCH, with chimera sequences identified and removed during the procedure. After removing singletons, OTU taxonomy was assigned using the Ribosomal Database Project classifier with a confidence threshold of 0.8 (Wang et al., 2007) against the SILVA database (release 132; Quast et al., 2013). OTUs classified as chloroplasts and mitochondria, and those unassigned at the domain level were removed.
Although the 515-F/806-R primer set also captures Archaea, the diversity and abundance of Archaea were found to be very low in this study. Specifically, only 171 OTUs were assigned to Archaea accounting for only 0.42% of the total number of 16S rRNA sequences, and thus they were removed from further analyses. A total of 23,714,066 sequences and 27,888 OTUs were retained. Proteobacteria (35.19%), Acidobacteria (23.48%), and Actinobacteria (12.45%) were the three most abundant bacterial phyla ( Figure S1). To account for unequal sequencing depth, samples were rarefied to 9,736 sequences.

Assignment of microbial functional groups
To explore the potential functional consequence of human activity, we assigned bacterial and fungal OTUs Art. XX, page 2 of 10 Chen et al: Associations between human impacts and forest soil microbial communities to functional groups, respectively. Functional annotation of bacterial OTUs was performed using the Functional Annotation of Prokaryotic Taxa, which uses a manually curated database of cultured microorganisms to extrapolate microbial functional profiles from taxonomic profiles (Louca et al., 2016). We focused on two bacterial functional groups (i.e., nitrifiers and cellulolytic bacteria) because they play crucial roles in nitrogen (Schimel and Bennett, 2004) and carbon cycling (Leschine, 1995), respectively. Fungal OTUs were assigned to three fungal functional guilds using FUNGuild (Nguyen et al., 2016). OTUs that were assigned to only one specific guild (i.e., saprotrophic fungi, ectomycorrhizal fungi, or plant pathogens) with "highly probable" or "probable" confidence ranking were retained (Wang et al., 2019).

Collection of explanatory variables
Explanatory variables in this study were grouped into four categories: human activity, climate, soil, and space (Table  S1). Human activity is multidimensional, and it includes demography, urbanization, economic growth, and pollution, among others (McGill et al., 2015). Based on the city or municipality a specific sample is located in, we compiled human population density, proportion of urban population, per capita gross domestic product, sulfur dioxide, nitric oxide, and smoke contents of the year 2012 from the China Statistical Year book (http://www. stats.gov.cn/tjsj/ndsj/2012/indexeh.htm; Yan et al., 2016). For each sample, soil properties including moisture, pH, electrical conductivity, total organic carbon, total nitrogen, total phosphorus, total potassium, available phosphorus, available potassium, available aluminum, available calcium, and available magnesium were measured as described in Wang et al. (2019). Moreover, we compiled 19 bioclimate variables from WorldClim (http://www.worldclim.org/bioclim) to describe the mean annual trends and seasonality of temperature and precipitation (Hijmans et al., 2005). Latitude and longitude were recorded in situ.

Statistical analyses
Statistical analyses were implemented in R (R Core Team, 2018). To avoid potential multicollinearity problems, we retained explanatory variables with a variation inflation factor below 10 (Torres-Romero and Olalla-Tárraga, 2015). All soil variables, all human impact variables, five climate variables including mean diurnal range of temperature, isothermality, mean temperature of wettest quarter, annual precipitation, and precipitation seasonality were retained. Latitude and longitude were retained to account for spatial autocorrelation (see the same approach in Delgado-Baquerizo et al., 2017). These multicollinearity free variables (x) were standardized to 0 to 1 range (x À x min )/(x max À x min ).
Variation partitioning was performed on microbial richness (i.e., number of rarefied OTUs) using partial linear regression (Legendre, 2008). The total variation in microbial richness was decomposed into fractions explained uniquely and jointly by human activity, climate, soil, and space. Negative values in the variance explained might be found and were interpreted as zeros (Legendre, 2008). To validate our results, spatial factors were also constructed using the Principal Components of Neighborhood Matrix (PCNM) approach (Dray et al., 2006). The PCNM approach decomposes the total spatial variation into a set of PCNM eigenvectors, where the first PCNM eigenvector represents the broadest spatial scale, and each successive eigenvector represents an increasingly finer scale (Dray et al., 2006). These analyses were conducted using the R package vegan.
Structural equation modeling (SEM) was employed to explore the association between human activity and microbial richness. An a priori model was first constructed ( Figure S2). In our SEM, human activity can uniquely impact the microbial richness, and it can be linked with microbial richness through its associations with climate and soil variables ( Figure S2). In addition, climate, soil, and spatial variables can shape microbial richness ( Figure  S2). We used random forest modeling (Breiman, 2001) as described in Delgado-Baquerizo et al. (2017) to select significant predictors of microbial richness. Different categories of explanatory variables were grouped into the same box for graphical simplicity, but they were not grouped into composite variables when constructing the SEM. We saturated the SEM by allowing all variables within a particular box to covary and then removed explanatory variables with the weakest correlation to release a degree of freedom in the SEM (Delgado-Baquerizo et al., 2017). The goodness of fit for SEM was diagnosed using chisquare test. Moran's I test was used to test for spatial autocorrelation in the residuals of the models (Moran, 1950). In addition, partial residual plots were used to show the associations between microbial richness and human population density when accounting for the effects of all the other explanatory variables. P values for multiple testing were corrected using the Benjamini-Hochberg method (Benjamini and Hochberg, 1995). These analyses were conducted using the R packages rfPermute and lavaan.

Associations between microbial richness and human activity
Results of variation partitioning showed that climate, soil, and spatial variables dominantly explained the geographic patterns of microbial richness (Figure 1). The importance of climate and soil variables is consistent with most studies showing that abiotic factors play dominant roles in structuring the geographic variation of microbial richness (Fierer and Jackson, 2006;Tedersoo et al., 2014). The effects of spatial variables can be attributed to the limited dispersal ability of microbes, the physical barriers that limit the dispersal of microbes, or the unmeasured but spatially structured environmental variables. Human activity explained a unique portion of the variations in soil microbial richness, and it also shared a part of explained variation with climate, soil, and space ( Figure 1). We found qualitatively similar results when using PCNM vectors to represent the spatial factors ( Figure S3).
We used SEM to explore the associations between human activity and microbial richness (Figure 2).  . Associations between microbial richness and human population density. Different categories of explanatory variables were grouped in the same box for graphical simplicity. Numbers adjacent to arrows are significant standardized path coefficients at P < 0.05. All human impact variables were used to construct SEM, but only human population density is shown. Only those climate, soil, and spatial variables that could mediate the associations between human population density and microbial richness are shown.The complete results of SEM are provided in Specifically, we focused on human population density because this variable has been commonly identified as a proxy measure of various activities associated with human settlements (Cincotta et al., 2000;McKee et al., 2004) and has been widely used previously in studies of plants and animals (Luck, 2007;Spear et al., 2013;D'agata et al., 2014;Paradis, 2018). We found a positive association between bacterial richness and human population density (r ¼ 0.09, P < 0.01; Figures 2 and 3). A recent study also found that soil bacterial diversity was higher in urban areas densely populated by humans (Wang et al., , 2018, suggesting that this observation in urban landscapes can be extended to natural ecosystems. A possible explanation might be that more non-native taxa are introduced in regions with higher human population density, but they do not outcompete the native taxa, resulting in a net increase of species richness. This has been shown to occur in assemblages of microorganisms (McKinney, 2002) and plants (Stohlgren et al., 2003). For fungi, we did not find a significant association between human population density and richness (r ¼ 0.03, P ¼ 0.60; Figures 2 and 3). These results suggest that the effect of human activity on richness differs for bacteria and fungi in forest soils. In addition, we found that different environmental variables drove the links between human population density and microbial richness for bacteria and fungi ( Figure 2). Soil pH drove the association between human population density and bacterial richness but not fungal richness (Figure 2). The possible explanations for this indirect association between human population density and bacterial richness are the soil pH changes related to the release of leachates from calcareous materials, decomposing organic waste and nitrogen deposition associated with increasing human activity (Jim, 1998;Boateng et al., 2006;Cofie et al., 2009). Moreover, previous studies have shown that there is a strong association between fungal richness and soil pH (Tedersoo et al., 2014;Glassman et al., 2017). This discrepancy between previous studies and our observation suggests the possibility that human activity changes other environmental variables that play important roles in structuring fungal communities.

Associations between richness of microbial functional groups and human activity
We then explored the potential functional consequence of human impacts by examining the associations between human population density and different microbial functional groups. While human population density was not related to the relative abundances of microbial functional groups ( Figure S4), we found that human population density was associated with the richness of certain microbial functional groups (Figures 2 and 3). One of the crucial Partial residual plots represent relationships between microbial richness and human population density when accounting for the effects of all the other explanatory variables. The solid red line represents the best fitted linear regression model when a significant association between human population density and microbial richness was found. DOI: https://doi.org/10.1525/elementa.005.f3 Chen et al: Associations between human impacts and forest soil microbial communities Art. XX, page 5 of 10 steps in nitrogen cycling is nitrification, in which ammonia is oxidized to nitrate (Schimel and Bennett, 2004). In this study, we found that human population density was positively associated with the richness of nitrifying bacteria (r ¼ 0.25, P < 0.001; Figures 2 and 3). A possible explanation is that a range of activities associated with human settlements such as fossil fuel combustion and industrial nitrogen fixation increase the inputs of available nitrogen, resulting in an acceleration of the nitrification process. For cellulolytic bacteria, we found that there was a negative association between human population density and the richness (r ¼ À0.34, P < 0.001; Figures 2  and 3). Previous studies have shown that reduced richness of cellulolytic bacteria is correlated with reduced plant litter decomposition (Knorr et al., 2005;Frey et al., 2014). Therefore, the lower richness of cellulolytic bacteria in more densely populated regions suggests the potential inhibition of plant litter decomposition due to human activity. Cellulose is one of the major organic carbon compounds in plant litter, and its decomposition is key for carbon cycling in soils (Leschine, 1995). Thus, the decline in richness of cellulolytic bacteria might influence the quality of soil organic matter and subsequently plant growth (Lladó et al., 2017). Yet, experimental evidence is needed to validate the causal relationship between the richness of cellulolytic bacteria and litter decomposition. Moreover, we found that both mean diurnal range of temperature and soil pH drove the associations between human population density and the richness of nitrifying bacteria and cellulolytic bacteria (Figure 2), implying that human activity is associated with these microbes through its linkages with climate and soil conditions. Overall, these results raise the potential concern that increasing human activity and its associated environmental changes may perturb soil nutrient cycling in forests by shaping the richness of key bacterial functional groups. We found that the richness of ectomycorrhizal fungi was negatively associated with human population density (r ¼ À0.11, P < 0.001; Figures 2 and 3), which is consistent with previous results from urban landscapes (Schmidt et al., 2017;Soudzilovskaia et al., 2019). A possible reason is that the richness of ectomycorrhizal fungi declines when ectomycorrhizal trees are negatively affected by human activity. A prominent avenue for future work is to explore what types of ectomycorrhizal trees and their associated ectomycorrhizal fungi are vulnerable to human impacts, ultimately leading to a better management strategy for forest ecosystems. In contrast, fungal plant pathogen richness was positively associated with human population density (r ¼ 0.12, P < 0.001; Figures 2 and  3), possibly because crop plants and crop trees are used more frequently in densely populated regions. Mean diurnal range of temperature and isothermality indirectly drove associations between human population density and the richness of ectomycorrhizal fungi and fungal plant pathogens in opposite directions (Figure 2). The contrasting associations between human population density and ectomycorrhizal fungi and fungal plant pathogen richness found in this study might represent the negative interaction between them. That is, it has long been recognized that ectomycorrhizal fungi can act as antibiotic agents to protect plant roots from infection of soil-borne plant pathogens (Marx, 1972;Perrin, 1990;Fisher et al., 2012;Laliberté et al., 2015). For saprotrophic fungi, we did not find a significant association between human population density and richness (r ¼ À0.01, P ¼ 0.94; Figures 2  and 3), but soil moisture indirectly mediated their association ( Figure 2).

Conclusions
Human activity was associated with soil microbial richness in forest ecosystems across eastern China. Soil microbial communities in more densely populated regions were characterized by higher richness of total bacteria, nitrifying bacteria, and fungal plant pathogens but lower richness of cellulolytic bacteria and ectomycorrhizal fungi. The different responses of ectomycorrhizal and pathogenic fungi to human activity might increase future stress on aboveground plants and ecosystem processes. Collectively, our results suggest that the associations between human activity and microbial richness vary among different microbial taxonomic and functional groups. These results call for a multi-taxa approach to the conservation of belowground communities.

Supplemental files
The supplemental files for this article can be found as follows: Text S1. Associations between human impacts and forest soil microbial communities. PDF.

Acknowledgments
We are thankful to members of the Barberán lab for constructive comments. We thank Marc William Cadotte for feedback on an early version of the manuscript. We thank two anonymous reviewers for providing constructive comments on the manuscript.

Funding information
This work was supported by National Natural Science Foundation of China (31800422 and 31600403).

Competing interests
The authors have no competing interests to declare.