Soil sampling and property analyses
The soil was sampled in November 2017 (autumn) from a long-term rice cultivation field 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 classified 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 cm3 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 first 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]. Briefly, 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 filtered 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/). Quantification 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 amplified by primers 515F (5′-GTGCCAGCMGCCGCGGTAA-3′) and 907R (5′-CCGTCAATTCCTTTGAGTTT-3′). Amplified 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 classified up to phylum level (Additional file: Figure S1). RDP training set v16 was used for taxonomy annotations with a confidence cut-off of 80% [28] after which aligned sequences of representative OTUs were used to construct a maximum-likelihood tree. The community matrix was rarefied 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 file: 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 first 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 R2 coefficients 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 first 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, significant 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 significant variables, Additional file: 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 (RCbray) 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 RCbray >+0.95 while the relative contribution of homogenizing dispersal was estimated as the percentage of pairwise comparisons with |βNTI| <2 and RCbray <−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 |RCbray| <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.