Analysis of publicly available 16S rRNA gene sequence datasets for taxonomic profiling

To summarize crustal bacterial and archaeal taxa for this review, we synthesized publicly-available sequence datasets of the 16S rRNA gene in environmental DNA extracted from seafloor and subseafloor basalts generated using 454, Illumina and Ion Torrent amplicon platforms. These include seafloor basalts from the Dorado Outcrop (Lee et al., 2015) and the Lō'ihi Seamount (Jacobsen Meyers et al., 2014) in the Pacific Ocean and subseafloor basalts from North Pond on the western flank of the Mid-Atlantic Ridge (Jørgensen & Zhao, 2016) and the Juan de Fuca Ridge flank in the northeastern Pacific Ocean (LaBonté et al., 2017). Datasets from rock colonization experiments conducted in the subseafloor at the Juan de Fuca Ridge flank site (Smith et al., 2016; Ramírez et al., 2019) were also included, as well as microbial community surveys of the subseafloor crustal fluids from the anoxic Juan de Fuca site (Jungbluth et al., 2016) and the oxic North Pond site (Tully et al., 2017; Meyer et al., 2016). For comparison, we included select reference datasets from oxic (Reese et al., 2018; Zinke et al., 2018) and anoxic sediment (LaBonté et al., 2017) and the overlying bottom seawater (Lee et al., 2015) from these same study sites.

Raw sequence data from the reviewed studies were downloaded from the NCBI Short Read Archive. Sequencing reads generated using Illumina and Ion Torrent platforms were quality filtered and processed to unique Amplicon Sequence Variants (ASVs) using DADA2 (Callahan et al, 2016), with taxonomy determined by the naïve Bayesian classifier in DADA2 using a training set from the SILVA v132 database (Quast et al., 2013; Yilmaz et al., 2014; Glöckner et al., 2017). For the 454 GS-FLX sequence datasets, operational taxonomic units (OTUs) constructed with 97% or greater sequence similarity in the original analyses were reprocessed in mothur V.1.37.6 (Schloss et al., 2009) against the same SILVA database. All short read datasets were merged and summarized to the relative abundance at phylum resolution (or to class level for Proteobacteria phyla) using Phyloseq v1.24.0 (McMurdie & Holmes, 2013). Figures were produced using ggplot2 R package version 2.2.1 (Wickham, 2016) in RStudio (RStudio Team, 2017). Taxonomic grouping in each sample separated taxa into common (>5% abundance in at least one sample) versus rare (never more than 5% in any sample). Supplemental Figure S1 shows the breakdown of Gammaproteobacteria families in the samples presented in Figure 4 of the main text, and Supplemental Figure S2 highlights the abundance of rare taxa (never >5% abundance in any sample). The Bray-Curtis distances between samples was calculated using the same dataset described above, summarized to relative abundance at the Family taxonomic level using Phyloseq and the Vegan package (Oksanen et al., 2018). A Non-Metric Multidimensional Scaling (NMDS) ordination was produced from this distance matrix. It should be noted that common rules for beta diversity comparisons, such as common library preparation/sequencing protocols and library-size normalization, were not performed in this analysis due to the diversity of the datasets being considered and the resulting NMDS ordination having high-stress (>20%). Therefore, the results should be viewed as broadly qualitative and not quantitative.

All data processing steps and markdown files are available via github: https://github.com/orcuttlab/ocean-crust-micro