Sequence Analysis of 16S rRNA Gene Amplicons

Sequences of 16S rRNA gene amplicons were analyzed using the Mothur 1.34 following the standard protocols (Schloss et al., 2009). Barcoded sequences were de-multiplexed and filtered to remove low quality reads (Phread score < 25). Reads that had more than two mismatches to the paired barcode sequences were removed by the noise reduction while retaining the information of the representative reads and the numbers of reads merged by pre-clustering (Huse et al., 2010). The unique reads were aligned to the Silva NR119 database¹. Reads not aligned in the same region were removed. The sequence regions beyond the primers were truncated. Potential chimeric sequences were detected and removed using the UCHIME program (Edgar et al., 2011). The number of sequences in each sample after quality filtering is shown in Supplementary Table S1.

The taxonomy of each unique sequence was assigned using the Silva SSU dataset of the NR 119 release as the reference. Taxonomy assignments with bootstrap values greater than 80% were considered to be valid. Sequences sharing more than 97% identity were further clustered into individual operational taxonomy units (OTUs) using the nearest neighbor algorithm (Schloss and Westcott, 2011). Based on the rarefied datasets (n = 10,008), alpha diversity indices, such as the number of observed OTUs, Chao-1, and Inverse Simpson (Hill, 1973; Chao et al., 1984; Faith, 1992), were computed.

The weighted Unifrac (Lozupone and Knight, 2007), which quantitatively incorporates the relatedness of community members with their abundances for the evaluation of difference between samples, was used to determine the community dissimilarity between samples from the sediment core. In addition, analyses of non-metric multidimensional scaling (NMDS) was conducted on the basis of the weighted Unifrac distance matrix using the software Mothur 1.34 (Schloss et al., 2009). The significance of environmental variables relative to the NMDS ordinations was computed using “envfit” and 999 permutations. Finally, the dissimilarity matrix between samples from the sediment core and incubations was constructed by the Bray–Curtis method (Bray and Curtis, 1957) and visualized by principal coordinate analysis (PCoA) using the R “Phyloseq” package.