Microbial Functional Responses in Marine Biofilms Exposed to Deepwater Horizon Spill Contaminants

Marine biofilms are essential biological components that transform built structures into artificial reefs. Anthropogenic contaminants released into the marine environment, such as crude oil and chemical dispersant from an oil spill, may disrupt the diversity and function of these foundational biofilms. To investigate the response of marine biofilm microbiomes from distinct environments to contaminants and to address microbial functional response, biofilm metagenomes were analyzed from two short-term microcosms, one using surface seawater (SSW) and the other using deep seawater (DSW). Following exposure to crude oil, chemical dispersant, and dispersed oil, taxonomically distinct communities were observed between microcosms from different source water challenged with the same contaminants and higher Shannon diversity was observed in SSW metagenomes. Marinobacter, Colwellia, Marinomonas, and Pseudoalteromonas phylotypes contributed to driving community differences between SSW and DSW. SSW metagenomes were dominated by Rhodobacteraceae, known biofilm-formers, and DSW metagenomes had the highest abundance of Marinobacter, associated with hydrocarbon degradation and biofilm formation. Association of source water metadata with treatment groups revealed that control biofilms (no contaminant) harbor the highest percentage of significant KEGG orthologs (KOs). While 70% functional similarity was observed among all metagenomes from both experiments, functional differences between SSW and DSW metagenomes were driven primarily by membrane transport KOs, while functional similarities were attributed to translation and signaling and cellular process KOs. Oil and dispersant metagenomes were 90% similar to each other in their respective experiments, which provides evidence of functional redundancy in these microbiomes. When interrogating microbial functional redundancy, it is crucial to consider how composition and function evolve in tandem when assessing functional responses to changing environmental conditions within marine biofilms. This study may have implications for future oil spill mitigation strategies at the surface and at depth and also provides information about the microbiome functional responses of biofilms on steel structures in the marine built environment.


INTRODUCTION
The marine built environment sustains the vast biodiversity of life present in the soft bottom habitat of the seafloor (Gage, 2004). Historic shipwrecks and other submerged structures, such as oil pipelines and drilling rigs, have the potential to transform into artificial reefs on the seabed, but marine biofilms are the biological foundation of these ecosystems. Pioneer microorganisms immediately initiate biofilm formation on submerged structures through surface adhesion, release of extracellular polymers, and recruitment of diverse microorganisms to the substrate. Once established, mature biofilms release chemical signals that serve as settlement cues for other micro-and macrofauna (Hadfield, 2011;Dobretsov and Rittschof, 2020). Exposure to aquatic pollutants, such as crude oil and dispersant, impact biofilm composition and function Mugge et al., 2019a,b), with potential downstream effects on the diversity of higher trophic level organisms recruiting to artificial reefs.
Following the 2010 Deepwater Horizon spill, oil plumes were detected at 1,000-1,300 m water depth at distances up to 16 km surrounding the 1,500-m wellhead until it was capped after 85 days (Valentine et al., 2014;Gao et al., 2015). An estimated 4-31% of total oil released from the wellhead was sequestered to the seafloor at estimated depths of 1,300-1,700 m on the continental slope, based on the detection of oil biomarkers in seafloor sediment (Valentine et al., 2014). To expedite bioremediation of the oil, an estimated total of 7 million liters of the chemical dispersant COREXIT EC9500A was applied at the surface and at the wellhead (Kleindienst et al., 2016). Within the acute spill footprint (∼16 km), where spill contaminants were deposited on the seafloor, a 5-cm-thick reddish-brown surface layer, high porosity (∼90%), and a microbiome dominated by bacteria with the functional potential for hydrocarbon degradation were observed in collected sediment cores . Sediment, biofilm, and water microbiomes from samples collected around historic shipwrecks within the acute spill footprint indicate that residual spill contaminants may induce changes to foundational biofilm microbiomes on artificial reefs, which may negatively impact preservation of these structures (Mugge et al., 2019a).
Access to the deep sea is limited, and therefore, biofilm formation on submerged structures and the impacts of environmental contaminants on biofilm composition and function are not well understood. Carbon steel disks (CSDs) are commonly used as biofilm recruitment surfaces in simulated experiments since the majority of marine infrastructure is constructed from carbon steel (McBeth et al., 2011). Salerno et al. (2018) examined biofilm response to spill contaminants by placing CSDs in microcosms filled with surface seawater (SSW) and amended with environmentally relevant concentrations of crude oil, dispersant, and dispersed oil to discover how submerged metal infrastructure may be compromised by oil spills. Mugge et al. (2019b) replicated this study using deep seawater (DSW), and the results from both studies reveal that exposure to spill contaminants immediately induces changes to the biodiversity and function of bacteria inhabiting marine biofilms. Investigating the metagenomes sequenced from the same biofilm samples may give insight into the functional redundancy of marine biofilms exposed to spill contaminants.
Functional redundancy is the ecological theory that taxonomically distinct species are capable of similar metabolic roles within communities and ecosystems (Rosenfeld, 2002); this theory has been extended to microbial assemblages (Allison and Martiny, 2008;Burke et al., 2011). From an evolutionary perspective, functional redundancy benefits ecosystems by allowing the diversity and abundance of species to be interchanged with little or no impact to the function of the ecosystem. However, the generality of functional redundancy may not hold true for microbial assemblages under changing environmental conditions and altering the composition of a microbiome may also alter its functions (Fetzer et al., 2015;Galand et al., 2018).
Marine microorganisms play vital roles in global carbon and nutrient cycling, feedback mechanisms, and organic matter fluxes to the seafloor (Bienhold et al., 2016). Sequencing technology has provided insight into the diversity, abundance, and functional capabilities of deep-sea bacterial communities, which are largely unclassified. In comparison to seawater and sediment habitats, marine biofilm microbiomes are distinct due to surface associations that increase the relative abundance of novel taxa and functional genes (Zhang et al., 2019). Additionally, environmental parameters such as pH, dissolved oxygen, light, and nutrient availability play major roles in determining the structure of marine biofilm microbiomes (Orland et al., 2019). While microbial communities are sensitive to disturbances, they can also be resistant and resilient to environmental changes (Allison and Martiny, 2008). The impacts of anthropogenic environmental disturbances on the diversity and functional redundancy of marine biofilms, however, is understudied.
Studying marine biofilm recruitment, and response to contaminants, over time may provide insight into how microbiomes on built structures respond to environmental disturbance. Accordingly, our study quantified changes in developing marine biofilm microbiomes following exposure to Deepwater Horizon oil spill contaminants, using metagenomics to determine if functional similarity occurs between developing biofilms sourced from SSW and DSW. We hypothesized that biofilms recruited to carbon steel will harbor distinct microbial taxa due to different seawater source populations but will exhibit similar functional responses to contaminant exposure.

Laboratory Microcosms
Two separate laboratory experiments were carried out using different source water to recruit microbiomes to C1020 mild CSDs to simulate biofilm growth on built marine infrastructure. CSDs (Metal Samples Company, Munford, AL, United States) measured 1.59 cm in diameter and 0.016 cm thick with a mill finish and were mounted using EpoThin2 epoxy (Buehler, Lake Bluff, IL, United States) with one face exposed, following the protocol of Lee et al. (2004Lee et al. ( , 2005. Inside microcosm tanks, CSDs were housed in short towers (16.51 cm outer diameter by 27.94 cm high) constructed from polyvinyl chloride (PVC) pipe (Commercial Industrial Supply, Rockhill, SC, United States) and secured using 100% waterproof silicon sealant (Silicone Max, all purpose). Mugge et al. (2019b) provides further description of the experimental design. To collect biofilm samples, CSDs were sacrificed biweekly (n = 2) and stored upright in sterile plastic collection cups at −80 • C. Week 2 samples were collected before contaminant additions of crude oil, dispersant, and dispersed oil to microcosms, and all other samples were collected following exposure to contaminants.
For the SSW experiment, water was collected from the US Naval Research Laboratory Key West pier (Key West, FL, United States) and stored for 1 week in a cold room maintained at 4 • C (±0.4 • C). Four, 25-L experimental tanks were each filled with 16 L of seawater and kept covered to exclude ambient light, and on stir plates for continuous water circulation. Salinity (37-38 PSU), temperature (3.4-4.4 • C), and dissolved oxygen (86.5-97.6%) were measured weekly for 16 weeks (Supplementary Table 1). After a 2-week acclimation period, each of the four tanks were amended with a single contaminant, as follows: 1 control tank (no contaminant), 1 oil tank (Louisiana sweet crude oil, 5 mg/L final concentration), 1 chemical dispersant tank (Corexit EC9500, 0.05 mg/L final concentration), and 1 dispersed oil tank (5 mg/L crude oil + 0.05 mg/L dispersant final concentration). Contaminant concentrations were calculated from a review of published data and from data obtained from the NOAA Office of Response and Restoration's Environmental Response Management Application (ERMA) database. 1 For the DSW experiment, water was collected from ∼10 m above the seafloor at the Anona shipwreck in the Viosca Knoll lease area in the Gulf of Mexico at an approximate depth of 1,200 m using twelve 12-L Niskin bottles attached to a conductivity, temperature, and depth (CTD) rosette on R/V Point Sur cruise PS17-26. Additional site information is provided in Hamdan et al. (2018). Seawater was collected into 13 autoclaved, 20-L polycarbonate carboys and stored at 4 • C (±0.5 • C) in a cold room. Each of the 12 experimental tanks (three replicate tanks for each of three contaminants, and three replicate control tanks) were filled with 16 L of seawater and maintained in the cold room, covered to exclude ambient light and on stir plates set to ∼620 rpm for continuous circulation. Salinity (35-39 PSU), temperature (4.8-5.6 • C), and dissolved oxygen (77.3-94.8%) were measured weekly during the 14-week experiment (Supplementary Table 1). Water-accommodated fractions (WAFs), which have been implemented in aquatic toxicity studies, were used to simulate oil spill conditions in the tanks (Forth et al., 2017). Following 2 weeks of acclimation, replicate (n = 3) tanks were amended with crude oil using a highenergy water-accommodated fraction (HEWAF, 80 mL/tank; 5 mg/L final concentration of oil), diluted chemical dispersant (Corexit EC9500, 0.05 mg/L final concentration), dispersed oil using a chemically enhanced water-accommodated fraction (CEWAF, 80 mL/tank; 5 mg/L final concentration of oil and 1 https://erma.noaa.gov/gulfofmexico/erma.html 0.05 mg/L final concentration of Corexit), or no treatment (control). HEWAF and CEWAF were made following the protocols as previously described (Mugge et al., 2019b).

DNA Extraction and Quantitation
Biofilms, which include corrosion products, were removed from CSDs and collected into pre-weighed lysing matrix E tubes for DNA extraction (MP Biomedicals, LLC, Santa Ana, CA, United States) following a modified version of the manufacturer's protocol of the BIO 101 FastDNA Spin Kit as previously described (Hamdan et al., 2013). Prior to DNA extraction from corrosion biofilms, bovine serum albumin (Promega, Madison, WI, United States) was added as a chelating agent to prevent inhibition of PCR amplification by iron ions Mugge et al., 2019b). Extracted DNA (150 µL) was stored frozen at −80 • C. Total extracted genomic DNA (10 µL) was quantitated on a Qubit 2.0 Fluorometric Quantitation system (Invitrogen, Carlsbad, CA, United States) and a NanoDrop spectrophotometer (Thermo Fisher Scientific, Waltham, MA, United States) was used to assess purity. 16S rRNA amplicon sequence taxonomy data and metagenome data were obtained from the same samples from each respective microcosm. The metagenome data are presented in this study while the 16S rRNA data are published elsewhere Mugge et al., 2019b).

Metagenome Analysis
Prior to metagenome sequencing at the Integrated Microbiome Resource (IMR) facility at Dalhousie University (Halifax, Nova Scotia, Canada), replicate samples from the DSW microcosm were pooled to minimize cost and for comparability to SSW metagenomes. DSW samples included one pooled sample from week 2 prior to contaminant addition, which served as the pretreatment sample, and pooled samples from replicate contaminant tanks for weeks 4, 8, 12, and 14 for a total of 17 DSW metagenomes. Sample pooling was based on statistical analysis of Shannon diversity and total assigned reads from replicate 16S rRNA gene amplicon sequences presented elsewhere (Mugge et al., 2019b;Supplementary Table 2). The SSW microcosm did not have replicate tanks and thus all samples were sequenced and only metagenomes from weeks 2, 4, 8, 12, and 14 were included for comparability with DSW metagenomes. Metagenome Libraries were prepared with the Nextera XT Library Preparation Kit (Illumina Inc.), a PCR-based library preparation approach, modified to use the Just-a-Plate TM 96 PCR Purification and Normalization Kit (Charm Biotech) as described in Salerno et al. (2018). Sequencing was performed on an Illumina NextSeq 550, generating 150-bp, paired-end sequences. Metagenome sequences were submitted to the NCBI SRA under BioProject PRJNA623475.

Metagenome Bioinformatics Analysis
Metagenomes were analyzed in accordance with workflows available in the Microbiome Helper repository (Comeau et al., 2017). Raw sequences were quality filtered (Andrews, 2010), assembled into paired-end reads (Zhang et al., 2013), trimmed (Bolger et al., 2014), and screened for contaminants (Langmead and Salzberg, 2012). MetaPHlAn2 was used to assign taxonomy based on unique clade-specific marker genes that unambiguously assign reads to microbial clades (Truong et al., 2015). OTU abundance was extracted from MetaPHlAn2 data, averaged by treatment within each microcosm, and relative abundance bar charts were created using Microsoft Excel, and excluded any taxonomic group that represented less than 2% of the total population per treatment. The relative abundance of all taxa at the OTU level for both experiments across all time points is included in Supplementary Table 3. HUMAnN2 (Franzosa et al., 2018) was used to profile gene families based on the UniRef90 database, and these features were renormalized to copies per million (CPM), and regrouped and renamed by using the Kyoto Encyclopedia of Genes and Genomes to assign KEGG orthologs (KOs; Kanehisa and Goto, 2000).

Statistical Analysis
Taxa and gene abundance data from all SSW week 2 samples (n = 4) were each pooled and averaged as a SSW pretreatment sample to create an analogous comparison with the DSW pretreatment sample. Species-level annotations were extracted from the MetaPHlAn2 results of SSW and DSW microcosms, and an operational taxonomic unit (OTU) table was imported into PRIMER (v.6.13, PRIMER-E Ltd., Plymouth, United Kingdom). Functional profiles of metagenomes were analyzed using HUMAnN2, which resulted in abundance tables of gene families. The default output unit from HUMAnN2, reads per kilobase (RPK), accounts for gene length but was renormalized to CPM to account for sequencing depth using the humann2_renorm_table script. Following renormalization, a combined gene abundance table for SSW and DSW microcosms was imported to PRIMER. For taxonomy and functional profiles, Bray-Curtis dissimilarities were calculated, hierarchical clustering (CLUSTER) illustrated similarity based on the group average linkage, and nonmetric multidimensional scaling (NMDS) plots were constructed to visualize community and functional differences in threedimensional space. Following a permutational analysis of variance (PERMANOVA) main test on each profile, separate pairwise PERMANOVA tests were run on the interaction term "microcosm × treatment" to test for significant differences between the same treatments within each microcosm. Shannon diversity of both taxonomy and functional profiles was calculated using Quantitative Insights Into Microbial Ecology (QIIME). Assigned reads were correlated with observed OTUs from all metagenomes, and Excel was used to calculate the correlation coefficient, t score, and p value.
To observe functional responses to treatments, pretreatment samples were excluded in downstream functional comparative analyses. Mapping of the merged HUMAnN2 gene families file of all metagenomes to KOs resulted in 57% unmapped reads and 30% ungrouped reads, leaving 13% of reads mapped to KOs. These KOs were renormalized by CPM, averaged by treatment for each microcosm, and then ranked by overall abundance. Averaged values were input to PRIMER to visualize major differences in functional profiles by running a CLUSTER analysis. A similarity percentage (SIMPER) analysis was run to determine functional categories driving similarities and differences between metagenomes. The top 20% of assigned reads (top 149 abundant KOs) were assigned to a functional hierarchy using the hierarchical structure of the KEGG Brite database. Each hierarchy classification was summed and a bar chart was created using Microsoft Excel. To interrogate the dataset for functional similarity, MaAsLin2 (Mallick et al., 2017) was used to build a generalized linear model using a minimum abundance threshold of 0.001 and a significance of p < 0.05. This model was used to associate source water metadata with control and treatment groups and yielded an output of significantly different genes per treatment between SSW and DSW experiments based on relative abundance. Percentages of significant genes were calculated from filtered data and outputs of the model within MaAsLin2. The output of the analysis is the percentage of genes with a statistically significant differential abundance pattern, and thus, a higher percentage illustrates dissimilarity and a lower percentage corresponds to greater similarity.

RESULTS
The SSW experiment had a larger number of observed OTUs (average 38) across all treatments compared to the DSW experiment (average 24) (Supplementary Table 2). Shannon diversity was greater, on average, in the SSW experiment (4.01), with the greatest average diversity in the oil treatment (4.17). Across all treatments, SSW metagenomes had a greater number of total assigned reads from gene families to KOs (900,000-1,000,000) as compared to DSW metagenomes (200,000-300,000) (Supplementary Table 2). Comparing taxonomic richness (observed OTUs) against functional richness (assigned reads to KOs) revealed a significant positive correlation (R 2 = 0.77, P < 0.001) in both DSW and SSW metagenomes (Figure 1).
An NMDS plot based on Bray-Curtis dissimilarities of taxa abundance data revealed separation between communities in the SSW and DSW experiments across all treatments, highlighting distinct microbial communities throughout the course of the experiments (Figure 2A). Although taxonomy in the DSW pretreatment sample displayed high separation from all other samples, the SSW pretreatment sample clustered at 70% similarity with one SSW dispersant sample. The cluster of DSW samples is visualized in an NMDS subset in Supplementary Figure 1A. A pairwise PERMANOVA comparing taxonomy data in the same treatments in SSW and DSW experiments returned statistically significant p values (P < 0.05), with control and dispersant communities exhibiting highest percent similarities (10.55%, 10.33%) (Supplementary Table 4).
An NMDS ordination revealed high separation between SSW and DSW metagenomic functional profiles based on Bray-Curtis dissimilarity of genes ( Figure 2B). Similar to taxonomy, the DSW pretreatment sample separated from all other samples. An NMDS subset of all other DSW samples can be found in Supplementary Figure 1B. The SSW pretreatment sample clustered with all other SSW samples at 70% similarity. A pairwise PERMANOVA comparing genes from the same treatments across experiments returned statistically significant p values  Table 5). Average similarity ranged from 61.73% between control samples and 64.33% between dispersed oil samples.
The SSW pretreatment community was dominated by an unclassified OTU in the Pseudoalteromonadaceae family while the DSW pretreatment community was dominated by Colwellia (Figure 3). Sulfitobacter were highly abundant (13-36%) in all SSW treatment communities in contrast to low abundance (0.3%) in the SSW pretreatment community. Two Marinobacter OTUs greatly increased in abundance in all DSW control and treatment communities as compared to pretreatment, while Colwellia decreased. Other OTUs that contributed to differences in pretreatment communities included Pseudoalteromonas haloplanktis and Pseudoalteromonas agarivorans in the SSW pretreatment and Colwellia psychrerythraea and Marinobacter in the DSW pretreatment. The DSW control and treatment communities were dominated by Marinobacter algicola, while Sulfitobacter was dominant in SSW control and treatment communities except for dispersant communities, which revealed P. haloplanktis to be most abundant. The DSW dispersed oil community revealed a high abundance of Pseudomonas pelagia compared to other DSW communities, and Arcobacter was more abundant in DSW oil and dispersant communities as compared to DSW control and dispersed oil communities. All SSW communities increased in diversity compared to pretreatment, with Marinomonas present only in SSW oil and dispersant communities, and a decreased abundance of Pseudomonadaceae and Pseudomonas agarivorans.
Using the KEGG Brite hierarchical structure to map individual KOs to functional categories (Supplementary  Table 7) in all metagenomes, which account for 20% of assigned reads (Figure 4). A CLUSTER analysis (Figure 5) was used to illustrate KEGG hierarchy similarity of the most abundant KOs from the different treatments and experiments. This shows that all treatments were 70% similar to each other. At the 80% similarity threshold, all metagenomes from the DSW experiment and two from the SSW experiment (oil and dispersant) clustered together, while the control and dispersed oil treatments from SSW grouped apart. In both experiments, oil and dispersant treatment metagenomes were greater than 90% similar to each other.
A SIMPER analysis revealed that membrane transport was the highest contributing category to overall differences between SSW and DSW, while translation drove similarities within each experiment. Translation was also the most abundant category overall and consisted of 30 KOs corresponding to large and small subunit ribosomal proteins, with higher reads in the DSW microcosm. Within this subset, nine categories of metabolism were annotated and the DSW microcosm had higher reads overall, with the DSW control being the highest. Membrane transport and cell motility KOs, which included chemotaxis and flagellar proteins, also had higher reads in the DSW microcosm compared to the SSW microcosm, with the highest abundance in the DSW control and the lowest abundance in the SSW dispersed oil. Lipid metabolism KOs had higher abundance in both SSW and DSW oil and dispersant metagenomes, which may have contributed to differences in SSW control and dispersed oil metagenomes compared to SSW oil and dispersant metagenomes. Signal transduction KOs had the highest abundance in SSW oil and dispersant metagenomes. Replication and repair KOs had similar abundances across both SSW and DSW oil and dispersant treatments, but had greater abundance in the DSW control.
Using MaAsLin2 to create a model to associate metadata with relative abundances revealed genes with a statistically significant (P < 0.05) differential occurrence pattern within the control group and each treatment group (Table 1). Output from analysis of filtered data revealed control metagenomes had the highest percentage (91.45%) of significant genes, indicating that this group had the most dissimilar response between SSW and DSW experiments. Oil metagenomes revealed that 49.03% significant genes and dispersant metagenomes had 49.44% significant genes. The most similar functional response was in dispersed oil metagenomes, with 45.99% significant genes. Lower percentages of significant genes observed in all treatment comparisons as compared to the control comparison indicates a more dissimilar functional response in the control groups from both experiments, and a more similar functional response in treatment groups from both experiments, specifically the dispersed oil metagenomes.

DISCUSSION
This study investigated the compositional and functional responses of marine biofilms, formed from two different source populations, on carbon steel following exposure to crude oil, chemical dispersant, and dispersed oil. Previous studies indicate that marine biofilm formation and recruitment are linked to environmental parameters, including pH, temperature, radiation, and dissolved oxygen (Lee et al., 2014;Toyofuku et al., 2016;Orland et al., 2019). Changes in these parameters can potentially influence the community composition of mature marine biofilms in shallow water, although few studies have examined marine biofilms in the deep ocean. For this study, biofilm microbiomes were hypothesized to form taxonomically distinct communities when exposed to the same experimental treatments, because the source of the microbiomes originated from different environmental settings (surface water vs. deep-sea water). It was further hypothesized that the microbiomes would have similar functional responses, especially in the formative phases of the biofilms on the same surfaces under the same conditions. This was hypothesized based on evidence from previous studies of marine biofilms that detail the necessary conditions and phases of initial formation (Kaplan, 2010;Grzegorczyk et al., 2018). SSW microorganisms isolated from Key West, FL, initially had access to sunlight, warm temperatures, and more organic matter prior to experimentation, as opposed to DSW microbial communities sourced from the Anona shipwreck site (∼1,200 m), where pressure, absence of light, and cold (4 • C) water temperatures were the prevailing in situ conditions. In support of testing these hypotheses, a side-by-side comparison of metagenome taxonomy profiles revealed SSW metagenomes harbored more observed OTUs and a greater average alpha diversity compared to DSW metagenomes. The microbial communities between SSW and DSW microcosms were distinct in pretreatment samples and remained distinct across time and different treatments.
Community analysis at a finer resolution revealed the OTUs contributing to community differences (Figure 3). Observations of different community structures between SSW and DSW experiments support the hypothesis that biofilms from two different source water populations will harbor distinct microbial taxa. The pretreatment biofilm communities for the SSW and DSW experiments contained some overlapping taxa, but as expected, given the different locations they were derived from, they contained distinct communities. A dramatic shift in community structure during the experiments (weeks 4-14) was evident in both SSW and DSW, and suggests that these early stage biofilms underwent a post-settlement maturation process (Salta et al., 2013). For the duration of the experiments, however, the DSW and SSW communities remained taxonomically distinct from each other, even though they were held under the same experimental conditions.
Rhodobacteraceae are known to be key members in initial marine biofilm formation (Elifantz et al., 2013;Kviatkovski and Minz, 2015;Dang and Lovell, 2016), and their prevalence in SSW control and treatment groups provides evidence of this process, even in the presence of contaminant addition. However, their lack of abundance in the DSW experiment indicates that other organisms must be involved in this process.
DSW metagenome communities also experienced a shift in community structure; while Colwellia dominated the pretreatment community, control and treatment groups during the experiment were dominated by M. algicola and another Marinobacter OTU, indicating a shift from free-floating to surface-attached bacterial communities (Grimaud et al., 2012;Mounier et al., 2014). While not as abundant, these same two OTUs appeared in the SSW experiment and combined were 2-6% of the total community. The genus Marinobacter is widespread in the deep ocean, and known for its capabilities of alkane degradation and EPS production (Handley and Lloyd, 2013;Gao et al., 2015;Arnosti et al., 2016;McKay et al., 2016;Hu et al., 2017). The latter capability (EPS production) potentially signals their involvement in initiating or maintaining marine biofilms. Others report differential abundances of Marinobacter and Rhodobacteraceae in comparative studies of biofilms formed offshore vs. inshore (Lema et al., 2019) and support the concept of taxonomically distinct communities with functionally redundant roles during this experiment.
Even though source populations were exposed to the same contaminants under the same environmental conditions, both experiments exhibited unique community responses during the experiments. Tremblay et al. (2017) showed that Marinobacter FIGURE 5 | Similarity CLUSTER dendrogram analyzed in PRIMER based on Bray-Curtis dissimilarity of all averaged metagenomes showing percent similarities of the group averages calculated from KEGG hierarchies in copies per million (CPM). In sample names, SSW is surface seawater metagenomes and DSW is deep seawater metagenomes. Samples are color coded by treatment: control metagenomes are red triangles, oil metagenomes are pink squares, dispersant metagenomes are green diamonds, and dispersed oil metagenomes are dark blue triangles.
abundance increased with exposure to dispersed oil. In the current study, the abundance of Marinobacter declined in the DSW dispersed oil treatment, and both the SSW oil and dispersed oil treatment. Another prior study showed that some Marinobacter also declined in the presence of chemical dispersants under laboratory conditions (Kleindienst et al., 2015b), suggesting that both oil and chemical dispersants impact the abundance and function of natural oil-degrading and EPSproducing bacteria in the marine environment.
The DSW dispersed oil community was also enriched with by P. pelagia-a psychrotolerant bacterium (Hwang et al., 2009) that has recently been suggested to degrade hydrocarbons in subarctic sediments (Gontikaki et al., 2018). P. haloplanktis, another cold adapted bacterium (Qi et al., 2019) with reported hydrocarbon degradation capabilities under laboratory conditions (Chronopoulou et al., 2015), was observed in elevated abundance in the SSW oil and dispersant treatments, relative to the control group. Colwellia OTUs, also endemic to cold and deep marine environments and known to degrade hydrocarbons (Techtmann et al., 2016), were dominant in the oil plume during the Deepwater Horizon spill (Valentine et al., 2014). Colwellia metagenomes, from three C. psychrerythraea strains isolated from a deep-sea basin, have the ability to adapt to local conditions, degrade chemical dispersants, and produce EPS for biofilm formation (Techtmann et al., 2016). The abundance of Colwellia OTUs in all DSW biofilms, specifically OTUs affiliated with C. psychrerythraea in oil and dispersant treatments, suggests that they too may have roles in responding to contaminant exposure in biofilm communities (Kleindienst et al., 2015a(Kleindienst et al., , 2016. Marinomonas and Pseudoalteromonas, taxa typically found in marine biofilms (Lawes et al., 2016), were only observed in SSW biofilms, the former in oil and dispersant treatments, and the latter in all treatments. Marinomonas are potential oil degraders (Brakstad et al., 2008) and Pseudoalteromonas have been implicated in hydrocarbon degradation in the water column (Gutierrez et al., 2013); however, the dominance of this taxon in both SSW control and treatment groups suggests that it may also have other roles in biofilm communities on carbon steel. Halomonas, which was abundant in all control and treatment groups in both experiments, has also been implicated in both EPS production and hydrocarbon degradation (Arnosti et al., 2016). The presence of these OTUs in metagenome taxonomy in this study corroborates results of biofilm 16S taxonomy data sourced from SSW in a previous study , which demonstrated the prevalence of taxa associated with biofilm formation. To fully understand the impacts of contaminant exposure on developing biofilms, functional responses must also be evaluated in tandem with compositional shifts (Galand et al., 2018), which was accomplished in the current study through analysis of the metagenome. Although community compositions between SSW and DSW biofilms had limited similarity (Supplementary Table 4), the functional profiles of all metagenomes were at least 60-70% similar to each other ( Figure 5 and Supplementary Table 5). It can be difficult to quantitate functional redundancy within microbial communities and currently there are no widely accepted methods. In this work, we used CLUSTER and NMDS plots to visualize major functional differences of highly abundant proteins mapped to functional hierarchies to help decipher functional redundancy within these communities. SIMPER analysis of these abundant KOs that determined similarity between SSW and DSW metagenomes was driven primarily by translation KOs and protein families. Signaling and cellular processes KOs were also the most abundant functional categories across all metagenomes (Figure 4). Although information about specific functions in marine biofilms is still emerging, protein translation is known to be an essential function during the transition from planktonic cells to surface-attached cells during biofilm formation (Qayyum et al., 2016). Cellular signaling within biofilms, particularly quorum sensing, is also known to play an important role in the biofilm formation process (Dang and Lovell, 2016;Wang et al., 2020). The prevalence of these functional groups in all metagenomes, including treatments, indicates that biofilms observed in this study were undergoing development, despite being exposed to contaminants, and can be considered functionally redundant in this regard.
The analysis of the most abundant KOs also showed that the SSW oil and dispersant metagenomes grouped with all DSW metagenomes while SSW control and dispersed oil metagenomes clustered separately (Figure 5). SIMPER (on abundant KOs) revealed that differences between SSW and DSW metagenomes were driven by membrane transport KOs, found in higher abundance in all DSW treatment groups relative to SSW. Membrane transport genes are important to particle-associated microbes in marine environments (Dang and Lovell, 2016), specifically for their role in transporting siderophores through the cellular membrane, which bind iron for energy (Schalk et al., 2012). This may be vital for marine biofilms growing on carbon steel, although it is unclear why DSW metagenomes have higher abundance of membrane transport KOs than SSW metagenomes. DSW metagenomes also had a higher abundance of cell motility proteins, including three chemotaxis proteins, which were highest in the DSW control and lowest in the SSW dispersed oil, suggesting that exposure to dispersed oil in SSW may impede bacterial movement, which is a crucial function in developing biofilms (Oliveira et al., 2016).
The MaAsLin analysis of all KOs provides the most direct comparison of metagenome functional data across the DSW and SSW experiments because it enables the treatment groups to be analyzed individually. The MaAsLin2 control comparison revealed that over 90% of gene associations were significantly different ( Table 1). This value declined below 50% in all other treatment comparisons, with the lowest dissimilarity (46%) between the SSW and DSW dispersed oil treatments. These data suggests that exposure to oil and dispersant in this experiment resulted in greater functional redundancy than in communities not challenged by these compounds. This is a significant outcome of this work, as it shows that while functional redundancy was generally evident across two experiments having different source communities, greater redundancy may emerge as a consequence of challenging environmental conditions.
Regardless of receiving the same treatment effects, under the same experimental conditions, taxonomically distinct metagenomes emerged between SSW and DSW experiments, likely as a result of different source water environmental conditions that dictated the phylogeny of communities present. This confirms prior amplicon sequencing-based studies of marine microbiome responses to the Deepwater Horizon spill Salerno et al., 2018;Mugge et al., 2019a). Advances in sequencing technology have permitted insight into the metabolic profiles of microorganisms and how taxonomically distinct groups can share common functions within an ecosystem (Rosenfeld, 2002). However, Galand et al. (2018) challenged functional redundancy within marine communities by demonstrating that altering microbial communities across space and time also alters their functional attributes, which further highlights the existing knowledge gap surrounding diversity and genetics of marine microorganisms. The current study and previous works (Lee et al., 2014;Lawes et al., 2016;Toyofuku et al., 2016) indicate that functional redundancy should be interpreted with caution (Fetzer et al., 2015), as multiple parameters can shape the biofilm microbiome, including interspecific interactions, inherent biodiversity, and ambient environmental conditions, and should be included in evaluations of microbial functional redundancy.
The goal of this study was to visit the topic of functional redundancy in biofilms recruited to carbon steel surfaces from different seawater source communities. Our prior 16S rRNA studies of these communities provided evidence of the taxonomically distinct nature of these communities, which was confirmed with the metagenome results presented here. However, the functional capabilities of the communities under oil spill conditions simulated in laboratory microcosms were unknown. This study demonstrates that partial functional redundancy is evident when these taxonomically distinct biofilm communities were held under the same experimental conditions. The degree of redundancy, however, was not uniform, providing evidence that exposure to crude oil and chemical dispersant may differentially impact biofilm functional redundancy. Specifically, there is less functional similarity in the absence of environmental challenges than under treatment effects. The study reinforces previous findings that the history of microbial communities will dictate their composition and functional potential (Fetzer et al., 2015;Galand et al., 2018;Louca et al., 2018). However, when marine microbiomes encounter anthropogenic contamination events, more uniform functional responses may emerge. This information may inform understanding of microbiome responses to oil spills in surface waters and in the deep sea especially when they occur proximate to steel structures, including historic shipwrecks and oil and gas infrastructure. The work may also provide better understanding of impacts oil spills may have on the microbiomes of artificial and natural reef and fouling communities. Although much remains to be discovered about functional responses in marine biofilm microbiomes, understanding how environmental effects, including exposure to anthropogenic contaminants, is a crucial step in shaping the broad understanding of observed biodiversity and function to address how human activities impact microbial ecosystems (Galand et al., 2018). Future studies incorporating the metatranscriptome of marine biofilms would also provide useful information to further explore the relationship, or lack thereof, between taxonomic composition and functional responses in biofilm assemblages.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material.

AUTHOR CONTRIBUTIONS
LH designed the study. RM and LH conducted the microbiome data analyses. LH and RM wrote the manuscript. All authors conducted the laboratory analyses and contributed to the manuscript.

FUNDING
Funding for this study was provided by the Bureau of Ocean Energy Management, Environmental Studies Program under Cooperative Agreement M13AC00015 and Interagency Agreement M13PG00020.