Variation of Oxygenation Conditions on a Hydrocarbonoclastic Microbial Community Reveals Alcanivorax and Cycloclasticus Ecotypes

Deciphering the ecology of marine obligate hydrocarbonoclastic bacteria (MOHCB) is of crucial importance for understanding their success in occupying distinct niches in hydrocarbon-contaminated marine environments after oil spills. In marine coastal sediments, MOHCB are particularly subjected to extreme fluctuating conditions due to redox oscillations several times a day as a result of mechanical (tide, waves and currents) and biological (bioturbation) reworking of the sediment. The adaptation of MOHCB to the redox oscillations was investigated by an experimental ecology approach, subjecting a hydrocarbon-degrading microbial community to contrasting oxygenation regimes including permanent anoxic conditions, anoxic/oxic oscillations and permanent oxic conditions. The most ubiquitous MOHCB, Alcanivorax and Cycloclasticus, showed different behaviors, especially under anoxic/oxic oscillation conditions, which were more favorable for Alcanivorax than for Cycloclasticus. The micro-diversity of 16S rRNA gene transcripts from these genera revealed specific ecotypes for different oxygenation conditions and their dynamics. It is likely that such ecotypes allow the colonization of distinct ecological niches that may explain the success of Alcanivorax and Cycloclasticus in hydrocarbon-contaminated coastal sediments during oil-spills.

Members of these ubiquitous genera occupy distinct trophic niches, where usually the aliphatic hydrocarbon degrader Alcanivorax blooms first, followed by the (poly-)aromatic hydrocarbon degrader Cycloclasticus, as polyaromatic hydrocarbons (PAH) are less amenable to degradation (Head et al., 2006).
Despite our improved knowledge of the ecology of MOHCB, their behavior when confronted with fluctuations in environmental parameters is far from understood. Coastal marine sediments are subjected to fluctuating oxygenation and thus redox conditions from tidal cycles, diurnal cycles (and thus photosynthetic oxygenation), and macrofaunal burrowing activities that in turn affect microbial degradation processes (McKew et al., 2013;Cravo-Laureau and Duran, 2014;Duran et al., 2015a,b). The time span of these oxygen intrusions often varies from few minutes to several hours (Wakeham and Canuel, 2006;Militon et al., 2015Militon et al., , 2016. Previous studies demonstrated that anoxic/oxic oscillations promote organic matter biodegradation (Abril et al., 2010), and more particularly hydrocarbon degradation Vitte et al., 2011Vitte et al., , 2013. However, the microbial processes underlying the biodegradation of hydrocarbons within the anoxic/oxic transitional zone is not well understood Cravo-Laureau and Duran, 2014), especially for MOHCB subjected to such extreme fluctuating conditions.
The interaction of microorganisms with their environment is a key question in microbial ecology. Microorganisms have developed several metabolic and behavioral strategies to survive under different environmental conditions (Shade et al., 2012), which include physiological versatility and plasticity (Dolla et al., 2006;Fourcans et al., 2008;Evans and Hofmann, 2012) and dormancy (Lennon and Jones, 2011). Alternatively, the adaptation could involve ecotypes, subpopulations with specialized adaptations to microenvironments (Hunt et al., 2008;Coleman and Chisholm, 2010). The ecotype formation is an adaptive process allowing subpopulations to occupy an ecological niche, a first step in the speciation process (Wiedenbeck and Cohan, 2011). An impressive example is provided by the multiple ecotypes of the cyanobacterium Prochlorococcus, which was found in different abundances according to environmental factors such as light, temperature, and phosphate and nitrate contents (Martiny et al., 2009). Such ecotype diversification supports the success of Prochlorococcus as a dominant phototroph in oligotrophic seawater (Martiny et al., 2006(Martiny et al., , 2009. Ecotypes are recognized as phylogenetic clusters occupying a specific habitat. The existence of ecotypes has be shown by the habitat specificity of 16S rRNA sequences, as demonstrated for Synechococcus adapted to different temperatures (Melendrez et al., 2011) and for the PAH-degrading Alteromonas adapted to different depth of seawater (Math et al., 2012). Based on the oligotyping approach, a computational method that reveals the micro-diversity within sequences clustering in a single OTU (Eren et al., 2013), Kleindienst et al. (2016) proposed ecotypes for the hydrocarbon-degraders Cycloclasticus, Colwellia and Oceanospirillaceae adapted to hydrocarbon gradients.
In order to test the hypothesis that the success of MOHCB, particularly the most ubiquitous Alcanivorax and Cycloclasticus, to colonize distinct niches is based on ecotype diversification, an experimental ecology approach that offers the possibility to test ecological hypotheses under controlled conditions (Cravo-Laureau and Duran, 2014) was implemented. A microbial hydrocarbon-degrading community was maintained in bioreactors and exposed to different oxygenation regimes including permanent anoxic, anoxic/oxic oscillation and permanent oxic conditions. The bacterial diversity dynamic was assessed by16S rRNA gene transcript sequences (using high throughput sequencing technology) and the micro-diversity of Alcanivorax-and Cycloclasticus-related sequences was examined. First, oligotypes were defined and correlated with the redox conditions, suggesting the presence of subpopulations. Then, single nucleotide difference analysis confirmed the presence of cohesive subpopulations for the different oxygenation conditions representing specific ecotypes. The distribution of the ecotypes throughout the incubations explained the dynamics of Alcanivorax and Cycloclasticus during the anoxic/oxic oscillations.

Bioreactors System and Experimental Set Up
The slurry (1.8 L) was distributed in nine bioreactors (with a working volume of 2 l) to apply the three following conditions in triplicate: anoxic/oxic oscillation, permanent oxic and permanent anoxic conditions as detailed (Terrisse et al., 2015). The incubations were carried out for 15 days in batch conditions, with stirring (250 rpm, Stuart SS20), in the dark and at room temperature (ranging from 17 to 24 • C, InPro6800 sensors, Mettler Toledo International Inc.). Bioreactors were maintained for 5 days, under oxic or anoxic conditions (for the permanent anoxic and the anoxic/oxic oscillation conditions) for a microbial community stabilization period, before addition of crude oil and starting anoxic/oxic oscillations for the oscillating conditions. At day 5, Russian export blend crude oil (REBCO) was added with a concentration of 21.2 ± 5.7 mg/g of dry weight sediments. The REBCO is an Ural type crude oil, distilled at 110 • C to eliminate the more volatile hydrocarbon compounds. This oil contained 59.9% of saturated and 24.8% of aromatic hydrocarbons, 10.2% of resins and 5.1% of asphaltenes. In order to be close as possible to conditions prevailing in the environment, where oxygen pulses often occur at timescales of minutes to hours (Wakeham and Canuel, 2006;Militon et al., 2015Militon et al., , 2016, anoxic/oxic oscillations consisted of an alternation of anoxic periods with two 1-day periods of aeration performed at days 7 and 10. Aeration periods were produced by injection of filtered air (Acro R 37 TF Vent Device with a 0.2 µm PTFE Membrane, Life Sciences) using air pumps (flow rate of 70 L/h, Rena air 50) into the gas and aqueous phases, generating bubbling in the latter. Permanent oxic conditions were produced in the same way. Periods of anoxic conditions were achieved by stopping aeration, sealing, closing the system and creating a slight overpressure with nitrogen gas. The same process was used to achieve the permanent anoxic conditions. Dissolved oxygen, temperature (InPro6800 sensors, Mettler Toledo International Inc.), pH and redox potential (InPro 4260i/SG/225 sensors, Mettler Toledo International Inc.) were measured twice a day.

Sample Collection
Samples were collected with a sterile syringe (TERUMO Corporation), connected to bioreactors by a Norprene tube, for chemical and biological analyses. The sampling system was purged before each sampling. 10 mL samples were collected for chemical analysis from each reactor along the incubation at 5.4 (i.e., 10 h after oil addition), 7, 8, 10, 11, and 15 days of incubation. Sampling was performed before switching condition. Samples were stored in amber glass bottles with polytetrafluoroethylene stoppers (WHEATON) at −20 • C. For molecular analyses 1.5 mL of slurry was collected and immediately mixed with 190 µL of RNA stabilization solution (the RNA stabilization buffer was prepared as follow: 5 mL of phenol were mixed with 5 mL of 1 M Na-acetate buffer pH 5.5, then the two phases were separated by centrifugation at 4000 × g for 3 min. The phenolic phase was finally added to 95 mL of pure ethanol) to preserve rRNA transcripts integrity. Samples with RNA stabilization solution were homogenized and centrifuged at 10000 × g for 5 min at 4 • C (Jouan MR 1812). The supernatant was then removed and tubes containing pellets were introduced immediately in liquid nitrogen and stored at −80 • C. Samples were collected in triplicate from each reactor at days 0, 5 (before oil addition), 5.4 (10 h after oil addition), 7, 8, 10, 11, and 15.
Aqueous samples were analyzed after addition of internal standards (perdeuterated PAHs), extracted according to the Stir Bar Sorptive Extraction protocol and analyzed by SBSE-GC/MS as previously described (Stauffert et al., 2013).
Chemstation software was used for determining concentrations of hydrocarbons. The total petroleum hydrocarbons and the target molecules of n-alkanes and PAHs were quantified relatively to the perdeuterated eicosane and PAHs (internal standards) using calibration curves of REBCO crude oil (from 0.05 to 5 mg/mL), of n-alkanes (TRPH Standard from 0.5 to 50 µg/mL, Ultra Scientific, Florida) and of PAHs (CUS-9306 from 0.25 to 25 µg/mL, LGS Standards, France), respectively.

Total DNA and RNA Co-extraction
Genomic DNA and RNA transcripts were co-extracted using the commercial RNA PowerSoil R Total RNA Isolation Kit (MoBio Laboratories). Extractions were performed according to the manufacturer's instructions with a slight modification: step 4 was amended by suspending nucleic acids in 100 µL of SR5 buffer (supplied in the kit) after precipitation with ethanol 70% v/v. The separation and purification of DNA and RNA were performed with the commercial Allprep DNA/RNA Mini Kit (QIAGEN) following the manufacturer's recommendations. RNA and DNA were separately eluted in 100 µL of free DNase/RNase sterilized MilliQ water. RNA extracts were treated with Turbo DNA free kit (Ambion, Applied Biosystems) according to the manufacturer's instructions to ensure total elimination of DNA. DNA total elimination from RNA extracts was checked by negative 16S rDNA PCR amplification on the extracts. The quality and the size of DNA and RNA obtained were verified by electrophoresis on a 1% w/v agarose gel in Tris-Borate-EDTA buffer. RNA quality and concentration were also investigated by micro-capillary electrophoresis using RNA 6000 Nano LabChips kit and an Agilent 2100 Bioanalyzer (Agilent Technologies). RNA Integrity Numbers (RIN) between 7 and 8.5 were obtained for the RNA extracts, indicating good quality. The DNA and RNA extracts were aliquoted and stored at −80 • C.

Reverse Transcription of RNA
The synthesis of a complementary DNA strand (cDNA) from a RNA template strand was performed by reverse transcriptase as previously described (Stauffert et al., 2014b). Equipment and solutions used were certified RNase and DNase free or previously treated with DEPC (diethylpyrocarbonate). The reaction (final volume of 20 µL) was carried out with 10-60 ng of RNA per sample, adding 40 U of RNase OUT (Invitrogen by Life Technology), 0.5 mM of dNTPs, 10 ng/µL of random hexamers (Roche), 0.01 M DTT, 1 × enzyme buffer and 200 U of the Moloney murine leukemia virus reverse transcriptase (M-MLV, Invitrogen by Life Technology), following suppliers' advice. The cDNA products were then used directly in a PCR reaction or stored at −80 • C.
Quantitative-PCR of 16SrDNA on DNA and cDNA DNA and cDNA were amplified using primers 338F (5 -CTCCTAC GGGAGGCAGCAGT-5 ) and 518R (5 -GTATTACCGCGGCTGCTG-3 ) targeting bacterial 16S rRNA as previously described (Giloteaux et al., 2010). LightCycler R 480 SYBR GREEN I Master was used to prepare 10 µL reactions containing 5 µL of 2 × MasterMix, 0.4 µM of each primer, 2.5 µL of 1/100 diluted DNA or cDNA extracts. Quantitative Polymerase Chain Reaction (Q-PCR) was performed with the Roche LightCycler 480 Real Time PCR system. The cycling program was as follows: 95 • C for 5 min followed by 40 cycles with denaturation step at 95 • C for 15 s, hybridization step at 60 • C for 15 s and elongation step at 72 • C for 20 s. At the end of the program, an increase of the temperature from 64 • C to 97 • C provided the melting curves, informing on the amplification quality. The LightCycler 480 Software was used to analyze fluorescence. 16S rRNA copy numbers per µL from DNA or cDNA amplifications were quantified relatively to calibration curves as previously described (Paissé et al., 2012).

Sequencing and Data Analysis
High-throughput sequencing analysis, targeting bacterial 16S rRNA, was performed from 0.1 µg cDNA at days 0 (mix from 3 extract samples for day 0), 5, 5.4, 7 (under oxic and anoxic conditions), and 8, 10, 11, 15. Amplicon 454-pyrosequencing was performed by the Molecular Research DNA laboratory (MR DNA, Shallowater, TX, United States) following process originally described by Dowd et al. (2008). The 16S universal bacterial primers 27Fmod (5 -AGRGTTTGATCMTGGCTCAG-3 ) and 519Rmodbio (5 -GTNTTACNGCGGCKGCTG-3 ) were used, targeting V1-V3 variable regions. A single-step PCR using HotStarTaq Plus Master Mix Kit (Qiagen, Valencia, CA, United States) was used with the following conditions: 94 • C for 3 min, followed by 28 cycles of 94 • C for 30 s, 53 • C for 40 s and 72 • C for 1 min, after which a final elongation step at 72 • C for 5 min was performed. The different PCR products were mixed in equal concentrations before purification using Agencourt Ampure beads (Agencourt Bioscience Corporation, Beverly, MA, United States). Roche 454 FLX titanium instruments and reagents were used to perform the sequencing following manufacturer's guidelines.
The open source software QIIME (Quantitative Insights Into Microbial Ecology) was used for sequence read analysis of the bacterial 16S rRNA gene sequence (Caporaso et al., 2010). From the raw sequencing output, data cleaning was initially performed (Quince et al., 2009) by denoising of the data and by eliminating chimera and sequences with length under 450 bp, with mistakes in primers sequences or with homopolymers. Sequences reads were aligned on a sequence length of 400 bp. The data OTU picking was carried out according to the Usearch method (Edgar, 2010) defining reference sequences at the similarity threshold of 0.97. The taxonomy assignment was performed comparing reference sequences to a reference database of known 16S rRNA genes, the Ribosomal Database Project (RDP) database. Results of the number of times an OTU was found in each sample were tabulated and the taxonomic predictions were added for each OTU. The OTU abundances for each sample were rarefied to the same number corresponding to the minimum read number observed for a sample (3,102 sequences). This normalization of the OTU abundance data per sample was performed with the package vegan (Oksanen et al., 2013; Version: 2.0-10) with R software. The complete dataset was deposited in the NCBI Sequence Read Archive (SRA) database and is available under the Bioproject ID PRJNA383383.

Statistical Analyses
The effects of the incubation conditions (anoxic/oxic oscillation, permanent oxic or permanent anoxic conditions), time and interactions among these factors, on the variability of gene expression rates, biodegradation ratios, coverage values and univariate diversity indexes (richness and Shannon's diversity estimated with R and Mothur software), were tested by two multifactorial analysis of variance. A linear mixed model was performed with the effects of the incubation condition and the time (and their interactions) as fixed factors and a random factor to take into account the effect of repeated measurements over time, corresponding to the non-independent nature of the various samples in bioreactors (lme function in R package nlme). A second linear model was conducted by testing the same factors as fixed factors but without taking into account the random factor (lm function in R). To test the effect of the random factor, analysis of variance (Kunin et al., 2008) has shown no significant difference between these two models, the linear model without repeated measurement was chosen. The logarithmic transformation of the numeric variable was sometimes necessary in order that the residuals follow a normal distribution. ANOVA and Tukey's HSD tests (Honest Significant Difference) were carried out following the linear models to compare between each other the different modalities of the various factors and highlight significant differences. Permutational multivariate analysis of variance (PerMANOVA) tested significance among incubation conditions and time as fixed factors and the random factor to take into account the effect of repeated measurements over time, with 16S rRNA sequencing data (based on dissimilarity matrix of the Bray-Curtis distances). Bacterial community composition (16S rRNA sequencing data) under the different incubation conditions was analyzed using non-metric multidimensional scaling (nMDS) implemented in Primer 6 (version 6.1.16). For nMDS ordinations, the Bray-Curtis distance was used to generate dissimilarity matrices. Confidence ellipses were based on cluster analyses (Primer 6 software; version 6.1.16). Schematic representations of OTU distribution according to incubation conditions and time were drawn with ade4 package in R software (Dray and Dufour, 2007).

Micro-Diversity Analyses
Alcanivorax and Cycloclasticus were perfect candidates for microdiversity analysis, since they were largely present within our samples playing a key role in oil-degrading microbial assemblage. First, subpopulations were examined by oligotyping (Eren et al., 2013). Oligotyping analysis generates oligotypes by systematically identifying nucleotide positions that represent information-rich variation among closely related sequences. The identification of nucleotide positions of interest was performed using the Shannon entropy (Shannon, 1948). A total of 53,842 and 26,120 sequences were extracted respectively from the dominant OTUs (represented by more than 100 sequences) of Alcanivorax ). An alignment was performed against the SILVA database (release 128) and sequences were trimmed to a consistent start and end position. Using the Shannon entropy analysis, a total of 35 and six information-rich positions contributing to the Alcanivorax and Cycloclasticus oligotypes were identified. To reduce the noise, only oligotypes that occurred in more than 1% of the reads for at least one sample in which the most abundant unique sequence represented more than 0.01% of all reads, were retained. The analysis of the distribution of oligotypes was conducted from the raw counts data without normalization, because the identified oligotypes presented low diversity in each sample. Then, in order to resolve in depth the distribution of each genus within the different oxygenation conditions, a new clustering of Alcanivorax and Cycloclasticus was performed using the swarm algorithm (single-linkage clustering method) implemented in Qiime (Caporaso et al., 2010;version 1.9.1), with a local clustering threshold of 1 (d = 1). Thus, clusters of sequence were obtained differing by one nucleotide defining swarm OTUs that correspond to ecotypes. Only swarm OTUs composed of more than 100 sequences were retained. The distribution of these swarm OTUs was represented in heatmaps carried out with DECIPHER (R package; Wright, 2016) using the log-transformed data in order to deal with the skewed and wide distribution of the raw data. Since swarm OTUs represented cohesive subpopulations specifically correlated with an oxygenation condition they were defined as ecotypes.

Redox Status and Hydrocarbon Degradation
The redox and oxygen saturation, followed during the 15 days bioreactor-incubation of the hydrocarbon-degrading bacterial community, indicated that the hydrocarbon-degrading microbial community was effectively subjected to three oxygenation regimes, namely permanent oxic, permanent anoxic and anoxic/oxic oscillation conditions (Supplementary Figure 1). Under anoxic/oxic oscillation condition, the oxygenation periods at days 7 and 10 were characterized by a liquid phase with oxygen saturation above 50% and redox potentials higher 15 mV while the anoxic periods were characterized by a liquid phase with oxygen saturation at 0% and redox potentials bellow −200 mV as observed for the permanent anoxic condition. Such redox oscillations modify the availability of oxidants and reductants, which in turn influences microbial community structure (Lipson et al., 2015) and metabolic activities (DeAngelis et al., 2010), supporting increased microbial activity (Fenchel and Finlay, 2008). Previous studies demonstrated that anoxic/oxic oscillations promote organic matter biodegradation (Abril et al., 2010), and more particularly hydrocarbon degradation Vitte et al., 2011Vitte et al., , 2013. In our study the n-C 17 /pristane ratio (McKenna and Kallio, 1971) and phenanthrene/dimethylphenanthrene ratio (Michel and Hayes, 1999), used as indexes reporting the biodegradation for n-alkane and PAH respectively, indicated that the most efficient biodegradation for n-alkanes was under anoxic/oxic oscillation condition while similar efficiency was observed for PAH biodegradation under permanent oxic and anoxic/oxic oscillation conditions at the end of incubation (Supplementary Figure 2). It was notable that n-alkanes were progressively depleted under anoxic/oxic oscillation conditions while PAH biodegradation was more efficient during the aeration periods. This observation suggests co-stimulation of aerobic and anaerobic metabolism for the biodegradation of n-alkanes as proposed for organic matter degradation (Abril et al., 2010) while PAH biodegradation involved exclusively aerobic metabolisms. The latter implies that PAH-degrading microorganisms developed strategies for their maintenance under anoxic conditions. Such an assumption was supported in our study by a significant correlation (Pearson correlation index r = 0.904, p-value < 0.05) between PAH biodegradation index (phenanthrene/dimethylphenanthrene ratio) and expression rate (rRNA transcripts copies number/rRNA genes copies number) of 16S rRNA gene (Supplementary Figure 3) indicating growth stimulation of PAH-degrading microorganisms after oxygen input. It is thus likely that PAH-degrading microorganisms were in dormancy during the anoxic periods of the anoxic/oxic oscillation condition. Dormancy strategies play crucial role in microbial community stability in fluctuating environmental conditions, particularly by providing seed bank maintaining metabolic diversity through various mechanisms including persistent-subpopulation ecotypes, niche complementation and functional redundancy (Lennon and Jones, 2011;Shade et al., 2012). Dormancy strategies have been evidenced for marine microbial communities in coastal areas (Campbell et al., 2011;Hugoni et al., 2013). However, in a previous study we demonstrated that rhd gene transcripts, encoding ring hydroxylating dioxygenase involved in the first step of PAH biodegradation, were continuously produced during anoxic/oxic oscillation conditions including the anoxic periods (Vitte et al., 2013) suggesting a metabolic activity of some PAHdegrading populations even in absence of oxygen. It is likely that some of the hydrocarbon-degrading populations may adopt a dormancy strategy while another part would remain metabolically active. Determining the dynamic of the hydrocarbon-degrading populations during fluctuating redox conditions at a fine taxonomic scale may help to unveil the microbial processes involved in the adaptation to fluctuating redox conditions.

Influence of Oxygenation Regimes on the Overall Bacterial Community Organization
In order to determine the mechanisms underlying the behavior of hydrocarbon-degrading bacterial communities under different oxygenation regimes, bacterial community compositions were characterized at different incubation times for the three oxygenation conditions by high throughput sequencing of 16S rRNA gene transcripts, thus assessing active bacterial communities. After trimming and rarefication 3,102 sequences were obtained per sample corresponding between 177 and 600 OTU 97 s at the species level (97 % similarity threshold) per sample. The characteristics of the revealed diversity are presented in Supplementary Table 1. Good's coverage were above 0.90 indicating that the number of sequences per sample was sufficient to describe the accessible diversity. The α-diversity indexes [Richness (R) and Shannon (H)] varied significantly with the condition (PerMANOVA, p-value < 0.001), the time (PerMANOVA, p-value < 0.005 and <0.05 respectively) and time variation was different between conditions (PerMANOVA, p-value < 0.001). Furthermore, R and H were higher for permanent anoxic and permanent oxic conditions (PerMANOVA, R > 400 OTUs; p-value < 0.001; H > 3.6 and 4.0 respectively) than for anoxic/oxic oscillating conditions after the first period of aeration (R < 250 OTUs; H < 2.5). These observations corroborated the effect of oxygen regime on bacterial community organization, which was further supporter by nMDS analysis showing three main clusters (i) anoxic/oxic oscillating communities, (ii) permanent oxic communities and (iii) anoxic/oxic oscillating + permanent anoxic communities (Figure 1). The Venn diagram comparing bacterial communities revealed that 1,736 OTU 97 s were shared between the three conditions, which represent 46, 49, and 56% of total OTU 97 s for the permanent anoxic, permanent oxic and anoxic/oxic oscillating conditions respectively (Figure 2), suggesting that a large proportion of microorganisms have the capacity to tolerate the different oxygenation regimes. It is notable that these OTU 97 s represented above 70% of the sequences retrieved in each condition corresponding thus to the most abundant OTU 97 s. Among these shared OTU 97 s, the most abundant were related to the MOHCB genera Alcanivorax (37%) and Cycloclasticus (28%), some species of which related species have been described as aerobic and microaerobic (Dyksterhouse et al., 1995;Yakimov et al., 1998;Lai et al., 2013). Regarding the specific OTU 97 s for each condition, low abundant OTU 97 s that may play an important role in the microbial assemblage organization and functioning were found. Among the less abundant OTU 97 s (Supplementary Figure 4), sequences related to genera know to present hydrocarbon degradation capacities were found, such as Thalassolituus (Yakimov et al., 2007), Marinobacter (Duran, 2010), Marinobacterium (Sherry et al., 2013), Pseudomonas (Paisse et al., 2011) and Desulfobacterium (Paissé et al., 2008). Interestingly, the specific sequences were also related to Cycloclasticus and Alcanivorax OTU 97 s, which dominated the less abundant OTU 97 s in each condition, except in the permanent anoxic condition where specific Alcanivorax OTU 97 s were observed at low abundance (<50 sequences). It is important to note that under permanent anoxic conditions the specific OTU 97 s represented 2% of the sequences corresponding thus to rare OTU 97 s (Figure 2), which were dominated by Cycloclasticus and Ilyobacter. Members related to Ilyobacter have been found specialized in the degradation of hydroaromatic compounds under anoxic conditions (Brune et al., 2002). Both Alcanivorax and Cycloclasticus are widely distributed, detected in various marine environments including surface water, hydrothermal vents, deep sea water bodies, coastal and mudflat sediments (Bordenave et al., 2004a(Bordenave et al., , 2008McKew et al., 2007;Yakimov et al., 2007;Staley, 2010;Paisse et al., 2011;Coulon et al., 2012). It is not surprising to observe these microorganisms together because they use different hydrocarbon substrates as carbon and energy sources (McKew et al., 2007;Coulon et al., 2012). More intriguing is their presence in anoxic conditions, particularly for Cycloclasticus, species of which related have been described as aerobic bacteria. However, Cycloclasticus species have been detected in bio-irrigated coastal marine sediments (Montgomery et al., 2008) and isolated from marine polychaete burrows in an intertidal mudflat (Chung and King, 2001), suggesting that related members of this genus have the capacities to withstand fluctuating environmental conditions. Additionally, Cycloclasticus relatives have been found active in methane-enriched microcosms (Sauter et al., 2012) and environments (Redmond and Valentine, 2012). Recently, the metabolic versatility of Cycloclasticus in the degradation of hydrocarbons was demonstrated (Rubin-Blum et al., 2017) and Cycloclasticus ecotypes, specific for the deep sea methane-rich hydrocarbon plume arose with the Deepwater Horizon oil spill, have been proposed (Kleindienst et al., 2016) based on the oligotyping approach (Eren et al., 2013).

Influence of Oxygenation Regimes on Alcanivorax and Cycloclasticus Populations
The distribution of Alcanivorax and Cycloclasticus OTU 97 s within the different oxygenation conditions was examined in order to determine whether specific OTU 97 s explain their survival under the distinct conditions. Alcanivorax and Cycloclasticus were represented by six and eleven OTU 97 s respectively, the relative abundance of which varied between the different conditions during the incubation period (Figure 3). It is important to note that the distribution of Alcanivorax and Cycloclasticus OTU 97 s did not show difference between biological replicates (PerMANOVA p-value > 0.05). Alcanivorax and Cycloclasticus genera adopted different strategies to face fluctuating conditions: Alcanivorax taking advantage of brief favorable condition (oxygenation) for growth under fluctuating conditions whereas Cycloclasticus requiring stable conditions to attain its maximun abundance (Figure 3).
Alcanivorax was the most abundant genus under anoxic/oxic conditions despite the fact that relatives were not detected in the metabolically active bacterial communities at the beginning of the incubations. They were below the detection limit under the permanent anoxic incubation. Alcanivorax related species have been described as nitrate reducing bacteria (Yakimov et al., 1998), but here the redox potential was obviously too low to allow their development. Two OTU 97 s (ID.A37 and ID.A2593) were found abundant under the permanent oxic condition and poorly represented under the other conditions ( Figure 3A). These OTU 97 s were closely related to Alcanivorax sp. NBRC 102021 (acc. no. AB681668) isolated from seawater for ID.A2593 and Alcanivorax sp. OM-2 (acc. no. AB053128) isolated from oiled marine sediments for ID.A37. Interestingly, four Alcanivorax OTU 97 s (ID.A0, ID.A9793, ID.A8753 and ID.A485), although they were either not detected or in very low abundance under the permanent oxic condition, were found with high relative abundances following periods of aeration under the anoxic/oxic oscillating condition (Figure 3A). This observation suggested that an episodic presence of oxygen was favorable for Alcanivorax growth. Episodic intrusion of oxygen in anoxic zones has been shown to support aerobic metabolism in typical anoxic environments (Ulloa et al., 2012). It is likely that such Alcanivorax related OTU 97 s lack the capacity to outcompete strict aerobes in the oxic condition and adopt a fast growing strategy (r-strategist) when oxygen becomes available in the anoxic condition, which is a typical lifestyle in unstable and unpredictable environments (Andrews and Harris, 1986). These OTU 97 s were closely related to a sequence retrieved in desalination plant (IDA484: clone UV-RV-025, acc. no. HQ326427) and a sequence detected in deep-sea petroleum contaminated sediment (ID.A0, ID.A9793, ID.A8753: clone TVG01-83, acc. no. KF545057). The presence of different Alcanivorax related OTU 97 s is in accordance with previous reports showing the functional redundancy of hydrocarbon degradation within the Alcanivorax genus . Alcanivorax species have been shown to exploit distinct hydrocarbon substrates as carbon and energy sources following different physiological strategies  and presenting different susceptibilities to hydrostatic pressure . Such metabolic diversity may explain Relative abundances of the specific OTU 97 s for each condition. The analysis was performed at the genus level applying a threshold similarity of 97% for OTU identification (OTU 97 s). The analysis includes all samples covering the whole incubation period. The analysis is based on biological triplicates. The other * group combine the genera related to rare OTUs that are represented by less than 1% of total sequences per sample. The relative abundances of the OTU 97 s belonging to the "other group" are presented in Supplementary Figure 4. niche differentiation as observed for Alcanivorax phylotypes, phylotype SK2 occupying floating biofilms while phylotype OM-2 living in the sediment (Coulon et al., 2012). In our study, the dissimilar behavior shown by OTU 97 s ID.A9793 and ID.A8783 between the two oxygenated phases suggested the presence of subpopulations. Cycloclasticus, which was also among the most abundant genus, was represented by OTU 97 s showing different behavior (Figure 3B). Eight OTU 97 s (ID.C9747, ID.C3, ID.C8389, ID.C4039, ID.C7228, ID.C1181, ID.C178, ID.C2805) were detected under all conditions, but they were most abundant under permanent oxic condition. These OTU 97 s were related to Cycloclasticus pugetii 15BN12L-10 (acc. no. KF470997) isolated from Artic Ocean deep-sea hydrocarbon-contaminated sediment (ID.C9747, ID.C3 and ID.C8389), sequences of clones (acc. nos. AM882527, EU438147 and KJ094255) detected in petroleumcontaminated coastal marine sediments (ID.C4039, ID.C178 and ID.C2805 respectively), sequences of clones (acc. nos. FJ981469 and FJ980914) obtained from deep-sea hydrothermal plume (ID.C7228 and ID.C1181 respectively). Two OTU 97 s (ID.C5592 and ID.C7005) were detected under both permanent anoxic and permanent oxic conditions while the OTU 97 s ID.C9052 was not detected under the permanent anoxic condition. These OTU 97 s were closely related to Cycloclasticus pugetii 15BN12L-10 (ID.C7005: acc. no. KF470997), extracellular symbiont BG-C1 of Benthomodiolus (ID.C5592: acc. no. AB679348) and a sequence obtained in petroleum-contaminated coastal sediment (ID.C9052: acc. no. FM242294). Cycloclasticus was more sensitive to anoxic/oxic oscillations than Alcanivorax (Figure 3). Although this condition was unfavorable for its optimal growth, Cycloclasticus showed almost constant low abundances with only slight increases during the oxygenated phases, behavior that is more typical of K-strategist with slow growth. Such a strategy may explain the metabolic activity of some PAH-degrading populations that were shown previously to express rhd gene transcripts even in absence of oxygen during anoxic/oxic oscillations (Vitte et al., 2013). Because Cycloclasticus species been described as obligate aerobes (Staley, 2010), the presence of a single OTU 97 in all conditions may FIGURE 3 | Relative abundances of Alcanivorax OTU 97 s (A) and Cycloclasticus OTU 97 s (B) according to oxygenation regimes. The heatmap was performed at the species level applying a threshold similarity at 97% for OTU identification (OTU 97 s). The OTU 97 s abundances were normalized by rarefication to the lowest read count of a sample. Phylogenetic trees are shown on the left (the bars represent 0.8% estimated sequence divergence). ID sequences and the highest BLAST hits are indicated on the right. Color legend shows OTU 97 s abundance. The analysis is based on biological triplicates. * indicates oxygenation period under anoxic/oxic oscillation condition.
suggest the existence of distinct ecotypes, which correspond to a subpopulation that has acquired the genetic capacity to inhabit a slightly different ecological niche (Konstantinidis and Tiedje, 2005).

Micro-Diversity
The micro-diversity of Alcanivorax and Cycloclasticus was examined in order to determine whether ecotypes explain their distribution in ecological niches distinguished by their oxygen availability. The resolution of 16S rRNA gene transcript sequences at the subpopulation level allows the identification of specific ecotypes inhabiting distinct ecological niches (Eren et al., 2013;Tikhonov et al., 2015). We acknowledge the inherent biases of 16S rRNA gene-based approaches to define ecotypes, particularly the limited capacity in defining phylogenetic cohesive populations (Berry et al., 2017) and functional traits (Martiny et al., 2013). But such 16S rRNA-based approaches provide useful information to draw hypotheses that could explain observed ecological behavior of hydrocarbon-degrading bacteria.
hydrocarbon polluted coastal sediment (acc. no. AM882527 associated to OTU97 ID.C4039). Interestingly, oligotype C4 emerged specifically after the second oxygenation phase under the anoxic/oxic oscillation condition (Figure 4B), which was observed within all three biological replicates. This oligotype was affiliated to Cycloclasticus pugetii (acc. no. KF470997 associated to OTU 97 s ID.C3 and ID.C9747). The oligotypes explained only partly the behavior of Alcanivorax and Cycloclasticus related species to withstand oxygenation conditions. For example, it is not clear whether specialized ecotypes drive the success of Alcanivorax under anoxic/oxic oscillation since oligotype A1 was found dominant during all incubation period for both permanent anoxic and anoxic/oxic oscillation conditions. Similarly, although the distribution of Cycloclasticus oligotypes showed specialized ecotypes for anoxic permanent condition, the distribution of oligotypes C1 and C3 found under both anoxic and oxic samples suggested non-homogeneous subpopulations. We thus explored the micro-diversity of Alcanivorax and Cycloclasticus more deeply (swarm analysis with a dissimilarity threshold = 1), which revealed 10 and 23 subpopulations respectively (Figure 5). We assume that these subpopulations represent ecotypes since they occupy distinct ecological niches as reported by Tikhonov et al. (2015), who defined ecological subpopulations differing by one nucleotide in their 16S rRNA gene sequence. Furthermore, swarm analysis revealed cohesive subpopulations specifically observed under an oxygenation condition (Figure 5). PerMANOVA did not reveal difference in the distribution of Alcanivorax and Cycloclasticus ecotypes within biological replicates (p-value > 0.5). Alcanivorax ecotypes were related to the OTU 97 s ID.A9793 and ID.A8753 whereas the OTU 97 s ID.A0 and ID.A37 (ecotypes A1 and A2 respectively) showed uniform and cohesive populations irrespective of the oxygenation condition ( Figure 5A). The ecotypes' succession reflected the distinct ecological niches occurring throughout the incubation under the anoxic/oxic oscillation condition characterized not only by oxygen availability but also by metabolites appearing during aliphatic hydrocarbon degradation as well as the presence of other organic substrates. Ecotypes A6 (91.5% of oligotype A6) and A8 (100% of oligotype A3) related to OTU 97 ID.A9793 emerged after the first oxygenation period while ecotype A5 (100% of oligotype A6) related to OTU 97 ID.A8753 was observed after the second oxygenation period (Figure 5A). Ecotypes A3 (100% of oligotype A1) and A7 (94.6% of oligotype A1) bloomed during the second anoxic period while ecotypes A4 FIGURE 5 | Relative abundances of Alcanivorax ecotypes (A) and Cycloclasticus ecotypes (B) according to oxygenation regimes. The heatmap was performed at the sub-species level (swarm dissimilarity threshold = 1 on log-transformed data) representing ecotypes, which correspond to specific cohesive subpopulations correlated to an oxygenation condition. Each ecotype represents more than 100 sequences. The color legend at the left indicates the ecotype abundance (sequence counts). ID sequences and their affiliation (highest BLAST hits) are indicated on the right. The analysis is based on biological triplicates. * indicates oxygenation period under anoxic/oxic oscillation condition.
(100% of oligotype A5), A9 (100% of oligotype A4) and A10 (100% of oligotype A6) appeared at the end of the incubation after the third anoxic phase. Exploring the micro-diversity at one nucleotide difference level allowed the identification of specific ecotypes occupying distinct ecological niches during the incubation under anoxic/oxic oscillation condition. It is likely that such Alcanivorax ecotypes may take advantages of oxygen-fluctuating hydrocarbon-contaminated environments allowing them to compete with other alkane-degraders such as Marinobacter and Thalassolituus detected in our study at lower abundances (0.27 and 0.51% of total OTU 97 s, respectively), and probably with other alkane-degraders such as Oleibacter, Oleispira and Colwellia (all below 0.005% in our study), which have been detected in hydrocarbon-contaminated environments (Coulon et al., 2012) including during Deepwater Horizon oil spill (Kleindienst et al., 2016;Yang et al., 2016a).

CONCLUSION
An investigation of the behavior of Alcanivorax and Cycloclasticus, the most widely distributed MOHCB, under well-controlled different oxygenation regimes, revealed that members of these genera adopted distinct strategies to develop under oxygen-fluctuating conditions in oiled sediments. Anoxic/oxic oscillations were more favorable for Alcanivorax, which was more abundant than under the other conditions, suggesting that Alcanivorax behaved as a typical r-strategist. In contrast, Cycloclasticus abundance was lower in such fluctuating conditions compared with both permanent anoxic and permanent oxic conditions. Oligotyping revealed that the distribution of Alcanivorax and Cycloclasticus subpopulations (oligotypes) correlated with the redox conditions, but the analysis was unable to fully explain how these hydrocarbon-degraders may withstand anoxic/oxic oscillations. Further micro-diversity analysis at one nucleotide difference level (swarm analysis) allowed the identification of specific subpopulations, which we assume correspond to ecotypes occupying distinct ecological niches during the incubation under anoxic/oxic oscillation condition characterized by substrates and oxygen availabilities. Such ecotypes allow the colonization of distinct ecological niches that may explain the success of these MOHCB genera during an oil-spill. However, further efforts are required to isolate and characterize Alcanivorax and Cycloclasticus ecotypes to gain new insights on their ecological role within microbial networks involved in oil degradation in marine environments.

AUTHOR CONTRIBUTIONS
RD, CC-L, and CC conceived and designed the study. RD, CC-L, CC, and FT ran the experiments. RD, CC-L, CC, FT, CN, AD, and TMG analyzed the resulting data. RD wrote the manuscript. CC-L, CN, and CC revised the manuscript.

FUNDING
We thank the French National Research Agency (ANR) for their support through the DECAPAGE (ANR 2011 CESA 006 01) project. We acknowledge the financial support by a Ph.D. grant from the French Ministry of Higher Education and Research to CN.