A compendium of multi-omic sequence information from the Saanich Inlet water column

Marine oxygen minimum zones (OMZs) are widespread regions of the ocean that are currently expanding due to global warming. While inhospitable to most metazoans, OMZs are hotspots for microbial mediated biogeochemical cycling of carbon, nitrogen and sulphur, contributing disproportionately to marine nitrogen loss and climate active trace gas production. Our current understanding of microbial community responses to OMZ expansion is limited by a lack of time-resolved data sets linking multi-omic sequence information (DNA, RNA, protein) to geochemical parameters and process rates. Here, we present six years of time-resolved multi-omic observations in Saanich Inlet, a seasonally anoxic fjord on the coast of Vancouver Island, British Columbia, Canada that undergoes recurring changes in water column oxygenation status. This compendium provides a unique multi-omic framework for studying microbial community responses to ocean deoxygenation along defined geochemical gradients in OMZ waters.


Background & Summary
Marine oxygen minimum zones (OMZs), areas of low dissolved oxygen (O 2 ) in subsurface waters, result from a combination of microbial respiration of organic matter raining down from surface waters and increased water column stratification [1][2][3] . As O 2 becomes limiting, microbial communities shift their energy metabolism to use alternative terminal electron receptors in a thermodynamically defined order resulting in increased nitrogen loss and the production of climate active trace gases including nitrous oxide (N 2 O) and methane (CH 4 ) [4][5][6][7][8][9] . Currently, OMZs are expanding throughout the global ocean 3,10-14 making it increasingly important to define the microbial metabolic networks driving nutrient and energy cycling under changing levels of water column O 2 -deficiency 2,10,15 .
Advances in sequencing technology are enabling the study of microbial communities at unprecedented scales 16 . Tag sequencing uses primers to amplify specific target genes and subsequently sequence them on high-throughput platforms, generating molecular barcodes that can be used to study microbial community structure and function depending on the marker. The small subunit ribosomal RNA (SSU rRNA) gene is a universally conserved marker commonly used to compare microbial communities among and between samples in a quantitative manner 17 . Metagenomic sequencing enables reconstruction of microbial community metabolic potential at the level of genes and pathways over relatively long time scales reflecting persistent information storage in the environment. However, metabolic potential does not necessarily indicate active processes as different genes may be expressed under changing environmental conditions 18 . Metatranscriptomic sequencing opens a window into microbial community gene expression patterns on relatively short time scales reflecting environmental response patterns 18 . Similarly, metaproteomics identifies proteins present in the microbial community on intermediate times scales, providing an alternative perspective on post-translational regulation and catalysis 18 . Collectively, multi-omic methods (DNA, RNA and protein) can be used to chart microbial metabolism at individual, population and community levels along defined geochemical gradients in OMZs [19][20][21][22][23][24][25][26] .
Saanich Inlet, a seasonally anoxic fjord on the coast of Vancouver Island, British Columbia is a model ecosystem for studying microbial community responses to changing levels of water column O 2deficiency 8,[20][21][22][27][28][29][30][31][32][33] . As microbial communities shift their energy metabolism to use alternative electron acceptors within the Saanich Inlet water column, differential modes of metabolic coupling can be observed including a modular denitrification pathway coupled to sulphide oxidation 20,21,28 . Such metabolic coupling is reminiscent of symbiotic associations and is likely widespread in OMZ ecosystems 2 . Although current research efforts are increasingly focused on defining co-metabolic innovations among and between ubiquitous OMZ microorganisms, many open questions remain regarding the regulatory and ecological dynamics modulating microbial community responses to OMZ expansion 2,5,[20][21][22]34 .
Here we present a compendium of multi-omic sequence information from the Saanich Inlet water column ( Fig. 1) encompassing 412 SSU rRNA pyrotag (V6-V8 region) samples (Data Citation 1), 82 SSU rRNA iTag (V4 region) samples (Data Citation 1) (Table 1), 90 metagenomes (Data Citation 1) (totalling 4.1 TB of cleaned reads or 16.2 GB of assembled data), 62 metatranscriptomes (Data Citation 1) (including 46 unique samples and 16 replicates, totalling 1.7 TB of cleaned reads or 2.88 GB of assembled data) and 68 metaproteomes (64 unique samples, totalling 5.2 million unique proteins) (Data Citation 2) ( Table 2). Together sequence read data is approximately 5.9 TB of data, comparable to nearly 2,000 human genome equivalents. These data sets, in combination with a cognate geochemical compendium 35 comprise a unique time-resolved framework for reconstructing microbial community metabolism along defined geochemical gradients and promoting the development of models to predict microbial community responses to changing levels of water column oxygen-deficiency 11,22,36 .
During sampling procedure, large volume waters were collected in 2 × 12 l Go-Flow bottles on a wire separated by less than one meter. Waters for metatranscriptomics and metaproteomics (2 l each) were collected consecutively from the Go-Flow into 2 l Nalgene bottles with sterile silicon tubing immediately following sampling for dissolved gases to minimise changes in microbial gene expression. Waters for metatranscriptomics were filtered on-board within 8 min of collection on deck, followed immediately by filtering waters for metaproteomics. For both metatranscriptomics and metaproteomics, a peristaltic pump was used to filter waters through a 0.22 μm Sterivex filter with an in-line 2.7 μm GDF pre-filter. Following removal of residual seawater by extrusion using a 10 cc or 60 cc syringe, 1.8 ml of RNAlater (Ambion) was added to metatranscriptomic sample filters and 1.8 ml of sucrose lysis buffer was added to metaproteomic sample filters. Filters were placed on dry ice, returned to lab and stored at −80°C until processing. A detailed small volume filtration protocol can be found online at http://www.jove.com/ video/1163/small-volume-1-3l-filtration-of-coastal-seawater-samples 37 . Waters for LV SSU rRNA gene tags and metagenomics were combined from two Go-Flows into a 20 l carboy using sterile silicon tubing. Carboys were then transported back to the lab for filtration approximately 6-10 h after collection. Approximately 10 l was filtered through 2 Sterivex filters per depth as described above. Following removal of residual seawater by extrusion, 1.8 ml of sucrose lysis buffer was added to LV sample filters. Filters were placed on dry ice, returned to lab and stored at −80°C until processing. A detailed large volume filtration protocol can be found online at http://www.jove.com/video/ 1161/large-volume-20l-filtration-of-coastal-seawater-samples 38 .
Waters for high resolution (HR) SSU rRNA gene tags were collected in 12 l Go-Flow or 8 l Niskin bottles and transferred into 1 l (February 2006-September 2009) or 2 l (September 2009 onward) Nalgene bottle with sterile silicon tubing and stored on ice. Bottles were then transported back to the lab for filtration approximately 12-16 h after collection. A peristaltic pump was used to filter waters through a 0.22 μm Sterivex filter without a pre-filter. Following removal of residual seawater by extrusion, 1.8 ml of sucrose lysis buffer was added to HV sample filters. Filters were placed on dry ice, returned to lab and stored at −80°C until processing. A detailed small volume filtration protocol can be found online at http://www.jove.com/video/1163/small-volume-1-3l-filtration-of-coastal-seawater-samples 37 .

Environmental DNA, RNA and protein extraction
Genomic DNA (HR and LV) was extracted from Sterivex filters as described in 27 . Video protocols describing the extraction process in detail can be found online at http://www.jove.com/video/1352/dnaextraction-from-022-m-sterivex-filters-cesium-chloride-density 39 . Briefly, after thawing Sterivex filters on ice, lysozyme (Sigma) was added and incubated at 37°C for 1 h with rotation followed by addition of Proteinase K (Sigma) and 20% SDS and incubation at 55°C for 2 h with rotation. Lysate was removed by extrusion and then filters were rinsed with sucrose lysis buffer. Combined lysate was extracted with phenol:chloroform followed by chloroform and the aqueous layer collected and loaded onto a 10 K Amicon filter (Millipore) cartridge, washed three times with TE buffer (pH 8.0) and concentrated to a final volume between 150-400 μl by centrifugation.
Total RNA was extracted from Sterivex filters using the mirVana miRNA Isolation kit (Ambion) 19,25 modified for Sterivex filters. Briefly, after thawing filters on ice, RNAlater was removed by extrusion, followed by rinsing with Ringer's solution (Sigma) and incubation at room temperature for 20 min with rotation. Ringer's solution was removed by extrusion, followed by addition of lysozyme and incubation at 37°C for 30 min with rotation. Lysate was removed by extrusion into 15 ml tube and subjected to organic extraction as described in the mirVana kit protocol, adjusting for total volume of lysate. Removal of DNA and purification of total RNA were conducted using the TURBO DNA-free kit (Thermofisher) and the RNeasy MinElute Cleanup kit (Qiagen) protocols respectively.
Total protein was extracted from Sterivex as described in Hawley et al. 40 . Briefly, after thawing Sterivex filters on ice, Bugbuster (Novagen) was added and incubated at room temperature for 20-30 min with rotation. Lysate was removed by extrusion and filters were rinsed with 1 ml lysis buffer. Buffer exchange was carried out on combined lysate using Amicon Ultra 10 K (Millipore) with 100 mM NH 4 HCO 3 a total of three times with a final volume between 200-500 μl. Protein concentration was determined using the 2-(4-carboxyquinolin-2-yl) quinoline-4-carboxylic acid (Bicinchoninic acid or BCA) assay. Urea was added to a final concentration of 8 M and dithiothreitol added to a final concentration of 5 mM and incubated at 60°C for 30 min, followed by 10-fold dilution with 100 mM NH 4 HCO 3 . Samples were then subject to trypsin digest at 37°C for 6 h followed by C18 solid phase extraction and strong cation exchange.

Environmental DNA and RNA sequencing
Extracted genomic DNA was used to generate small subunit ribosomal RNA (SSU or 16S/18S rRNA) gene pyrotag libraries 20,41 . Pyrotag datasets from HR and LV samples were generated by PCR amplification using universal three-domain forward and reverse bar-coded primers targeting the V6-V8 hypervariable region of the SSU rRNA gene 42 : 926F (5′-AAA CTY AAA KGA ATT GRC GG-3′) and 1392R (5′-ACG GGC GGT GTG TRC-3′). Samples were purified using the QIAquick PCR Purification Kit (Qiagen), and sequenced by 454-pyrosequencing 43 at the DOE Joint Genome Institute (JGI), or Génome Québec Innovation Centre at McGill University (Table 3 and Supplementary Table 1). iTag datasets from HR and LV samples were generated by PCR amplification using forward and reverse barcoded primers targeting the V4-V5 hypervariable region of the bacterial and archaeal SSU rRNA gene: 515F (5′-Y GTG YCA GCM GCC GCG GTAA-3′) and 806R (5′-CCG YCA ATT YMT TTR AGT TT -3′) 44,45 . Samples were sequenced according to the standard operating protocol on an Illumina MiSeq platform at the JGI 46 .
Extracted environmental total RNA was used to generate paired end sequenced Illumina metatranscriptome libraries 53 at the JGI on the HiSeq and (Table 4 and Supplementary Table 2).

Environmental protein sequencing
Extracted environmental protein was sequenced using tandem mass spectrometry (MS/MS) as described previously 40 . Samples were analysed by capillary liquid chromatography-tandem mass spectrometry (Thermo, LTQ ion trap mass spectrometer or Thermo LTQ-Orbitrap mass spectrometer) in datadependent mode. Spectra were matched to peptide sequences using the search tool MSGFDBPlus 54  amino acid sequence database used for matching spectra to protein sequences was constructed from metagenomic information from Saanich Inlet and the Line P transect in the Northeastern subarctic Pacific Ocean with additional full-length fosmid libraries 20 , and single cell genomes from Saanich Inlet 55 totalling over 23 million protein sequences.

Data Records
A Table unifying all multi-omic samples and Saanich Inlet geochemical samples from Torres-Beltrán et al. 35 is summarised in Table 2 and detailed in Supplementary Table 2.
Small subunit rRNA gene tag sequences Small subunit rRNA gene pyrotag and iTag sequences are available from the National Center for Biotechnology Information (NCBI) Sequence Read Archive (Data Citation 1). A summary of Pyrotag and iTag samples can be located in Table 1, key to the Pyrotag and iTag data files is located in Table 3

Metaproteome
Metaproteomic datasets have been deposited in the ProteomeXchange Consortium via the PRIDE 56 partner repository with the dataset identifier PXD004433 (Data Citation 2). The key detailing the file types available on the PRIDE repository is located in Table 5. Previously published metaproteomic data sets 20 were archived at Mass Spectrometry Interactive Virtual Environment.

Technical Validation
Small subunit ribosomal RNA gene tag sequences For SSU rRNA gene pyrotags and iTags the quality of extracted genomic DNA was verified on 0.8% agarose gels stained with ethidium bromide or SYBR Safe DNA dye (Thermofisher). Samples were run at 16 V overnight with molecular ladders (10 μl of 50 ng ml −1 of 1 kb+; 2, 5 and 10 μl of 50 ng ml −1 HindIII ladder) to determine size and estimate concentration of extracted sample DNA (5 μl of extract run on gel). After running, the gel was observed under a UV gel documentation system to check for approximate molecular weight (>36 kb) and evidence of shearing or degradation. In addition, DNA concentration was quantified and corrected to volume filtered using PicoGreen (Thermofisher) following the vendor's protocol.

Data field Description
Sample ID Identifier of unique time-series time point and depth in which seawater sample for dataset was obtained, links to geochemical time series data (Torres Beltrán et al.)

Cruise ID Numerical identifier of individual cruises
Year Year of cruise  Quality control for 454-pyrosequencing entails accurate quantification of prepared library fragments to optimize the sequencing run output, therefore JGI has developed custom qPCR methods to quantify 454 libraries 43 . In addition, an optimized emulsion PCR protocol was developed to significantly improve the coverage in high GC regions that otherwise would be biased 49 . Génome Québec Innovation Centre quality control protocol for 454 pyrosequencing entails the use of the default parameters assigned to the signal processing software for Long Amplicon #3 pipeline as indicated in the Roche '454 Sequencing System Software Manual' V3.0 Section 1.3 53 . For produced 454 pyrotag datasets for both high resolution (HR) and large volume (LV) samples a histogram of raw read counts versus read length (Fig. 2a) was used to determine the success of a run e.g. a majority of reads exceeding 450 base pairs. A plot of read counts versus read length for all HR and LV samples is provided in Fig. 2a.
Quality control protocol for iTag sequencing entails QC amplification with V4-V5 primers prior to sample submission to ensure that samples will be successful in JGI sample prep and that there are no contaminants present that will inhibit amplification. Amplified products are run on agarose gel or Bioanalyzer to ensure sample amplification occurred with the expected band size,~450 bp for bacterial and archaeal SSU rRNA gene V4-V5 region (including primers). Samples showing proper amplification and sizing are passed for Amplification QC and approved to ship to the JGI for processing. All data from the sequencer were demultiplexed and stored in JGI's archiving and metadata organizer system (JAMO). Read data was processed through JGI's centralized rolling quality control system verifying that there were no sequencing issues and removing known contaminant reads using a kmer filter in bbduk. Quality controlled reads were processed by iTagger 45,46 .

Metagenomic and metatranscriptomic data validation
The quality of DNA used for metagenomic sequencing was determined at the same time as verification for SSU rRNA gene pyrotag and iTag amplification as described above. JGI quality control protocol for

Cruise ID Numerical identifier of individual cruises
Year Year of cruise

Station
Indicates the sampling station from which seawater sample was obtained

Depth
Depth at which seawater sample was obtained

Tag Data
Indicates if SSU rRNA tag data exists for that sample and what type of tag (see Table 4)

Search Files
Parameter and settings files used for the database search of spectra to peptide

Peak Files
De-isotoped values of mass, observed charged states, and chromatographic elution times from the mass spectrometry runs

RAW Files
Mass spectrometry run files, in original format

FASTA Files
Amino Acid sequence file for all detected proteins from all Saanich Inlet metaproteome samples.

SBI_Metagenome2015_AllPeptides
Tabular lists of identified peptides, associated confidence scores, and protein reference names

AllProteinsAllExperiments
Protein lists from all samples, including redundant peptide to protein matches  metagenomic sequences prior to assembly entails rolling QC, an in-house sequence QC pipeline that performs a set collection of analyses and produces a summary report for each lane of Illumina data produced by the sequencing group. This set of analyses calculates read quality, measures sequence uniqueness, and detects abnormal sequence motifs. An assembly, using Velvet, was used to measure coverage and detect contamination 48 . For individual sample assemblies the average fold coverage versus the contig length (Fig. 2b) was plotted and should have a distinct shape for different samples where peaks in contig length representing at a specific coverage represent a given closely related microbial population. Additionally, the percent GC versus average coverage can be plotted, again with distinct shapes for different samples and clusters representing divergent microbial populations. The quality of purified RNA was verified on the Bioanalyzer using a RNA nano Analysis Kit (Agilent Technologies) in order to check on the RNA integrity and sample quantification before cDNA library production and sequencing at JGI. Due to variation inherent in environmental samples the RNA integrity number (RIN) varied between 5.5 and 9 for the sample sand averaged 7.3. JGI quality control protocol for metatranscriptomic sequencing preparation follows the 'TruSeq Stranded Total RNA Sample Preparation Guide' (Illumina). Briefly this protocol removed ribosomal rRNA with RiboZero, followed by RNA fragmentation for first strand cDNA synthesis. This was followed by second strand synthesis and subsequent ligation of adapters. After PCR amplification library quality was checked using Bioanalyzer for fragment size (260 bp) and purity. Indexed (barcoded) libraries were normalized to 10 nM and pooled in equal volumes. Transcriptomes were assembled de novo and or mapped to a corresponding metagenome. For Additional quality assessment of sequencing run for each sample, histograms of percentage of reads verse average read quality and of reads per percent GC were generated (Fig. 3a). Information on run quality for individual samples is available via the JGI genome portal.

Metaproteomic data validation
The quality of LC-MS runs were monitored visually by the instrument operator by individually inspecting each analytical run for instrument response, retention time characteristics, and background interference. QC was further monitored using a whole cell digest of Shewanella oneidensis that is prepared in bulk and routinely used across all LC-MS systems in the EMSL production labs. Analysis of this QC standard is fully automated once uploaded to an in-house data management system and reports back unique peptide identifications, chromatographic peak width, and mass error within~1 h. Additionally, these data are subjected to dozens of other analytical metrics that can be further assessed as needed (i.e., when the first pass visual and 1 h results are in question). This QC standard was run at least once per week but generally more often, between sample batches (a single project), and sometimes between sample blocks (when a batch is relatively large and a blocking scheme has been utilized) due to the large diversity of samples analysed in-house. New LC columns were conditioned and tested prior to use on project samples by running a minimum of three QC standards. Blanks were always run between QC standards and between sample batches, but not necessarily between sample blocks. For peptide mapping to full length protein sequences the False Discovery Rate (FDR) was calculated using the spectra to peptide matches that resulted in reversed hits from the on-the-fly reversed database search and a filter on the MSGF value. Number of peptides and proteins detected varies between samples (Fig. 3b). Due to the large size of metagenomic dataset used and redundancy in protein sequences because of multiple sampling of the same environment in the Saanich Inlet time series most peptides mapped to multiple identical proteins, resulting in a greater number of proteins identified then peptides.

Usage Notes Suggested modes of downstream data analysis
In brief we describe the main workflow used in the analysis of SSU rRNA gene pyrotag and iTag datasets. Tag sequences were clustered into operational taxonomic units (OTUs) using the Quantitative Insights Into Microbial Ecology (QIIME) software package 57 and annotated using the SILVA or GreenGenes databases 58,59 . Community structure was determined using relative abundance and distribution of obtained microbial OTUs along with statistical analyses such as hierarchical clustering and indicator species analysis to identify characteristic groups of OTUs occurring under specific water column conditions such as different water column O 2 concentrations.
Illumina metagenomes, metatranscriptomes and metaproteomes were analysed using MetaPathways (version 2.0 and 2.5), an open source pipeline for functional annotation and prediction of reactions and pathways in metagenomes and metatranscriptomes 60 . Direct link to software download and specifications can be found online on the Hallam Lab Github repository (https://github.com/hallamlab/ metapathways2).
With respect to the metaproteomic datasets there is redundancy in the protein database used for peptide mapping in that amino acid sequences will have different names but identical sequences, and thus a given peptide may map to multiple sequences. To manage this one could cluster the proteins by sequence identity at the amino acid level prior to mapping. In previous publications we have calculated a normalised spectral abundance factor (NSAF) as a pseudo-quantative metric to compare abundance of the same proteins between samples 40 .