Succession of microbial populations and nitrogen-fixation associated with the biodegradation of sediment-oil-agglomerates buried in a Florida sandy beach

The Deepwater Horizon (DWH) oil spill contaminated coastlines from Louisiana to Florida, burying oil up to 70 cm depth in sandy beaches, posing a potential threat to environmental and human health. The dry and nutrient-poor beach sand presents a taxing environment for microbial growth, raising the question how the biodegradation of the buried oil would proceed. Here we report the results of an in-situ experiment that (i) characterized the dominant microbial communities contained in sediment oil agglomerates (SOAs) of DWH oil buried in a North Florida sandy beach, (ii) elucidated the long-term succession of the microbial populations that developed in the SOAs, and (iii) revealed the coupling of SOA degradation to nitrogen fixation. Orders of magnitude higher bacterial abundances in SOAs compared to surrounding sands distinguished SOAs as hotspots of microbial growth. Blooms of bacterial taxa with a demonstrated potential for hydrocarbon degradation (Gammaproteobacteria, Alphaproteobacteria, Actinobacteria) developed in the SOAs, initiating a succession of microbial populations that mirrored the evolution of the petroleum hydrocarbons. Growth of nitrogen-fixing prokaryotes or diazotrophs (Rhizobiales and Frankiales), reflected in increased abundances of nitrogenase genes (nifH), catalyzed biodegradation of the nitrogen-poor petroleum hydrocarbons, emphasizing nitrogen fixation as a central mechanism facilitating the recovery of sandy beaches after oil contamination.


Results and Discussion
In this study, an in situ experiment was conducted whereby standardized sediment-oil-agglomerates (sSOA) in stainless steel meshballs (3.8 cm diameter) were attached to PVC pipe and buried in the Pensacola Beach supratidal zone from 10 cm to 50 cm sediment depth (Fig. 1). Ten of these arrays were buried on 22 October, 2010 and retrieved 41,89,131,181,235,279,327,445,735, and 1152 days after burial. SOAs along with sandy sediment surrounding each of the sSOAs were collected to characterize the impacts to indigenous microbial communities. Control sands that showed no oil contamination were also collected from the surface of a nearby sand dune after 41,89,279,445,735,1152 since the onset of the experiment. sSOAs, sSOA-surrounding sands, and control sands were then collected over a 3 year time series to investigate in situ biodegradation of SOAs and microbial community dynamics. Total enviromental DNA was extracted from the source material (sSOA at day 0), incubated SOAs, SOA-surrounding sands, and control sands and used for downstream analysis.
Diversity and composition of microbial communities in SOAs. No significant trends in alpha diversity were observed over the time course of the experiment for microbial communities associated with SOAs and SOA-surrounding sands. Alpha diversity as determined by Shannon indices showed substantial variation (Fig. 2a) in SOAs and SOA-surrounding sands from both sediment depths, while diversity in control sands was more consistent with time. Whereas a previous study of buried small oil particles and oil films in the supratidal zone of Pensacola Beach (PB) showed a >50% reduction of Shannon indices during the initial 6 months after oil came ashore 17 , here we detected little to no reduction of Shannon indices in the SOAs. We attribute the lack of diversity change to the fact that the source material used in this study was already colonized by an established community of hydrocarbon-degrading bacteria. However, statistical analysis showed that Shannon indices in SOAs were distinct from those of control sands. The mean number of observed OTUs declined by 23% in SOAs and SOA-surrounding sands (827 ± 147 observed OTUs) in comparison to control sands (1069 ± 148 observed OTUs) (Fig. 2b).
In this study, the source material used for the experiment was produced from SOAs that had freshly formed from MC252-oil mousse washing onto the shore within a few months after the DWH oil spill (end of June 2010), and thus had a similar hydrocarbon composition as the weathered mousse 32 . SOAs formed as the mousse was stranded and mixed with beach sand. Microbial communities in the SOAs were distinct from those present in SOA-surrounding sands or in control sands, as determined by beta diversity of SSU rRNA genes according to the Bray Curtis distance metric (Fig. 3a). SOA community composition more closely resembled the source material in the beginning of the experiment and then strongly diverged from the source material and surrounding control sands across the time series, showing no evidence of recovery (Fig. 3a). To further examine the relationships between beta diversity and oil contamination, a canonical analysis of principal coordinates (CAP) analysis was performed exclusively on SOAs over time. For the first 6 months of the time series, community structure was strongly linked to total petroleum hydrocarbon concentrations, indicating that the input of labile hydrocarbon compounds was driving microbial community dynamics more than time or depth. Across intermediate time scales, after labile hydrocarbons were depleted, beta diversity appeared to vary with the concentration of polycyclic aromatic hydrocarbons. Subsequently, the influence of petroleum hydrocarbons diminished from the end of year 1 to year 3, and communities were structured by depth and time (Fig. 3b).
The impact of oil contamination on microbial community composition was further examined in buried SOAs by investigating the dynamics of specific taxa over the 3-year time series. The mean of relative abundances from the two different depths were pooled for further analysis. A clear succession of microbial populations was observed in SOAs with time, whereas these same taxa were much lower in relative abundance and showed no noticeable pattern in SOA-surrounding sands and control sands (Fig. 4). At the phylum and class level, microbial communities in SOAs were dominated by Gammaproteobacteria (42 to 58% relative abundance), Alphaproteobacteria (38 to 58%), and Actinobacteria (up to 10%) (Fig. 4). Few studies are available from microbial communities of oil emulsions or SOAs/ tar balls and taxonomic characterization in previous work was often limited to the phylum or class level 9,24,33 . In an investigation of 3 mousse samples collected off the coast Gammaproteobacteria, in tar balls collected within 13 months of the Texas city "Y" spill, although a phylum to class level classification of microbial communities was not provided in their study 24 . SOA communities appear to more closely resemble those of mousse rather than pristine sands or sands containing more diffuse oiled particles 17,34,35 . In particular, SOAs and mousse foster much higher levels of Alphaproteobacteria and Actinobacteria.
Similar to previous work in environments impacted by the DWH spill, the Gammaproteobacteria showed a maximum relative abundance of 58% early in the time course at approximately 100 days, whereas the Alphaproteobacteria peaked at 58% later at ~ 300 days, followed by a bloom of up to 10% Actinobacteria after 400 days [34][35][36] (Fig. 4). Both Gammaproteobacteria (Alcanivorax, Marinobacter) and Alphaproteobacteria (Rhodobacteraceae) were enriched in Pensacola Beach sands after the Deepwater Horizon oil spill and were   www.nature.com/scientificreports www.nature.com/scientificreports/ identified as key players in oil degradation in previous studies that monitored microbial community shifts after weathered oil contamination over a relatively shorter time period 34,35 . As in these studies of more diffuse oil contamination on shorelines 34,37 , groups of known hydrocarbon-degrading bacteria within the Gammaproteobacteria (e.g. Alcanivorax) responded first to oil in SOAs followed by members of the Alphaproteobacteria at later stages when recalcitrant oil hydrocarbons predominated. Major differences between the SOA time series presented here and time series of more diffuse oil contamination in buried sands 17,33,34 are that microbial succession occurred over much longer time scales (>3 years) and no recovery occurred in the SOAs. In previous work at Pensacola Beach 17,34,35 , PHCs returned to background levels one year after oil came ashore, and a typical beach sand microbial community had reestablished that showed little to no evidence of oil hydrocarbon degradation potential. At the order to genus level, three distinct maxima in microbial populations were observed with time in the SOAs, whereas these same taxa showed a much lower relative abundance in control samples ( . Early responders included the Caulobacterales within the class Alphaproteobacteria and Oceanospirillales within the class Gammaproteobacteria that increased in relative abundance during the first 100 days, from 0.5% to 12% and from 7% to 13%, respectively (Supplemental Fig. 2). At the genus level, the relative abundance of known alkane degraders such as Alcanivorax, Hyphomonas, Phenylobacterium, and Mycoplana increased early from 0.13-3% to 3-10% in SOAs and not in control sands (Fig. 5). Alcanivorax is known to degrade relatively short-chain alkanes and not being capable of degrading aromatic hydrocarbons 38,39 . In previous studies of oiled sands and tar balls, Alcanivorax spp. were the most abundant OTU in the oil-contaminated samples 24,34 . A high relative abundance of Alcanivorax (up to 20%) was also observed in a gravel beach after the Xingang oil spill in China 40 . The genus Hyphomonas, known to be able to utilize aromatic hydrocarbons, was enriched in oiled sand and microcosms with crude oil 35,41,42 . A DNA-based and stable isotope probing (SIP) study with [U-13 C]anthraquinone from PAH-contaminated soil showed that the genus Phenylobacterium is responsible for anthraquinone degradation 43 . The genus Mycoplana is also known to degrade aromatic hydrocarbons 44,45 . Our observations are corroborated by petroleum hydrocarbon analysis, which revealed that short-chain alkanes (C 15 ) and relatively low-molecular weight PAHs such as naphthalene were rapidly depleted over the first ~100 days when populations of early responders were enriched 32 .
Distinct microbial groups were enriched during the mid phase, between 131 and 235 days, of the time series. The genus Pseudomonas within the order Pseudomonadales was enriched from 5 to 15.2% after 131 days in SOAs, whereas Pseudomonas remained generally below 1% abundance in control sands. Pseudomonas has been reported to produce biosurfactant during PAH degradation and was also enriched in previous studies of oil-contaminated sands 34,35,46,47 . The relative abundance of Rhizobiales, Actinomycetales, and Rhodospirillales increased in the later stages of the time series after 235 days. Parvibaculum within the order Rhizobiales is known to degrade both aliphatic and aromatic hydrocarbons, and this genus was highly enriched in SOAs and not in control sands late in our time series, in agreement with previous work in oiled sands 35,48,49 . Relative abundance of Mycobacterium, a known PAH degrader, increased from <0.01% to 9% during the time series (Fig. 5). Evidence of Mycobacterium was also detected using a DNA fingerprinting approach from SOAs collected from supratidal www.nature.com/scientificreports www.nature.com/scientificreports/ zone of coastal headland beach in Louisiana 9 . Mycobacterium is the only genus shown to degrade 4-ring PAHs such as the chrysene observed in SOAs 9,50 . Mycobacterium is also known to be adapted to low moisture content and periods of desiccation 36 , which resembles the characteristics of dry SOAs 16 . Another known PAH-degrader Stenotrophomonas was a major microbial group associated with tar balls in the intertidal zone of a headland beach 9 , but its relative abundance in this study was very low at 0-0.22% throughout our time course. Members of the genus Marinobacter, a generalist group known to degrade both alkane and PAHs, were enriched in many oil-impacted environments 34,36,40,51,52 but their abundances in SOAs decreased during the experiment from 6% to near 0% (Fig. 5). Petroleum hydrocarbon analysis showed that mid-chain alkanes (C 16 -C 30 ), phenanthrene, and dibenzothiophene were degraded from ~100 days to ~300 days after the SOA burial. This implies that secondary responders e.g. Pseudomonas and some of early the responders such as Hyphomonas and Mycoplana were capable of degrading these hydrocarbon compounds. After approximately 300 days of incubation, primarily long-chain alkanes (C 30 -C 40 ) remained to be degraded 32 , which coincided with the rapid increase of the relative abundance of Parvibaculum. A few isolates of the genus Parvibaculum such as P. lavamentivorans and P. hydrocarboniclasticum are known to be capable of utilizing n-alkanes or linear alkylbenzenesulfonates 49,53 , which implies degradation of long-chain alkanes by Parvibaculum at a later stage of incubation.
The abundance of overall bacteria and diazotrophic communities in SOAs. The abundance of SSU rRNA genes on average was three orders of magnitude higher in SOAs in comparison to the surrounding sands or control sands. This bacterial bloom indicates that SOAs are hotspots of microbial growth. In the SOAs, overall bacterial abundance increased from 3.4 × 10 7 to 4.4 × 10 8 copies g −1 during the first 89 days, whereas overall bacterial abundance in SOA-surrounding sands and control sands remained at 2.0 × 10 5 to 2.4 × 10 7 copies g −1 . Characterization of the microbial communities associated with the SOAs by SSU rRNA gene sequencing, as discussed above, revealed a number of abundant groups that contain members known to be capable of nitrogen fixation (diazotrophy) including the Rhizobiales, Frankiales, Rhodobacterales, and Rhodospirillales, and we hypothesized that microbial nitrogen fixation would be enhanced with time. In order to estimate nitrogen fixation potential during in situ SOA incubation in Pensacola Beach sand, the abundance of genes encoding nitrogenase enzyme (nifH), the best studied molecular marker for nitrogen fixation 54 , was quantified. Results revealed that both SSU rRNA gene and nifH abundance in SOAs differed from those of SOA-surrounding sands and control sands (Fig. 6). Initially, the ratio of nifH to SSU rRNA gene abundance remained 0-0.02 for the first 300 days and then showed a large increase to 0.44 towards the latter stages of the time course, suggesting that diazotrophs bloomed after 1 year. The abundance of nifH genes in SOAs then was one to three orders of magnitude higher in comparison to the nifH abundance in the surrounding sands and control sands (7.6 × 10 5 to 8.4 × 10 7 copies g −1 in SOAs; 2.3 × 10 4 to 7 × 10 4 in surrounding; 2.6 × 10 4 to 5.1 × 10 4 in control). Previous studies have shown that SOAs contain extremely high C/N ratios indicative of nitrogen limitation 16 . In addition, the results from this SOA experiment are corroborated by previous studies of more diffuse oil in sands at Pensacola Beach, which also showed elevated nitrogenase abundance in oil-contaminated sand layers 35 . Results indicate that a paucity of nitrogen as evidenced by a high C/N ratio triggered microbial nitrogen fixation to produce bio-available nitrogen from atmospheric nitrogen.
Sequencing of nifH amplicons was conducted to further elucidate the microbial groups responsible for diazotrophy in the oiled dry beach environment. Diazotroph microbial diversity as determined by Shannon indices decreased over the time series (Fig. 7). At the phylum to class level, results showed that the Alphaproteobacteria and Actinobacteria were the most dominant diazotroph groups in SOAs (relative abundances of 64-71% and 15-21%, respectively) across the time series (Supplemental Fig. 3a). At the order level, members of the Rhizobiales (36-44% relative abundance) of the Alphaproteobacteria and the Frankiales of the Actinobacteria (15-21%) were the most abundant throughout time series. The late bloom of diazotroph abundance is concurrent with the maximum relative abundance of Rhizobiales and Frankiales in the times series, determined by sequencing of SSU rRNA as well as nifH amplicons, which peaked after 400 days. At the genus level, Methylobacterium within the order Rhizobiales was the most abundant diazotroph group, which constituted from 25 to 33% of nifH gene relative abundance (Supplemental Fig. 3b). Methylobacterium increased rapidly from 25% relative abundance after 445 days of incubation to 31% abundance at the end of the time course. The second most abundant diazotroph group was Frankiales which comprised up to 15-21% relative abundance.
Members of the Rhizobiales as well as Frankiales are well known as nitrogen-fixing symbionts associated with plant roots 55 , and free-living members of the Rhizobales are also thought to catalyze nitrogen fixation in a variety of ecosystems 56 . Methylobacterium was shown to grow on PAHs as well as produce biosurfactants in oil-contaminated systems 57,58 . Recently, it was shown that Frankia grows with PAHs as the sole carbon source and contains genes for alkane degradation 59 . Methylobacterium was also abundant in oil mousse collected from the sea surface and in salt marshes in the northern Gulf impacted by the DWH oil spill 33,60 . Therefore, multiple lines of evidence indicate a close coupling of petroleum hydrocarbon degradation to nitrogen fixation in SOAs. In our time series, nitrogen appears to become limiting after the first year, resulting in selection for microbial populations capable of coupling nitrogen fixation to hydrocarbon degradation. inferred metagenomic analysis. Based on microbial community composition in SOAs as determined by SSU rRNA gene amplicon sequencing, inferred metagenomic analysis was performed to assess the metabolic potential of the communities across the time series. Given the abundance of hydrocarbon-degraders and diazotrophs in the time course, our analysis focused on functional genes for hydrocarbon degradation and nitrogen fixation (Supplemental Fig. 4). The predicted relative abundance of alkane-1-monooxygenase (alkB) genes peaked within the first approximately 90 days post burial and then decreased rapidly, suggesting that relatively simple hydrocarbon substrates such as alkanes were utilized by bacteria at this early stage of the time series. Microbial groups that were predicted to contribute to alkane degradation include members of the www.nature.com/scientificreports www.nature.com/scientificreports/ Rhodobacteraceae, Pseudomonadaceae, Alcanivoraceae, and Alteromonadaceae. The predicted relative abundance of naphthalene 1,2-dioxygenase genes (nahAc, ndoB, nbzAc, dntAc) and other PAH dioxygenase genes (nidA, nidB) reached a maximum later in the time series, at approximately 200 and 300 days post initiation, respectively. Pseudomonas was predicted to contribute to the degradation of more recalcitrant, aromatic hydrocarbons later in the time series. Lastly, predicted relative abundance of nitrogenase genes increased rapidly after 400 days and peaked at approximately 750 days, indicating enhanced bacterial nitrogen fixation at a later stage of SOA  www.nature.com/scientificreports www.nature.com/scientificreports/ incubation. Microbial groups that are predicted to contribute to nitrogenase gene abundance include Rhizobiales, Rhodobacterales, and Rhodospirillales (Supplemental Fig. 5).
Predictions from PICRUSt were corroborated by the chemical evolution of petroleum hydrocarbons, as determined in our companion study 32 , along with temporal trends in the abundance of overall bacteria and diazotrophs. The predicted abundance of hydrocarbon degradation genes peaked during the first 400 days in parallel with overall bacterial abundance as well as the degradation of alkanes and aromatic compounds. Predictions of alkane monooxygenase (alkB) abundance showed good agreement with the consumption of short chain (C 15 ) alkanes, with both showing maximum changes during the first 100 days. Between 100 and 400 days, the consumption of longer chain alkanes (C 18 -C 22 ) was not concurrent with the predicted alkB abundance, indicating that other enzyme pathways are responsible for the degradation of these compounds. Finally, the predicted abundance of nitrogenase shows good agreement with the observed nitrogenase abundance, with both showing the largest increases between 400 and 750 days. Thus, multiple evidence support the coupling of hydrocarbon degradation to diazotrophy in SOAs.

conclusions
Macroscopic sediment-oil agglomerates (SOAs) were formed when MC252 oil from the Deepwater Horizon disaster reached the shores of the northern Gulf of Mexico and interacted with the sediment. Hydrocarbon-degrading bacteria were enriched and a succession of microbial populations was observed that paralleled the chemical evolution of the petroleum hydrocarbons over longer time scales (>3 years) in comparison to previous work on more diffuse oil contamination in beach sands. We provide evidence of bacterial blooms in SOAs, underlining that these large aggregates are hotspots of microbial growth. Our quantification of diazotrophs in large aggregates shows that nitrogen-fixing taxa predominate in oil-degrading microbial communities during the late stages of the time course when nutrients likely become depleted. The coupling of nitrogen fixation to hydrocarbon degradation thus represents a key process for the microbial decomposition of macroscopic oil aggregates.

Materials and Methods
Sample collection and experimental design. Sediment oil agglomerates (SOAs) were collected at Pensacola Beach, FL, USA (30.3261 N, 87.1744 W) on 30 June, 2010. SOAs were homogenized, and then filled into 3.8 cm diameter-stainless steel meshballs producing standardized SOAs. After determining the initial masses of each standardized SOA, 10 filled meshballs were attached in pairs to a PVC pipe (1.3 cm diameter) at 10 cm intervals (Fig. 1). Ten of these meshball arrays were buried in the supratidal zone at Pensacola Beach on 22 October, 2010 such that the meshballs were located at 10, 20, 30, 40 and 50 cm sediment depth ( Fig. 1 and Supplemental  Fig. 6). The arrays were retrieved 41,89,131,181,235,279,327,445,735, and 1152 days after burial, and the mass of each SOA was determined again before freezing at −20 °C in clean glass jars for further microbial community analysis. Together with the arrays, sandy sediments from the region surrounding each deployed SOA were also collected in order to identify possible impacts to indigenous microbial communities. This study analyzed SOAs and associated sand buried at 10 and 50 cm sediment depth. Control sand without oil contamination was collected from the surface of a nearby pristine sand dune at the same study site after 41, 89, 279, 445, 735, 1152 days of incubation.
Nucleic acid extraction and microbial community analysis. Total genomic DNA was extracted from SOAs using a MoBio PowerSoil DNA isolation kit (MoBio Laboratories, Carlsbad, CA) with slight modifications from the manufacturer's protocol. Briefly, 0.25 g of thawed SOA sample was placed into a 2 ml bead tube and homogenized for 1 min using a Talboys High Throughput Homogenizer (Troemner, Thorofare, NJ). Overall microbial communities and nitrogen-fixing prokaryotes were characterized by targeting SSU rRNA and nitrogenase (nifH) genes, respectively. PCR amplification of SSU rRNA genes was performed using 515 F (5′-GTGCCAGCMGCCGCGGTAA-3′) and 806 R (5′-GGACTACHVGGGTWTCTAAT-3′) primers as described by the Earth Microbiome Project (http://www.earthmicrobiome.org/empstandard-protocols/dna-extraction-protocol/) for Illumina sequencing 61 . For nifH, the primer set IGK3 (5′-GCIWTHTAYGGIAARGGIGGIATHGGIAA-3′) and DVV (5′-ATIGCRAAICCICCRCAIACIACRTC-3′) was used for PCR amplification as described previously 62 . PCR products were barcoded using an Access Array Barcode Library (Fluidigm, South San Francisco, CA), purified using the E.Z.N.A Cycle Pure Kit (Omega Bio-tek, Norcross, GA), and pooled based on DNA concentration. Purified and pooled PCR amplicons were sequenced on an Illumina MiSeq platform at the DNA services facility at the University of Illinois at Chicago (https://rrc.uic.edu/). Sequence analysis was accomplished using the software QIIME ver. 1.9.1 63 and Mothur ver. 1.38.0 64 . Sequences with quality score below 20 were removed using Mothur ver. 1.38.0 and clustered into operational taxonomic units (OTUs) by 97% and 92% sequence identity for SSU rRNA and nifH genes 62 , respectively using UCLUST 65 implemented in QIIME ver. 1.9.1. Representative sequences were aligned against the SILVA ver. 123 database (https://www.arb-silva.de/) or nifH reference alignments database 62 , and chimeric sequences were removed using UCHIME 65 implemented in Mothur ver. 1.38.0. These high-quality sequences were taxonomy assignment to the SILVA SSU rRNA or nfH reference alignments database with the RDP classification algorithm with a minimum confidence threshold of 50%. (https://www.arb-silva.de/, Gaby et al., 2018). The resultant OTU table was scaled using the CSS algorithm implemented in QIIME ver. 1.9.1 66 . The microbial diversity calculations and statistical analyses were performed with default R functions or with "phyloseq" and "vegan" R packages 67-69 . For alpha diversity analysis, Shannon indices were calculated with QIIME ver. 1.9.1. To assess shifts in the diversity and community composition over time, a Bray-Curtis distance matrix was calculated from the rarefied OTU table and used for a principal coordinate analysis (PCoA) and canonical analysis of principal coordinates (CAP) analyses. The oil effect on community similarity and dispersion was estimated with a PERMANOVA and BETADISP statistical tests with 1,000 permutations. Additionally, CAP analysis was performed to assess the correlation between microbial community structure and the following variables: Incubation.Time + Depth + Total.Oil + Total.PAH + Total.Alkens using "capscale" function in R package vegan. The significance of the CAP models was tested using the "permutest" function in a vegan package with 999 permutations. Finally, the Mantel correlation test with 10,000 permutations was applied to determine the similarity between the patterns of the chemical components and microbial communities.
Based on microbial community composition in SOAs as determined by SSU rRNA gene amplicon sequencing, Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt) was employed in order to predict metagenome functional content of differentially abundant OTUs 70 and predict Kyoto Encyclopedia of Genes and Genomes (KEGG) Ortholog functional profiles. The OTU table was normalized using the software QIIME ver. 1.9.1 63 , and each OTU was divided by SSU rRNA gene abundance. The resultant OTU table was used to create the final metagenome functional predictions. All raw sequences have been uploaded to NCBI under Bioproject PRJNA450618.