Station ALOHA (22°45’N, 158°00’W) is located in the NPSG 100 km north of Oahu, Hawaii. The Hawaii Ocean Time-series (HOT) program¹ (Karl and Lukas, 1996) has made repeated observations of the hydrography, chemistry, and biology of the water column at Station ALOHA approximately once per month since October 1988. Sampling and analytical protocols are described elsewhere (Dore and Karl, 1996; Karl and Lukas, 1996; Letelier et al., 1996; Karl et al., 2001; Church et al., 2005). Metadata used for multivariate statistical analyses were downloaded from the HOT Data Organization and Graphical System (HOT-DOGS) website². Environmental variables used in this study are listed in Table A2 in Appendix.

Seawater samples were collected from discrete depths using 12 L polyvinyl chloride (PVC) bottles attached to a CTD rosette sampler, and sub-sampled into acid-rinsed 4 L polypropylene bottles. Three different sampling protocols were used throughout the study period (November 2004 to July 2008). For samples from cruises conducted between 2004–2006 (HOT cruises 165–185), 1–2 L of each sub-sample was filtered through 25 mm diameter, 0.22 μm pore-sized hydrophilic polyvinylidene difluoride (PVDF) membrane filters (Durapore; Millipore Corp., Billerica, MA, USA). The filters were placed in 2.0 mL microcentrifuge tubes containing 500 μL of Tris–EDTA (TE) buffer, immediately quick frozen in liquid nitrogen, and stored at −80°C until further processing. DNA was extracted and stored at −20°C following the protocol described by Church et al. (2009). For samples from cruises 186–193, 0.5–4 L of each sub-sample was filtered onto 13 mm diameter, 0.2 μm pore-sized polyethersulfone membrane filters (Supor-200; Pall Corp., Port Washington, NY, USA). The membrane filters were subsequently stored in 2 mL screw-cap tubes containing 250 μL of DNA lysis buffer [20 mM Tris HCl (pH 8.0), 2 mM EDTA (pH 8.0), 1.2% (vol/vol) Triton X100, and 20 μg/mL lysozyme] at −80°C until extraction. For HOT cruises 197–200, ca. 4 L sub-samples were filtered onto 25 mm diameter, 0.2 μm pore-sized polyethersulfone membranes (Pall Corp.), and subsequently stored at −80°C in 500 μL of DNA lysis buffer. DNA was extracted from cruises 185–200 using the Qiagen DNeasy 96 Tissue Kit (Qiagen Inc., Valencia, CA, USA) following the manufacturer’s recommended protocol for bacteria. Extracts were quantified using the PicoGreen fluorescent assay (Invitrogen Corp., Carlsbad, CA, USA) on a SpectraMax M2 plate reader (Molecular Devices Corp., Sunnyvale, CA, USA).

Two microliter of each DNA extract (equivalent to ca. 1–100 ng of genomic DNA) was used as template for the polymerase chain reaction (PCR) amplification of bacterial small subunit (SSU) ribosomal RNA (rRNA) genes. Individual 50 μL reactions contained 1× PicoMaxx reaction buffer (Stratagene, La Jolla, CA, USA), 200 μM of each deoxynucleotide, 200 nM each of the fluorescently labeled general bacterial SSU rRNA gene oligonucleotide primer 27F-B-6FAM (5′-[6∼FAM]AGRGTTYGATYMTGGCTCAG-3′) and universal SSU rRNA gene oligonucleotide primer 1492R (5′-GGYTACCTTGTTACGACTT-3′; Lane, 1991) and 0.625 units PicoMaxx thermostable DNA polymerase (Stratagene). Cycling conditions consisted of an initial denaturation step of 95°C for 5 min, followed by 35 cycles of 95°C denaturation for 30 sec, 55°C annealing for 1 min, and 72°C extension for 2 min, with a final extension step at 72°C for 20 min. PCR products were purified using the QIAquick Multiwell PCR Purification system (Qiagen Inc.), and digested in 20 μL reactions containing 4 μg acetylated bovine serum albumin (BSA), 1X MULTI-CORE enzymatic reaction buffer (Promega, Madison, WI, USA), and 10 units of HaeIII restriction endonuclease (Promega). After a 6-h incubation at 37°C, the digests were terminated by incubation at 70°C for 20 min, and purified using the Millipore MultiScreen Assay System (Millipore Corp.) in conjunction with Sephadex G-50 Superfine (GE Healthcare, Piscataway, NJ, USA). Approximately 100 ng of the purified products were electrophoresed on an ABI 3100 Genetic Analyzer (Applied Biosystems, Foster City, CA, USA). Operational taxonomic units (OTUs) were identified as peaks between 30 and 600 bp in length, using GeneMapper v4.0 software (Applied Biosystems, Foster City, CA, USA). Terminal restriction fragment length polymorphism (T-RFLP) profiles were first run through an in-house MATLAB (The MathWorks, Natick, MA, USA) script that corrects for instrumental drift while rounding fractional peaks to whole numbers. In brief, this script identified peaks with ambiguous decimal increments (i.e., x.3–x.6) and ensured that they were rounded to the same integer throughout all profiles by using a majority rule (i.e., rounded up if most were > x.5 and vice versa). Non-ambiguous decimal increments (i.e., <x.3 and >x.6) were simply rounded to their nearest integer. The profiles then underwent a limited manual curation in order to further compensate for instrumental drift, and finally normalized using a variable threshold (Osborne et al., 2006). Peaks were transformed into relative units by dividing integrated peak areas by the total peak area for that sample. Diversity (evenness) was estimated from the relative units using the Gini coefficient (Wittebolle et al., 2009). The Gini coefficient provides a proxy for the equality of taxa distributions: the higher the Gini coefficient, the more uneven a community is.

Putative identification of abundant OTUs was performed using clone library data obtained from both coastal and open ocean sites in the vicinity of Hawaii, including Station ALOHA. Representative clones were characterized by Sanger sequencing and fragment sizing after restriction digests to correlate phylogenetic identity with terminal restriction fragment length. Sequences were deposited in Genebank under accession number JN166096 to JN166367.

Statistical analyses. Statistical analyses were performed using the software packages PRIMER 6 (PRIMER-E Ltd., Ivybridge, UK), EcoSim (Gotelli and McCabe, 2002), PAST (Hammer et al., 2001), and MATLAB (The Mathworks). Sorenson (Bray–Curtis) distances were calculated from rank-transformed T-RFLP data using a natural weighting as recommended by McCune and Grace (2002), and visualized using non-metric multidimensional scaling (NMS; McCune and Grace, 2002).

Statistical analyses required the exclusion of samples with missing values. In order to account for such gaps in the availability of metadata, environmental variables were grouped into biological, chemical, and physical data classes in order to increase sample size for analyses (Table A2 in Appendix). Within PRIMER 6, BIOENV procedures were used to explore the correspondence between microbial community structure and environmental metadata (Clarke and Warwick, 2001). In this analysis, the environmental variables were first Log(X + 1) transformed and normalized by subtracting the mean from each entry of a single variable and dividing by the standard deviation for that variable. This allowed the identification of variables or combinations of variables that best corresponded with the community data. As this analysis was based on a limited number of samples due to missing data, the best fitting combinations were re-analyzed after adding samples with non-missing values (in the identified variable combination) using the RELATE test (Clarke and Warwick, 2001) in PRIMER 6 (PRIMER-E Ltd.). The Analysis of Similarity (ANOSIM in PRIMER 6; Clarke, 1993; Clarke and Warwick, 2001) test was used to examine for differences between bacterial communities from different depth horizons in the water column.

Two taxa co-occurrence indices were calculated using the EcoSim 700 software package (Gotelli and McCabe, 2002). The Checkboard index (C-board) calculates the number of taxa that never co-occur, which can be interpreted as evidence for competitive exclusion (segregation) or differences in habitat preference. The second index (C-score) is calculated by: CU = (Ri − S)(Rj − S) where Ri and Rj are the number of sites where taxon i and j occur, and S is the number of sites containing both taxa. The C-score is then calculated from the mean number of all pair-wise CUs. Monte Carlo simulations were used to test if the observed co-occurrence measures were random following the most conservative permutation procedure. Following Gotelli and McCabe (2002), we calculated a standardized effect size (SES) that allowed comparing the degree of co-occurrence across temporal and spatial scales.

Spearman’s rank correlation coefficient (ρs) was calculated in MATLAB in a pair-wise fashion between OTUs, using the relative contribution of each OTU to the total community for each sample. Only OTUs that were detected in more than two independent samples were included in the analysis (n = 397). A significance level of p < 0.05 was used to calculate the percentage of significant positive and negative correlations from the total number of pair-wise comparisons. In addition, Bonferroni correction was used to address the problem of multiple comparisons. In the same way, ρs between relative contributions of each OTU and depth were calculated.

To test if the number and distribution of significant pair-wise correlations was different from random, Monte Carlo simulations were performed. Following a permutation procedure where the rank of OTUs within samples was randomized, a total of 1000 randomized matrices were created and subsequently used in pair-wise correlations. The maximum and minimum ρs values of correlations within each 0.01 interval (from 1 to −1) were determined from all 1000 randomized matrices.

A total of 461 different OTUs, each represented by a specific terminal restriction fragment length, were detected in 296 samples collected over 3.5 years (November 2004 to February 2008) from depths ranging from 5 to 4000 m at Station ALOHA. Of these, 397 OTUs were present in at least two T-RFLP profiles, and so were used for further analyses. Co-occurrence indices calculated from this dataset differed significantly from random simulations: more OTU pairs exhibited larger C-scores and C-boards than expected by chance (Table 1). A higher C-score than expected by chance indicates low co-occurrence, whereas a lower C-score than expected by chance indicates high co-occurrence. Low co-occurrence can be interpreted as evidence for competitive exclusion (segregation), or differences in habitat preference. Of the Spearman’s rank correlation coefficients (ρs) from all pair-wise comparisons of OTUs (n = 78,607), 27.3% were significant at a 5% significance level, of which 6.6% were negatively correlated (ρs < 0) and 20.7% were positively correlated (ρs > 0; Figure 1A). Both significantly positive and negative correlations occurred more often than expected by chance when compared to simulations where the rank order of species were randomized within a community. In addition, maximum ρs values were much higher for the observed matrix compared to the simulated matrices (Figure 1A). Using Bonferroni correction for addressing the problem of multiple comparisons, the number of significant correlations decreased from 21,485 to 3,276 (3.5% significant positive and 0.5% significant negative correlations).

On average, 108 significant correlations were observed per OTU; 22 OTUs were correlated with more than 50% of the 397 OTUs included in the dataset. This subgroup of 22 included the most abundant OTUs, which are putatively identified as Prochlorococcus and SAR11 subclades IB and IA. Other abundant OTUs putatively identified as SAR116 subclade II and SAR11 subclade IIB were significantly correlated with >150 OTUs. However, many OTUs contributing <1% in relative abundance (on average) to the total bacterioplankton community exhibited a high number of significantly correlations with other OTUs in our dataset (Figure 1B). Hence, we did not observe a significant correlation between the average relative abundance of an OTU and the number of correlations it exhibited.

Variation in the C-score along the depth gradient was investigated in detail by calculating the SES from C-score calculations grouped by depth (SEScscore; Figure 2). A positive SEScscore significantly different from random (p < 0.01) was observed at all depths, indicating low co-occurrence. The SEScscore was not correlated with depth, but exhibited intriguing variability along the vertical axis (Figure 2). The highest value was estimated at 400 m, and decreased again toward the oxygen minimum zone at ca. 800 m. The SEScscore did not vary with the size of the data set, as measured by either the number of samples or the number of different OTUs included in a given incidence matrix (ρ = 0.103, p = 0.696, and ρ = 0.359, p = 0.137, respectively). It could also be argued that the number of OTUs represents an estimate of community richness. Hence, this apparent lack of relationship between richness and SEScscore, together with evenness and SEScscore also not being related (ρ = −0.279, p = 0.276) indicates that co-occurrence patterns are not related to alpha diversity at least in our study.

The Gini coefficient, an estimate of community evenness, ranged from 0.007 to 0.067 in the 296-sample dataset. These low values indicate a high degree of equality, as a value of 0 indicates complete equality while a value of 1 indicates maximal inequality among OTUs. Spearman’s rank correlation analysis revealed a significant negative relationship between the Gini coefficient and depth (ρ = −0.27, p < 0.001) indicating an increase in community evenness with depth (Figure 2).

In order to visualize potential spatial and temporal patterns in community structure, communities were plotted using NMS to two dimensions based on a Bray–Curtis distance matrix (Figure 3). ANOSIM and RELATE analyses corroborated initial visual observations of depth-specific clustering of bacterioplankton communities at Station ALOHA (R = 0.613, p < 0.001, and ρ = 0.67, p < 0.001, respectively). An R-value of 0 resulting from the ANOSIM analysis indicates that there were no differences among the communities from different depths, while an R-value of close to 1 would indicate that the dissimilarities between communities from different depths were larger than those from each individual depth. However, pair-wise ANOSIM revealed that community structure from different depth horizons in the upper 100 m were not significantly different whereas, below 100 m, all depth horizons were significantly different from each other (Table A1 in Appendix). Exceptions were adjunct depths where, in most cases, differences were not significant (Table A1 in Appendix).

Average Bray–Curtis dissimilarities between communities from different depths increased with increasing distance along the vertical (i.e., depth) profile. The significant ρ value of the RELATE analysis indicated that communities from more distant depth horizons were more different than communities from adjunct depths. When Bray–Curtis similarities were averaged over the sampling period for each depth, values stayed rather constant above the mixed layer (Figure 2). In the depth layers where the mixed layer and chlorophyll maximum temporally oscillate, Bray–Curtis similarity averages were low. The relative abundance of nearly early half (46%) of all OTUs was significantly correlated with depth (p < 0.05), of which 32% were positive and 14% negative (Figure 4).

Unlike the clear depth-specific zonation of bacterioplankton communities, temporal patterns were not readily apparent in NMS plots (data not shown). However, we were able to statistically test for temporal patterns in bacterioplankton community structure, since the sampling schema spanned over more than three consecutive years. When months and seasons were used as discriminating factors in ANOSIM analyses, and when months, seasons, and Julian days were used in cyclic RELATE analyses, significant temporal patterns were observed at various depths. However, significant annual cycles in bacterioplankton community composition were only found in the upper 200 m of the water column (Table 2) and were not detected in the deeper horizons.

BIOENV procedures were employed to search for relationships (correspondence) between community and environmental variability, and to identify the strongest relationships. Since the environmental dataset was incomplete, variables were divided into three categories: physiochemical, biogeochemical, and biological (Table A2 in Appendix). Data from each category was used separately in the BIOENV procedure, and the variable combinations that best matched the patterns in community composition were identified and then analyzed by RELATE after inflating the number of samples (Table 3). Parameters from all three general categories revealed significant correspondence along the sampled spatial and temporal scales (Table 3). This included depth, which represented the single variable that corresponded the best with community composition. Variation in bacterioplankton community structure in the upper 200 m of the water column was highly correlated with concentrations in heterotrophic prokaryotes, Prochlorococcus, Synechococcus, and pigments such as chlorophyll a and phaeopigments. In addition, inorganic nutrients such as nitrite, nitrate, and particulate carbon corresponded well with community composition. Most of these variables were also correlated with depth (Table A3 in Appendix), and annual periodicity in these variables was either weak or not present (data not shown).