Deciphering Centimeter-Scale Microbial Biogeographic Patterns and Community Assembly Processes From A 3D Perspective

Background: Revealing the effects of multi-dimensional spatial distribution on microorganisms is crucial for the further understanding of microbial diversity, turnover and ecological processes. However, microbial community assembly and the factors that shape it are still unknown from a three dimensional (3D) perspective. Here, a 3D model was created by performing an exhaustive sampling strategy to a 4x4x4 soil matrix. We examined the dynamics of microbial diversity, biogeographic patterns and microbial assembly processes when transfroming sampling scheme from 2D to 3D. Results: Our results indicated that dispersion of microbial community and signicance of distance decay relationship was higher in the 3D compared with 2D sampling scheme, suggesting increased microbial turnover when transforming the model from 2D to 3D. Only a small fraction of community variation can be explained by environmental, spatial factors and spatial canonical axes, possibly due to unmeasured environmental variables. The assembly of microbial community was dominated by deterministic processes that shifted from homogeneous selection to variable selection as we transformed the model from 2D to 3D. The importance of stochasticity increased when homogeneous and variable selection processes were well balanced. However, heterogeneity of existing environmental and spatial variables failed to explain the dynamics of community assembly. Conclusions: Our study revealed signicant dynamics of microbial diversity and assembly processes when assessed from 2D and 3D perspectives. As microorganisms are spatially distributed in soil, this spatial dependent diversity and assembly suggests that microbial ecological questions need to be considered in more dimensions than they usually are. Further, new models that integrate all data sets are still needed to disengle the microbial processes in multiple dimensions.


Background
Biogeography is the study of the distribution of organisms across space and time and has long been applied to higher organisms such as plants and animals [1,2]. In recent decades, as the development of genetic methods has allowed researchers to explore microbial information to greater depths, microbial biogeographic patterns have also been examined in a wide variety of environments such as soils, oceans and the phyllosphere [3][4][5]. However, characterizing the underlying processes that shape microbial biogeographic patterns continues to be a challenge [5,6].
The Distance Decay Relationship (DDR) is a fundamental biodiversity phenomenon that has been regarded as evidence of biogeographic patterns [7,8]. The relationship is observed as a decreasing relatedness in microbial community similarity with increasing spatial distance and can be in uenced by both deterministic and stochastic processes. Deterministic processes are niche based and include abiotic selection pressures such as climate and nutrient conditions, and biotic selection pressures such as competition and trade-offs. Stochastic processes are neutral based and include dispersal, drift and mutations independent of species traits [9]. In general, selection and drift strengthen DDR, dispersal weakens DDR, and mutations modify the variance and height of DDR slope [7]. It is, however, di cult to quantify the relative importance of these ecological processes. Investigations into the relative importance of deterministic and stochastic processes in shaping microbial communities have shed some light on the mechanisms controlling diversity, function, succession and biogeography [10][11][12]; null model approaches have been used to disentangle the in uence of these two processes and further quantify the relative importance of selection, dispersal and drift on community assembly in soil systems [13][14][15].
Ecological processes operate and overlap at all scales of life and result in complex, multiscale patterns ranging from the microscopic particle to the superregional biome. Therefore, identifying pattern scales and relating them to ecological processes are paramount goals in modern ecology [16]. Spatial scaling plays an important role in microbial biogeographic patterns of community diversity with DDR and microbial assembly processes reported for a wide range of habitats at local, regional and continental scales [5,17,18]. However, as most studies focus on large scale models in two dimensions (2D; horizontally or vertically), there is still limited knowledge of microbial turnover and assembly processes from a three-dimensional (3D) perspective. This research gap prevents a thorough understanding of community ecological processes and biogeography because it underestimates the importance of microenvironmental variables that shape microbial composition and activity [19]. This limitation also masks the spatial heterogeneity of soil physical, chemical and microbial properties that result in a variety of 3D niches and habitats [20].
Microorganisms consume little space compared to plants and animals, and are mainly in uenced by the abiotic and biotic factors of their microenvironment [21,22]. Thus, microbial composition and assembly varies at small scales. However, typical soil sampling strategies for microbial research generate samples vastly larger than the organisms of interest [23]. This limitation is exacerbated by a microbe's relatively short life span and poor mobility which means that microbial community habitats may be far smaller than expected. The 3D nature of microbial distributions can be observed in microbial abundance and species richness in and between aggregates [24] and in different pore size classes in undisturbed soil [20]. Although DDR has been detected within several meters [22,25], there is limited literature on microbial biogeographic patterns and assembly processes at the centimeter scale. Understanding small scale distributions of soil microbial communities will help narrow down the factors driving community structure and may link structure to observed heterogeneity in biogeochemical processes.
Our study addresses this research gap through an intensive sampling strategy to determine microbial diversity, biogeographic patterns and assembly processes at a centimeter scale from a 3D perspective.
This strategy used an 8 cm 3 soil cube from a paddy eld divided into a 4 × 4 × 4 cubic grid (64 sub-cubes) to construct a 3D soil matrix by introducing a third axis into the 2D system. Soil basic properties and 16S rRNA genes of each sub-cube were analyzed to answer the following ecological questions: (1) do microbial biogeographic patterns differ between 2D and 3D sampling schemes; and (2) which factors are related to any differences in microbial assembly processes in 2D compared with 3D sampling schemes?

Materials
Soil sampling and property analyses The soil was sampled in November 2017 (autumn) from a long-term rice cultivation eld after harvest in Wenling, Zhejiang province, China (28°21'N, 121°15'E). This region has a subtropical monsoon climate with a mean annual rainfall of 1650 mm and a mean annual temperature of 17.3°C. The soil is a paddy soil and classi ed as a loamy mixed active Thermic Endoaqualf. Undisturbed surface soil (0-10m) was collected with a soil corer (internal diameter of 50.46 mm) and stored at 4°C until micro-computed tomography (mCT) scanning. Then, an 8 cm 3 cube (2.0 cm side length) of the sample was selected and cut into 64 sub-cubes (4x4x4 with 0.5 cm side length, Fig. 1). To compare the differences in microbial community composition and assembly processes from a 3D view, three sampling schemes were designed to gradually transition the data from two to three dimensions. The rst sampling scheme (s1; 2D) used x and y axes only, thus dividing the soil matrix into 16 samples (4 sub-cubes of each column united as one sample, Fig. 1-z1). Next, the second sampling scheme (s2; 2-3D transition) introduced a single z axis into s1, thus separating the 16 4x1 columns into 32 2x1 columns ( Fig. 1-z2). Finally, the third sampling scheme (s3; 3D) introduced two further z axis divisions of the columns, separating the soil cube into 64 equally proportioned samples ( Fig. 1-z3). This sampling strategy was performed from three orientations to increase the robustness of the data; there are therefore three sampling schemes repeated three times: s1 (x1, y1 and z1), s2 (x2, y2 and z2) and s3 (same samples for all orientations; see Figure 1). Soil structure was analyzed by mCT scanning as described by Yu and Lu [26]. Brie y, the soil matrix was scanned with an industrial Phoenix v|tome|x m X-ray CT (General Electric, GmbH, Wunstorf, Germany).
After which a CT reconstruction of the selected soil matrix was performed by Datos | × 2.0 software (Wunstorf, Germany) with a ltered back-projection (FBP) reconstruction algorithm. Digital image processing and analysis was performed by Image J 1.51W (National Institute of Health, Bethesda, MD, USA; http://rsb.info.nih.gov/ij/). Quanti cation of soil porosity was completed by AVIZO 9.2 (Hillsboro, OR, USA). The soil carbon (C), nitrogen (N), hydrogen (H) and sulfur (S) content in each sub-cube was determined by Vario EL Cube elemental analyzer (Elementar, Germany).

DNA extraction and amplicon sequencing
After mCT scanning, sub-cube DNA was extracted using a MP FastDNA SPIN Kit for soil (MP Biomedicals, Solon, OH, USA) as per the manufacturer's instructions. 16S rRNA genes were ampli ed by primers 515F (5′-GTGCCAGCMGCCGCGGTAA-3′) and 907R (5′-CCGTCAATTCCTTTGAGTTT-3′). Ampli ed PCR products were sequenced using an Illumina HiSeq PE250 sequencing platform (Illumina, San Diego, CA, USA). Sequences were clustered into operational taxonomic units (OTUs) with USEARCH11 using a sequence similarity threshold of 0.97 [27]. Across all samples, we obtained 761 137 quality controlled unique sequences with a total of 10 898 OTUs classi ed up to phylum level (Additional le: Figure S1). RDP training set v16 was used for taxonomy annotations with a con dence cut-off of 80% [28] after which aligned sequences of representative OTUs were used to construct a maximum-likelihood tree. The community matrix was rare ed to the same depth (41752 OTUs per sample) in subsequent analysis. The relative abundance of the top 10 phyla was selected to give an overview of microbial structure (Additional le: Figure S1). To keep a consistent total OTU richness, we use the average OTU abundance of all the sub-cubes in each sample of S1 and S2 in the following data analyses.

Statistical analyses
All statistical analyses were performed using the R environment (v3.6.1; http://www.r-project.org/). Microbial α-diversity across the three sampling schemes was determined using Shannon index, Faith's diversity and observed OTUs as calculated using 'vegan' package. Principal co-ordinates analysis (PCoA) using the weighted-Unifrac dissimilarity matrix and a PERMDISP test (999 permutations) were performed to compare microbial structure of the three sampling schemes. To assess the relationship of community dissimilarity with spatial distance and environmental distance in the 3D soil matrix, we rst located the geometric center of each sample of the three sampling schemes using 3D coordinates. Then, Euclidean distances of each pair of samples were determined to generate a distance matrix. For the analysis of community similarity, Bray-Curtis distance was calculated using the 'vegdist' function in the R package 'vegan' [29]. Environmental distance was built using the Euclidean distance of environmental variables between each pair of samples. Mantel and Partial Mantel tests (Pearson, permutations=999) were performed to test the effect of environmental and spatial distance in shaping microbial community.
We further performed a redundancy analysis (RDA) based variation partitioning analysis (VPA) with adjusted R 2 coe cients to quantify the effect of environmental distance, spatial distance and spatial canonical axes (x, y or z) on community similarity in each sampling scheme [30]. For this analysis, we rst standardized OTU abundance using the Hellinger method and then calculated community similarity matrix as a response variable. Environmental variables including soil C, N, H, S and porosity were used in the calculation. For spatial variables, we extracted a set of spatial variables generated using principal coordinates of neighbor matrices (PCNM) analysis based on the 3D coordinates of each sample. For spatial canonical axes, we used the coordinates of each sample at each axis [30]. Before the VPA, the importance of environmental variables, PCNM variables and spatial canonical axes in explaining community similarity was determined by a separate RDA analysis using Monte Carlo permutation tests (999 unrestricted permutations). Then, signi cant environmental variables, spatial variables and coordinates were selected following a forward selection by permutation of residuals under a reduced model from each of the explanatory sets (no VPA was performed on x1 as no environmental and canonical axes were selected as signi cant variables, Additional le: Table S2).
The potential importance of stochastic and deterministic processes on community assembly of the three sampling schemes was determined using a null model [14]. To infer community assembly processes, we calculated β-nearest taxon index (βNTI) which is the difference between observed β mean nearest taxon distance (βMNTD) and the mean of the null distribution of βMNTD normalized by its standard deviation following the description of Stegen et al [13]. Generally, |βNTI| >2 means the turnover in phylogenetic community composition is governed by deterministic process. The relative contribution of stochastic process was estimated as the percentage of pairwise comparisons with |βNTI| <2. The relative contributions of variable and homogeneous selection were estimated as the percentage of pairwise βNTI values that fell above +2 and below −2, respectively [15]. Further, a Bray-Curtis based Raup-Crick matrix (RC bray ) was calculated to partition the pairwise comparisons that were assigned to stochastic processes [14]. The relative contribution of dispersal limitation was estimated as the percentage of pairwise comparisons with |βNTI| <2 and RC bray >+0.95 while the relative contribution of homogenizing dispersal was estimated as the percentage of pairwise comparisons with |βNTI| <2 and RC bray <−0.95 [14]. Pairwise comparisons that did not fall into any of these fractions indicate that no single process dominated community assembly. This undominated fraction was therefore estimated as the percentage of pairwise comparisons with |βNTI| <2 and |RC bray | <0.95.
To assess the effect of environmental variation and spatial distance on microbial assembly processes in each sampling scheme, βNTI values were regressed against spatial distance and environmental distance using a Mantel test with 999 permutations. Similarly, to assess the relationship between βNTI and environmental variables or spatial distance after controlling for spatial or soil environmental distance, we also performed a partial Mantel test with 999 permutations while controlling for the other distance variable.

Results
Microbial diversity in 2D and 3D sampling schemes To understand how the bacterial community structure varied across different sampling schemes, we examined microbial taxonomic diversity. Microbial α-and β-diversity indicated signi cantly different structures among the three sampling schemes. α-diversity (Shannon index, Faith's diversity and observed OTUs) signi cantly decreased from s1 to s3 in all three orientations (one-way ANOVA, p<0.001, Fig. 2a).
The community structures of the three sampling schemes were not clearly separated in two-dimensional plot ( Fig. 2b-d). However, we found increased β-diversity dispersion from s1 to s3 in all three orientations (PERMDISP, p<0.001), indicating that the microbial compositions were more clustered in 2D compared with 3D sampling schemes.
Microbial biogeographic patterns in 2D and 3D sampling schemes Distance decay relationship was examined in both 2D (s1) and 3D (s2 and s3) sampling schemes to detect the spatial effect on microbial biogeographic patterns. We found signi cantly negative correlations (p<0.05) between community similarity and spatial distance both before and after controlling for environmental distance (except x1, Fig. 3, Additional le: Table S1). Further, the signi cance levels (p value) strengthened from s1 to s3 in all three orientations. However, most of the t values were low (R 2 < 0.1), indicating weak DDR of community similarity with spatial distance. Community similarity signi cantly decreased from s1 to s3 which indicating microbial turnover rate increased as the third axis was introduced into the model. Environmental distance signi cantly (p<0.05) increased with spatial distance in the z orientation but also with low t value (R 2 < 0.1, Fig. 3). Environmental distance increased from s1 to s3 but not in signi cant level. Environmental distance only signi cantly negatively correlated with community similarity before and after controlling for spatial distance in z1 (p<0.05) with low t values (R 2 <0.1, Additional le: Figure S2, Table S1).
Environmental factors, spatial factors and canonical axes shaped the microbial community in 2D and 3D sampling schemes We performed a VPA to quantify the community variation explained by environmental factors, spatial factors and spatial canonical axes in the 2D and 3D sampling schemes. Only a small proportion of microbial community variation can be explained by environmental factors, spatial factors, spatial canonical axes and their shared effect in the three sampling schemes (10-53%, Fig. 4). In general, pure environmental factors, spatial factors and spatial canonical axes explained up to 18%, 7% and 7% of the total variation, respectively. The community turnover was mostly in uenced by the spatial factors from y orientation and z2 (p<0.05, based on 999 Monte Carlo permutations). Although a pure canonical axes effect was not evidenced in most cases, interactions with spatial factors could explain up to 19% of community variation. The forward selection procedure selected different subsets of environmental variables, PCNM and spatial axes which had signi cant effects on microbial composition in the three sampling schemes (Additional le: Table S2). A large amount of microbial variation (up to 90%) remained unexplained, possibly resulting from stochastic processes and unmeasured environmental variables.
Relative importance of deterministic and stochastic processed in shaping microbial assembly in 2D and 3D sampling schemes A null model was performed to determine the relative importance of deterministic and stochastic processes that substantially changed in the three sampling schemes (Fig. 5). Generally, deterministic process dominated community assembly processes in all three sampling schemes (more than 68%). βNTI values signi cantly increased from s1 to s3 in all three orientations, suggesting that microbial assembly processes shifted from homogeneous selection (βNTI<-2) to stochastic processes and then to selective selection (βNTI>2) (Fig. 5a). Homogeneous selection dominated the assembly processes in s1 (85%, average of the three orientations), then decreased in s2 (12%) before reaching a low of 2% in s3. In contrast, variable selection increased from s1 (3%) to s2 (64%), reaching a high of 92% in s3. Consequently, the proportion of stochastic processes (βNTI>|2|) increased from s1 (11%) to s2 (24%) and then decreased in s3 (6%). Homogeneous dispersal dominated stochastic processes (more than 70%) of the total assembly processes in all the three sampling schemes. Pairwise comparisons of βNTI values was signi cantly correlated with spatial distance both before and after controlling for environmental distance in x2 and z2 (Additional le: Table S2). However, no relationship was found between βNTI values and environmental distance or between βNTI values and the differences of individual environmental variables before and after controlling for spatial distance (Additional le: Table S3).

Discussion
In this study, microbial biogeographic patterns and the ecological assembly of microbial communities in paddy soil was examined from both 2D and 3D perspectives at the centimeter scale. We found increased variation of microbial composition and stronger distance decay relationships in the 3D sampling schemes compared with the 2D sampling scheme. Microbial assembly signi cantly depended on the sampling schemes with ecological processes shifting from homogeneous selection to stochasticity and variable selection when moving from 2D to 3D models. These results suggest that more dimensions should be taken into consideration when assessing microbial diversity and assembly processes.
Distance decay relationship is widely used to observe microbial biogeography and microbial community turnover rates and is known to signi cantly differ at different sampling scales [31]. The microbial community exhibited a stronger DDR in our results than large scale studies (when transferred to the same units), which is potentially a result of undersampling [17,32,33]. Undersampling causes decreased detection of low abundance taxa, leading to the community appearing to have less taxonomic turnover across space. This problem is more apparent in large scale studies due to the impracticality of sampling soil microbes to saturation. However, it should not be the case for the different DDR within our study between 2D and 3D sampling schemes as we retained the same OTU richness and abundance in all three sampling schemes. It is the case though that as the third axis was introduced into the soil cube, a greater number of ecological niches (both environmental and spatial) are likely simultaneously introduced. Thus, the increasing signi cance of DDR from s1 to s3 could have resulted from the increased spatial heterogeneity and environmental dissimilarity [21]. The decreasing community similarity from s1 to s3 also answers our rst ecological question of the effect of dimensional sampling on microbial biogeographic patterns.
Although environmental factors, spatial factors and spatial canonical axes play minor roles in shaping microbial community, microbial compositions were more correlated with spatial factors than environmental factors according to VPA. This result is consistent with Franklin and Mills [21] where they demonstrated that spatial processes play a stronger role in structuring ectomycorrhizal fungi communities than environmental factors. However, the contribution of environmental factors to bacterial community turnover was lower in our study when compared to larger scale studies [34]. This lower contribution may be due to: 1) the limited environmental variables engaged in our calculations; 2) the inherent stability of paddy soils under long term wet-dry cycling [35], resulting in less environment ltering (consistent with our result of no relationship between microbial community similarity and environmental distance after controlling for spatial distance; Tab. 1); or 3) undersampling of environmental factors, resulting in increased environmental dissimilarity and causing an overestimation of environmental effects on taxonomic turnover. Although the variation explained by spatial canonical axes was also relatively low, the shared effect with spatial or environmental factors or both was higher (up to 26%), revealing an important effect of the 3D spatial scale when assessing microbial diversity and turnover. However, even though several environmental and spatial factors signi cantly explained the variation in the RDA, a large proportion of community variation remained unaccounted. This nding is in line with previous studies demonstrating that measured factors play a minor role in shaping the microbial community [36,37]. Therefore, a major contribution of the unaccounted variation could be the underestimation of environmental factors resulting from unmeasured environmental variables such as pH, EC, oxygen content and soil mineral composition.
Ecological processes are spatially dependent and are imposed by different environmental conditions [15]. In our study, βNTI increased as we introduced the third axis, indicating that phylogenetic turnover was increasingly larger than expected from s1 to s3. The dependency on this sampling strategy may result from: 1) the use of average OTU abundance of all sub-cubes in s1 and s2, leading to a homogenization in microbial community differences; or 2) the increase in environment heterogeneity between pairwise samples as the third axis was introduced to the model. Homogeneous selection decreased as variable selection increased from s1 to s3, which can be explained by our results that β-diversity dispersion and environmental distance increased from s1 to s3 (Fig. 2 and Fig. 4). Consequently, stochasticity increased from s1 to s2 and decreased in s3, indicating that stochasticity increased when selective environment was balanced between heterogeneous and homogeneous selection. These results are consistent with the ndings of Tripathi et al [38], which show that stochastic processes are more important when environmental ltering is weak. Apart from the offset of homogeneous selection and variable selection, the increasing stochasticity form s1 to s2 may also resulted from a "Anna Karenina principle", in which the microbiological changes induced by many perturbations are stochastic, and therefore lead to transitions from stable to unstable community states [39]. The PCoA results in our study also support this assumption as microbial community dispersion signi cantly increased from s1 to s2. However, as more stress was imposed onto the system from s2 to s3, stochasticity decreased as variable selection dominated the assembly processes. These results suggest that sampling schemes (or sampling grain) should be determined according to the scale of ecological questions considered as spatial heterogeneity and distribution of soil microbial communities are scale-dependent [21,40,41].
A limitation of this study that should be considered is that the null model result is largely dependent on microbial α-diversity when quantifying the relative importance of ecological processes. The null model framework is dependent on OTU richness and abundance during the permutational calculations.
However, this dependency can be optimized by controlling the observed α-diversity constant [10]. In our study, we used the same sequencing depth for all three sampling schemes to minimize the abundance effect. Nevertheless, OTU richness was different in samples of the three sampling schemes (Fig. 2a). While OTU richness may have a considerable effect when analyzing community turnover rate taxonomically and phylogenetically, we consider that this is a reasonable outcome as OTU richness logically increases with sampling grain in natural soil systems [42]. Furthermore, as microbes are spatially distributed and microbial diversity is largely determined by habitat environmental and spatial factors [43], we regard the differences induced by OTU richness as part of the spatial and environmental effects.
It is generally accepted that both environmental ltering and spatial heterogeneity help maintain microbial assembly and composition with the caveat that the relative importance of these two factors in explaining community variation depends on the model used. In the present study, VPA indicated that spatial factors explained a larger proportion of community turnover, while the null model demonstrated that community turnover was mostly due to environmental ltering. These contrasting results may be due to the unmeasured environmental variables as mentioned above. Thus, although no signi cant correlation between community turnover and environmental factors was detected, we hypothesise that enviromental factors still play an important role in shaping microbial composition. Further, the underestimation of environmental effects can also lead to an overestimated importance of spatial factors and undominated processes [44], which can also explain the different results of VPA and the null model. Finally, while VPA is mainly used for partitioning β-diversity of microbial communities in a region [45], the βNTI null model approach is designed for microbial communities that are mainly in uenced by selection in an optimal organism environment [38]. In this case, the null model approach may also underestimate spatial processes when community composition is strongly in uenced by dispersal based community assembly. Thus, new models that integrate all relevant data (taxonomic, phylogenetic, environmental and spatial) are still needed, and more dimensions (for example time, a fourth dimension) should be taken into consideration when relating the relative in uence of ecological processes to community assembly.

Conclusions
This study showed microbial diversity and assembly processes dynamically varied with 2D and 3D models at a sub centimeter scale. Although deterministic processes dominated the community assembly in both 2D and 3D sampling schemes, the importance of homogenous selection to variable selection changed. These results indicated that microbial community should be estimated in more dimentions to improve our understanding of microbial biodiversity and biogeographic patterns. However, unmeasured enviromental variables and the different emphasis between VPA and null model approaches lead to a contrasting explaination of community turnover. Thus, we suggest the use of integrated models when deciphering microbial diversity and processes in multiple dimensions.

Declarations
Ethics approval and consent to participate Not applicable.

Consent for publication
Not applicable.

Availability of data and material
The raw reads of 16S rRNA sequencing have been deposited in the NCBI Sequence Read Archive (SRA) database (accession number SUB7287031).

Competing interests
The authors declare that they have no competing interests.   Bray-Curtis similarity of microbial community against spatial distance and sampling schemes from three orientations a) Distance-decay curves showing the relationship between Bray-Curtis similarity of microbial community and spatial distance from three orientations. Linear trend is the ordinary leastsquares linear regressions; shaded area indicates con dence interval (0.95); b) Patterns between community similarity and sampling schemes from three orientations. Environmental distance against spatial distance and sampling schemes from three orientations a) Relationship between environmental distance and spatial distance from three orientations. Linear trend is the ordinary least-squares linear regressions; shaded area indicates con dence interval (0.95); b) Patterns between environmental distance and sampling schemes from three orientations.
Page 18/19 Figure 5 Variation partitioning of microbial abundance and diversity into soil environment (magenta), spatial distance (green) and spatial canonical axes (blue) components from three orientations. Variation (%) explained by the interactions of two or all factors are indicated on the sides or middle of the triangles, respectively. Black rectangles below indicate residual variation (%).