Study area description and field soil sampling
The present study was conducted on the south coast of Hangzhou Bay (121°05′-121°35′′E, 30°0′-30°25′N) in Cixi County, Zhejiang Province, China. The study area is a typical northern subtropical monsoon climate zone, with a mean annual precipitation of 1273 mm and temperature of 16.0°C. Salt marshes formed from the Qiantang Estuary saltwater sediments and diked lands were exposed to tidal effects at an elevation range from 2.6 to 5.7 m [23]. The tidal salty marsh was eventually reclaimed for crop growth, beginning 1000 years ago with the formation of dikes at different historical stages. Based on the construction times of dikes, the years of soil reclamation were calculated, and a detailed description is available online in Chinese at (www.cixi.gov.cn). In this study, around 1000 years of ancient spanning of soil chronosequence was identified, including undisturbed coastal salty marsh (T0) and soil reclaimed 5-years (T5), 20-years (T20), 50-years (T50), 60-years (T60), 120-years (T120), 200-years (T200), 220-years (T220), 280-years (T280), 500-years (T500), and 1000-years (T1000) ago.
Distances between these soils were no more han 30 km in a similar topography. Soil samples were collected and mixed from six random locations, and then six mixed soils sampled from the same reclaimed year were considered replicates. A stainless-steel auger (3-cm diameter) was used for soil collection. Briefly, soil samples were taken from the plough layer (0–20 cm), and a total of 66 samples were obtained representing specific reclamation years (Fig. S1). The taxonomic class of soil samples was fluvisol as per Food and Agriculture Organization/United Nations Educational, Scientific and Cultural Organization (FAO/UNESCO), which was formed from fluvial deposits around the Qiantang River, and the soil texture was silty loam to clay-loam [24, 25]. In December 2015, reclamation site samples (T5-T1000) were collected from a traditional soybean-broccoli rotation field during active broccoli growth from crop free sites to avoid the possible impact of plant roots. Compound fertilisers were applied @ 900 kg·ha-1, of which nitrogen (N):phosphorus (P2O5):potassium (K2O) was estimated to be 1:1:1. The salt marsh soils (T0) were sampled in the high-tide area under Phragmites australis. Samples were immediately taken to the laboratory on ice packs and then sieved through a 2 mm mesh to remove roots and another residue. After that, sieved soils were divided into three parts; two parts were kept at -80°C and 4°C for total DNA extraction, mineralised N content, and nitrification analysis. The third part was air-dried for the determination of physicochemical parameters of the soil.
Total DNA extraction
The soil total genomic DNA was extracted using the Soil DNA kit (Omega Bio-Tek, Inc., Norcross, GA, USA) according to the manufacturer's instructions from 0.5 g of soil. The purity and quality of DNA were measured by 1% agarose gel electrophoresis and a spectrophotometer, respectively (Nanodrop, PeqLab, Germany).
Amplification and sequencing of the amoA gene for AOA and AOB
The amoA gene of archaea and bacteria was amplified by polymerase chain reaction (PCR) using Arch-amoAF (5'- STAATGGTCTGGCTTAGACG-3'), Arch-amoAR (5'-GCGGCCA TCCATCTGTATGT-3') and amoA-1F (5′-GGGGTTTCTACTGGTGGT -3′) and amoA-2R (5′- CCCCTCKGSAAAGCCTTCTTC -3′) respectively. PCR was performed with 15 μL of Phusion® High-Fidelity PCR Master Mix (New England Biolabs), 0.2 μM/μL forward and reverse primers, and 10 ng/μL DNA template. The following cycling parameter were used: incubation at 95℃ for 5 min, followed by 30 cycles of 95℃ for 30 s, 56°C and 58°C for 60 s, 72°C for 60 s and 35 s (archaea and bacteria, respectively), and a final 10 min extension at 72℃. The PCR product was recovered and purified by 2% agarose gel electrophoresis and GeneJET Gel Extraction kit (Thermo Fisher Scientific USA). According to the manufacturer's instructions, sequencing libraries were constructed using the Ion Plus Fragment Library Kit 48 rxns (Thermo Scientific). The library quality was estimated with a [email protected] 2.0 Fluorometer (Thermo Scientific), and the library was sequenced on an Ion S5TM XL platform.
Quantification of the amoA gene in AOA and AOB
To estimate the amoA functional gene for AOA and AOB, quantitative PCR (qPCR) was performed. The ammonia-oxidising archaeal amoA gene was quantified by fluorescence qPCR using the primers shown above (section 2.3). Absolute quantification was carried out on the Applied Biosystems QuantStudioTM 6 Flex Real-Time PCR System (Life Technologies Corporation, Carls-bad, CA, USA) amplifier using Hieff® qPCR SYBE® Green Master Mix (YEASEN). The fluorescence qPCR reaction volume was 20 μl, containing 10 μl 2 × SYBR Green Mix, 0.8 μl upstream and downstream primers (10 pmol / μl), 1 μl diluted total DNA template, and 7.4 μl double distilled water. The protocol for the PCR reaction was as follows: pre-denaturation at 95°C for 5 min; 40 cycles at 95°C for 30 s, annealing at 56°C and 58°C for 60 s, and extension at 72°C for 60 s and 35 s (archaea and bacteria, respectively); the melting curve temperature range was 65–95°C.
Sequencing data processing and analysis
Single-end reads were assigned to samples based on their unique barcode, and Cutadapt (V1.9.1) was used for quality filtering to obtain high-quality clean reads [26]. Reads were analysed by the Silva database using the UCHIME algorithm to classify and delete chimera sequences [27, 28]. The Uparse program grouped the purified sequences (Uparse v7.0.1001). The same operational taxonomic units (OTUs) were assigned sequences with 97% identity [29]. The Mothur algorithm [30] was used to perform annotated OTU analysis to obtain representative species and count the number of OTUs per sample in taxonomic information (Kingdom, Phylum, Class, Order, Family, Genus, and Species level). Archaeal and bacterial OTUs were identified and assigned in the archaeal and bacterial Silva Database [27]. For multiple sequence alignment, MUSCLE software (Version 3.8.31) was used to examine the phylogenetic relationship between different OTUs and the classification of the dominant species in different samples [31].
Phylogenetic diversity matrices analysis
The Mean Phylogenetic Distance (MPD) determined by Net Relatedness Index (NRI) for all co-occurring populations revealed ‘basal dispersal’ within the population, while the Net Taxon Index (NTI) measures the Mean Nearest Taxon Distance (MNTD) between populations, thus calculating the population 'terminal' phylogenetic dispersion [32]. The MPD and MNTD standardised effect size was measured using the ses.MPD and ses.MNTD commands in the ‘Picante’ package of R [33]. The distribution of 999 null values, computed by shuffling the tip labels in the tree, was used to account for temporal changes of the species pool. ses.MPD. observed ses.MNTD. The observed values were multiplied by -1 to be equivalent to the NRI and NTI [33].
NRI= -(MPDobserved–MPDrandomised)/(sdMPDrandomised) [32]
NTI = -(MNTDObserved − MNTDRandomised)/ (sdMNTDRandomised) [33]
The beta Means Nearest Taxon Distance (βMNTD) was computed as previously described [34]. The beta Nearest Taxon Index (βNTI) and standard deviation of βMNTD were calculated in R with the package MicEco by using command ses.comdistnt (abundance.weighted= True) [35]. The significance was calculated using 9999 permutations with the ‘ade4’ package in R. Overall environmental and geographical differences were calculated as standardised Euclidean distances between samples.
Statistical analyses
Microbial diversity and a comparison of gene copy numbers were evaluated by one-way analysis of variance (ANOVA), followed by Duncan's multiple range test at SPSS 25.0 for Windows. Discrepancies were considered statistically significant at P < 0.05. Soil ecosystem function (potential nitrification and nitrogen mineralisation rates) correlation with environmental parameters and amoA gene abundance, as well as reclaimed soil age, were calculated by Spearman’s correlation analysis. Alpha diversity analysis namely observed-species, Chao1, Shannon, Simpson, ACE, and good-coverage, were calculated with QIIME (Version1.7.0). Beta diversity analysis was performed on QIIME software (Version1.7.0) to evaluate the differences in community complexity among samples and to calculate the Unique Fraction (UniFrac) distance. The Unweighted pair group method with arithmetic mean (UPGMA) sample clustering was conducted in R software (Version 2.15.3). NMDS (non-metric multi-dimensional scaling) analysis and non-parametric Analysis of Similarity (ANOSIM) based on the Bray-Curtis distance for analysing compositional differences between ammonia oxidiser prokaryotic communities were performed in R software (Version 2.15.3). Differential abundance of ammonia oxidiser between individual samples was calculated by MetaStat analysis [36].
Temporal variations in phylogenetic beta diversity were analysed by Permutational Multivariate Analysis of Variation (PERMANOVA) using R function ADONIS based on the Bray-Curtis distance. The effect of physical and chemical factors on ammonia-oxidising microbial communities was analysed by canonical correspondence analysis (CCA). The cca-envfit function was used to test the impact of each environmental factor on the distribution of species. Spearman rank correlation coefficient analysis was used to correlate the relationship of environmental factors and reclaimed soil stages with ammonia oxidiser microbial abundance (alpha diversity) and species distribution. We conducted correlation analyses to assess the relationship of environmental variables with NTI and NRI to examine the change in community structure by Spearman’s correlation. The mechanisms of ammonia oxidiser community assembly were determined by calculating βNTI. Briefly, if βNTI > 2 or βNTI < -2, deterministic processes may be important in shaping the community composition across all sites, whereas stochastic processes may play a significant role in community assembly processes when the values of βNTI are between -2 and 2. The βNTI > 2 revealed significantly more phylogenetic turnover than predicted, which is often interpreted by chance as variable selections, while βNTI < 2 referred to less phylogenetic turnover than expected, i.e., homogeneous selection (Stegen et al. 2013). If the measured βMNTD values did not bring clarity to significant differences from the null distribution of βMNTD, e.g. |βNTI|<2, the observed phylogenetic variability was not the consequence of selection [37]. To overcome this problem, the Bray-Curtis-based Raup-Crick metric (RCbray) was further determined as described by [38]. Briefly, the values of (RCbray) ranged between -1 and 1, and we compared |βNTI|<2 and (RCbray) values. The relative contribution of dispersal limitation was estimated as the percentage of pairwise comparison between |βNTI| < 2 and (RCbray) values > 0.95, whereas |βNTI| < 2 and (RCbray) values > -0.95 indicated homogenous dispersal. The undominated process was calculated as |βNTI| < 2 and (RCbray) values > 0.95. The undominated concept described a state wherein the primary cause of variations between population compositions was neither dispersion nor selection, namely ecological drift (population sizes fluctuating due to stochastic birth and death events) [39]. We applied ordinary least-squares regression analysis to determine the slope of the association among phylogenetic relatedness with environmental factors. The Mantel tests assessed the association of phylogenetic distance and environmental parameters with geographic distance. Using partial Mantel tests with Pearson's correlation coefficient and 999 permutations, the relationships between βNTI, geographical distance, and Euclidean distances in environmental parameters were analysed.