Soil Mesofauna Community Changes in Response to the Environmental Gradients of Urbanization in Guangzhou City

There has been a recent increase in interest on how urbanization affects soil fauna communities. However, previous studies primarily focused on some limited land use types or line transects of urban-rural gradients. At family and higher taxonomic levels, we investigated the changes of soil mesofauna communities (abundance, species richness, and community structure) with urbanization intensity along different disturbance features in 47 sites evenly located in downtown Guangzhou and adjacent regions. The 47 research sites were classified into four ecosystem types mainly according to the location (rural/urban), vegetation cover, and management intensity. In turn, the four types with increasing urbanization intensity were rural forest, urban forest, urban woodland, and urban park. Firstly, the role of urban soil property (soil physicochemical characteristic and soil heavy metal content) in regulating soil mesofauna community was investigated. The results showed that soil mesofauna abundance and diversity decreased with increasing soil pH, total nitrogen content (TN), and heavy metal comprehensive index (CPI). Soil Pb decreased soil mesofauna species richness (taxa number) and regulated soil mesofauna community structure. Secondly, we examined the effects of landscape changes on the soil mesofauna community. We found impervious surface (IS) ratio did not predict changes in soil mesofauna abundance, species richness, or community structure. Instead, IS ratio was positively correlated with soil pH, soil TN, and CPI. After excluding sites that belonged to rural forests and urban parks, site area was positively correlated with soil mesofauna abundance. Thirdly, our results revealed significant differences in soil property, landscape trait, and soil mesofauna community among the four ecosystem types. Interestingly, urban forest, the one lightly disturbed by urbanization, but not rural forest, had the highest soil mesofauna abundance. Soil mesofauna abundance in urban woodlands was similar to that in urban parks, which was about half of that in urban forests. Species richness in urban parks was 21% lower than that in rural forests. Our results also showed that urban woodland and urban parks had distinct mesofauna community structures compared to those in rural forests and urban forests. In conclusion, the present study suggested that (1) soil property changes due to urbanization, such as increased pH and heavy metal enrichment in urban soil, decreased soil mesofauna abundance and species richness, changed community structure, and mediated the effect of landscape change on soil mesofauna community; (2) however, soil and landscape changes could not explain the increase of abundance in urban forests, which supported the intermediate disturbance hypothesis.

There has been a recent increase in interest on how urbanization affects soil fauna communities. However, previous studies primarily focused on some limited land use types or line transects of urban-rural gradients. At family and higher taxonomic levels, we investigated the changes of soil mesofauna communities (abundance, species richness, and community structure) with urbanization intensity along different disturbance features in 47 sites evenly located in downtown Guangzhou and adjacent regions. The 47 research sites were classified into four ecosystem types mainly according to the location (rural/urban), vegetation cover, and management intensity. In turn, the four types with increasing urbanization intensity were rural forest, urban forest, urban woodland, and urban park. Firstly, the role of urban soil property (soil physicochemical characteristic and soil heavy metal content) in regulating soil mesofauna community was investigated. The results showed that soil mesofauna abundance and diversity decreased with increasing soil pH, total nitrogen content (TN), and heavy metal comprehensive index (CPI). Soil Pb decreased soil mesofauna species richness (taxa number) and regulated soil mesofauna community structure. Secondly, we examined the effects of landscape changes on the soil mesofauna community. We found impervious surface (IS) ratio did not predict changes in soil mesofauna abundance, species richness, or community structure. Instead, IS ratio was positively correlated with soil pH, soil TN, and CPI. After excluding sites that belonged to rural forests and urban parks, site area was positively correlated with soil mesofauna abundance. Thirdly, our results revealed significant differences in soil property, landscape trait, and soil mesofauna community among the four ecosystem types. Interestingly, urban forest, the one lightly disturbed by urbanization, but not rural forest, had the highest soil mesofauna abundance. Soil mesofauna abundance in urban woodlands was similar to that in urban parks, which was about half of that in urban forests. Species richness in urban parks was 21% lower than that in rural forests. Our results also showed that urban woodland and urban parks had distinct mesofauna community structures compared to those in rural forests and urban forests.
In conclusion, the present study suggested that (1) soil property changes due to urbanization, such as increased pH and heavy metal enrichment in urban soil, decreased

INTRODUCTION
Urbanization has a tremendous effect on the economy and environment in the modern world. Habit loss and fragmentation, land use change, environmental pollution, and anthropogenic interruption during the urbanization process has caused biodiversity declines in mammals, birds, and amphibians (McKinney, 2006(McKinney, , 2008Diego Ibanez-Alamo et al., 2017). Soil mesofauna, such as collembola and mite, can exert substantial effects on soil fundamental functions (Wolters, 2001;Bardgett and van der Putten, 2014). As a result, it is important to know how a soil mesofauna community would respond to the environmental changes induced by urbanization.
Urbanization always proceeds with enormous land use type change, which may be the strongest factor influencing soil fauna communities. Many previous studies investigated the urbanization effects on soil fauna by comparing communities among urban ecosystems of different land use types. In Phoenix, the soil arthropod community is significantly different among residential, industrial, agricultural, and desert remnants, especially when springtails, ants, and mites are excluded from the analysis (McIntyre et al., 2001). In a semi-arid city, Ge et al. (2019) showed that urban lawn soils have lower soil invertebrate evenness and diversity, but higher density than remnant patches of native coastal sage scrub. Joimel et al. (2018) investigated collembola communities in both extensive and productive green roofs, and found no differences in collembola taxonomic structures. However, it is not easy to know exactly which environmental factors shape a soil fauna community from such studies, because most environmental conditions vary dramatically and inconsistently among land use types.
More and more studies have begun to investigate the soil fauna response to a specific agent of urban environmental change. Soil in a city is sometime called "man-made soil, " because it is someway re-created during city construction and residential daily lives (Gilbert, 1989). During the re-creation, man-made materials, such as broken bricks, glass and china, kitchen refuse, fertilizer, and so on, are mixed into urban soils, which result in a widely reported elevated soil pH (Jim, 1998;Pouyat et al., 2015;Asabere et al., 2018), increased soil nutrients (Pouyat et al., 2002;Trammell et al., 2020), and heavy metal accumulation (Wei and Yang, 2010). In urban green spaces under management, removal of aboveground litter, application of fertilizer and irrigation, and so on may also substantially change soil characteristics (Norton, 2011;Mao et al., 2014;Tresch et al., 2018Tresch et al., , 2019. Many efforts have been made to reveal the relationship between heavy metal and soil fauna. It is reasonable that soil fauna will be inhibited by heavy metal, which has been widely reported in ecotoxicological tests in lab incubations (Crommentuijn et al., 1993;Didden and Rombke, 2001;Herbert et al., 2004). However, there is no consensus in field studies. While some field studies found heavy metal in urban soils decrease soil fauna abundance and richness (Fiera, 2009;Santorufo et al., 2012), other studies showed contrary results (Fountain and Hopkin, 2004;Joimel et al., 2017;Milano et al., 2018;Sterzynska et al., 2018). In contrast, there are only a few studies exploring the relationships between soil fauna communities and other soil characteristics in urban ecosystems. In the city of Warsaw, Rzeszowski et al. (2017) found that factors contributing to significant collembola communities' variation are soil potassium and phosphorus concentrations and soil pH, whereas soil nitrogen and carbon do not have significant effects. In the Beijing Olympic Park, Song et al. (2015) also showed that there are significant correlations between soil physicochemical property and the soil mesofauna community. However, due to the limited number of studies, it is not clear that whether such relationships widely exist in different urban soils.
Most urban ecosystems are embedded within human construction and are of a small area with a large ratio of girth to area (Bradshaw, 2003;Hruska, 2006), which is typical of fragmented habitats. Such landscape characteristics could affect individual and gene exchanges, thus shaping the process of community construction (LaPoint et al., 2015;Lepczyk et al., 2017). Studies on aboveground organisms suggested that habitat area, connectivity, and landscape diversity are important in determining organism community composition and structure (Faeth et al., 2011;Braaker et al., 2014, Burrow, 2018. In contrast, studies on urban soil fauna seldom take landscape features into account (but see Bolger et al., 2000;Milano et al., 2018;Xie et al., 2018). A recent systematic review suggests that the underlying mechanisms regulating soil organisms may be inconsistent with those generated from studies on aboveground organisms (Thakur et al., 2019). As a result, it remains unclear whether landscape features could have significant effects on soil fauna community.
Some previous studies found species within a family or even a genus could have different capacities to resist environmental stress as a result of their different functions and physical traits (Magura et al., 2010;Santorufo et al., 2014). However, taxonomic analyses at a species level in soil fauna is expensive, both in terms of equipment and time. In recent years, many studies found it is realizable to indicate soil environmental changes at family or higher taxonomic levels (Parisi et al., 2005;Santorufo et al., 2012;Menta et al., 2018). These studies hypothesize that soil fauna within a higher taxonomic (family or order, etc.) level have similar capabilities to adapt to soil environmental changes, which was supported by many studies, especially in natural grassland (Madej and Kozub, 2014), farmland (Rudisser et al., 2015), and degraded soils Madej and Kozub (2014).
In the present study, we made soil and soil mesofauna samples in 47 sites evenly located in the downtown of Guangzhou city and adjacent regions, representing urban soil ecosystems of different soil physicochemical properties, soil heavy metal pollution degree, and landscape trait. Our main objective was to advance the knowledge on the effect of specific urbanization agents on soil mesofauna abundance, diversity, and community structure at family and higher taxonomic levels.

Study Area and Sample Design
The city of Guangzhou (22 • 26 -23 • 56 N and 112 • 57 -114 • 03 E) is the capital of Guangdong province in South China. The soil type in natural ecosystems in Guangzhou is latosolic red soil. The climate in Guangzhou is a subtropical, marine monsoon climate with an annual average temperature of 21.5-22.2 • C and an annual average precipitation of 1,623.6-1,899.8 mm.
Guangzhou has an annual average number of 149.2 rainy days, with most rainy days occurring during April to October. Typically, July is the hottest month with a mean temperature of 28.6 • C, and January is the coldest with a mean temperature of 13.6 • C.
To establish representative sites distributed evenly in the downtown of Guangzhou city and adjacent regions, we made preliminary selections of urban ecosystems in online satellite images along (but not closed to) the arterial roads across the study area. At first, 55 sites were chosen on the satellite images. After that, we went to each potential site to discriminate whether the site met our standards. In this procedure, sites that could not be reached or with newly translocated soils or newly planted plants were excluded. At last, 47 sites were established in the study area (Figure 1).
The 47 sites were classified into four ecosystem types according to the location, vegetation cover, amount of anthropogenic solid waste, dominant plant type, and management intensity ( Table 1). The four typical ecosystem types with increasing urban disturbance intensity are rural forest (8 sites), urban forest (16 sites), urban woodland (13 sites), and urban park (10 sites). A rural forest site was always located in a consecutive forest and far from the city's infrastructure (such as roads and building), while an urban forest site was surrounded by human construction and isolated from other soil ecosystems. However, the urban forests were usually located on hills, thus only suffering mild human disturbances. Both the rural and urban forests were broadleaf evergreen forests and had high vegetation covers (>80%). The dominant plants in the rural forest and urban forest were native or/and alien trees. An urban woodland site was usually small and close to city construction. As a result, the soils in urban woodland were heavily polluted by solid anthropic wastes (such as plastic, wire, paper, brick, clothes, and so on), indicating a strong human disturbance. The dominant plants in urban woodlands were street trees or/and garden shrubs. There were only sparse and young plants in the urban woodlands, and thus with lower vegetation cover. The dominant plants in urban parks were similar to those in the urban woodlands. However, urban parks were managed ecosystems. The vegetation in the urban parks were pruned, watered and fertilized, and the aboveground plant litter was removed for sanitization. The soils may be turned over, and solid wastes polluted.

Soil Mesofauna Collection and Identification
From October 2018 to December 2018, we sampled soil mesofauna in the study area once for each site. At each site, a transect including three 5 × 5 m 2 sub-plots separated by 10 m was established. In each sub-plot, three soil cores were randomly taken from surface soil (0-10 m) with a steel corer (5 cm inner diameter) to generate a soil sample for soil mesofauna extraction. Therefore, we collected three samples per site. Immediately after collection, the samples were transported to our laboratory, and soil collembola and mites were extracted using Tullgren dry funnels for 48 h. All specimens were sorted and counted and examined with an Olympus BX41 research microscope (Olympus, Tokyo, Japan). The mites were allocated to three groups, namely Mesostigmata, Oribatida, and Prostigmata, while the collembola were grouped according to family and others to order or class, mainly according to Yin (1998). At each site, soil mesofauna abundance (number of soil mesofauna per m 2 ) and species richness (total number of taxa, S) were recorded and used to calculate Shannon diversity index as (Shannon, 1948): where P i is the ratio between the abundance of species i and the total soil mesofauna abundance. Evenness (J) was evaluated according to the Pielou index (Pielou, 1969):

Soil Properties Analysis
At each transect, another three soil samples were taken using a steel corer (3 cm inner diameter) to generate a soil sample to test for soil physicochemical properties and heavy metal levels. Apart from five polluted samples, we assessed soil pH, soil organic matter (SOM), and soil total nitrogen (TN) for each site. Soil pH was measured using a 1: 2.5 soil: water suspension with the potentiometric method. SOM was determined using H 2 SO 4 -K 2 Cr 2 O 7 oxidation method. Soil TN was quantified by the Kjeldahl acid digestion method. Heavy metal concentration (Zn, Cu, Cd, and Pb) was analyzed after digestion in a mixture of nitric, perchloric acid, and hydrogen peroxide. Soil heavy metal comprehensive index (CPI) was calculated as (Li et al., 2008): where M i is the pollution index of heavy metal i; C i (mg kg −1 ) is the measured heavy metal i; S i (mg kg −1 ) is the environmental background value of Guangzhou City (Zhou et al., 2009); and CPI is the comprehensive pollution index.

Assessments of Impervious Surface and Site Area
Land cover information was generated from a set of global 10-m resolution land cover map from Sentinel 2 images via the Google Earth Engine. Buildings and pavements were classified into impervious surface (IS). Circular buffers at radiuses of 200, 500, and 1000 m were established for each site with the site as circle center. After that, IS ratios of each circle buffer were calculated. The IS ratio at a radius of 500 m was chosen for the next analyses because it had higher correlation coefficients with soil property and soil mesofauna parameter in Pearson's correlation analysis (Supplementary Table 1). In addition, to show the overall view of the IS ratio of the study area, the original map of land cover was transformed to a map of geographic grid cells at 500 m. The IS ratio in each grid cell was calculated, respectively, and visualized in Figure 1. The area of each site was measured in the LocalSpaceViewer 4.0 using Google Earth Map. The area reported in this study did not include those of consecutive rural forest sites, because they were part of consecutive forests. In addition, the areas of urban park sites were not reported in the present study, because internal impervious pavements in urban parks could not be distinguished in the map.

Statistical Analysis
For all analyses, values from three soil samples (sub-plots) in a site were averaged and treated as a single replicate per site (N = 47 for soil mesofauna parameters, N = 42 for soil properties). Data were square-rooted or log transformed to fulfill the hypothesis of normality and homogeneity. We plotted scatter plots and conducted Pearson's correlation to examine the relationship among landscape trait (IS ratio and site area), soil property (pH, SOM, TN and soil total Zn, Cu, Cd, Pb, CPI), and soil mesofauna community parameter (abundance, richness, Shannon's index, and Pielou's index). After that, according to the result of Pearson's correlation, linear regression was used to determine the variation of each soil mesofauna parameter along each soil property and landscape trait. We used oneway ANOVA to examine the differences in landscape trait, soil property, and soil mesofauna parameter among the four ecosystem types with different urbanization intensities. Least significance difference (LSD) was used after the ANOVA to test the differences between any two ecosystem types. Pearson's correlation, one-way ANOVA, and regression analysis were conducted in SPSS 19.0 (SPSS, Chicago, IL, United States).
To visualize the soil mesofauna community pattern, we performed ordination on community composition with nonmetric multidimensional scaling (NMDS). Chao method was chosen because of the unbalanced replicates among ecosystem types (Anderson and Millar, 2004). To this end, environmental variables were fitted on this NMDS space. A permutational multivariate analysis of variance (PERMANOVA) was conducted to test the significance of each environmental variable. We also used PERMANOVA to quantitatively examine the difference in soil mesofauna community among ecosystem types using a Chao dissimilarity matrix. The PERMANOVA and NMDS were conducted using Vegan in R 3.6.1. Significance level was set at α = 0.050.

Effects of Urban Landscape Changes on Soil Mesofauna Abundance, Diversity, and Soil Property
Impervious surface ratio was not significantly correlated with any soil mesofauna parameters (Supplementary Table 3). In contrast, both soil pH (R 2 = 0.255, p < 0.001) and soil total nitrogen (R 2 = 0.180, p = 0.005) were positively correlated with IS ratio (Figure 3). Heavy metal CPI also increased with IS ratio (R 2 = 0.119, p = 0.026). In addition, as the Pearson's correlation analysis showed, soil total Zn, Cu, and Cd were significantly positively correlated with IS ratio (Supplementary Table 3).
After excluding the sites belonging to rural forests and urban parks, soil mesofauna abundance was positively correlated with site area (R 2 = 0.259, p = 0.002) (Figure 4). Both soil pH (R 2 = 0.288, p = 0.005) and heavy metal CPI (R 2 = 0.358, p = 0.002) were negatively correlated with site area (Figure 4). In addition, soil mite and collembola abundance significantly increased with site area, while total soil Zn, Cu, and Cd were significantly negatively correlated with site area (Supplementary Table 3). There were no significant correlations between soil mesofauna diversity and site area (Supplementary Table 3).

Soil Property, Landscape Trait, and Soil
Mesofauna Abundance and Diversity in Rural Forest, Urban Forest, Urban Woodland, and Urban Park All measured landscape and soil property variables varied significantly among the four ecosystem types (all p < 0.031), except the SOM and total soil Pb. Impervious surface ratio, site area, soil pH, TN, total soil Zn, Cu, Cd, and CPI increased from rural forest, urban forest, urban woodland, to urban park ( Table 2).
Soil mesofauna abundance in urban forests was the largest in the four ecosystem types (Figure 5). Rural forest had the second largest soil mesofauna abundance. Soil mesofauna abundance in urban woodland were similar to those in urban parks, which were only about half of that in urban forests. Specifically, total soil mesofauna abundance in the four ecosystem types were 76, 076 ± 10,848, 98,156 ± 9,299, 51,341 ± 5,748, and 50,768 ± 4,985 ind. m −2 , respectively. The total soil mesofauna abundance in urban forests was significantly higher than those in urban woodland (p < 0.001) and urban parks (p < 0.001), and 20% higher than that in rural forests (p = 0.080).

Effects of Urbanization on Soil Mesofauna Community Structure
Differences in mesofauna community structure among sampling sites were visualized with a non-metric multidimensional scaling (NMDS) (Figure 6). In addition, soil mesofauna community structure among the four ecosystem types were significantly different (p = 0.036) ( Table 3). Except for the difference between rural forest and urban forest, differences between any two ecosystem types were significant (all p ≤ 0.033). The results of permutational multivariate analysis of variance showed that only soil total Pb concentration was significant in explaining the variation of soil mesofauna community structure (R 2 = 0.174, p = 0.036) (Supplementary Table 4).

Soil Property Regulated Soil Mesofauna Abundance and Diversity in Urban Soils
Most soil pH values in the study area were higher than those of natural forests in this region (Supplementary Table 2; Liu et al., 2010). It is well known that soil alkalization can induce destruction in soil physical structure, which may decrease soil mesofauna abundance and diversity. Less clear is why soil mesofauna abundance and diversity decreased with increased soil pH at a lower range in urban soils as in the present study. In natural grassland, decrease in soil pH from 7.57 to 4.76 depresses soil nematode abundance by decreasing soil base cation availability and soil microbial biomass (Chen et al., 2015). In a study on mineral-poor dry heathlands, an increase in soil pH from 3.50 to 4.13 under liming treatments alleviates P limitation, thus increasing micro-arthropod abundance (mite and collembola) (Siepel et al., 2018). In the city of Warsaw, an increase in soil pH associated with lower concentrations of potassium and phosphorus depressed some collembola genus, thus changing soil collembola community structure (Rzeszowski et al., 2017). These studies may indicate that, besides physical structure, the effects of changed soil pH on an urban soil mesofauna community could be mediated by soil nutrient availability.  Compared to geogenic heavy metal concentrations in the city of Guangzhou (Zhou et al., 2009), most of the sampling sites had higher heavy metal contents (Supplementary Table 2). Toxicological evidence has been widely reported showing that high soil heavy metal concentration can induce soil fauna death in laboratory conditions, especially for collembola and enchytraeid (Crommentuijn et al., 1993;Didden and Rombke, 2001;Herbert et al., 2004). In urban field studies, there is still much controversy on whether accumulated heavy metal in soils could significantly inhibit soil mesofauna or change the community structure (Fountain and Hopkin, 2004;Fiera, 2009;Santorufo et al., 2012;Sterzynska et al., 2018). In ecosystems other than urban ecosystems, most of the controversy could result from the taxon-specific resistance to heavy metal toxicity. For example, in a grassland contaminated by Zn, while abundance of most soil fauna decreased, abundance of coleoptera and arachnida increased (Nahmani and Lavelle, 2002). In a shooting range with Pb and Sb 20 times higher than the surrounding soils, while total abundance and richness of soil mesofauna does not significantly change in the polluted soils, diplura and protura seem to benefit from heavy metal pollution, as they only occurred in soils of the shooting range (Migliorini et al., 2004). In the present field study, at family and higher taxonomic levels, our results suggested that heavy metals that accumulated in urban soils may be an important factor regulating soil mesofauna abundance, richness, diversity, and community structure. Moreover, our result indicated that comprehensive heavy metal pollution degree was more important in regulating soil mesofauna abundance and species richness than any single heavy metal, while soil Pb content seemed to be a key factor influencing soil mesofauna biodiversity FIGURE 4 | Variations of soil mesofauna abundance and soil property with site area. R 2 and p were the fitness and significance of the linear regression models. Note that sites that belonged to rural forest and urban park were excluded in the regressions. The "×" indicates the data point not included in the line regression analysis. The letters indicate significant differences at the level of α = 0.050 among different ecosystem types. IS ratio, impervious surface ratio; SOM, soil organic matter content; TN, soil total nitrogen content; Zn, Cu, Pb, Cd, soil total Zn, Cu, Pb and Cd concentration, respectively; CPI, heavy metal comprehensive index. and community structure. That is, soil mesofauna may respond differently to different soil heavy metal species/heavy metal index. Consistent with our result, Santorufo et al. (2012) found collembola abundance is negatively correlated with soil Pb and Zn concentration, but is not sensitive to Cu. A study on oribatida inhabiting leaf litter suggested that the oribatida community structure is regulated by litter Cd but not Zn, Pb, or Cu (Khalil et al., 2009).

Landscape Change Effects on Soil Mesofauna Community Were Mostly Mediated by Soil Property
Impervious surface (IS) ratio is a direct and important index indicating urbanization intensity (Seress et al., 2014). In stream ecology studies, total impervious area is thought to be an excellent surrogate variable summarizing the effect of variables associated with urban stream syndrome (Walsh et al., 2005;King and Baker, 2010). However, the present study did not find a significant relationship between soil mesofauna community and IS in urban soils, as some previous studies did (Fogaca et al., 2013;Cordonnier et al., 2019). A study conducted in French cities showed that the dependence of ants on IS ratio varies with species identity and spatial scale studied (Fogaca et al., 2013). However, the relationships between soil mesofauna community and IS ratio at 200, 500, and 1000 m radiuses was not significant. Decrease of habitat area is a consequence of habitat fragmentation due to urbanization. In this study, site area was positively correlated with soil mesofauna abundance. Contrasted to this result, soil collembola abundance in dry grassland does not show a significant relationship with site area (Querner et al., 2018). In a systematic literature review, Thakur et al. (2019) found only seven out of 12 studies on soil mesofauna provided support for the theory of island biogeography. Therefore, the question of whether the theory of island biogeography is applicable for soil mesofauna should be further investigated. Nevertheless, this study showed that higher site area could support higher soil mesofauna abundance in urban forests and urban woodland.
On the other hand, we found soil pH, TN, and CPI increased with IS ratio, indicating that IS ratio was a good indicator for these urban soil properties. Similar to IS ratio, site area was negatively correlated with soil pH and CPI. We concluded that habitat fragmentation due to urbanization may depress soil mesofauna abundance, however, to a larger extent, through influencing local soil property. This is in line with previous studies on aboveground arthropods in urban ecosystems (Smith and Schmitz, 2016;Kyro et al., 2018), which suggested local soil property is more important in regulating soil mesofauna communities than landscape trait. FIGURE 6 | Non-metric multidimensional scaling (NMDS) plot of soil mesofauna community structure (Stress = 0.250; method = Chao) and the environmental variables fitted on the NMDS spaces. IS ratio, impervious surface ratio; SOM, soil organic matter content; TN, soil total nitrogen content; Zn, Cu, Pb, Cd, soil total Zn, Cu, Pb, and Cd concentration, respectively; CPI, heavy metal comprehensive index. Circles: rural forests, triangles: urban forests, pluses: urban woodlands,corsses: urban parks.

Soil Mesofauna Community Significantly
Varied Among Rural Forest, Urban Forest, Urban Woodland, and Urban Park Apart from the landscape and soil property changes, urbanization may affect soil mesofauna communities through species  (Bai et al., 2017). In the present study, 47 sites across Guangzhou City were classified into four ecosystem types with increasing urbanization intensity. Located at hills, protected by laws, and far from urban buildings, the rural forest sites underwent the least anthropogenic disturbance. Consistent with the present study, many previous studies found a higher soil mesofauna richness in ecosystems with lower disturbance (Santorufo et al., 2012;Nagy et al., 2018;Lovei et al., 2019). Interestingly, urban forests had the highest soil mesofauna abundance. Urban forest sites were more isolated than rural forests, with an average IS ratio 290% higher. In addition, both soil pH and CPI, significantly and negatively correlated with soil mesofauna abundance, were higher in urban forests than those in rural forests. It seemed that we failed to define the specific disturbance features which overcame the negative effects of landscape and soil property changes and supported a more abundant soil mesofauna community. Nevertheless, consistent with some studies (e.g., Bogyo et al., 2015), our study suggested that light urban disturbance could increase soil mesofauna abundance, which supported the intermediate disturbance hypothesis (IDH). However, some other studies found that soil fauna abundance and diversity decrease with increasing urbanization intensity (Gray, 1989;Nagy et al., 2018;Lovei et al., 2019), thus not supporting IDH. The contradiction could be attributed to taxa-or parameterspecified sensitivity to urban disturbance (Gonzalez et al., 2016;Yoccoz et al., 2018;Nielsen et al., 2019). As in the present study, the IDH predicted the soil mesofauna abundance but not the Shannon's diversity index or Pielou's evenness index in the four urban ecosystem types. In addition, most previous studies studied the urbanization effect on soil fauna by comparing several ecosystems along an urban -rural urbanization gradient (e.g., Santorufo et al., 2014;Nagy et al., 2018;Lovei et al., 2019), which might fail to capture a complete spectrum of environmental variation due to urbanization.
Typically, urban woodland and urban park sites were of small area (fragmented), adjacent to impervious pavements or buildings. As a result, they were subjected to serious anthropogenic disturbances. The stronger disturbances could cause substantial resource availability changes for soil mesofauna and, as the present study suggested, increases in soil pH and soil heavy metal concentration. These changes may account for the decrease of soil mesofauna abundance and change in soil mesofauna community structure in these two ecosystem types.
Besides abundance and community, soil mesofauna species richness substantially decreased in urban parks compared to rural forests. In addition, though urban woodland and urban parks had similar soil mesofauna abundance, they had distinct soil mesofauna community compositions, revealing that garden managements had significant and different effects on soil mesofauna from those in urban woodlands. Garden management practices, such as sanitation, tillage, and application of pesticide, can substantially alter resource availability and quality for soil mesofauna (Norton, 2011) and change soil physiochemical properties (Norton, 2011;Tresch et al., 2018Tresch et al., , 2019, thus affecting soil mesofauna. In line with the present study, some previous studies found garden management decreased soil mesofauna abundance and diversity, thereby changing soil fauna community composition (Komarek et al., 2010;Dequiedt et al., 2011;Norton, 2011).

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this manuscript will be made available by the authors, without undue reservation, to any qualified researcher.