Samples were collected during March 2010 during R/V Cape Hatteras cruise GC-5 as described by Tolar et al., in revision; see also Figure A1 in Appendix for station locations and Table A1 in Appendix for selected sample data). Samples were obtained using Niskin bottles and a General Oceanics rosette sampling system equipped with a SBE25 CTD data package. They were pressure filtered through 0.22 μm Durapore filters (about 1 l at ∼60 kPa; Millipore, Inc., New Bedford, MA, USA). Filters were frozen at −20°C in 2 mL lysis buffer (0.75 M sucrose, 40 mM EDTA, 50 mM Tris; pH 8.3). After enzymatic hydrolysis (lysozyme and proteinase K with sodium dodecyl sulfate), DNA in 800 μL of lysate was purified by phenol-chloroform extraction. Purified DNA was stored in 50 μL of Tris-EDTA buffer (pH 8) at −80°C (Bano and Hollibaugh, 2002).

Triplicate PCR reactions for each sample consisted of 12 μL of PCR grade water, 2.5 μL 10× High Fidelity PCR Buffer (Invitrogen), 0.75 μL 25 mM dNTP mix, 1 μL 50 mM MgSO₄, 5 μL 5× Bovine Serum Albumin (Promega), 1.5 μL each of forward and reverse primers (10 mM stocks), 0.2 μL Platinum Taq DNA Polymerase High Fidelity (Invitrogen), and 0.5 μL template DNA. Primers 515F and 806R were used to amplify Bacteria and Archaea 16S rRNA genes (Liu et al., 2007; Jones et al., 2009; Bergmann et al., 2011). Primer 515F included the Roche 454-B pyrosequencing adapter and a GT linker. Primer 806R included the Roche 454-A sequencing adapter, a 12-bp unique barcode (Hamady et al., 2010), and a GG linker. The PCR program consisted of an initial denaturation step (94°C; 3 min), followed by 26 cycles of 94°C (1 min), 54°C (1 min), 68°C (2 min), with a final extension of 10 min at 68°C. Amplicons were visualized by electrophoresis on a 1% agarose gel. Products from triplicate reactions were pooled and purified with UltraClean PCR Clean-up kits (MoBio; Folsom, CA, USA) according to the manufacturer’s instructions. After quantifying DNA concentrations of the cleaned reactions using a Nanodrop spectrophotometer, equal masses of PCR product for each sample were combined and shipped to the Environmental Genomics Core Facility at the University of South Carolina where pyrosequencing was performed using a Roche 454 automated sequencer using titanium chemistry.

Sequences were analyzed using PANGEA (Giongo et al., 2010) and Mothur (Schloss et al., 2009) pipelines. After trimming the initial set of reads (435,290 sequences), PANGEA produced 239,983 sequences (>150 bp and ≥20 quality score) with identifiable barcodes distributed among the 44 samples (minimum sequence number per sample = 764, maximum = 9154; mean = 5454, SE = 362). A total of 38,045 reads could not be associated with a barcode. PANGEA identified the phylogenetic affiliations of the reads using MEGABLAST with a database of 170,273 Bacteria and Archaea isolate 16S rRNA sequences obtained from the Ribosomal Database Project. PANGEA clustered unidentified sequences for each sample using CD-HIT and threshold values (D) of 0.80, 0.90, 0.95, 0.97, and 0.99. OTUs comprised of sequences that could be identified with MEGABLAST were clustered similarly for each sample, and these clusters were combined with unidentified OTUs to produce the final composition for each sample. Spatial patterns in the resulting sample compositions, including OTU abundances, were analyzed using multivariate statistics [e.g., principal components analysis (PCA), non-metric multidimensional scaling, and canonical correlation analysis (CCorA)] with XLSTAT after excluding cyanobacterial and mitochondrial sequences, and removing OTUs represented by singletons within the 44-sample set.

Sequences were also processed using the Mothur platform (Schloss et al., 2009). After trimming barcodes and primers, and filtering for quality, sequences ≥150 bp were aligned to the Greengenes ProkMSA aligned database. Chimeric sequences were identified with the Chimera Slayer algorithm from the Broad Institute and removed. The pre.cluster command was used to minimize errors introduced during pyrosequencing. Cyanobacterial sequences were excluded from analysis. A total of 109,867 sequences comprised of 15,071 unique sequences were included in the analysis. For OTU-based analyses of sample composition, the “average neighbor” clustering algorithm was used to group sequences at a similarity level of 97%. The taxonomy function in Mothur was used with the Ribsosomal Database Project (RDP) training set to identify each OTU for phylotype-based. Because the RDP database used for analysis lacked Thaumarchaeota sequences, and because this phylum is highly represented in the Gulf of Mexico (Tolar et al., in revision), representative sequences for all OTUs that were identified as Archaea at a distance of 0.03 were analyzed with BLAST against the Greengenes prokMSA database to verify their identity. We also used the NCBI database to conduct a manual BLAST analysis of the representative sequences for the most abundant OTUs at a distance of 0.03 to verify identities, and to assess membership within SAR11, SAR86 and SAR92, and SAR324 clades for α-, γ-, and δ-Proteobacteria, respectively. These results were further verified by analyzing representative sequences in the ARB-SILVA database using the classify function. Mothur was also used to generate diversity indices based on samples with normalized numbers of sequences, to select specific phylogenetic groups for further analysis and to provide input for fast UniFrac (Hamady et al., 2010), which was used to assess spatial patterns in bacterioplankton communities.

The phylogenetic structure of nGoM bacterioplankton was also analyzed using Phylocom (Webb et al., 2002) to calculate the “net relatedness index” (NRI) and “nearest taxon index” (NTI). NRI is a standardized estimate of the mean pairwise phylogenetic distances for all pairs of OTUs in a sample compared to the mean pairwise distances of a random or null set of OTUs; NRI provides a measure of the extent of tree-wide clustering, including terminal and deep branches. NTI is a similarly standardized measure, but assesses the phylogenetic distance of each OTU in a sample to its nearest neighbor OTU; NTI provides a measure of terminal clustering independent of deeper clustering (Webb et al., 2002). Both indices were generated with 9999 randomized runs and included OTU frequencies.

Sequence data have been deposited with MG-RAST (metagenomics.anl.gov) at accession numbers 4509220.3–4509263.3. Metadata are available via the project page, “Analysis of composition and structure of coastal to mesopelagic bacterioplankton communities in the nGoM.”

This first detailed assessment of nGoM bacterioplankton spatial variability, diversity, and phylogenetic community structure reveals that composition varies substantially between surface (≤100 m) and sub-surface (>100 m) samples (e.g., Figures 1A,B and 6). Although changes in community composition with depth typify marine systems (e.g., Lin et al., 2006; De Corte et al., 2009; Church et al., 2010; Agogue et al., 2011; Eiler et al., 2011; Zinger et al., 2011; Friedline et al., 2012), the observed distribution of nGoM Proteobacteria differs from patterns reported by others. For example, in most other systems, α-Proteobacteria abundance typically exceeds that of γ-Proteobacteria regardless of depth (e.g., Cottrell and Kirchman, 2000; Pommier et al., 2007; Barberán and Casamayor, 2010; Galand et al., 2010; Zinger et al., 2011). While this trend holds for nGoM surface samples (≤100 m), γ-Proteobacteria abundance exceeds that for α-Proteobacteria in most nGoM deep-water sites, in some cases by >2-fold (Figure 2). Similar relationships have also been described for bathypelagic bacterioplankton in the Sea of Marmara (Quaiser et al., 2011) and the deep Arctic Ocean (Galand et al., 2010), but factors that account for this similarity are unclear.

The controls of nGoM α-Proteobacteria distribution are uncertain, as are reasons for the increased relative abundance of γ-Proteobacteria with depth. Multiple linear regression analyses show that neither depth, temperature, dissolved oxygen, pH, fluorescence (a measure of chlorophyll concentration), nor particle concentration (beam attenuation) individually accounts for α-Proteobacteria variability. However, the interaction of temperature and oxygen explains >70% (r² = 0.760, p < 0.0001). While this suggests that α-Proteobacteria might be sensitive to the combined effects of several variables that affect metabolic activity, factors other than or in addition to those included here likely play important roles.

Both the individual variables noted above and interactions among variables have limited explanatory power for γ-Proteobacteria distributions (e.g., maximum r² = 0.244, p = 0.003). This indicates that nGoM α- and γ-Proteobacteria respond to different ecological determinants. The availability of aliphatic and aromatic hydrocarbons emanating from extensive hydrocarbon deposits on the nGoM continental slope (Sassen et al., 2001; Milkov and Sassen, 2003; Liu et al., 2009) might account for the increased relative abundance of γ-Proteobacteria in deeper waters of the nGoM compared to other systems. Several prior studies in other systems have shown that γ-Proteobacteria respond strongly to hydrocarbon availability (Gerdes et al., 2005; Al-Awadhi et al., 2007; Berthe-Corti and Nachtkamp, 2010), and Hazen et al. (2010), Valentine et al. (2010), Kessler et al. (2011) have shown that various γ-Proteobacteria formed blooms at depths of about 1000–1200 m in response to the DWH hydrocarbon plume. In support of a role for hydrocarbon availability as a structuring factor for deep γ-Proteobacteria, members of the hydrocarbon-oxidizing genus, Pseudoalteromonas (Melcher et al., 2002), are among the most abundant phylotopes in deep nGoM waters (Figure 4B).

The increase in relative abundance of both γ- and δ-Proteobacteria with depth (Figures 3A and 4B) might also be supported by chemolithoautotrophic growth as proposed by Swan et al. (2011). However, Swan et al. (2011) characterized clades of γ-Proteobacteria (ARCTIC96BD-19 and Agg47) that are poorly represented in the nGoM, and it is not yet known if some of the better represented groups, e.g., SAR86 and SAR92, grow chemolithoautotrophically. Initial genomic analyses suggest that SAR86 lacks genes essential for CO₂ fixation (Dupont et al., 2012). In contrast, Swan et al. (2011) provided evidence for chemolithoautotrophic metabolism in a δ-Proteobacteria clade that is important in deep-water nGoM communities. This clade, SAR324, accounts for about 66% of nGoM δ-Proteobacteria. Thus, chemoautolithotrophy might account for changes in δ-Proteobacteria, while other factors determine the distribution of γ-Proteobacteria.

The abrupt increase in Thaumarchaeota relative abundance below 100 m as reported (Figure 3C) here agrees with results of a separate qPCR study of Archaea and Bacteria 16S rRNA genes at the same sites (Tolar et al., in revision). Tolar et al. (in revision) show that compared to Bacteria, absolute and relative Thaumarchaeota abundances increase in deep samples. An earlier study of nGoM Archaea at a single station (samples from 10, 400, 900 m) also observed that Marine Group I Crenarchaeota (i.e., Thaumarchaeota) dominate at depth, while Group II-β Euryarchaeota dominate in surface waters (Liu et al., 2009).

Although multiple reports have documented increased Thaumarchaeota abundance with increasing water column depth, the controls of this pattern remain unclear. Some have proposed that Thaumarchaeota are better adapted to environments characterized by low metabolic energy fluxes (see Valentine, 2007), and thus are better able than other lineages to couple growth to the limited resources (e.g., ammonium and perhaps some heterotrophic substrates) available below the epipelagic zone (Pester et al., 2011).

This hypothesis has implications for the spatial partitioning of thaumarchaeal clades reported by others (Santoro et al., 2010; Hu et al., 2011a,b; Yakimov et al., 2011) and also observed here (Figure A3 in Appendix). In particular, a “shallow water” clade consisting of phylotypes from surface waters and sediments has been distinguished from a “deep-water” clade found largely in meso- to bathypelagic samples (Beman et al., 2008; Hu et al., 2011a,b; Tolar et al., in revision). For nGoM samples, deep-clade phylotypes increase in relative abundance abruptly below 100 m, and then continue to increase linearly with depth (Figure A3 in Appendix). This might indicate that the surface-clade phylotypes are less well adapted than deep-clade phylotypes to low metabolic energy fluxes. Differences in gene expression for members of each clade in surface and deep samples could provide an indication of adaptations and competitiveness.

Although nGoM bacterioplankton composition changes substantially with depth, neither richness (e.g., Chao 1), nor evenness (e.g., Shannon index) and dominance (e.g., inverse Simpson’s index) measures of diversity show any distinct geographic trend (Figure 7A; Table A2 and Figure A5 in Appendix). The absence of trends within nGoM bacterioplankton at a regional scale is consistent with patterns reported for samples at a global scale. In particular, no distinct general trends as a function of depth or location have been reported for Chao 1 and the Shannon and inverse Simpson’s indices for assemblages ranging from epipelagic to bathypelagic depths in the sub-tropical Pacific to the Arctic Ocean (summarized in Table 2, Stevens and Ulloa, 2008), although site-specific trends have been described for Archaea and Bacteria (e.g., Brown et al., 2009). This similarity suggests common patterns for community assembly regardless of location, depth, temperature, or other variables. Comparisons among data sets should be treated cautiously, however, since diversity indices have well known limitations (Hill et al., 2003). These include sensitivity to sequencing effort, numbers of OTUs, and patterns of rank abundance within a given community. Nonetheless, Hill et al. (2003) have shown that Chao 1 and the Shannon index can provide robust measures of changes among samples or systems.

To promote comparisons for the nGoM, we have used the Mothur pipeline (Schloss et al., 2009) to produce normalized samples comprised of equal numbers of sequences randomly selected from the pool of sequences available for each sample. We have then compared theoretical minimum and maximum Shannon indices with observed values (Figure 7A). Minimum values assume a single dominant OTU, and that all others occur as singletons. Maximum values assume that all OTUs are equally abundant. With one notable exception, a sample from 760 m, most of the assemblages are characterized by Shannon indices that range between about 50 and 70% of theoretical maxima. This provides further support for the notion that the structure of nGoM bacterioplankton communities develops independently of major physical-chemical variables (e.g., depth, temperature, oxygen) and even some biological variables (e.g., community composition, chlorophyll concentration).

Two measures of phylogenetic community structure, NRI and NTI, provide additional insights about nGoM bacterioplankton communities. For this study, NRI values have been calculated using the full set of nGoM OTUs for null models; however, the choice of null model (see http:/phylodiversity.net/phylocom/phylocom_manual.pdf) did not affect the outcome.

Like nGoM diversity indices, neither NRI nor NTI vary consistently spatially (Figure 7B). This suggests that substantial changes in numerous physical-chemical and biological variables have little effect on phylogenetic community structure. Negative NRI values for all samples indicate that nGoM bacterioplankton communities are phylogenetically overdispersed throughout the water column. Overdispersion arises in a community when its OTUs are more distantly related to each other across terminal and deep branches of a community phylogenetic tree than expected by chance. Multiple processes have been proposed to account for overdispersion, including competitive exclusion when traits that determine community membership are phylogenetically conserved, and habitat filtering when such traits arise in different lineages by convergence or horizontal gene transfer (Webb et al., 2002; Vamosi et al., 2009).

Distinguishing between these or other processes that might lead to overdispersion requires a deeper understanding of the factors that constrain bacterioplankton richness than is currently available. However, since oligotrophs dominate the marine microbiota (Lauro et al., 2009), heterotrophic substrate limitation might promote competitive exclusion, particularly if different groups of bacteria specialize on different carbon sources. Lauro et al. (2009) have recently identified genomic/metagenomic markers for copiotrophs and oligotrophs. Phylogenetic analysis of these markers in nGoM and other bacterioplankton could provide a means to test relationships between substrate utilization and phylogenetic community structure.

In contrast to uniformly negative NRI values, uniformly positive NTI values (Figure 7B) indicate that nGoM OTUs occurring at the tips of the community trees form clusters more closely related than expected by chance. OTU clustering has been attributed to “habitat filtering,” which occurs when one or more environmental variables determine the patterns of community assembly. However, since clustering characterizes all nGoM samples, the environmental variables responsible for habitat filtering are uncertain, and likely vary among phylogenetic groups and with depth.

Habitat filtering has also been suggested as a structuring agent for surface bacterioplankton communities based on analyses of 16S rRNA gene sequences obtained from globally distributed samples (Barberán and Casamayor, 2010; Pontarp et al., 2012). For many of these communities, positive values have been obtained for both NRI and NTI, although in some cases similar to the nGoM, negative NRI, and positive NTI values have also been reported. The more consistently positive NRI values in global scale studies (Barberán and Casamayor, 2010; Pontarp et al., 2012) might have arisen in part from differences in the null communities used to calculate NRI (Webb et al., 2002). Null communities used in global scale studies have been derived from a global OTU pool, while the nGoM analysis has been based on a regional pool. Nonetheless, NTI values from all samples are consistent with environmental selection as a driver for assembly of both surface and deep communities.