Oceanographic habitat and the coral microbiomes of urban-impacted reefs

Study region

The Florida reef tract (FRT) is a coral reef system in the USA, that stretches from the Dry Tortugas (83.0W, 24.5N) east and north to the lower edge of the South Atlantic Bight (80.0W, 27.3N). The northern third of the FRT lies along the southeast Florida shelf, offshore of Miami-Dade, Broward, and two other counties. For this study, four reefs were sampled (Fig. 1): Barracuda and Oakland Ridge reefs in Broward, and Emerald and Pillars reefs in Miami-Dade. Three distinct sites within each reef were selected for replicate coral tissue collections, resulting in a total of 12 sampling sites. These sites are in waters between 6 and 12 m deep, and lie between 1.1 and 5.5 km offshore of mainland Florida.

Sample collection

The data for this study were leveraged from previous projects. Details about the sites and collection methods have previously been presented (Sinigalliano et al., 2018; Staley et al., 2017). Briefly, samples for this study were collected during March, June, and September of 2015 (see Table 1). Visually healthy corals were sampled at random from two coral species, Porites astreoides and S. siderea, with SCUBA and a syringe biopsy method. This entailed loosening the polyps from the calyx by gently scraping/digging at the polyp tissue before adhering a five mL syringe against the coral polyp to create a suction and extract that polyp from the coral head. This was repeated twice per coral head and combined into one sample in the lab. When divers reached the boat, samples were transferred to a 50 mL tube with 45 mL of 95% molecular grade ethanol, placed on ice, then preserved at −80 °C at the NOAA AOML laboratory. In addition, divers measured sea temperature with a dive-watch (see Physical oceanographic habitat analysis), and both overall and per-species coral cover and coral mean size (Sinigalliano et al., 2018) at each of the reef sites. A total of 24 Porites astreoides and 31 S. siderea samples were collected across four reefs (Barracuda, Oakland Ridge, Emerald, and Pillars) during the 3 months of sampling (March, June, and September).

DNA extractions and microbiome amplicon sequencing

Coral polyps, which included the tissue and remnant mucus, were processed for DNA extraction and sequencing, briefly as follows (Staley et al., 2017): samples were vortexed to resuspend coral tissue, and the entire contents were placed on a 0.2 µm polycarbonate filter under a vacuum. Both the filter and any residual large coral tissue were transferred to a bead beating Lysing Matrix “E” tube. The samples were then homogenized by bead beating with a FastPrep 24 instrument (MP, Biomedicals, Irvine, CA, USA) by two separate rounds of bead beating at an impact speed of 6.0 m/s for 60 s each (for a total of 120 s at 6 m/s). Each lysed sample was then purified with the FastDNA SPIN kit for soil DNA Extraction (MP, Biomedical) following the manufacturer’s instructions and eluted in 100 μL of TE buffer. DNA was amplified using dual index primers, 515F and 806R, that target the variable-4 (V4) region of the 16S rRNA gene. Amplified DNA was sequenced on the MiSeq platform, generating paired-end reads of 300 base pairs (bps).

Bioinformatic analysis

Demultiplexed sequences from the 55 samples were retrieved from the National Center for Biotechnology Information (NCBI) Sequence Read Archive database (#SRP076111). Once downloaded, the sequence primers were trimmed using the program cutadapt. The sequences were analyzed with qiime2-2018.6 using the DADA2 pipeline for classification of amplicon sequence variants (ASVs), and were run with a max expected error of 2 (Callahan et al., 2016). To select sequence trimming positions, the quality score of a random set of sequences were analyzed. Where sequences showed a drop-in quality score (<Phred = 30), those bp positions where chosen for trimming. All forward sequences were truncated to the same length of 220 bps and trimmed at the first eight bps. The reserve reads were also truncated to be of the same length at 120 bps, and trimmed at the first 20 bps. The DADA2 pipeline was used for ASV selection, dereplication, merging paired-end reads (min overlap = 20 bps) and removing chimeras (with the “consensus” option). Taxonomy was assigned with the Silva database version132 training set, which was trained on the V4 region by the fit-classifier classify-sklearn function. Only sequences that were taxonomically identified as Bacteria or Archaea were analyzed. The biome data, ASV table, fasta file, phylogenetic tree, and metadata for this analysis can be found on figshare (DOI 10.6084/m9.figshare.7388672). The specific NCBI accession number for each sequence file is listed under “sampleID” in the metadata file.

Statistical analysis

Microbial alpha-diversity was evaluated with qiime2-2018.6, using the Shannon diversity index and species evenness metrics. To assess alpha-diversity, ASV counts were rarefied to (randomly filtered to an equal sequencing depth of) 403 counts, so that all 55 samples could be evaluated. Since 403 is a low number of sequences, the ASV table was also rarefied to 2,055 (N = 53). We compared these two methods and since both methods resulted in the same conclusions, a rarefaction of 403 was used in order to include all 55 samples in alpha-diversity analysis. Categorical habitat data (see Table 1) were tested for significant relationships with Shannon diversity index and species evenness using the alpha-group-significance function, which applies a pairwise/all-group Kruskal–Wallis test (Gregr & Bodtker, 2007). Continuous habitat data (see Table 1) were tested vs. alpha-diversity with a simple linear regression, using the function lm on R.3.4.3 (R Core Team, 2017).

For beta-diversity tests, the data were not rarefied, but instead were filtered by retaining only those taxa that were seen four or more times in at least 10% of the samples. The counts were then transformed into relative abundances. Continuous habitat data were evaluated for pairwise covariance; from each pair with a correlation >80%, one variable was selected for further analysis. Principal coordinate analysis (PCoA) ordination was applied to the relative abundance of ASVs with a weighted UniFrac distance metric (Lozupone & Knight, 2005). Discrete variables (i.e., reef site, species, and month) were tested using Analysis of similarities (ANOSIM) and a Permutational multivariate analysis of variance (PERMANOVA) with 999 permutations to identify which groups were different (Clarke & Warwick, 1994). The counts matrix was also used to build a model with forward and backward model selection, using the vegan function ordistep, which selects significant variables that best explain the counts data. Prior to inputting habitat variables into ordistep, their respective variance-inflation factors (Salmerón, García & García, 2018) (VIF) were tested; variables with a VIF score >20 were removed from further analysis. The variables selected by the model were then used to construct a canonical or constrained correspondence analysis (CCA) (Henriques et al., 2006). These statistical analyses were conducted with the packages vegan 2.4.6 and phyloseq 1.22.3.

To further evaluate patterns seen in the beta-diversity results, we used a classification method to evaluate which microbial taxa were important to the physical habitat variables that were found to be significant using ordistep (above). The classification method used on each of the significant physical habitat variables was the sample-classifier regress-samples function from qiime2. This function performs random forest (a supervised machine learning algorithm) and was used with parameters --p-parameter-tuning, --p-estimator RandomForestRegressor (RFR), and --p-n-estimators 500. The model outputs importance values and the highest importance values have the greatest prediction power. Thus, the five taxa with the highest importance values from each random forest test were selected for simple linear regression analysis vs. the corresponding physical habitat variable. Only the taxa with a significant correlation (p < 0.05) and a coefficient of determination R² > 0.2 were plotted. The code for the statistical analysis for the microbial data and the corresponding figures are found on figshare (DOI 10.6084/m9.figshare.7925495).

Physical oceanographic habitat analysis

The minimum near-bottom sea temperatures were collected by divers (see Sample collection). Replicates for each measurement are listed in Table 1. Seafloor depths (bathymetry) were derived from NOAA’s 10 m horizontal-resolution port-bathymetry project (Carignan et al., 2015). From bathymetry, we estimated local seafloor depth ([m]), seafloor slope (rise/run), isobath orientation (alongshore direction, deg. True), and distance from shore. We also estimated distance from each sampling site to the nearest major inlet running through the Intracoastal Waterway (Fig. 1; Table 1).

Satellite multispectral ocean-color was used to produce a proxy for in-water turbidity, a product called relative turbidity or color index (CI). Since 2000, the MODerate-resolution Imaging Spectrophotometer instruments aboard the polar-orbiting Aqua and Terra satellites have produced multispectral ocean color data at 250 m horizontal resolution. CI datasets were produced from these data, using an algorithm developed by the University of South Florida’s Optical Oceanography Laboratory (Barnes et al., 2015).

The CI product is subject to bias in shallow waters where the seafloor may be visible from space. In this study, CI values were normalized using the median and 7th and 93rd percentiles for 15 years of CI data (2000–2015) from each pixel. Percentile values were used because CI for coastal pixels is often not normally distributed (Gramer & Hendee, 2018). On days when clouds obscured a pixel during one or both daytime satellite overpasses, no CI value was available for that site. Therefore, for the purposes of this study, we used a 14 days arithmetic mean of CI values centered on each sample date, resulting in relative turbidity time series with no gaps (Figs. 2A–2C). Finally, time series for satellite CI were validated by linear regression vs. in-water turbidity measurements (NTU) done with boat-based casts using a factory-calibrated nephelometer throughout the sampling period (2013–2015) (Fig. S1; Staley et al., 2017).

Ocean currents were extracted as daily snapshots from the Gulf of Mexico HYbrid-Coordinate Ocean Model (GoM HYCOM; Gierach, Subrahmanyam & Thoppil, 2009) at the vertical layer closest to the surface. Northward and eastward vector-components were bilinearly interpolated from native model resolution (1/25th degree) to individual sampling location. Current vectors were then rotated so that positive cross-shore currents were directed perpendicular to the local isobath (i.e., “offshore”) at each site. Daily snapshots of currents from GoM HYCOM alias the diurnal and semidiurnal tidal cycles, showing the distinct signature of the 27.55 days M_m tidal constituent. Furthermore, although we expect GoM HYCOM to provide useful statistics for the amplitude of lower-frequency (multi-day) current variability, the model may not always correctly model the phase of that variability (Gramer, 2013). Thus, centered moving 14 days averages of each current component (cross-shore and alongshore) were analyzed (Figs. 2D–2F).

To estimate surface wave action, we used significant wave height from NOAA’s WaveWatch III multi-mesh grid operational model for the western Atlantic, with four nautical-mile nominal resolution (Lee et al., 2009; Tolman, 2014). A wave attenuation algorithm was applied to wave height, bilinearly interpolated to the horizontal resolution of the bathymetry (10 m) (Fig. S2; Gramer & Hendee, 2018). This algorithm linearly reduces significant wave height at bottom-depths between 0 and 20 m, resulting in high-resolution wave fields that reach all the way inshore, with zero wave-height at the beach (Fig. S3).

Ocean model outputs and diver-watch sea temperatures were validated using in situ hourly observations from acoustic Doppler current profilers (ADCP) (Fig. S4). These ADCPs were moored at two sites, one at 7 and one at 26 m depth, six km north of Barracuda Reef (Sinigalliano et al., 2018). They measured near-bottom sea temperatures (nominal accuracy 0.01 K) and three-dimensional ocean current profiles (bin sizes 0.5 m at 7 m depth and 1 m at 26 m, nominal accuracy 10⁻⁴ m s⁻¹) every 20 min throughout all three sampling periods of March, June, and September of 2015.

Tide heights for each site were modeled using the TPXO 7.2 global tide solution (Stammer et al., 2014). Hourly in situ measurements of near-bottom sea-pressure from the FACE mooring at seven m depth were also used with mean depth removed (“sea-level anomaly”), both to validate model tide-stage predictions and to consider the additional effects of setup from wind, waves, and larger-scale sea-level variation offshore. The time of day of collection was not logged for the tissue samples analyzed in this study, nor can the TPXO 7.2 global solution model fine-scale variability in tide heights and currents on the topographically complex southeast Florida shelf (~1 km, diurnal or faster). For these reasons, a 14 days centered average of both the modeled tide height and the sea-level anomaly (Fig. S5) were analyzed in the present study. The Ecoforecasts toolkit, a set of MATLAB functions developed by Gramer as a part of NOAA CHAMP, was used to produce some time series analyses and maps presented here. Due to the prevalence of non-Gaussian distributions in habitat data (Fig. S6), results for habitat variables are cited as median ± half of interquartile range values. Elsewhere, for example, for relative abundances, mean and standard deviation are reported. The code used to generate physical ocean habitat variables and the corresponding figures are found on figshare (DOI 10.6084/m9.figshare.7925546).

The physical oceanographic habitat of four southeast Florida coral reefs in 2015

To extract all other physical habitat variables for each of the four reefs, we used longitude and latitude coordinates from each of the 12 sample-collection sites (three sites per reef, Fig. 1) together with the dates of sampling, as described in the Methods section. March was a period of geographic patchiness in CI (Fig. 2A). However, patchiness was largely confined during March sampling dates to waters away from the sampling sites. On the other hand, both June and September (Figs. 2B and 2C) showed less intense regional patchiness, but patches of higher relative turbidity did reach some of the sites during sampling. Satellite relative turbidity (CI) showed a complex seasonal pattern than the other variables considered (Fig. 3C). The highest extreme and the greatest range of values occurred during June sampling (−0.43–0.23, median −0.16), with sampling in September scattered over a mainly lower range (−0.50 to −0.08, median −0.23), and that from March clustered within a very narrow range (−0.13 to −0.04, median −0.09). Barracuda was the reef that experienced the greatest positive and negative extremes of normalized CI during the 2015 sampling, with a peak at the three Barracuda sites of between +0.06 and +0.23 on the sampling day of June 12th, and lows between −0.30 and −0.50 on September 11th. Interestingly, June was also the month with the highest CI value at Oakland Ridge (0.11), despite lower modeled waves. (Results of regressing satellite CI values vs. in-water turbidity measurements are shown in Fig. S1)

Alongshore currents from the Gulf of Mexico HYCOM model (Figs. 2D–2F and 3D) showed a relative peak at all sampling sites during June, with weekly medians of +0.30 to +0.65 m s⁻¹ (positive being a northward current). Further peaks in July and August were followed by a rapid decline during sampling weeks in September, to medians of +0.10 to +0.45, which was also similar to medians during the March sampling. At the time of coral tissue sampling, in situ ocean current profilers were deployed near the sampling sites. The midsummer peak could also be seen in these instruments (Fig. 3D, gray and black lines). According to climatology models (Fig. 2E), this summer peak corresponded to a westward meander of the Florida Current; this June meander however, affected the southernmost sampling sites, those at Emerald Reef (median +0.60 ± 0.00 m s⁻¹), more than at the other three reefs (medians +0.30 to +0.48 ± 0.01 m s⁻¹). Emerald Reef actually experienced higher alongshore currents than the other sites during all three sampling months, although these differences were less extreme in March and September than they were in June.

The minimum diver-watch sea temperature was recorded during each tissue-sampling dive (Table 1). For validation, dive-watch temperatures were compared to sea temperatures measured at two moorings (that at 7 and that at 26 m depth). Dive-watch sea temperature followed closely the coincident mooring temperatures at seven m (Fig. 3A; black line): regressions were significant (p < 0.0001, Fig. S4), with root mean-squared error (RMSE) 0.28 at Barracuda, 0.78 at Oakland Ridge, 0.31 at Pillars, and 0.29 K at Emerald. Regressions vs. the deeper mooring at 26 m were also significant (p < 0.05) with RMSE of 0.8 at Barracuda, 1.0 at Oakland, 0.34 at Pillars, and 0.25 K at Emerald.

Like the moorings, dive-watch sea temperature across the four reefs followed a seasonal pattern (Fig. 3A). For March, there was a median 23.2 and interquartile range ± 0.4, for June 28.0 ± 0.2, and for September 29.2 ± 0.5 °C. However, temperatures were somewhat cooler than might be expected from historical in situ measurements in this region, especially during the warmer months (Gramer et al., 2018). In particular, the mooring at seven m depth (Fig. 3A) recorded near-bottom temperatures of 24.7 ± 0.5 in March, 29.2 ± 0.9 in June, and 31.2 ± 0.4 °C in September of 2015. Similarly, near-surface temperature measurements in two m of water at the nearby Virginia Key pier (NOAA National Ocean Service station “VAKF1”) showed warmer temperature ranges in June and September of 29.7 ± 1.1 and 30.1 ± 0.6 °C, respectively (Fig. 3A; red line).

Wave heights from the attenuated NOAA model showed a relatively simple seasonal pattern for all sites (Fig. 3B), with the highest waves at all sites during sampling in March (0.25–0.35 m), and intermediate waves during sampling in late September (0.2–0.3 m). Waves in June were generally lower (0.15–0.25 m), due to the absence of storm systems in the nearby Atlantic and Gulf of Mexico, and to generally light local summer winds. The lowest wave estimates, however, were those for samples taken in early September at Oakland Ridge and Barracuda (0.08 m); these two reefs were generally the most sheltered (Fig. 1; Fig. S3). Furthermore, these samples were gathered prior to a series of wind events in mid-September, and perhaps more importantly, prior to the development of a major tropical weather system in the Atlantic in late September (Fig. 3B; see the Discussion section).

Tide models used for this study do not provide sufficient spatial or phase resolution to estimate diurnal or semidiurnal variability in water heights or tidal currents. However, these models do resolve variability between sites associated with the M_m tidal constituent, which has a 27.55 days period. Tide height variations for this constituent were at or near their minimum (−0.008 to −0.005 m) during all March sampling dates. In June, sampling was done near the M_m peak (+0.002 m) at all reefs except Oakland Ridge, where samples were gathered near the time of the minimum or ebb M_m (−0.005 m). Similarly, September samples were collected near the ebb at both Oakland Ridge and Barracuda (−0.006 to −0.004 m), and near the peak at Pillars and Emerald Reef (+0.001 m).

Unlike with the tide model, sea-level anomaly as estimated from the seven m mooring shows the effect of sea-level setup from wind and possibly other forcings, in addition to the effect of the barotropic tide (Cushman-Roisin & Beckers, 2011). A sea-level anomaly in March at all reefs clustered between −0.04 and −0.08, and in September, it clustered even more tightly at −0.09 for most reefs, with Pillars the sole September outlier at +0.16 m. Finally, in June, sea-level anomaly ranged more widely, from a high of +0.08 at Pillars, to +0.06 at Emerald, +0.01 at Barracuda, and −0.10 m at Oakland Ridge.

Fictibacillus and Endozoicomonas were dominant bacterial members in two hard coral species in four southeast Florida reefs

The 55 16S rRNA amplicon sequences from coral tissue samples (24 Porites astreoides and 31 S. siderea) were leveraged from a previous study (Staley et al., 2017). The program DADA2 identified a total of 22,451 ASVs and after filtering, 15,201 ASVs remained. The minimum frequency of an ASV per sample was 404 and the maximum frequency was 196,536. The median ASV frequency per sample was 102,667 with a mean of 98,118.

The microbial taxa of the samples are summarized in Fig. 4. On average, across the four Southeast Florida reefs in 2015, the bacterial genus Fictibacillus, family Bacillaceae made up the majority of the microbial community (N = 55; mean relative abundance, meanRA = 13.9%, ±standard deviation 1.8%). The bacterial genus Endozoicomonas, family Endozoicomonadaceae also had a high meanRA (12.2% ± 9.7%). When analyzing the two host-species separately, Porites astreoides was dominated by Endozoicomonas (N = 24; 10.7% ± 8.0%); S. siderea on the other hand, had a higher abundance of the genus Fictibacillus (N = 31; 10.4% ± 11.1%), while Endozoicomonas was only this host’s 9th most abundant taxon (1.5% ± 3.0%). During the month of September (N = 24) specifically, across samples from both host species, the family Bacillaceae was found at higher abundances, and was composed of the genera Paenibacillus (8.9% ± 12.2%) and Fictibacillus (7.7% ± 9.6%). However, for the months of March (N = 20) and June (N = 11), across all sites and both host species, Endozoicomonas showed the highest meanRA of 4.5% ± 5.2% (March) and 5.2% ± 9.0% (June). Two reef sites had the highest overall average abundances for Fictibacillus across all months and both hosts: Barracuda (N = 16; 4.3% ± 10.0%) and Emerald (N = 10; 3.8% ± 9.6%). The other two reefs showed their highest overall average abundances for Endozoicomonas across month and host: Oakland Ridge (N = 13; 4.9% ± 8.0%) and Pillars (N = 16; 4.3% ± 6.8%).

Bacterial alpha diversity correlated with the physical oceanographic habitat

Table 1 lists three (3) categorical and 14 numerical variables. A total of 12 physical habitat variables were selected for analysis to identify univariate correlations between these variables and coral microbiomes (Overall coral cover and overall mean size were not considered; only coral cover and mean size for the two host species were considered - Porites astreoides and S. siderea.). Latitude and longitude were removed from this analysis, as they correlated (ρ > 80%) with other independent habitat variables.

Microbial Shannon diversity was negatively correlated with alongshore current speed (R² = 0.05, p = 0.044), tide height (R² = 0.12, p = 0.004), and temperature (R² = 0.13, p = 0.002; Figs. 5A–5C). Similarly, microbial evenness was negatively correlated with sea-level (R² = 0.09, p = 0.014), tide height (R² = 0.12, p = 0.005), and temperature (R² = 0.15, p = 0.002; Figs. 5D–5F). Among independent categorical variables, host species was significant to Shannon diversity (p = 0.04). With regard to interactions between sampling months, pairwise comparisons were significant for both Shannon and evenness (p < 0.05) for March vs. June and March vs. September, but not for June vs. September. Finally, there were no significant pairwise differences between reefs, in either Shannon diversity or evenness.

Coral bacterial community composition correlated with physical oceanographic habitat

In the beta-diversity analysis, ordination axis 1 of the PCoA explained 89.5% of the variance, and PCoA 2 explained 5.8% of the variance (Fig. 6A). Of correlations between the 12 numeric habitat variables considered, and axis 1 or axis 2 of the ordination, only relative turbidity was significant to axis 2 (R² = 0.27, p < 0.05). In addition, for categorical variables, an ANOSIM showed that the microbial community correlated with host species (R² = 0.38, p = 0.001) and month (R² = 0.21, p = 0.001), but not with reef location (R² = 0.04, p = 0.094). A PERMANOVA also showed the same results for host-species (p > 0.001), month (p > 0.001), and reef location (p = 0.86). Four of the 12 habitat variables had explanatory power for the microbiome count matrix: alongshore currents (p = 0.005), relative turbidity (p = 0.015), significant wave height (p = 0.005), and temperature (p = 0.005). All four of these variables had a VIF (diagonal of the inverse of the correlation coefficient matrix) <2, indicating that they represent effects which are independent of one another (Salmerón, García & García, 2018).

The canonical or CCA ordination was constrained by these four habitat variables, with CCA1 and CCA2 verified as significant by one-way ANOVA. Categorical data were analyzed vs. CCA1 and CCA2: month (R² = 0.52, p = 0.001) and reef (R² = 0.15, p = 0.009) were found to be significant, but the host species was not. Given that host species was not a significant factor to CCA, we analyzed samples from both species together. Accordingly, the CCA ordination plot annotated by month and reef (Fig. 6B) shows that microbial communities cluster primarily by month, with samples collected in March (N = 20) correlating with significant wave height and relative turbidity, and samples collected in September (N = 24) correlating with temperature; samples collected in June (N = 11) correlated with alongshore current speed.

Microbial taxa correlated with changes in waves and temperature

To identify individual microbial taxa that correlated with changes in the physical oceanographic habitat, we used a RFR for variable selection, and then applied a simple linear regression on individual variables. Based on ordistep model outputs, we analyzed four variables with RFR—significant wave height, relative turbidity, alongshore currents, and sea temperature. The RFR classification model was significant for significant wave height (p < 0.001, R² = 0.85), alongshore currents (p < 0.01, R² = 0.60), and sea temperature (p < 0.001, R² = 0.85). The model was not significant for relative turbidity (p = 0.92, R² = 0.001). The five microbial taxa with the highest importance values from each significant RFR model are listed in Table S1. That table also notes which of these “important” taxa were significantly correlated (p < 0.05, R² > 0.2) with their respective habitat variables. For significant wave height, two taxa were significant from the order Flavobacteriales, and one from Marine Group II. For temperature, one taxon was significantly correlated from the order SAR86 (Fig. 7). There were no microbial taxa that were found to be correlated significantly with alongshore currents (i.e., none with p < 0.05 and R² > 0.2).