Surface seawater samples were collected bi-weekly for 1 year (09/2015 to 09/2016; n = 26) from the end of the Pensacola Beach pier (30° 19.640′ N, 87° 08.514′W) which extends approximately 0.3 km into the Northwestern Gulf of Mexico. Samples were transported at in situ temperatures in the dark to the laboratory. Samples were classified into seasons based on the astronomical calendar as follows: fall was defined as September 23rd – December 20th, winter was defined as December 21st – March 18th, spring was defined as March 19th – June 19th, and summer was defined as June 20th – September 21st.

In situ seawater temperature, salinity, NO₃^– + NO₂^– concentrations, NH₃ concentrations, total Kjeldahl nitrogen (TKN) concentrations, orthophosphate concentrations, and total phosphorus (TP) concentrations were measured for each water sample (Table 1). NO₃^– + NO₂^– and orthophosphate concentrations were measured on a Lachat QuickChem 8500 using EPA standard method 353.2 (1993) for NO₃^– + NO₂^– concentrations and 365.1 (Kopp, 1979) for orthophosphate concentrations. NH₃, TKN, and TP concentrations were measured using EPA standard method 350.1 (1993), EPA standard method 351.2 (1993), and EPA standard method 365.4 (1974), respectively.

Samples for bacterial counts were preserved with 0.2 μm filtered, buffered formalin. Preserved samples were stained with 4′,6-diamidino-2-phenylindole, dihydrochloride (DAPI) using the method of Porter and Feig (1980) and counted using an epifluorescence microscope. Samples for chlorophyll a concentrations were filtered (200 mL) onto 25 mm GF/F filters in triplicate, extracted in 90% acetone overnight, and measured fluorometrically (Turner Trilogy Laboratory Fluorometer) using a standard curve (Welschmeyer, 1994).

Aged Gulf of Mexico seawater (collected ∼40 km offshore Pensacola, FL and >3 years old) was filtered using a 0.2 μm pore-size polycarbonate filter (Millipore). Twenty-five milliliter of filtered seawater was aliquoted into 35 mL Teflon bottles (Nalgene FEP). Bottles containing filtered seawater were pasteurized at 70°C for 2–4 h. Once cooled, pasteurized seawater was amended with surrogate crude oil (Pelz et al., 2011) to contain a final volume of 2% crude oil. Bottles were incubated in a 20°C temperature-controlled water table (Fisherbrand, Isotemp 4100) under full solar exposure for 5 days during the summer on the roof of the Environmental Sciences building at the University of West Florida (Vaughan et al., 2016). Bottles were shaken twice per day and returned to the water table. Samples were pooled into a separatory funnel, and the aqueous fraction (i.e., the WAF) was collected and transferred to 20 mL scintillation vials in 10 mL aliquots and stored frozen at −20°C. Here, we make the assumption that because we used filtered, pasteurized aged seawater, any organic carbons in the seawater were dissolved hydrocarbons. Therefore, total organic carbon (TOC) was used as a proxy for hydrocarbon concentrations. TOC concentrations of the WAF were determined as non-purgeable organic carbon with a Shimadzu TPC-VVSN Analyzer using standard method 5310 (2018). WAFs contained an average TOC concentration of 68.9 ppm (SD = 1.1, n = 4). Fractions were thawed prior to each sensitivity assay. In this way, each bi-weekly water sample was exposed to the same WAF for the duration of the project.

Bacterial production was estimated through incorporation of ³H-leucine. The standard WAF exposure consisted of a dose response curve of 0, 0.5, 1, 2.5, 5, and 10% v/v WAF (final concentration amended to the seawater sample) for each time point. WAF was aliquoted into 4 replicate 5 mL polystyrene snap-cap tubes for each treatment. Four replicate controls received 100 μL of filtered seawater (0.2 μm pore size syringe filter). Seawater was amended with ³H-leucine (52.9 Ci mmol^–1 PerkinElmer, Bridgeport, CT, USA) to a final concentration of 10 nM; 3.1 mL of labeled seawater was added to each tube. Samples were capped, mixed, and incubated in the dark at in situ temperature for 4 h. To terminate leucine incorporation, triplicate 1.0 mL subsamples were removed from each snap cap tube and placed into 2 mL microfuge tubes containing 50 μL of 100% trichloroacetic acid. Samples were processed following the procedure of Smith and Azam (1992). Liquid scintillation counting using a Packard Tri-Carb 2900 was performed to determine ³H-leucine incorporation in samples.

Primary production was determined through the fixation of ¹⁴C-bicarbonate under increasing light exposures (PI curves) (Matrai et al., 1995) using a modified photosynthetron. Triplicate dose response curves of 0, 0.5, 1, 2.5, 5, and 10% WAF treatment were generated for each time point. WAF and 33 mL of seawater were aliquoted into 50 mL conical centrifuge tubes for each treatment which were then amended with ¹⁴C-bicarbonate (2 μCi/mL). Triplicate controls received 1 mL of filtered seawater (0.2 μm pore size syringe filter) and 33 mL of seawater that were then amended with ¹⁴C-bicarbonate. Each treatment was aliquoted into eight, 4 mL snap cap tubes that were incubated at in situ temperature for 4 h in a photosynthesis irradiance incubator. The photosynthetically active radiation of each position was measured with a quantum scalar laboratory radiometer (Biospherical Instruments Inc.) and recorded. Irradiances for the eight positions were ∼180, 90, 50, 30, 20, 15, 10, and 7 μE cm^–2 for each treatment and time point. Lower irradiances were used to maximize the number of points within the linear part of the curve for calculating inhibition of primary production (see below). Fixed ¹⁴C was determined after overnight acidification of the samples via liquid scintillation counting using a Packard Tri-Carb 2900. Photosynthetic efficiency was determined as described by Matrai et al. (1995).

Bacterial production inhibition and primary production inhibition were determined using dose response curves (0, 0.5, 1, 2.5, 5, and 10% WAF). For bacterial production inhibition, the average disintegrations per minute (DPM) of each treatment replicate were expressed as a percent of the average DPM of the control (0% WAF). The percent control was log-transformed, and a linear regression of the log-transformed percent of the control versus percent WAF concentration was performed in Excel. The positive slope of the regression for each experiment was then used as a quantitative value for degree of inhibition and is referred to as ‘‘bacterial production inhibition’’ throughout. For primary production inhibition, PI curves were constructed for each control and each treatment replicate. Linear regressions were performed on the linear part of the curve in Kaleidagraph. The slope of each treatment replicate was expressed as a percent of the average slope of the control. The percent control was log-transformed, and a linear regression of the log-transformed percent of the control versus percent WAF concentration was performed in Excel. The positive slope of the regression for each experiment was then used as a quantitative value for degree of inhibition and is referred to as ‘‘primary production inhibition’’ throughout. The standard error of the regression slopes was obtained from Kaleidagraph and was used to assess the uncertainty in bacterial production inhibition and primary production inhibition. Raw data and example calculations are available on GitHub.¹

To determine how temporal variation in the background bacterial community contributed to temporal changes in the inhibition of bacterial production, samples for DNA extraction were collected bi-weekly. Samples were collected by filtering 2 L of seawater through three 0.22 μm GPWP (MilliporeSigma) filters. Samples were stored at −80°C until further processing. Half of each previously frozen filter and 250 μL of extraction buffer (Milli-Q water, 5 mM EDTA, 25 mM Tris, and 50 mM glucose) were added to tubes containing a mixture of 0.1 mm silica and 0.5 mm glass beads. Filters were ground using a sterile plastic pestle and were homogenized for two 1-min cycles at 2,000 rpm. Samples were cooled to −80°C, heated to 80°C for 10 min, cooled to −80°C again, and then brought to room temperature. Lysozyme (final concentration = 1.5 mg mL^–1) was added to each sample, and samples were incubated for 90 min at 37°C. Proteinase K (final concentration = 3 mg mL^–1) was added to each sample, and samples were incubated for 90 min at 50°C. Sodium chloride (final concentration = 0.5 M) and M1 buffer from the Omega E.Z.N.A. Mollusc DNA Kit were added to each sample. The Omega E.Z.N.A. Mollusc DNA Kit was used to wash and collect purified DNA. Extracted DNA was quantified using a NanoDrop spectrophotometer and were stored at −20°C. DNA extracts were sent to the University of Illinois at Chicago’s Sequencing Core for amplification and sequencing. The V4–V5 region of the 16S rRNA gene was amplified using the 515F/926R universal primers (Needham and Fuhrman, 2016). Amplicons were pair-end sequenced (2 ×300) with the MiSeq Illumina platform. Sequence files are available at the NCBI Sequence Read Archive under BioProject ID: PRJNA894536. Accession numbers for each sample are reported in Supplementary Table 1.

Forward and reverse primers were removed using cutadapt (Martin, 2011) in QIIME2 (Bolyen et al., 2019). Forward and reverse reads were quality filtered with fastq-mcf (Aronesty, 2013). A window-size of 10 was used to calculate mean quality score. Reads were truncated when the mean quality score was less than 20. After trimming, reads that were shorter than the minimum length threshold of 150 bp were removed. Reads that contained N-calls were also removed. Forward and reverse reads were merged based on a minimum overlap threshold of 10 bp, minimum merge length threshold of 350 bp, and number of maximum differences of 5 bp allowed in the overlapping region using usearch (Edgar, 2010). Final trimming, quality filtering, clustering of amplicons, and removal of chimeras was performed using DADA2 (Callahan et al., 2016) in QIIME2. The merged reads were trimmed to a length threshold of 365 bp to maintain alignment. Reads that matched to the PhiX genome or that contained more than 3 expected errors were removed. The error model was trained using a minimum of 800,000 reads. Samples were then dereplicated, reads were clustered into amplicon sequence variants (ASVs), and chimeric ASVs were removed using a consensus procedure. QIIME2 artifacts generated from this bioinformatics workflow are available on GitHub² and the contents of each artifact are described in Supplementary Table 2.

Taxonomic, diversity, and statistical analyses were performed in R (version 3.5.1) (R Core Team, 2018). All R code is available on GitHub.³ All colors used in figures were checked for accessibility using Adobe Color’s “Color Blind Safe” Accessibility Tool. Taxonomy was assigned (“assignTaxonomy” function; dada2 package; Callahan et al., 2016) using RDP’s Naïve Bayesian classifier (Wang et al., 2007) and the SILVA 138 reference database (Quast et al., 2012). ASVs matching to eukaryotes, archaea, or that were unassigned at the bacterial Kingdom or Phylum level were removed from subsequent taxonomy and diversity analyses. Temporal trends in the 25 most abundant genera were examined across the time series. Counts of all genera in each sample are available in Supplementary Table 3.

For alpha-diversity analyses, sequencing depth was normalized by rarefying each sample to 16000 sequences (“rrarefy” function; vegan package; Oksanen et al., 2022). Alpha-diversity was calculated using the Shannon index (“diversity” function; vegan package) which is a composite metric of richness and evenness. Therefore, alpha-diversity was also calculated as richness (i.e., the number of ASVs in each sample; “richness” function; microbiome package) and evenness using Pielou’s index. Pielou’s index was calculated as H/log(S) where H is the Shannon index and S is richness. A one-way analysis of variance (ANOVA) was performed to determine if the three metrics of alpha-diversity varied by season. Prior to performing ANOVAs, assumptions of normality and homogeneity were checked using the Shapiro–Wilk test and Levene’s test, respectively. Tukey’s HSD was then performed to identify differences in alpha-diversity between seasons. To visualize similarities and differences in ASV composition by season, a Euler diagram was constructed using the rarefied ASV count table (“ps_euler” function; MicEco package; Russel, 2021). Additionally, the nestedness and turnover of ASVs between seasons were calculated (“beta.temp” function; betapart package; Baselga et al., 2022) to better understand ecological succession throughout the year. Lastly, correlation analysis of the three metrics of alpha-diversity with temperature and TKN were performed (“cor.test” function, method = “pearson”).

For beta-diversity analyses, sequences were normalized using the variance stabilizing transformation (“varianceStabilizingTransformation” function; DESeq2 package; Love et al., 2014). Principal component analysis (PCA) was conducted on the normalized ASV table to visualize differences in microbial community structure by season (“ordinate” function; phyloseq package; McMurdie and Holmes, 2013). PC1 was extracted and correlation tests of temperature and TKN with PC1 were conducted (“cor.test” function; method = “pearson”). A PCA was also conducted on the environmental variables (“rda” function; vegan package), loadings were extracted, and the environmental loadings were overlaid as vectors onto the community ordination plot to create a biplot. Permutational ANOVA was performed to determine if community structure varied by season (“adonis” function; vegan package). Prior to performing the permutational ANOVA, the assumption of homogeneity of dispersion among groups was checked (“betadisper” function; vegan package). To identify which seasons significantly differed in their community structures, pairwise permutation multivariate ANOVA was performed (“pairwise.perm.manova” function; RVAideMemoire package; Hervé, 2020). To prevent inflation of Type I error rate due to multiple comparisons, the Hochberg method was applied to calculate adjusted p-values.

To determine if changes in the inhibition of primary and bacterial production varied by season, ANOVAs were performed, as described above. Simple linear regressions were performed to identify significant linear relationships between inhibition and biological/environmental variables. A mantel test was performed to determine if variation in the inhibition of bacterial production was correlated with changes in microbial community structure (“mantel” function; vegan package). The mantel test was performed using 10,000 permutations on two normalized dissimilarity matrices: (1) a Euclidean distance matrix of inhibition of bacterial production and (2) a Bray–Curtis distance matrix constructed from the rarefied ASV count table (“vegdist” function; vegan package).

Surface coastal waters exhibited a strong temporal temperature trend but lacked temporal variability in salinity and nutrient concentrations (Figure 1 and Table 1). Sea surface temperatures during the summer were significantly higher than all other seasons (p < 0.05) with a maximum value of 28.2°C, while temperatures during the winter were significantly colder than all other seasons (p < 0.01) with a minimum value of 13.6°C (Figure 1A). Salinity ranged from 27 to 37 with no significant differences between seasons but with higher variability observed at the end of spring and during summer due to increased rainfall (Figure 1B). NH₃ and TP concentrations were below the minimum detection limits (MDL). NO₃^– + NO₂^–, TKN, and orthophosphate concentrations were above the MDL but exhibited minimal variation and had no temporal trends (p > 0.05) (Figures 1C-E).

Patterns in bacterial diversity, community structure, and taxonomic composition were examined by sequencing the V4–V5 region of the 16S rRNA gene. Alpha-diversity was calculated using three different metrics: the Shannon index, ASV richness, and Pielou’s evenness index (Supplementary Table 1). The Shannon index and ASV richness did not exhibit any temporal trends (Supplementary Figures 1A, B), but evenness did significantly vary between winter and summer with winter having higher evenness than summer (p < 0.05) (Figure 2A and Supplementary Figure 1C). While alpha-diversity varied minimally across seasons, there were significant differences in community structure by season. Although PCA indicates that there is some overlap in community structure (Figure 2B), PERMANOVA confirmed that they are statistically distinct from each other according to season (p < 0.05). Partitioning beta-diversity into turnover and nestedness components revealed that the differences in community structure are largely due to turnover, with the turnover component accounting for 84–93% of the dissimilarity between seasons. This analysis is supported by a Euler diagram which shows that there are larger proportions of ASVs unique to each season compared to the proportion of ASVs shared between subsequent seasons (Figure 2C). There are also noticeable transitions between seasons in community composition. Candidatus Actinomarina, SAR11 clade Ia, HIMB11, NS2b marine group, NS4 marine group, NS5 marine group, OM60 (NOR5) clade, and Synechococcus dominated the bacterial community throughout the year (>1% average relative abundance) (Figure 3A), but their relative abundances shifted across seasons (Figure 3B). Within each season, particular genera became dominant (>1% average relative abundance). During the spring, Blastopirellula, Candidatus Aquiluna, and Cyanobium were dominant genera, while during the summer, Balneola and Cyanobium were dominant. The largest shifts in dominant genera occurred during the fall and winter with Cyanobium, Formosa, MB11C04 marine group, and SAR11 clade Ib becoming dominant in the fall, and Ascidiaceihabitans, Blastopirellula, Formosa, and the OM43 clade being dominant in the winter. Thus, we observed that these bacterial communities have distinct structures and compositional shifts according to season.

To determine the temporal effect that exposure to water-soluble crude oil components had on the inhibition of primary and bacterial production, triplicate dose response curves (0, 0.5, 1, 2.5, 5, and 10% WAFs) were generated bi-weekly for 1 year. Exposure to WAFs led to inhibition of primary production across all seasons (Figure 4A). Inhibition of primary production was significantly lower during the winter compared to the spring and summer (p < 0.05) (Figure 4A). Inhibition of primary production exhibited strong negative relationships with chlorophyll a concentrations (p < 0.0001; Adj. R² = 0.765) (Figure 4B). Additionally, inhibition of primary production had a positive relationship with temperature (p = 0.0002; Adj. R² = 0.444) (Figure 4C) but no relationship with NO₃^– + NO₂^– concentrations (p = 0.08) or with orthophosphate concentrations (p = 0.51). However, inhibition of primary production had a negative relationship with TKN (p = 0.002; Adj. R² = 0.332) (Figure 4D). These results indicate that water soluble crude oil constituents have the strongest inhibitory effect on primary production when surface coastal waters are warm, when TKN concentrations are low, or when phytoplankton abundance is low.

Exposure to WAFs also led to inhibition of bacterial production across all seasons. Inhibition of bacterial production was significantly lower during winter compared to all other seasons (p < 0.05) (Figure 5A). There was no relationship of the inhibition of bacterial production with absolute bacterial abundance (p = 0.34), but there were weak, negative relationships of bacterial production inhibition with bacterial alpha-diversity (Shannon index: p = 0.028, Adj. R² = 0.152; Pielou’s evenness index: p = 0.012, Adj. R² = 0.205) (Figure 5B). Additionally, there was a weak, positive correlation of bacterial production inhibition with bacterial community structure (p = 0.002; mantel r = 0.28), indicating that the severity of inhibition is partially dependent on the background bacterial community. Although inhibition of bacterial production had no relationship with inorganic nutrient concentrations (NO₃^– + NO₂^– and orthophosphate) (p > 0.10), it did have a weak, negative relationship with TKN (p = 0.032; Adj. R² = 0.143) (Figure 5D) as well as a strong, positive relationship with temperature (p < 0.0001; Adj. R² = 0.64) (Figure 5C). These results are similar to what was seen for inhibition of primary production and indicate that water soluble crude oil constituents are most inhibitory to bacterial production when temperatures are high, when TKN concentrations are low, or when alpha-diversity is low.

Additionally, there were strong relationships between environmental conditions and the bacterial community which corresponded to the degree of bacterial production inhibition. Temperature had a strong positive correlation with PC1 (p < 0.0001, r = 0.756) (Figure 2B), a moderate negative correlation with the Shannon index (p = 0.028, r = −0.431), and a moderate negative correlation with Pielou’s evenness index (p = 0.01, r = −0.498), while TKN had a moderate negative correlation with PC1 (p = 0.025, r = −0.438), a moderate positive correlation with the Shannon index (p = 0.048, r = 0.392), and a moderate positive correlation with ASV richness (p = 0.035, r = 0.415). Also, winter, which was the season with the coldest temperatures, had more unique ASVs compared to warmer seasons (Figure 2C). This demonstrates that changes in temperature strongly corresponded with changes in bacterial community structure and with variability in alpha-diversity through changes in evenness, while TKN had a weaker correspondence with bacterial community structure and with variability in alpha-diversity through changes in richness. Combined, the data showed that lower temperatures corresponded to a bacterial community with more unique ASVs, higher evenness, and less inhibition of bacterial production, while higher concentrations of TKN corresponded to a bacterial community with higher richness and lower inhibition of bacterial production.

These results partially supported our hypothesis that there is a temporal response in inhibition due to acute exposure to water-soluble components of crude oil and that this response would be primarily driven by high temperature and low inorganic nutrient availability. We observed that exposure led to inhibition of primary and bacterial production across all seasons. The impact of exposure was found to be highest during warm months. Interestingly, we did not observe any effects of inorganic nutrient concentrations (NO₃^– + NO₂^– and orthophosphate), likely because there was low variability throughout the year. Surprisingly, there was a negative relationship between inhibition and TKN concentrations for both primary and bacterial production. This relationship may suggest links between human inputs of nitrogen, microbial community response, and the inhibition of production due to water soluble crude oil constituents. TKN primarily comes from human inputs and is a composite measurement of organic nitrogen, ammonia, and ammonium. When TKN concentrations increase in coastal waters, some heterotrophs and phytoplankton can use these compounds for their growth. Thus, increases in TKN may alter the community and buffer the impacts of water-soluble crude oil components. This is further supported by the relationships we observed between TKN, microbial community structure, and microbial community diversity. Additionally, we observed that changes in temperature strongly corresponded with changes in the bacterial community which then corresponded to the severity of bacterial production inhibition. For example, we observed that the impact of exposure to WAF was partially dependent on the background phytoplankton and bacterial communities as seen through the positive relationship of primary production inhibition with chlorophyll a concentrations (i.e., phytoplankton abundance), the negative relationship of bacterial production inhibition with alpha-diversity, and the positive relationship of bacterial production inhibition with bacterial community structure. Therefore, it appears that there are tight linkages between environmental factors, the microbial community, and the response to WAF. Combined, these results suggest that changes in environmental factors, such as temperature and TKN, corresponded with changes in the microbial community, which thus resulted in varying degrees of production inhibition when exposed to WAF.