We collected 115 mineral soil samples from typical forests with different vegetation types across north and south China with the latitudes ranging from 18.70°N to 51.53°N (Figure 1). These soil samples harbor a wide range of soil types and edaphic and environmental characteristics (Supplementary Table S1). Mean annual temperature (MAT) and mean annual precipitation (MAP) data in sampling locations were obtained from WorldClim¹ All samples were taken in July and August, 2014. At each site, eight to ten randomly selected soil cores (0–10 cm, 5 cm in diameter) were collected within an area of about 400 m². Soil samples were combined into one composite sample for each site, and then transported at 4°C to the Institute of Applied Ecology, Chinese Academy of Sciences at Shenyang, China. Soil samples were sieved through 2-mm mesh to thoroughly homogenize and remove roots, plant detritus and stones. A portion of each soil sample was stored at -20°C until DNA extraction. The remaining soils were used to determine extractable and contents, soil microbial biomass carbon (MBC) content and soil physicochemical properties.

Soil pH was measured using a pH meter in the supernatant after shaking soil – water (1:5 w/v) mixture for 30 min. Soil total carbon, total nitrogen were determined using an Elemental analyzer (VarioEL III, Germany), while soil available phosphorus and total phosphorus were determined as previously described methods (Kuo, 1996). Soil and were extracted with 2 M KCL solution for 1 h on a shaker, and their contents were determined using a flow injection analyzer (Futura, Alliance, France). The soil MBC was estimated using the chloroform fumigation-extraction method (Vance et al., 1987; Joergensen, 1996). The soil exchangeable K⁺, Na⁺, Ca²⁺, and Mg²⁺ were determined by extracting the soils with ammonium acetate (Thomas, 1982). Amounts of Ca²⁺ and Mg²⁺ in the extracts were analyzed by atomic absorption spectrometry (AAS) and K⁺ and Na⁺ were analyzed by flame photometry.

Each soil DNA was extracted from the 0.25 g freeze-dried soil after sampling using a Mobio PowerSoil DNA Isolation Kit (MoBio Laboratories, Carlsbad, CA, USA) according to the manufacturer’s instructions. DNA was eluted with 100 μl Tris buffer (10 mM), quantified by spectrophotometer at 260 nm and stored at -20°C until use.

Soil DNA samples were sent to Novogene Company (Beijing, China) for high-throughput sequencing. The amplicon targeting V4-V5 hypervariable region of bacterial 16S rRNA was amplified with primer set 515F/806R, which contained sample specific 6-bp barcodes in the 5′ ends of them. The sequences of 515F and 806R were 5′- NNN NNN (barcode) GTG CCA GCM GCC GCG GTA A -3′ and 5′- NNN NNN (barcode) GGA CTA CHV GGG TWT CTA AT -3′, respectively. All PCR reactions were carried out in a volume of 30 μl mixture containing 15 μl of Phusion^® High-Fidelity PCR Master Mix (New England Biolabs), 0.2 μM of each primer, about 10 ng template DNA, and ddH₂O filled to 30 μl. Thermal cycling included an initial denaturation at 98°C for 1 min, followed by 30 cycles of denaturation at 98°C for 10 s, annealing at 50°C for 30 s, and elongation at 72°C for 1 min, with a final extension at 72°C for 5 min. PCR products were detected by 2% agarose gel electrophoresis, and those with bright main strip between 400 and 450 bp were chosen for further experiments. Equal amounts of the PCR product from each sample were pooled and then purified with GeneJET Gel Extraction Kit (Thermo Scientific). The sequencing library was generated using NEB Next^® Ultra^TM DNA Library Prep Kit for Illumina (NEB, USA) following the manufacturer’s instructions and thus sequencing adapters were added to 5′ ends of amplicon. The library quality was assessed on the Qubit^@ 2.0 Fluorometer (Thermo Scientific) and Agilent Bioanalyzer 2100 system. At last, the qualified library was sequenced on the Illumina MiSeq platform, producing 250 bp/300 bp paired-end reads.

Paired-end reads from the original amplicon were merged using FLASH (Magoč and Salzberg, 2011) which is designed to merge paired-end reads when there are overlaps between reads1 and reads2. Paired-end reads was assigned to each sample according to the unique barcodes which were removed together with primers subsequently. Sequences were analyzed using QIIME software package (Quantitative Insights Into Microbial Ecology²) (Caporaso et al., 2010), and in-house Perl scripts were used to analyze alpha-(within samples) and beta-(among samples) diversity. First, merged reads were filtered by QIIME quality filters. Then the clean sequences obtained with ≥97% similarity level were assigned to the same operational taxonomic units (OTUs). A representative sequence from each OTU was picked and annotated using the RDP classifier for taxonomic information (Wang et al., 2007) and aligned with “Core Set” in the GreenGene database for phylogenetic information (DeSantis et al., 2006). In order to unify the survey (Shaw et al., 2008), a subsample of randomly selected 7300 sequences in each sample was used for bacterial alpha-diversity (phylotype richness and phylogenetic diversity) and beta-diversity (community dissimilarity index) analyses. Observed Species and Phylogenetic Diversity Whole Tree (PD for short) indexes formed during alpha-diversity analysis were used to indicate phylotype richness and phylogenetic diversity in samples, respectively. Unifrac metric was used to compare the difference of overall community composition between each pair of samples (Lozupone and Knight, 2005), and thus generating unweighted and weighted pairwise unifrac distance matrixes. All sequences in this study are available in Sequence Read Achieve (SRA) database of NCBI under accession number SRP070864.

Correlation (Pearson’s rank correlation) or regression analysis between soil/site characteristics and individual phyla or between soil/site characteristics and indexes of community diversity and composition were performed in SPSS 17.0 for Windows. The other statistical analyses were conducted using the program R v.3.2.0 (R Development Core Team, 2006). The “Bray–Curtis” dissimilarity matrix for the bacterial community composition and the “Euclidean” dissimilarity matrices for geographic distance and environmental variables were constructed with the “vegdis” “function in the “vegan” package (Oksanen et al., 2016). The non-metric multidimensional scaling (NMDS) and cluster analysis of soil samples in the bacterial community composition was conducted with the “metaMDS” (Minchin, 1987) and “hclust” functions (Murtagh, 1985) based on the “Bray–Curtis” dissimilarity matrix within the package “vegan,” respectively. Additionally, we conducted 1 minus Unifrac distance (unweighted or weighted) in the total community structure to estimate the bacterial community similarity. BioEnv procedure (Clarke and Ainsworth, 1993) was performed to select the environmental variables which were further used to construct environmental distance matrix with the “vegdist” function. Using principle coordinates of neighbor matrices (PCNM) method (Borcard and Legendre, 2002), the geographic coordinates of the sites were transformed to significant vectors that could be used to construct geographic distance matrix across sites. Mantel tests with 999 permutations (Legendre and Legendre, 2012) were used to examine the correlation (Pearson’s rank correlation) between geographic or environmental distance and bacterial community distance within the vegan package. The canonical correspondence analysis (CCA) (Legendre and Legendre, 2012) was employed to identify the most important soil environmental factors shaping bacterial community structure. Monte Carlo permutation test (permutest) and “envfit” functions (Legendre et al., 2011) were used to test the significant environmental variables during CCA analysis. These significant PCNM vectors and environmental variables were used as explanatory variables in constrained ordinations (CCA) for variation partition analysis.

The latitude of each sampling site was highly correlated with the site’s MAT (r = -0.982, P < 0.001) and MAP (r = -0.873, P < 0.001) (Table 1). Soil pH showed significant correlations with concentrations of three types of exchangeable cations, i.e., K⁺ (r = 0.685, P < 0.001), Ca²⁺ (r = 0.843, P < 0.001), Mg²⁺ (r = 0.776, P < 0.001). Soil pH and these three cations increased with increasing latitude (Table 1). Soil total C (TC) was significantly positively correlated with soil total N (TN) (r = 0.733, P < 0.001), and they both were correlated with soil pH and these three cations. The latitudes of sampling sites were not significantly correlated with TC, TN, soil or contents. Soil and contents were significantly correlated with soil pH and TC, TN, exchangeable Ca²⁺ and exchangeable Mg²⁺ contents. Soil MBC had no significant relationship with other soil and site characteristics except for soil (Table 1).

We obtained 4,667,656 sequences from all 115 samples, with an average of 40,588 sequences per sample. The range of sequences per sample in the whole dataset was from 7317 to 190250, and most samples (80%) had sequences between 16000 and 65000 (12 samples had less than 16000 sequences and 11 samples had more than 65000 sequences). The length of these sequences ranged of 191–351 bp, with a mean of 253 bp. Among these sequences, 98.6% could be classified. At the 97% similarity level, the sequences in all soils could be grouped into 325,433 phylotypes, with an average of 2,830 phylotypes per sample. Actinobacteria, Acidobacteria, Alphaproteobacteria, Verrucomicrobia and Planctomycetes (relative abundance >5%) were dominant groups across all sequence data, and they accounted for more than 73% of the bacterial sequences (Supplementary Table S2). Moreover, groups of Chloroflexi, Betaproteobacteria, Deltaproteobacteria, Gammaproteobacteria, Gemmatimonadetes, Nitrospirae, Bacteroidetes and AD3 (relative abundance >1%) were less abundant (accounting for 22% of the bacterial sequences), but still existed in all soils. The rest of sequences could be classified into 56 groups, and 40 groups were rare (relative abundance <0.01%) (Supplementary Table S2).

The abundance of some dominant bacterial groups was significantly correlated with soil and site characteristics (Table 2). The abundance of Verrucomicrobia (r = 0.472, P < 0.001; r = 0.511, P < 0.001), Gemmatimonadetes (r = 0.628, P < 0.001; r = 0.328, P < 0.001) and Armatimonadetes (r = 0.409, P < 0.001; r = 0.464, P < 0.001) increased with increasing geographic latitude; while the abundance of Alphaproteobacteria (r = -0.429, P < 0.001; r = -0.331, P < 0.001) and Gammaproteobacteria (r = -0.601, P < 0.001; r = -0.523, P < 0.001) decreased with increasing geographic latitude (Table 2). These five bacterial groups were also closely related to MAT and MAP of their locations (Table 2).

Soil pH and exchangeable K⁺, Ca²⁺, and Mg²⁺ contents affected many bacterial groups. For example, the relative abundance of Chloroflexi, Betaproteobacteria, Deltaproteobacteria, Gemmatimonadetes, Nitrospirae and Bacteroidetes was positively correlated with soil pH and exchangeable K⁺, Ca²⁺, and Mg²⁺ contents, while the relative abundance of Alphaproteobacteria was negatively correlated with these parameters (Table 2; Figures 2 and 3). Although the Acidobacteria group did not show correlation with soil pH and exchangeable K⁺, Ca²⁺, and Mg²⁺ contents, most of Acidobacteria subgroups had significant correlations with these parameters (Supplementary Table S3; Supplementary Figures S1 and S2).

Soil total carbon, nitrogen and phosphorus were important factors for some bacterial groups. For example, these three parameters were all positively correlated with the relative abundance of Nitrospirae and Bacteroidetes (Table 2; Supplementary Figures S3 and S4). The relative abundance of Deltaproteobacteria and Gammaproteobacteria showed positive correlations with TC and TN, while the abundance of Actinobacteria and Firmicutes showed negative correlations with TC and TN. In addition, the relative abundance of Alphaproteobacteria was only correlated with TP but not with TC and TN.

Soil available phosphorus (AP) and extractable and contents are nutrients directly used by microbes and therefore were related to many bacterial groups. Soil AP was positively correlated with Bacteroidetes but negatively correlated with Alphaproteobacteria. The relative abundance of Actinobacterial, Betaproteobacteria, Deltaproteobacteria, Gammaproteobacteria, Nitrospirae and Bacteroidetes had significantly positive relationship with soil extractable and contents. Soil MBC was only correlated with the relative abundance of Gemmatimonadetes group.

The pairwise correlation analysis between the indices of the soil bacterial community diversity and soil location, or soil physical and chemical properties showed that soil pH, TC and exchangeable K⁺, Ca²⁺, and Mg²⁺ and extractable and contents were significantly positively correlated with both phylotype richness and phylogenetic diversity (Table 3; Figure 4). Soil TN was positively correlated with phylogenetic diversity. A parabolic relationship was found between MAT and bacterial phylotype richness, and between MAT and phylogenetic diversity (Figure 4). Latitude presented similar effects on these indices of the soil bacterial community diversity with the tipping point at between 33.5°N and 40°N.

Mantel test showed that “Bray–Curtis” distances of bacterial community composition was positively correlated with geographic distances (r = 0.259, P < 0.001, Table 4). Except for soil available phosphorus and MBC, all other examined environmental variables presented significant correlations with soil microbial community composition (Table 4). Soil pH was most strongly correlated with bacterial community composition (r = 0.700, P < 0.001). Utilization of both soil pH and soil TN predicted bacterial community composition better (r = 0.708, P < 0.001), while the addition of the other factors did not improve the regression’s efficiency. Moreover, NMDS visualization showed that variation in bacterial community composition was associated with variation in soil pH and differences in geographic region (Figure 5). Similarly, the significantly linear relationship between NMDS1 of NMDS scores and soil pH confirmed the importance of soil pH (Figure 6A). Soil exchangeable Ca²⁺ content had similar effects on soil bacterial community composition to soil pH (Figure 6B). Additionally, both NMDS1 and NMDS2 scores were closely correlated with MAT and MAP (Figures 6C–F). CCA analysis showed that soil pH, exchangeable Ca²⁺ and Mg²⁺, soil , and MAT and MAP of soil locations were important factors on soil bacterial community composition (Supplementary Figure S6). The directions of pH, exchangeable Ca²⁺ and Mg²⁺ were closely correlated with CCA1, while the directions of soil , MAT and MAP were correlated with both CCA1 and CCA2.

Based on the “Bray–Curtis” dissimilarity matrix, the bacterial communities of the 115 soils were roughly clustered into two big groups (Figure 7). Group I consisted of 32 samples which were mainly from southern forests of China with low latitudes (ranging from 18.70°N to 29.65°N). Group II was composed of 83 soils which were located in forests of northern China with middle latitudes ranging from 31.30°N to 51.53°N. Group II could be further divided into two subgroups (A and B). Subgroup A was consisted of 56 soils, most of which were sampled from northeastern China. Subgroup B contained 27 soils, which were sampled from Beijing, Qinling and Shennongjia with the latitude ranging from 31.30°N to 39.96°N (intermediate zone). These results corresponded to the results of NMDS analysis (Figure 5), which showed that bacterial community structure differed greatly between forest soils with acidic pH at lower latitude sites (less than 30°N) and forest soils with near-neutral, neutral or weakly alkaline pH at mid latitude (31°N – 40 °N) and high latitude sites (41°N – 52°N) in China. Additionally, the relationship between geographic distances, environmental distance and bacterial community similarity in community composition (Figure 8) indicated that more distinct bacterial communities could be found in two soils far from each other than in two soils with a near distance and also in two soils with more different soil properties.

The variance partitioning analysis showed the relative contributions of the geographic distance and environmental variables to the bacterial community structure (Figure 9). CCA analysis selected a subset of environmental variables (MAT, MAP, pH, TP, EK, ECa, EMg, TC, TN, and ) which together explained 21.11% of the bacterial community variation, more than the geographic distance (15.88%) (Figure 9). Therefore, the soil characteristics and environmental factors were more important than the geographic dispersal limitation in determining the bacterial community structure in Chinese forest soils. These selected environmental variables, i.e., MAT, MAP, pH, TP, EK, ECa, EMg, TC, TN, and , explained 1.86, 2.29, 3.05, 1.33, 1.91, 2.85, 2.48, 1.19, 1.04, 1.04, and 2.07% of the bacterial community variation, respectively. These environmental variables combining geographic distance explained 30.16% of the bacterial community variation, leaving 69.84% of unexplained variation, indicating that the overlapping effect of environmental variables and geographic distance on the bacterial community variation was 6.83% (Figure 9), and that there were many unmeasured or unknown factors that contributed to the large portion of unexplained variation in this study.

We found Actinobacteria phylum was the most dominant group (22%) in our studied forest soils. Acidobacteria accounted for 18% of all bacterial communities, while the relative abundance of Bacteroidetes was only about 1.4%. These results only partly agree with findings in the 88 soils across North and South America (Lauber et al., 2009) and the 26 black soils in northeastern China (Liu et al., 2014). They both found Acidobacteria was the most abundant phylum and Bacteroidetes was 11.2 and 5.6% as reported by Lauber et al. (2009) and Liu et al. (2014), respectively. Moreover, the relative abundance of Verrucomicrobia (8.68%) and Planctomycetes (6.75%) phyla in our study was much higher than results in Lauber et al. (2009) (0.9 and 0.09%, respectively) and Liu et al. (2014) (3.22 and 4.85%, respectively). Previous studies may have underestimated the abundance of Verrucomicrobia due to the bias of primers (Bergmann et al., 2011). However, using the same primers as this study, Fierer et al. (2012) found a huge variability in the relative abundance of the major bacterial taxa among different biomes including tropical forest, temperate forest, and boreal forest soils collected from different sites (for example, 5.22–40.29% for Verrucomicrobia and 1.61–5.56% for Bacteroidetes). Our results agreed with their findings and suggested that at large spatial scales, the dominant bacterial groups may be quite different among different regions.

We further studied the controlling factors of the dominant bacterial groups and found different bacterial groups responded differently to soil properties and local climate (MAT and MAP) gradients (Table 2). For example, Alphaproteobacteria and Gammaproteobacteria were more abundant in southern tropical zone than in northern temperate zone while Verrucomicrobia, Gemmatimonadetes, and Armatimonadetes presented the opposite trend. MAT and MAP may have played an important role in influencing some bacterial taxa, especially for Verrucomicrobia and Armatimonadetes, which had no or weak relevance to soil properties but strong relationship with local climate (Table 2). A significant negative correlation was observed between MAP and the relative abundance of Verrucomicrobia in forest soils in this study (Table 2), while the opposite trend was found in grassland soils of the arid and semiarid areas in China (Wang et al., 2015). Therefore, the responses of Verrucomicrobia to climatic conditions were different between arid/semiarid areas and semi-humid/humid areas, and that the intermediate amount of precipitation (e.g., 400–500 mm MAP) was probably the most beneficial for Verrucomicrobia. In general, the relative abundance of most phyla was positively correlated to soil parameters such as pH, exchangeable cations, C, N, P, and (Table 2; Figure 2), which indicated that most bacterial taxa exhibited copiotrophic attributes and seemed to be favored by neutral pH and high carbon availability. Prior studies also found that most bacteria benefits from optimum living conditions (McCaig et al., 1999; Axelrood et al., 2002; Padmanabhan et al., 2003; Fierer et al., 2007). However, Actinobacteria, Alphaproteobacteria, Firmicutes and the dominant subgroups of Acidobacteria (GP1 to GP3) showed negative relationships with those soil parameters (Table 2, Supplementary Table S3; Figure 3, Supplementary Figures S1 and S2). Therefore, they possibly had oligotrophic lifestyle and were adapted to low-nutrient and low pH soils, which were consistent with previous reports (Fierer et al., 2007; Naether et al., 2012).

We found the relationship between soil bacterial diversity and latitude was a parabola shape with the tipping point (maximum) at between 33.50°N and 40°N (Figure 4), which falls into the warm-temperate zone of China with MAT of about 7–9°C (Supplementary Table S1; Figure 4). The changes of animal and plant diversity with latitudinal gradients have been well documented and studied for centuries with the well-established conclusion that plants and animals generally exhibit an increase of diversity with decreasing latitude (Lomolino et al., 2006). However, whether microbial diversity also exhibits a latitudinal gradient has not reached a general conclusion. Some studies found that bacterial diversity increased or decreased with latitude (Buckley et al., 2003; Fuhrman et al., 2008; Liu et al., 2014), while other studies found no relationship between bacterial diversity and latitude (Fierer and Jackson, 2006; Corby-Harris et al., 2007; Lauber et al., 2009; Chu et al., 2010). The latitudinal trend found in our study was based on the large scale data with a latitude range between 18.70°N and 51.53°N. The parabola shaped relationship between bacterial diversity and latitude indicated that optimum conditions for highest bacterial diversity in our studied areas were located at the warm-temperate zone and bacteria did not show a simply latitudinal diversity gradient as previously reported (Buckley et al., 2003; Fuhrman et al., 2008; Liu et al., 2014). In the areas with higher bacterial diversity, most forests belong to the deciduous broad-leaved forest type, which have higher substrate availability (such as C or N sources) for bacterial growth compared to coniferous forest (Huang et al., 2004; Geng et al., 2009), or have a more comfortable soil physical environment for bacteria (Wallenstein et al., 2007). The near-neutral pH (6–7) and temperate climate conditions (Figure 4) may also contribute to the high bacterial diversity in this area.

According to the dissimilarity matrix in the whole bacterial community composition, all 115 soil samples in this study were clustered into two main groups (Figure 7), which were located in the low latitude zone (18.70°N to 29.65°N) and the middle latitude zone (31.30°N to 52.53°N) (Supplementary Table S1). Moreover, a significant correlation was observed between geographic distances and bacterial community dissimilarities (Figures 8A,B). These findings suggested that the soil bacterial communities in the Chinese forest soils zone were distributed geographically. It is clear that soil bacterial communities in southern forests are distinct from those in northern forests and 30°N could be coarsely considered as the dividing line between them. However, some soil samples in tropical or subtropical forests (for examples JFL07, JFL08, JFL09, HNZZ01 and HNZZ03) were clustered into group A (mainly consist of northern temperate and boreal forest soils); while some bacterial community structure in temperate forest soils (BSLHJL and QL03) were more similar to that in southern forest soils (group I). This result suggested that soil properties and climatic factors (Figures 6 and 8C,D) also play a very important role in determining bacterial community composition.

The variation of soil bacterial phylotype and phylogenetic diversity along latitude gradients was similar to that along soil pH, ECa²⁺, EMg²⁺ and TC gradients (Figure 4; Supplementary Figure S5). Therefore, these soil properties and climatic conditions (MAT and MAP) together determined the biogeographic distribution of soil bacteria in our studied areas although these parameters were inter-correlated with each other (Table 1). Multiple variable analysis suggested that soil pH was the most important determinant of soil bacterial community structure (Figures 5 and 6; Supplementary Figure S6), which has been broadly documented in soils from a broad range of ecosystems (Fierer and Jackson, 2006; Baker et al., 2009; Liu et al., 2014). In fact, the cluster analysis results which divided bacterial composition into southern (Group I) and northern forest (Group II) groups (Figure 7) was also related to changes of soil pH because most soils (28 out of 32 samples) from Group I were acidic with pH < 5, while the soils classified into group II mainly (74 out of 83 samples) had pH > 5 (Figure 5, Supplementary Table S1). Our results (Figure 4) agreed with previous findings that acidic soils usually showed lower phylogenetic diversity than neutral soils (Fierer and Jackson, 2006; Lauber et al., 2009; Chu et al., 2010) and were mostly dominated by particular taxa (Griffiths et al., 2011). Therefore, soil pH = 5 probably can be used as a dividing line between northern and southern China regarding soil bacterial community composition and also a threshold below which soil bacterial diversity may decline and soil bacterial community structure may change significantly.

It is noteworthy that some soil exchangeable cations (Ca²⁺, Mg²⁺, K⁺) were also correlated with soil bacterial composition and diversity and some taxonomic groups (Table 3, Figures 3 and 4; Supplementary Figure S2). Soil pH strongly influenced these cations (Table 1), which agreed with previous findings conducted in tropical soils (Sanchez, 1977; Fearnside, 1984) and ferrosol soils (Lacey and Wilson, 2001). Therefore, it is understandable that the effects of these cations on microbial diversity were similar to the effects of soil pH. Likewise, the threshold of soil ECa²⁺ at about 32 mmol kg^-1 can also be used for partitioning bacterial community structures into group I (31 in 32 samples with ECa²⁺s < 32 mmol kg^-1) and group II (82 in 83 samples with ECa²⁺s > 32 mmol kg^-1) (Supplementary Table S1), suggesting that soil ECa²⁺ may be another suitable marker for cluster-dividing of soil bacterial community composition besides soil pH.

We found environmental factors played a more important role in driving bacterial community pattern than geographic distance (explained 21.11 and 15.88% of the variation in bacterial community structure, respectively) at the large spatial scale of this study (Figure 9). This result was similar to the results obtained at a smaller scale in the black soils of northeast China (Liu et al., 2014), but different from a study conducted along a transect of arid and semi-arid grasslands in northern China, which showed geographic distance (36.02%) explained more of the variation in bacterial community structure than environmental variables (24.06%) (Wang et al., 2015). This was probably because the latter study was conducted along a latitudinal transect and the variations of climatic conditions and geographic distances were highly constrained. In addition, there was a significant correlation between geographic and environmental distance with a weak strength (‘Mantel test’, r = 0.16 for pearson’s rank correlation and r = 0.36 for spearman’s rank correlation, respectively, p = 0.001, data not shown) across sites in this study, and this indicates that variation in bacterial community composition may be associated with both geographic distance and environmental dissimilarity between sites.