Comparison of Hydrocarbon-Degrading Consortia from Surface and Deep Waters of the Eastern Mediterranean Sea: Characterization and Degradation Potential

: The diversity and degradation capacity of hydrocarbon-degrading consortia from surface and deep waters of the Eastern Mediterranean Sea were studied in time-series experiments. Microcosms were set up in ONR7a medium at in situ temperatures of 25 ◦ C and 14 ◦ C for the Surface and Deep consortia, respectively, and crude oil as the sole source of carbon. The Deep consortium was additionally investigated at 25 ◦ C to allow the direct comparison of the degradation rates to the Surface consortium. In total, ~50% of the alkanes and ~15% of the polycyclic aromatic hydrocarbons were degraded in all treatments by Day 24. Approximately ~95% of the total biodegradation by the Deep consortium took place within 6 days regardless of temperature, whereas comparable levels of degradation were reached on Day 12 by the Surface consortium. Both consortia were dominated by well-known hydrocarbon-degrading taxa. Temperature played a signiﬁcant role in shaping the Deep consortia communities with Pseudomonas and Pseudoalteromonas dominating at 25 ◦ C and Alcanivorax at 14 ◦ C. Overall, the Deep consortium showed a higher efﬁciency for hydrocarbon degradation within the ﬁrst week following contamination, which is critical in the case of oil spills, and thus merits further investigation for its exploitation in bioremediation technologies tailored to the Eastern Mediterranean Sea.


Introduction
The Eastern Mediterranean Sea (EMS) presents unique physical-chemical characteristics, with high salinity levels of approximately 39 psu (compared to 35 psu in the open ocean), elevated bottom seawater temperatures of 14 °C (as opposed to 4 °C on the average at similar depths and latitudes in the Atlantic Ocean) and ultra-oligotrophic conditions owing to extreme phosphorus limitation throughout the water column [1].The waters of the EMS are stratified, with distinct water masses of unique physical and chemical properties present at different depths, each harboring distinct microbial communities in terms of properties, activity levels and metabolic capacities [2].
Oil pollution in the Mediterranean Sea has been attributed mainly to intense shipping activities resulting in a number of surface oil spills [3].Oil contamination from shipping has been greatly reduced in recent years; however, the EMS is now facing a greater challenge: the intensification of offshore oil and gas exploration activities in deep and ultra-deep waters [4].In the past two decades, exploration for oil and gas is moving towards increasingly deeper waters due to the depletion of easily accessible hydrocarbon reserves on land, in combination with the increasing demand for fossil fuels.An empirical analysis of company-reported incidents (e.g., well blowouts, injuries and oil spills) in the Gulf of Mexico (GoM) has identified a positive correlation between platform incidents and water depth [5].An oil spill in the semi-enclosed basin of the Mediterranean Sea, similar to the Deepwater Horizon (DWH) in the GoM, would have severe consequences for the marine ecosystem and the economic life of the surrounding countries [6].
The DWH blowout is considered the world's largest offshore oil spill so far, with unprecedented quantities of hydrocarbons (oil and gas) being released at 1500 m below sea level in the northern GoM [7].Approximately half of the total oil released, along with all the gas, remained in the water column, creating a deep hydrocarbon plume [8].Another 40% reached the surface with an important percentage of it (14%) sinking back at depth in the form of marine oil snow (MOSSFA event) [9,10].Several studies have been conducted in the aftermath of DWH, which focused on the microbial response to oil contamination [11][12][13][14].Gammaproteobacteria dominated both surface and subsurface waters but selected for different microbial populations [15].In offshore surface waters, the microbial community shifted drastically and became enriched in Cycloclasticus with Alteromonas, Pseudoalteromonas and Halomonas increasing in abundance to a lesser extent [16].Oceanospirillales, Cycloclasticus and Colwellia were the three most dominant taxa present in the deep plume during the different phases of the spill [17].
The extended research performed in the GoM following DWH revealed several unexpected aspects of hydrocarbon dynamics in the event of a deep-sea release of live crude oil.One of the major lessons learned from the DWH accident was the necessity of area-specific data on the microbial composition and activity, pre and post the event of an accidental oil release, in order to predict the fate of the oil spill in simulation models and assess the recovery of the ecosystem [18].In the EMS, the future exploitation of recently discovered oil and gas reserves in deep waters may be translated into an increased risk of accidental oil spills.Nevertheless, our current knowledge regarding the response of microbial life to oil contamination in the EMS is extremely limited [19,20].This study aims to contribute towards that direction by assessing the biodegradation capacity of two hydrocarbon-degrading consortia enriched from surface and deep waters of the EMS.Key hydrocarbon degraders and succession patterns during crude oil degradation were identified in time-series microcosm experiments containing crude oil as the sole source of carbon and energy.To our knowledge, this is the first study to provide quantitative data on the capacity of microbial communities from the Eastern Mediterranean for oil-spill remediation in surface and deep-water layers and sets the baseline for site-specific bioremediation solutions and oil-spill monitoring.

Sample Collection
The circulation of the Eastern Mediterranean can be described as a three-layer system: Modified Atlantic Water (MAW) flows in along the surface through the Straits of Sicily (~upper 150 m).This water flows in an easterly direction getting more saline under the hot and dry conditions extant in the region, particularly in summer.The Levantine Intermediate Water (LIW), formed in the Levantine basin at intermediate depths (150-400 m), is characterized by temperatures around 15 °C, higher nutrient concentrations and high salinity (39 psu), and forms the major part of the return flow of water out of the basin through the Straits of Sicily.The Eastern Mediterranean Deep Water (EMDW), occupying depths below 400 m, is characterized by particularly warm temperatures (compared to similar depths and latitudes in the Atlantic), stabilizing at approximately 13.5 °C [1,21].
Seawater samples were collected on board the R/V Aegaeo (Hellenic Centre for Marine Research) on 23-28 August 2019.The sampling station was located off Southern Crete, Greece (Koufonisi station: 26.131288E, 34.815447N).Surface (10 m) and deep (1040 m) seawater samples were collected in 12 L Niskin bottles mounted on a CTD rosette and stored (within 1 h) in acid-washed plastic containers in the dark until use in laboratory experiments.Seawater temperature and salinity at each sampling depth were 27 °C/39.6psu at 10 m and 14 °C/38.8psu at 1040 m.Two liters of seawater from each sampling depth were filtered onboard through 0.2 μM polyethersulfone (PES) membrane filters (Rephile, USA) and stored frozen until DNA extraction.

Enrichment of Hydrocarbon-Degrading Microbial Consortia
Hydrocarbon-degrading microbial consortia originating from Eastern Mediterranean surface and deep seawater samples, hereafter referred to as Surface and Deep consortia, respectively, were produced through sequential enrichments in ONR7a medium (DSMZ medium 950 (https://www.dsmz.de/microorganisms/medium/pdf/DSMZ_M/edium950.pdf(accessed on 18 October 2019) and 0.5% v/v filter-sterilized Iranian light crude oil (CO, density ~0.7821 g mL −1 ).The enrichment cultures were initially inoculated with 20% v/v aged seawater (maintained in the dark at 14 °C for 4 months) and subsequently by transferring 10% v/v of culture to new ONR7a+CO medium.Microbial biomass was measured daily as optical density (OD) at 600 nm.Transfers took place when the culture growth entered the early exponential phase (OD 0.2-0.3).The cultures were incubated at room temperature (23-25 °C) for Surface and at 14 °C for Deep enrichments.Following each transfer, part of the remaining volume of the culture was stocked in 25% glycerol aliquots and the rest was filtered through 0.2 μM PES filters (Rephile, USA) for DNA extraction.Glycerol stocks and filters were stored at −80 °C.The community composition of background and aged seawater, as well as that of each transfer step, is given in Appendix A (Figure A1).The 2nd transfer (S2, D2) from both enrichments was selected for the time-series biodegradation experiment.

Time-Series Biodegradation Experiment
Batch microcosms were set up in 250 mL Erlenmeyer flasks containing 125 mL ONR7a medium and 0.5% v/v filter-sterilized CO.The microcosms were inoculated with thawed glycerol stocks of the enriched consortia (S2 and D2 for Surface and Deep microcosms, respectively).Three treatments were examined in triplicate microcosms: the Surface consortium at 25 °C (Surface25) and the Deep consortium at 14 °C (Deep14) and 25 °C (Deep25).Sub-samples for microbial community analysis were collected after 6, 12, 18 and 24 days on 0.2 μM PES filters, as described above.A parallel set of triplicate flasks per treatment and sampling points at 6, 12 and 24 days were prepared for destructive sampling for gas chromatography-mass spectrometry (GC-MS) analysis of hydrocarbon concentrations.Abiotic hydrocarbon losses were assessed in triplicate flasks prepared without inoculum and sampled on Day 24.

Hydrocarbon Extraction and GC-MS Analysis
Liquid-liquid extraction was performed to obtain the microbial activity extract, free from the aqueous culture medium (ONR7a).For the extraction of the organic compounds, an equal volume of dichloromethane (DCM Suprasolv ® , Merck KGaA, Darmstadt, Germany) was used (3 × 20 mL for each extraction), in 100 mL Erlenmeyer flasks with a glass spout at the bottom.The flasks were shaken manually to assist the dissolution of the organic compounds in the solvent and the extract was filtered through columns of anhydrous sodium sulfate and fiberglass to remove any water residues.Following solvent removal on a rotary evaporator, the dried samples were transferred to 4 mL vials with a small amount of DCM Suprasolv ® and concentrated by evaporation at low heat on a hot plate (~40 °C) with simultaneous nitrogen blow.The solid extracts that occurred were then separated in saturated and aromatic hydrocarbon fractions by elution through SPE columns (Bond Elute TPH, Agilent Technologies, Inc., Santa Clara, CA, USA) with nhexane Suprasolv ® (Merck KGaA, Darmstadt, Germany) and DCM Suprasolv ® , respectively.GC-MS analysis was performed on an Agilent GC-MS HP 7890/5975C system, with an Agilent HP-5MS 5% phenyl methyl siloxane column (60 m × 250 μm × 0.25 μm).The hydrocarbon mixture consisted of an Oil Analysis Standard (Absolute Standards Inc. ® , Hamden,CT, USA) containing 44 compounds, and a 17a(H),21b(H)hopane (Chiron AS ® , Trondheim, Norway).The standard composition of the hydrocarbon mixture was normal alkanes from C10 to C35, pristane and phytane and 16 polycyclic aromatic hydrocarbons (PAHs) (naphthalene, phenanthrene, anthracene, fluorene, fluoranthene, benzo(k)fluoranthene, benzo(e)pyrene, benzo(a)pyrene, perylene, indeno(g,h,i)pyrene, dibenzo(a,h)anthracene and benzo(1,2,3-cd)perylene).

Kinetic Evaluation
Kinetic analysis was used to calculate the degradation rate of the hydrocarbon compounds.Biodegradation of crude oil is typically assumed to follow first-order kinetics, when hydrocarbons are consumed early in the process [22,23].In our experiments, background nutrients (such as N and P) were apparently depleted by the microbial consortia within the first 6 days, after which biodegradation slowed down.In order to compare the different consortia, we determined the initial degradation rate of each treatment during the first few days when the initial bacterial concentrations were nearly identical.During the initial period, we can assume first-order kinetics apply and hence the concentration is given by where C(t) denotes the concentration of hydrocarbons (ppm) at time t, t is the time expressed in days, C0 is the initial concentration of hydrocarbons (ppm) and k is the apparent first-order degradation rate constant (day −1 ), which strongly depends on the average concentration of the hydrocarbon-degraders during the period of interest.The half-life (t1/2), which represents the time needed for the C0 to be reduced to half of its amount, was calculated by the following equation:

DNA Extraction
DNA was extracted according to [24], with modifications.Briefly, each PES membrane filter was cut in smaller pieces and divided into two Eppendorf tubes containing 670 μL of CTAB extraction buffer (1 M Tris-HCl (pH 8), 0.5 M EDTA (pH 8), 1 M NaH2PO4 (pH 8), 5 M NaCl and 5% cetyltrimethylammonium bromide) and mixed using a high-speed vortex for 2 min.The samples were then subjected to 3 freeze-thaw cycles in liquid nitrogen (1-2 min) to 65 °C (10 min) followed by addition of 100 μL (10 mg mL −1 ) lysozyme and 10 μL (20 mg mL −1 ) Proteinase K and incubation at 37 °C for 30 min.After incubation, 60 μL of 20% sodium dodecyl sulfate (SDS) were added and the samples were left at 65 °C for 2 h with gentle mixing every 10 min.The supernatant was collected after centrifugation at 10,000× g for 10 min and transferred into a new 2 mL microcentrifuge tube where it was mixed with an equal amount of phenol:chloroform:isoamylalcohol (25:24:1).The aqueous phase obtained after centrifugation (15,000× g for 15 min) was subjected to a second purification step with 700 μL chloroform: isoamylalcohol solution (24:1).DNA was precipitated overnight at 4 °C by addition of a x0.7 volume of ice-cold isopropanol to the aqueous phase.The precipitated DNA was centrifuged at room temperature (13,000× g, 45 min) and the supernatant was discarded.The DNA pellet was washed with ice-cold 70% moleculargrade ethanol and centrifuged again for 30 min at 13,000× g.Finally, the ethanol was removed, and the pellet was left to dry for 15 min in a laminar flow hood.Genomic DNA was resuspended in 50 μL TE buffer.DNA concentration was measured on a Qubit 4.0 fluorometer (Invitrogen, Thermo Fisher Scientific, Waltham, MA, USA) using the Qubit dsDNA high sensitivity assay.All chemical reagents were purchased from Sigma-Aldrich Inc. (Merck KGaA, Darmstadt, Germany).

16S rRNA Amplicon Sequencing and Bioinformatic Analysis
The V3-V4 regions of the 16S rRNA gene was sequenced on the Illumina MiSeq (2 × 300 bp) platform by Biosearch Technologies, LGC Genomics GmbH (Berlin, Germany) using the universal primers 341F and 806R [25].Bioinformatic analysis was performed in the "R" programming environment (version 4.0.2) [26] using the DADA2 R package and pipeline for microbiome data analysis [27].The SILVA SSU taxonomic training data for DADA2 (Silva version 138) [28] was employed for taxonomy assignment.A neighborjoining phylogenetic tree was constructed in the phangorn R package (version 2.5.5)[29] following multiple alignment of ASVs using the DECIPHER R package [30].The DADA2 pipeline resulted in a single object containing the sample-by-sequence feature table, the experimental metadata, the sequence taxonomies and the phylogenetic tree, which was subsequently imported into the phyloseq R package for further analysis [31].A total of 1,563,382 reads were obtained from 49 samples after Illumina paired-end sequencing of 16S rRNA gene amplicons (V3-V4 region) with 1,082,378 total reads remaining after quality filtering and chimera removal.The number of sequences per sample ranged from 4993 to 43,902, with an average of 22,083 sequences.One feature for which a taxonomic rank could not be assigned at the phylum level was removed.The remaining dataset consisted of 1415 unique ASVs.Supervised prevalence filtering was further applied, with five more phyla removed from the dataset.Following prevalence filtering, 1404 unique ASVs remained, distributed among 19 Phyla.

Statistical Analysis
Statistical analysis was performed in the "R" programming environment [26].All plots were produced using the R package ggplot2 unless otherwise specified [32].Linear regression was used to test the effect of treatment (Deep14, Deep25 and Surface25) and time (categorical variable, levels: 0, 6, 12, 18 and 24) on alpha diversity (richness, Shannon-Weaver, Simpson and Chao1 indices).The optimal models were identified based on AIC (step function) by forward and backward selection.In case the optimal model included non-significant regression parameters, further model selection steps were undertaken using the F-test (drop1 function) [33].Model validation was applied to check that the underlying statistical assumptions were not violated.Post-hoc multiple comparisons were carried out using Tukey's honest significant difference test.
Principal coordinate analysis (PCoA) and distance-based redundancy analysis (dbRDA) on weighed Unifrac distance were used to explore similarities between microbial communities and assess the effect of depth, temperature and time using the vegan R package [34].Count data were normalized by the cumulative sum scaling (CSS) method prior to multivariate analysis [35].Negative binomial regression on DESeq2transformed count data was used to identify differentially abundant taxa related to explanatory variables using the DESeq2 package [36,37].DESeq2 was run using the Wald test, with automatic filtering of low abundance ASVs, automatic calculation of adjusted pvalues and an alpha of 0.01.Levin's niche breadth and niche overlap indices of ASVs were measured with the MicroNiche R package [38].Prior to niche breadth analysis, the dataset was rarefied to 4993 sequences per sample.The MicroNiche functions were applied to the ASV count table after genus taxa assignment.A limit of quantification (LOQ) was applied, defined as a threshold calculated from the distribution of microbial taxa within the dataset and a 95% certainty that these taxa will fall within a null distribution, where the mean taxon abundance is zero.In our case study this was set to 1.3× standard deviations from zero in order to exclude any false-positive taxa that are misclassified as specialists/generalists. Levin's Overlap values (LO12 and LO21) are plotted as heatmaps using the gplots R package [39].

Hydrocarbon Degradation
The results from the GC-MS analysis revealed an effective degradation of CO by Day 24 in all three treatments (Surface25, Deep14 and Deep25), achieving a total of 51%, 53% and 55% reduction in the total petroleum hydrocarbon (TPH) concentration, respectively (Figure 1A).More than 85% of the TPH degradation by the Deep consortium at both temperatures (Deep14 and Deep25) was achieved by Day 6.The Surface consortium presented a lag phase until Day 6 but reached a comparable TPH reduction by Day 12.The crude oil fractions, examined by GC-MS analysis, consisted of about 99% saturated components, out of which 89% were short-chain (C14-C25) and 11% were long-chain alkanes (C26-C35), and 1% PAHs.The light alkanes fraction was rapidly degraded by the Deep consortium regardless of temperature; 50% of light alkanes were degraded by Day 24, of which 86-88% were removed by Day 6 (43-44% of total light alkanes) (Figure 1B).Light alkanes degradation reached 47% and 50% by Day 12 in the Deep14 and Deep25 treatments, respectively, but was negligible between 12 and 24 days.Almost no degradation of light alkanes occurred within the first 6 days in the Surface25 treatment, yet a 48% reduction in the concentration of light alkanes was achieved by Day 24 (85% of total reduction between 6 and 12 days) (Figure 1B).
Long-chain alkanes degradation was proportionally higher than that of light alkanes and ranged between 58 and 63% of the total heavy alkanes concentration on Day 24 in all treatments (Figure 1C).The maximum degradation rate of this hydrocarbon group was observed in the 0-6 time interval for both Deep treatments and 6-12 for Surface25.Unlike light alkanes, the Surface consortium immediately started degrading the heavier compounds of the oil.Deep25 presented the highest degradation of the heavier compounds (17% reduction) followed by Deep14 and Surface25 (13% and 12%, respectively).
PAH degradation was initiated immediately in all three consortia, reaching a total PAH reduction of ~15% by Day 24 (Figure 1D).Both Deep consortia performed better than Surface25 at the early stages, similar to the pattern observed for saturated hydrocarbons.The majority of PAHs were consumed by Day 6 in Deep14 and Deep25, after which the degradation slowed down until Day 24.PAH degradation in Surface25 occurred in the 0-6 and 12-24 time intervals; between 6 and 12 days, when the maximum degradation of alkanes is observed, PAH degradation in Surface25 stalled.Supporting data for the hydrocarbon degradation and the different components can be found in the Supplementary Materials section (Figures S1-S3, Table S1).
The results of the kinetic evaluation are presented in Table 1.As seen in Figure 1, no significant degradation occurred after Day 6 for the Deep consortia and after Day 12 for the Surface, due to the depletion of background nutrients; hence, in Table 1, the initial rates of biodegradation (up to Day 6 for the Deep and up to Day 12 for the Surface consortium) are provided.

Community
Crude Oil Components kapp (day The Deep consortium at 14 °C (Deep14 treatment) was highly enriched in Vibrio up to Day 6. Alcanivorax replaced Vibrio as the most abundant taxon between Days 12 and 24, reaching a relative abundance of 40-60%.Alteromonas and Pseudoalteromonas were also present in the Deep14 treatment but maintained constant relative abundance throughout the experiment (Figure 2).The Deep consortium at 25 °C (Deep25) was similar to that of the Deep14 samples up to Day 6. Changes in the dominant taxa due to increased temperature became prevalent on Day 12 onwards, with Pseudomonas and Pseudoalteromonas collectively reaching 80-90% in relative abundance (Figure 2).Alphaproteobacteria of the order Rhodospirillales (genus Thalassospira) were significantly enriched in Surface samples, reaching up to 40% of the total abundance (Figure 2, Figure S4).Halomonas and Alteromonas were also highly abundant in Surface samples.Idiomarina became enriched through time, reaching approximately 15% of the microbial abundance on Day 24 (Figure 2).

Alpha Diversity
A steep drop in alpha diversity was observed with seawater storage, i.e., between the background (fresh) and aged seawater (Figure 3, Figure A1).Alpha diversity decreased further with oil contamination (starting consortium Day 0) but did not change further with time during the incubation experiment.Statistical analysis to test for differences between treatments and time effects on alpha diversity metrics was performed for Days 6, 12, 18 and 24, in which triplicate samples were collected.Alpha diversity, as measured by the Shannon index, differed significantly between treatments (F = 42.09,d.f.2, p < 0.001) but did not change with time (F = 1.41, d.f.3, p = 0.257).The same pattern was observed using other diversity metrics (ASV richness, Simpson index and Chao1 index).Pairwise comparisons with Tukey's post-hoc test showed consistently higher ASV richness and diversity in the Deep14 treatment (Figure 3, Table S2), with the exception of the Chao1 index; the diversity estimator Chao1 index, which calculates the expected ASVs based on the low abundance species (singletons and doubletons), was not significantly different between the Deep14 and Deep25 treatments (t-value = −1.81,p = 0.181).

Beta Diversity
Microbial communities clustered into three groups according to the treatment (Figure 4A).Samples of the same treatment clustered together with the exception of Deep25 samples on Day 6, which clustered with Deep14, reflecting the gradual shift of the Deep community due to increased temperature.Treatments were significantly different to each other based on PERMANOVA multilevel pairwise comparisons (R2: Surface25-Deep25: R2 = 0.88, p = 0.003; Surface25-Deep14: R2 = 0.91, p = 0.003; Deep25-Deep14: R2 = 0.15, p = 0.012).Depth, temperature and time explained 69% of the variation in the structure of microbial communities based on dbRDA analysis using weighed Unifrac distances (Figure 4B).The distribution of the samples on the dbRDA triplot indicated clustering based on depth along the first axis, explaining 60% of the variance in the data, and on temperature and time along the second axis (Figure 4B).The ANOVA-like permutation testing for dbRDA (anova.ccafunction in vegan R package) indicated that all explanatory variables were significant at the 5% level (depth: F = 58.12,d.f.1,34, p = 0.001; temperature: F = 4.99, d.f.1,34, p = 0.007; time: F = 4.88, d.f.1,34, p = 0.011).

Identification of Influential Taxa Based on DESeq2 Analysis
Differential abundance analysis was conducted to identify ASVs that were strongly influenced by depth, temperature and time.DESeq2 analysis by depth identified 20 influential ASVs that were strongly enriched in the Surface25 treatment (Figure 5A, Table S3).Characteristic taxa of the Surface25 treatment belong to the Gammaproteobacteria families Halomonadaceae (genus Halomonas) and Idiomarinaceae (genus Idiomarina), and the Alphaproteobacteria Thalassospiraceae (genus Thalassospira).Halomonas and Thalassospira became enriched within the first 6 days while Idiomarina became influential after Day 12 (Table S4).DESeq2 analysis by temperature revealed influential ASVs in the Deep14 treatment (negative log2FC) (Figure 5B, Table S3), namely, Rhodobacteraceae (Sulfitobacter), Alcanivoracaceae (Alcanivorax), Vibrionaceae (Vibrio) and Marinobacteaceae (Marinobacter).DESeq2 analysis per treatment at each time interval (Table S4) identified common aspects in the behavior of Alcanivorax and Sulfitobacter to temperature, which both became depleted in Surface25 at 0-6 days and in Deep25 between 6 and 18 days.On the other hand, Thalassospira was favored by higher temperature, as evidenced by enriched ASVs in the Deep25 treatment between 12 and 24 days (Table S4).

Levin's Niche Analysis
Levin's niche analysis can provide valuable quantitative information as to which genera are considered specialists or generalists and assess taxon-taxon relationships [38].This is achieved by using Levin's niche breadth (BN) and Levin's overlap (LO), respectively [38].We performed niche analysis of bacterial genera identified in the three treatments of the experiment (Deep14, Deep25 and Surface25).A limit of quantification (LOQ) threshold was employed on rarefied data to exclude any taxa with low, random abundance from being falsely characterized as specialists/generalists (Figure S5).Levin's BN values indicate which taxa are present under any environmental conditions (generalists) and which are restricted in a unique environment, in this case, a treatment (specialists).Generalists have BN values closer to 1 while specialists have values closer to 0 (Figure S6, Table 2).Idiomarina, Marinobacter, Vibrio, Thalassospira and Sulfitobacter were the strongest specialists classified in our samples (p < 0.001, Table2).In addition, Pseudoalteromonas and Pseudomonas were also classified as significant specialists (p < 0.05), resulting in BN values below the 0.05 quantile threshold (BN < 0.57) (Table 2, Figure S6).Alcanivorax, Alteromonas and Halomonas presented BN indices that are neither proportionally equal (generalists), neither associated with a unique environment (specialists).Identified specialists are shown in boxplots indicating the treatment they were enriched in (Figure S7).Levin's Overlap (LO) takes into account the species abundance across environments and indicates which taxa coexist on a heatmap plot (Figure 6).LO12 and LO21 axis values do not necessarily have the same indices as that depends on the distribution of the compared taxa and their classification as specialists or generalists [38].For example, a strong overlap between Alcanivorax and Sulfitobacter was observed on the LO12 axis but not on LO21 (Figure 6).This indicates that Sulfitobacter co-occurred with Alcanivorax but the absence of Sulfitobacter did not entail absence of Alcanivorax.In addition, Idiomarina depended on the presence of Halomonas but not the opposite (LO12 overlap only).On the other hand, Halomonas and Thalassospira, as well as Pseudomonas and Pseudoalteromonas, overlap on both LO12 and LO21 axes, meaning that both species coexist in specific environments.Finally, Marinobacter depended on the presence of Vibrio (LO12 overlap only) (Figure 6).

Discussion
Hydrocarbon-degrading bacterial consortia were enriched from surface and deep waters of the EMS with the aim to investigate the degradation capacity of indigenous microorganisms and the potential use of enriched consortia in bioremediation.We did not observe a significant temperature effect on the degradative capability of the Deep consortium, which responded rapidly to hydrocarbon contamination, removing 47% of the alkanes and 12% of the PAHs within 6 days at both tested temperatures (25 °C and 14 °C).In comparison, the Surface consortium had a lower degradative capability, removing only 8% of the total alkanes and 8% of the PAHs during the same period with close to zero degradation of the low molecular weight alkanes.The response lag of the Surface consortium is reflected in the first-order degradation rate constants calculated over the active degradation period during the incubation (6 days for the Deep consortium and 12 days for Surface); 0.1048, 0.1064 and 0.0422 d −1 for TPH in Deep14, Deep25 and Surface25, respectively.The half-life of TPH was 16 days for Surface25 and 6-7 days for the Deep14 and Deep25 treatments.Despite the response lag, both the Deep and Surface consortia removed similar amounts of TPH by Day 24, ranging from 51% (Surface25) to 55% (Deep25).
Between-study heterogeneity in experimental conditions hampers direct comparison of the degradative characteristics of consortia.Temperature, nutrient levels, the geographic origin of microbial communities and previous exposure to hydrocarbon sources as well as the composition and concentration of the hydrocarbon substrate are part of a non-exhaustive list of factors affecting hydrocarbon degradation rates [11,22,40].In the GoM, a surface-water consortium degraded alkanes at a rate of 0.15 d −1 at 25 °C and removed almost 90% of total the alkanes and 77% of the PAHs after 8 days, while the deep-water consortium at 5 °C was equally capable in hydrocarbon removal albeit with a considerable lag phase [41].The high degradation rates by GoM consortia in the aforementioned study may be attributed to the low concentration of crude oil used in the experiments that is known to result in more efficient hydrocarbon removal [11,22,23].At comparable, to this study, crude oil concentrations (400 mg L −1 ), Bacosa et al. measured 65% degradation of alkanes after 50 days, and also similar to this study, degradation rate constants (0.05 d −1 ) in surface waters of the GoM [23].Microbial consortia originating from shallow and deep sub-arctic sediments of the NE Atlantic were less efficient in diesel and model oil degradation (1% v/v concentration; 50-58% and 33-42% TPH removal by the shallow and deep consortia, respectively) at 20 °C [42].
Temperature often emerges as the major environmental parameter determining the composition of oil-degrading microbial communities and the rate of hydrocarbon degradation [43].Deep-sea microbial communities exhibit a lower hydrocarbon degradation potential than shallower communities when incubations are performed at in situ conditions primarily due to the low temperatures prevailing in the deep sea [44][45][46].At equivalent incubation temperatures, there is not enough evidence to support an intrinsically lower potential of deep-water microbial communities for hydrocarbon degradation compared to shallow-water counterparts although the latter seem to perform better in most cases [40,42].The Deep consortium, in this study, performed better in terms of response time and overall performance even at 14 °C although differences with the Surface consortium were negligible by Day 24.Remarkably, raising the incubation temperature of the Deep consortium by 11 °C (from 14 to 25 °C) did not produce the expected increase in biodegradation rate based on a Q10 (change in metabolic rates with a 10 °C increase) range of 2-3, which is widely accepted for temperature compensation in biodegradation rate calculations (e.g., a Q10 of 2 is applied in the OSCAR Oil Spill Contingency and Response model) [47].Based on the Q10 approach, the hydrocarbon biodegradation rate in Deep25 should be at least double that in the Deep14 treatment.Yet, Q10 remains close to 1 for all hydrocarbon fractions, demonstrating the adaptation of deep-water microbes in the Eastern Mediterranean to ambient conditions and their ability to degrade hydrocarbons at comparable, if not greater, rates to surface communities.Our results confirm findings of faster-than-expected hydrocarbon degradation rates in deep waters of the GoM following DWH [48].
The narrow temperature gradient between surface and deep waters in the Mediterranean, selecting for mesophilic microorganisms throughout the water column, may partly explain why temperature is of secondary importance in this environment compared to other oceanic locations where the temperature of deep waters is typically 4 °C or less [40,46].The origin of the consortium from surface or deep waters was the most important variable in shaping the structure of the bacterial community, followed by temperature.These two factors selected for different dominant oil-degrading taxa within each treatment.Thalassospira and Halomonas were characteristic of the Surface consortium throughout the incubation period.Idiomarina, identified as a strong specialist based on the Levin's BN index, appeared solely in the surface consortium after Day 12 and increased in relative abundance thereafter.Several members of the Thalassospira and Idiomarina genera are known PAH degraders [49,50].In the Eastern Mediterranean, Idiomarina emerged as a potential PAH degrader in surface sediments of contaminated sites following the Agia Zoni II oil spill in Greece [51].Although highly abundant in the original S2 consortium, Alcanivorax was not present in the Surface25 treatment at any time point.It is not clear why Alcanivorax was not competitive at the given experimental conditions, but its absence may well explain the response lag in the degradation of light alkanes by the Surface25 community.In the Deep consortium incubations, Alcanivorax and Vibrio dominated at 14 °C but were replaced by the cosmopolitan Pseudoalteromonas and Pseudomonas at 25 °C after Day 6.The role of Alcanivorax species in deep-sea oil-spill remediation is ambiguous; several species of Alcanivorax have been isolated from the deep sea [52][53][54] yet Alcanivorax was not reported among the dominant taxa in the deep hydrocarbon plume during DWH and, based on ex situ experiments, this genus is not considered as characteristic of the microbial community response to oil contamination in the deep sea [11,45,55,56].Nevertheless, Alcanivoraceae were identified as key players in crude oil degradation in deep waters of the EMS [20], in agreement with the results of this study.However, the importance of Alcanivorax as a key taxon in hydrocarbon degradation in the deep Eastern Mediterranean needs to be confirmed under high hydrostatic pressure conditions.
The genus Vibrio was characteristic of the deep consortium at 14 °C.Members of the genus are capable of alkane and PAH degradation [57,58], and Vibrio has been found at high concentrations in oil-contaminated sites, particularly in association with biofilms [51,[59][60][61].Nevertheless, Vibrio strains isolated from oil-contaminated beach sands in the Gulf of Mexico following DWH showed little-to-no oil consumption [62].In, similar to ours, oil-enriched microcosms from GoM, Vibrio represented ~80% of the total community (DNA) abundance in surface water consortia, but only 25% of the active community (RNA) [41].Here, monitoring the community structure in each enrichment step allowed us to observe that Vibrio was not consistently present in all transfers of the deep consortium.The presence of Vibrio in high concentration was restricted to the second transfer consortium (D2), which introduces a random element in the outbursts of this genus.Further biologically-relevant conclusions can be extracted from hierarchical testing (DESeq2), in conjunction with Levin's niche analysis, both of which are useful for gaining insight on influential low-abundance taxa.Of interest in this study were Marinobacter and Sulfitobacter (family Rhodobacteraceae), both of which were associated with the Deep consortium at 14 °C and are common encounters in oil-contaminated sites under various environmental conditions [62,63].Interestingly, the presence of Sulfitobacter was tightly coupled to that of Alcanivorax both spatially and temporally.The preference of Sulfitobacter for lower temperatures has been previously reported [40].Little is known about the degradative capability of Sulfitobacter although utilisation of crude oil as a carbon source has been previously reported [56] and draft genome sequences of Sulfitobacter have revealed the presence of genes involved in aromatic hydrocarbon degradation [64].

Conclusions
Over the last few years, the area south of Crete has drawn attention concerning its oil and gas reservoirs that, if proven financially exploitable, will result in the development of a number of oil and gas production platforms in waters exceeding 1000 m in depth and in close proximity to land.This investment would increase the risk of oil-spill accidents in the area.The unique physical-chemical properties and the distinct microbial communities of the EMS require site-specific data to improve oil-spill preparedness and develop appropriate mitigation measures in the event of an accident.This study generated enriched hydrocarbon-degrading bacterial consortia from deep waters of the EMS that proved to be readily effective in the biodegradation of crude oil constituents within the first week following contamination, which is even more critical in the case of oil spills.This fundamental finding sets the foundations for further study on the efficacy of the Deep consortium under field environmental conditions, i.e., high pressure and low nutrient levels, with possible applications in bioremediation technologies tailored to the EMS.

Supplementary Materials:
The following are available online at www.mdpi.com/1996-1073/14/8/2246/s1.   Figure S4: Relative abundance of bacterial taxa at the order level plotted against time (days) and grouped by treatment (Deep14, Deep25, Surface25).Brackets on the x-axis correspond to the triplicate samples for the given timepoint.Figure S5: Rank distribution of genus taxa.The limit of quantification (LOQ), shown in blue circles, is 1.3 standard deviation from zero in order to address false positive results.Taxa are displayed in black while the lognormal rank distribution model is given in red circles.Figure S6: Null model distribution of Levin's BN niche.On the x-axis, the Levin's BN distribution across the three different treatments (Deep14, Deep25, Surface25) versus the frequency of random BN values on the y-axis.The red dotted lines indicate the 0.05 and 0.95 quantiles.Samples below the 5th quantile are considered specialists while those above the 95th are generalists.Figure S7: Direct counts from rarefied data of specialist taxa, identified by Levin's BN niche, by treatment.Table S1: % Degradation by Day24 for each of the 16 EPA priority PAHs examined in this study in treatments Deep14, Deep25, and Surface25.Empty values correspond to non-detectable peaks.Table S2: Alpha diversity indices of microbial communities in treatments Deep14, Deep25, and Surface25 at each time-point.Table S3: DESeq2 analysis results.Table S4: Successional patterns in the microbial community structure as shown in DESeq2 analysis by time and treatment.

Figure 1 .
Figure 1.Hydrocarbon biodegradation over time for the Deep14, Deep25 and Surface25 treatments for (A) total petroleum hydrocarbons (TPH), (B) light alkanes (C14-C25), (C) heavy alkanes (C26-C35) and (D) PAHs.Values on the y-axis correspond to the concentration range of each hydrocarbon group to aid visualisation.Error bars represent the standard deviation (n = 3).

Figure 2 .
Figure 2. Relative abundance of bacterial taxa at the genus level plotted against time (days) and grouped by treatment (Deep14, Deep25 and Surface25).Taxa that are less than 1% in abundance are grouped as "Other".Brackets on the x-axis correspond to the triplicate samples for the given timepoint.Day zero corresponds to the original D2 and S2 consortia and is common for the Deep14 and Deep25 treatments.

Figure 3 .
Figure 3. Measures of alpha diversity: (A) richness, (B) Shannon's diversity index, (C) Chao1 index and (D) Inverse Simpson index."Backgr" corresponds to the background microbial community in seawater at a 10 and 1040 m water depth."Aged" is the seawater that was used as inoculum for the enrichment of consortia.Day 0 corresponds to consortia S2 and D2, which were selected for the microbial succession experiment.

Figure 4 .
Figure 4. (A) Principal coordinate analysis (PCoA) of the weighted UniFrac distances.(B) Distance-based redundancy analysis (dbRDA) of the weighted UniFrac distances quantifying the impacts of depth, temperature and time on microbial community composition.Cumulative sum scaling (CSS) normalization was applied to the data prior to ordination analysis to account for differences in library size.

Figure 5 .
Figure 5. ASVs strongly influenced by depth (A) and temperature (B) based on DESeq2 analysis.Log2FC is the logarithmic fold change of taxa abundance between two conditions (Surface-Deep and 14-25 °C).

Figure 6 .
Figure 6.Heatmap of Levin's Overlap (LO) indicating coexistence of species taking into account their environmental distribution in the three different treatments (Deep14, Deep25 and Surface25).Only taxa that pass the LOQ line are included.The scale bar shows values between 0 for poor overlap (dark red) and 1 for complete overlap (yellow).
Figure S1: Concentration changes of the light alkanes over time in treatments Deep14, Deep25, and Surface25.Error bars represent standard deviation (n = 3).
Figure S2: Concentration changes of the heavy alkanes over time in treatments Deep14, Deep25, and Surface25.Error bars represent standard deviation (n = 3).
Figure S3: Concentration changes of the PAHs over time in treatments Deep14, Deep25, and Surface25.Error bars represent standard deviation (n = 3).

Author
Contributions: Conceptualization, methodology, and supervision E.G. and E.A.; validation, and formal analysis, E.G., G.C., E.F. and E.A.; microbial cultures, bioinformatic analysis and data curation G.C., A.B.D.M., K.A.K. and E.G.; chemical analysis, and data curation E.F. and E.A.; resources, N.K., N.P., P.N.P.; writing-original draft preparation, G.C. and E.F.; writing-review and editing, E.G., E.A., N.K., P.N.P., A.B.D.M., K.A.K., N.P.; funding acquisition, E.G. and E.A.All authors have read and agreed to the published version of the manuscript.Funding: This research was funded by the General Secretariat for Research and Technology (GSRT) and the Hellenic Foundation for Research and Innovation (HFRI) projects HEALMED with grant number 1874 (awarded to E.G.) and DEEPSEA with grant number 1510 (awarded to E.A.).Institutional Review Board Statement: Not applicable.Informed Consent Statement: Not applicable.

Table 2 .
Taxa that pass the limit of quantification (LOQ) threshold and can be characterized with confidence as specialists or generalists by their Levin's BN indices.Taxa with values closer to 0 are considered specialists while those closer to 1 are considered generalists.