Impact of Oil on Bacterial Community Structure in Bioturbated Sediments

Oil spills threaten coastlines where biological processes supply essential ecosystem services. Therefore, it is crucial to understand how oil influences the microbial communities in sediments that play key roles in ecosystem functioning. Ecosystems such as sediments are characterized by intensive bioturbation due to burrowing macrofauna that may modify the microbial metabolisms. It is thus essential to consider the bioturbation when determining the impact of oil on microbial communities. In this study, an experimental laboratory device maintaining pristine collected mudflat sediments in microcosms closer to true environmental conditions – with tidal cycles and natural seawater – was used to simulate an oil spill under bioturbation conditions. Different conditions were applied to the microcosms including an addition of: standardized oil (Blend Arabian Light crude oil, 25.6 mg.g−1 wet sediment), the common burrowing organism Hediste (Nereis) diversicolor and both the oil and H. diversicolor. The addition of H. diversicolor and its associated bioturbation did not affect the removal of petroleum hydrocarbons. After 270 days, 60% of hydrocarbons had been removed in all microcosms irrespective of the H. diversicolor addition. However, 16S-rRNA gene and 16S-cDNA T-RFLP and RT-PCR-amplicon libraries analysis showed an effect of the condition on the bacterial community structure, composition, and dynamics, supported by PerMANOVA analysis. The 16S-cDNA libraries from microcosms where H. diversicolor was added (oiled and un-oiled) showed a marked dominance of sequences related to Gammaproteobacteria. However, in the oiled-library sequences associated to Deltaproteobacteria and Bacteroidetes were also highly represented. The 16S-cDNA libraries from oiled-microcosms (with and without H. diversicolor addition) revealed two distinct microbial communities characterized by different phylotypes associated to known hydrocarbonoclastic bacteria and dominated by Gammaproteobacteria and Deltaproteobacteria. In the oiled-microcosms, the addition of H. diversicolor reduced the phylotype-richness, sequences associated to Actinobacteria, Firmicutes and Plantomycetes were not detected. These observations highlight the influence of the bioturbation on the bacterial community structure without affecting the biodegradation capacities.


Introduction
Crude oils are complex mixtures of hydrocarbons, asphaltenes, resins, and associated metals, mainly vanadium and nickel [1], which may have serious environmental toxic effects. The biodegradation of hydrocarbons by microorganisms is a crucial process by which petroleum and other hydrocarbon pollutants are eliminated from the environment [2].
The fate of hydrocarbons in the environment is not only dependent on microbial activities, but also on sediment characteristics, physical-chemical factors (electron acceptors, temperature), nutrients, co-metabolite availability, and bioturbation [3]. In particular, bioturbation is created and aided by the benthic macrofauna which physically mixes sediments and introduces oxygen into sediments by burrow ventilation may greatly affects the biogeochemical cycles of nutrients [4], the fate of contaminants and the microbial metabolisms [4]. Several studies have demonstrated that bioturbation processes play an important role in the burial and the degradation of aliphatic hydrocarbons [5], PAHs [6] and acyclic hydrocarbons [7]. In addition, polychaetes, a class of marine annelid worms, have the enzymatic potential to solubilize aromatic compounds [8,9]. Burrowing activity can also affect the overall microbial metabolic activities such as nitrifying [10] and sulfate-reducing activities [11,12]. It has been shown that the increase in microbial abundances and activities is due to complex biogeochemical interactions induced by bioturbation [13][14][15]. Previous reports suggest that bioturbation by macrofauna drives microbial community structure and activities, which can also be affected by the presence of hydrocarbons. To the best of our knowledge, only one study reported specific bacterial patterns due to bioturbation activities in the presence of hydrocarbons [16]. The lack of studies is most likely due to the difficulty in creating experiments for studying natural bioturbated microbial communities subjected to oil pollution. To overcome these difficulties, different strategies have been developed by performing either in situ studies [7] or microcosm [16] studies using sediments from which macrofauna have been eliminated and then reintroduced. The former is not a suitable strategy to investigate the effect of a pollutant since in situ oil spiking requires special authorization and special technical means for oil containment, while the latter would modify the sediment structure. Our strategy was to increase bioturbation activity by adding polychaetes to bioturbated mudflat sediment in order to investigate how oil affects the bacterial community structure and composition. We hypothesized that the addition of polychaetes stimulates the bioturbation activity which could result in the selection of a particular bacterial community with an increased biodegradation capability. Thus, an original microcosm system maintaining the structure of muddy sediments under tidal cycles was set up in order to ensure conditions close to those in the natural environment. A set of chemical and microbial analyses was performed over a nine-month period to monitor i) the removal of petroleum by estimating hydrocarbon contents ii) the macrofaunal reworking activity and iii) the structures and composition of bacterial communities. This multidisciplinary approach, combined with multivariate analyses, showed that a specific bacterial community was obtained in response to oil addition in bioturbated sediments. The results presented here showed that polychaetes strongly affected bacterial community structures without changing oil degradation capacities.

Sampling Area
Mud sediments for microcosm experiments were collected in autumn 2007 at the tidal basin Aber-Benoît, located in the Brittany region (Treglonou, France, 48u 33912.4099N; 4u 3298.6999W).
This sediment showed a pH of 7.8 with an organic carbon content of 6.8% and a salinity of 35%. The sediment moisture content was 60% and chemical analysis detected about 4 mg of biogenic hydrocarbons per gram of dry sediment in the upper layer of sediment (2 cm). Sediments were composed of 3% clay (,3.9 microns), 31% silt (3.9-62.5 microns) and 66% sand (62.5-2000 microns). The density and characterization of autochthonous macrofauna were determined from three cores collected from the sampling site. The macrobenthic community was composed of four species, the marine annelid worms Nephtys hombergi and Hediste (Nereis) diversicolor, the mollusk Hydrobia ulvae and the bivalve Mya arenaria, for an average total density of 241966.08 ind.m 22 . More specifically, the gallery-forming polychaete Hediste diversicolor represented 61 g.m 22 (509 ind.m 22 ).

Microcosm Setup
Sediments were sampled with a core collector (20-25 cm sediment depth), and transferred while maintaining their integrity, over geotextile membrane into microcosm boxes (65 cm650 cm641 cm; approximately 30 L of sediment). Twelve microcosm boxes were set up and connected to a device supplying natural, sand-filtered, and UV-treated seawater from Oceanopolis aquarium (Brest, France), thus ensuring constant sea water quality. Tidal cycles (12 hours) were applied using a water level control consisting of an up and down drainage system. Seawater was renewed with each tidal cycle (20 L of water per microcosm box).
At high tide, the water-level was maintained about 5 cm above the sediment surface using a faucet ballcock (Fig. 1). The microcosms were set up in a hall at room temperature (ranging from 10 to 20uC). They were not directly exposed to sunlight limiting evaporation and photodegradation of hydrocarbons. After 7 days of sediment stabilization, four conditions were applied in triplicate (N = 12) as follows: (i) CTRL: control condition, (ii) BAL: oil addition, (iii) NEREIS: addition of H. diversicolor and (iv) NEREIS+BAL, addition of oil and H. diversicolor. In order to standardize sediments, the 2 cm upper layer of each microcosm box was collected, pooled and mixed for homogenization. Half of the homogenized sediments were mixed with oil in order to obtain 25.6 mg of oil per gram of wet sediment. This oiled sediment was then spread in a 2 cm layer over the contaminated conditions (BAL and NEREIS+BAL) sediment surface. Thus, each oiled microcosms began the experiments with the same contamination level. We used Blend Arabian Light crude oil (BAL) which corresponds to a standardized Arabian crude oil, often used as model crude oil. The BAL oil used was distilled at 110uC to eliminate the more volatile hydrocarbon compounds; such oil is named BAL 110 oil (Arabian Light crude oil topped at 110uC). For the non-contaminated conditions (CTRL and NEREIS) the half un-oiled homogenized sediments was spread over the sediment surface (2 cm layer).
Hediste (Nereis) diversicolor, common polychaete worm found in mudflat sediments able to live within the first 10 cm of sediments, was added to the microcosms in order to increase burrowing activity and, thus, sediment reworking (NEREIS and NEREIS +BAL). Thus the H. diversicolor density was doubled (final density reaching 982 ind.m 22 ) but the biomass increased 10-fold (addition of 614 g.m 22 ), due to the addition of larger polychaete individuals.
Two hours after the experimental conditions were established, the tidal cycles were initiated. Subsamples of sediment were collected after 2, 7, 30, 90, 180 and 270 days of microcosm incubation (November 2007-July 2008). For chemical analyses, one core sample (Falcon tube, 2.5 cm68 cm) was collected for each microcosm (triplicates per condition). For molecular analyses, superficial samples (0-2 cm, Eppendorf tube, 2 mL) were collected in duplicate for each microcosm (sampling replicate). Thus, six replicates per condition at each sampling time were obtained. The sub-samples were collected 10 min. after reaching ''low tide'' randomly in order to consider the sediment's heterogeneity; the objective being to analyze the bulk sediments. The sampling hole was plugged with a plastic tube in order to maintain the sediment's structure. Samples were quickly frozen in liquid Nitrogen and stored at 280uC until analyzed at the end of the experiment when all samples were collected. It is important to note that the pH was stable around 7.5 during the experiments; there was no significant variation between high tides and low tides.
TPH contents were analyzed according to the standard NF EN ISO 9377-2 protocol using high temperature gas chromatography with a flame ionization detector (HTGC/FID) HP 6890N equipped with an autosampler. 1 mL of sample was injected in purged splitless mode (splitless time: 2 minutes, flow: 50 mL/min) at 425uC. The oven temperature was initially 30uC, then was increased at 20uC per minute 1 minute after injection to a final temperature of 415uC, then held for 30 minutes. The metallic chromatographic column was a DBHT-SIMD (Agilent J&W): 5 m60.53 mm ID60.15 mm film thickness. The helium carrier gas was maintained at a constant flow of 15 mL/min. Different concentrations of BAL 110 diluted with methylene chloride were used to calibrate the method and then to quantify the TPH extracted from the sediment samples.
Alkanes and aromatics were quantified using gas chromatography coupled with mass spectrometry (GC/MS). The GC/MS was an HP 6890N equipped with a split/splitless injector (Splitless time: 1 min, flow 50 mL/min) coupled to an HP 5973 Mass Selective Detector (MSD) (Electronic Impact: 70 eV, voltage: 1200 V). The injector temperature was maintained at 300uC. The interface temperature was 300uC. The GC/MS temperature gradient was from 50uC (1 min) to 300uC (20 min) at 5uC/min. The carrier gas was Helium at a constant flow of 1 ml/min. The capillary column used was an HP-5 MS: 30 m60.25 mm ID60.25 mm film thickness. n-Alkane and PAH semi-quantifications were performed using Single Ion Monitoring (SIM) mode with the most representative fragment (saturates) or the molecular ion (PAHs) of each compound at a minimum of 1.4 cycles/s. The TPH concentrations were tested by a three-way analysis of variance (ANOVA) in order to investigate the effects of time (T), the addition or non-addition of H. diversicolor (N), the depth (D) and all their interaction with the TPH contents of the microcosm boxes. A prior analysis was run to test the dependence of temporal observations (i.e repeated measures) with a linear mixed model with T,N,D and their interactions as fixed factors and the temporal pseudo-replication specified as random effect (using the package nlme in R, [18]). However, no temporal autocorrelation was detected and only the fixed factors were retained in the threeway ANOVA. Prior to statistical analysis, homoscedasticity was assessed by Hartley's Fmax test and data was log transformed in order to meet the test assumptions. Significant differences between means of groups were determined using a Duncan's pairwise multiple comparison test. The significance level was set at a = 0.05. Statistical data analyses were performed with STATIS-TICA software, version 5.5 (StatSoft, USA) and the R software (version 2.11.1).

Reworking Activity
The sediment reworking activity in the different experimental microcosm boxes was estimated using fluorescent particulate tracers (luminophores; [19]). One hour after the introduction of H. diversicolor into the NEREIS and NEREIS-BAL microcosm boxes, a suspension of 100 g of orange luminophores (63-125 mm) was homogeneously spread over the sediment surface of each experimental microcosm box. At the end of the experiment, one sediment core (I.D.: 10 cm, length: 15 cm) was sampled from each experimental microcosm box. Each core was then sectioned in 0.5cm thick layers from the surface down to 2 cm depth, then in 1-cm thick layers down to 10 cm. The sediment from each layer was removed separately and sieved through a 250 micrometer mesh to remove macrofauna from sediment particles. Sediments were then freeze-dried, gently crushed to powder and homogenized. Sediment subsamples were taken for luminophore counting under UV-light (digital camera Olympus C-2500L; image analysis software Image-Pro Plus). This resulted in a particulate tracer distribution profile over depth, for each sediment core. Then, sediment reworking (i.e. production of a biodiffusion-like coefficient Db describing particle transport in the whole sedimentary column) was quantified by applying a one-dimensional diffusion model (script developed in MatlabH) to the data [20].  kit (MoBio Laboratories) according to the manufacturer's recommendations. The protocol was modified as follows: sediment was centrifuged at 50006g for 15 min, to remove the water phase; the initial step of horizontal vortexing was performed for 30 min. DNA was eluted in 50 mL water and stored at 220uC.

RNA Isolation and Reverse Transcription
Total RNA extractions were performed from the 136 samples collected in the upper layer of sediment (2 cm, 2 g) using the RNA Power Soil kit (MoBio Laboratories) according to the manufacturer's recommendations. Only RNase-free certified plastic-ware was used and all solutions were prepared with sterile diethyl pyrocarbonate-treated (DEPC) water. DNA was removed from the total nucleic acid extraction by DNase digestion, using Turbo DNA-free TM protocols (Ambion). The complete removal of DNA was verified by a control PCR using primers of bacterial 16S rRNA directly on RNA extract (PCR conditions described below) and the RNA quality was controlled by capillary electrophoresis (Agilent Bioanalyzer). A concentration of 0.5 mg of DNase-treated RNA (Turbo DNase; Ambion) was reverse-transcribed with Moloney murine leukemia virus reverse transcriptase (M-MLV RT; USB Corporation) as previously described [21] with minor modifications including the use of random hexamers (5 mL at 0.1 U A260; Roche) and an incubation time of 2 h at 42uC. The RT products were used immediately for PCR amplifications and the remaining products were stored at 220uC. Possible DNA contamination of RNA templates was monitored by PCR amplification of RNA aliquots without a reverse transcription step. No DNA was detected in these reactions.

Amplification of Bacterial 16S rRNA Gene
16S rRNA genes were amplified from the 136 nucleic acid samples using the primers 63f (59-CAGGCCTAACACATG-CAAGTC-39) and 1387r (59-GGGCGGWGTGTAACAAGGC-39) specific for Bacteria [22]. The forward primer was labeled at the 59 end with phosphoramidite fluorochrome carboxyfluorescein (FAM). PCR conditions were as follows: initial denaturation (94uC for 3 min) followed by 35 cycles of denaturation (94uC for 30 s), annealing (58uC for 30 s), and extension (72uC for 60 s) and a terminal extension (72uC for 10 min). The reaction mix (50 mL final volume) contained 50 mM buffer, 0.4 mM dNTP, 0.2 mM of each primer, 1.25 U of Taq polymerase (Ozyme) and 10 ng of DNA template. PCR products were visualized by agarose gel electrophoresis and purified with the PCR purification kit (GE Healthcare).

T-RFLP Analysis
Purified 16S rRNA amplicons were separately digested by 3 U of HaeIII and HinfI restriction enzymes at 37uC (New Englands Biolabs) in a final volume of 10 mL for 3 h. A desalting step by dialysis was performed using a 0.022 mm pore filter (Millipore). Then, digested products (1 mL) were mixed with deionised formamide (10 mL) and 0.25 mL of internal size standard (TAMRA GS-500, Applied Biosystems). The samples were denatured by heating at 95uC for 5 min and then immediately transferred onto ice. Fluorescently labeled fragments were separated and detected with an ABI PRISM 310 capillary sequencer (Applied Biosystems) run in GenScan mode. Injection was performed electrokinetically at 15 kV for 30 s and the runs at 15 kV were completed within 30 min. T-RFLP profiles were analyzed using GENSCAN version 3.1 software (Applied Biosystems). T-RFLP profiles were normalized by calculating relative abundances of each peack from height fluorescence intensity. Only terminal fragments whose size ranged from 35 bp to 500 bp and whose height was greater than 30 fluorescence units were considered for analysis.
The comparison of DNA2 and RNA-based T-RFLP profiles was used to determine the active changes in the bacterial communities. For each sample, we calculated the percentage of OTUs that were present in both fractions (DNA+ RNA+), those were present only in DNA-based fractions (DNA+ RNA2) and those only present in RNA-based fractions (DNA2 RNA+). Relationships between the different fractions (RNA+DNA+, RNA2DNA+, RNA+DNA2) and hydrocarbon concentrations (TPH, n-alkanes and PAHs) across time were explored with simple linear regressions (lm function in R software).

Statistical Comparison of T-RFLP Profiles
The T-RFLP profiles were standardized by Dunbar's method [23]. Then, raw data were normalized by dividing each peak height by the total signal intensity of the corresponding profile. To deal with pseudo-replication within microcosm box (two sampling replicates within each microcosm box), the average over the pseudo-replicates were calculated in order to obtain three profile replicates per condition (biological replicates, BRs). The dissimilarities among BRs were estimated with a Bray-Curtis distance (Table S1 in the supplemental material). The effects of time, BAL addition and NEREIS addition on relative abundance of DNA Operational Taxonomic Units (OTUs) and cDNA OTUs in communities were tested with permutational multivariate analysis of variance (PerMANOVA; [24]). We performed a four-way PerMANOVA with time, BAL effect (BAL addition or not), and NEREIS effect (NEREIS addition or not) and BRs as additional factor included to account for repeated units. To assess this model, permutations were constrained within BRs, so that the effects of time and its interactions with other factors were taken into account (model I in Table S3 in the supplemental material). We also ran a three-way factorial PerMANOVA with BAL effect, NEREIS effect and their interaction and with time as an additional factor. In this model, permutations occurred within each sampling time in order to obtain a proper estimation of the BAL and NEREIS effects (model II in Table S3 in the supplemental material). These analysis were carried out with the adonis function in the vegan R package, [25] and based on untransformed abundance data, using a Bray Curtis distance and 999 random permutations (a = 0.05). Since most of terms (main factors and their interactions) of PerMANOVA model I and II were significant, two dimension non-metric multidimensional scaling (nMDS) ordinations (based on Bray Curtis distance) were used to visually represent dissimilarities between conditions at each sampling time and between sampling time for each conditions. The robustness of nM-MDS analysis was evaluated by calculating Kruskal Stress. Instead of performing Duncan post hoc pairwise comparisons, we determined the most similar communities across time and conditions with a hierarchical analysis performed with the average profiles from the three BR profiles. The groups were formed according to the similarity obtained between the BRs and the significant differences between groups were tested with an analysis of similarities (ANOSIM). Hierarchical clustering, nMDS analysis and ANOSIM were performed with Primer 6 software (Primer E, Plymouth, UK).

16S cDNA Library Construction and Analysis
16S cDNA libraries were obtained with unlabeled PCR products from 270 days for the four conditions (CONTROL, NEREIS, BAL, NEREIS+BAL). PCR products were obtained from the six replicates per condition and pooled before cloning. The pooled PCR products were cloned in Escherichia coli TOP10F' using the pCR2.1 Topo TA cloning kit (Invitrogen Inc.). The presence of insert DNA was checked by PCR using M13 primers (Eurogentec) surrounding the cloning site. For each cDNA library, 100 clones were taken randomly for sequence determination. The clones were sequenced by GATC Biotech SARL (Konstanz, Germany). The sequences were compared with the GenBank nucleotide database library by BLAST on-line searches [26]. The presence of chimeric sequences was checked with Mallard and Pintail tools. Multiple sequence alignment was performed using the CLUSTALX [27]. Phylogenies were constructed with the MOLECULAR EVOLUTIONARY GENETICS ANALYSIS v3.0 program [28] using the Kimura two-parameter model and the neighbor-joining algorithm. The significance of branching order was determined using bootstrap analysis with 1000 resampled data sets. The DOTUR program was used to determine Operational Taxonomic Units (OTUs) defined as sequence groups in which sequences differed by 3% for 16S rRNA sequences [29]. PAST (Paleontological Statistics v1.60) software was used to perform rarefaction analysis and calculate diversity indexes for each 16S cDNA library. The coverage value is given as C = 12(n1/N) where n1 is the number of clones which occurred only once in the library [30] and species richness was calculated using the Chao Estimator website (http://www2.biology.ualberta. ca/jbrzusto/rarefact.php) [31]. The compositions of bacterial communities were compared by clustering analysis using Pearson's correlation as a similarity distance and a heatmap was drawn with gplots package in R software [32]. In silico restrictions of clone sequences were performed using the Mobyle website (http:// mobyle.pasteur.fr/cgi-bin/portal.py) [33].

Nucleotide Sequence Accession Numbers
The sequences determined in this study have been submitted to the Genbank database under accession Nos. JF774416 to JF774826.

Characterization of the Biological Reworking Activity
The abundance of H. diversicolor was maintained throughout the experiment, and was not significantly affected by oil addition (NEREIS: 17826382 ind.m 22 ; NEREIS+BAL: 23346989 ind.m 22 ; p-value = 0.26). Biodiffusion-like coefficients (Db) varied according to the experimental conditions. The microcosms to which H. diversicolor were added showed higher Db (NEREIS: 0.706 0.08 cm 2 .y 21 ; Mean 6SD; n = 3) compared to the other microcosms (Db from 0.18 to 0.33 cm 2 .y 21 ) indicating a marked increase in burrowing with the addition of both H. diversicolor and oil (NEREIS+BAL: Db = 1.0760.12 cm 2 .y 21 ; Mean 6SD; n = 3). This observation indicated that the reworking activity was enhanced by the oil addition. The luminophore distribution indicated that most of the reworking activity took place within the top 2 cm (see Fig. S1 in the supplemental material). We thus analyzed microbial community in the top 0-2 cm layer, given this luminophore distribution, and because it is the most oxygenated zone where efficient hydrocarbon degradation is expected.

Fate of Petroleum Compounds
In order to assess the fate of BAL 110 oil, the total petroleum hydrocarbon (TPH) content was monitored during the 270 days of the experiment ( Fig. 2A, 2B). The distribution of TPH content according to depth highlighted the burrowing activity of H. diversicolor, particularly at an early stage when statistically significant differences between BAL and NEREIS+BAL condi-tions were observed for the 2-4 cm layer (p-value,0.05) ( Fig. 2A). Indeed, in the polychaete bioaugmented condition (NEREIS+-BAL), a mean of 23% of TPH was extracted from the 2-4 cm layer while the TPH content represented only 10% in the microcosm without H. diversicolor addition (BAL). The increase in the macrobenthic population and its related reworking activity clearly allowed deeper burying of TPH into the sediments.
The TPH concentration significantly varied with the time, the NEREIS addition, the depth and the different interactions among these factors (Table S2 in the supplemental material). According to the Duncan post-hoc test, the TPH concentration was significantly reduced irrespective of conditions after 90 days of incubation (BAL p-value = 0.0062; NEREIS+BAL p-value = 0.0048). By the end of incubation, 60% of the initial TPH had been removed (0-8 cm depth, see Fig. S2 in the supplemental material). Considering the 0-2 cm layer (Fig. 2B), no statistically significant difference was observed between the BAL and NEREIS+BAL conditions. 57% of the initial TPH was removed in the BAL condition and 48% was removed in the NEREIS+BAL condition. This result indicated that, throughout the 0-8 cm layer, an average of 87% of TPH removal was located within the 0-2 cm layer in both oilcontaminated conditions.
Focusing on the chemical composition of the biodegradable fraction (n-alkanes and PAHs; Fig. 2C), the n-alkane and the PAH contents significantly varied with time and the NEREIS addition (Table S2 in the supplemental material). The n-alkane content significantly reduced from day 180 (according to the Duncan test, p-value,0.05) in the BAL (68% at 180 days, 92% at 270 days) and NEREIS+BAL (57% at 180 days, 80% at 270 days) conditions, while PAHs were depleted at the end of the experiment only in the BAL condition (75% removal). In our experimental system, the removal of hydrocarbons could be explained by physical factors such as washing out due to the constant input of seawater for tidal simulation. However, we assume that washing out was a limited phenomenon because hydrocarbons were present until the end of the experiment (40% TPH at T 270d ), and that PAHs (16 PAH of the EPA list) were detected only at T 90d at 0.1 mg/L on the weekly survey of the discarded water. Furthermore, it is important to note that the n-C 18 /phytane ratio (Fig. S3 in the supplemental material) decreased progressively from T 90d (2.260.1 for BAL; 2.160.2 for NEREIS+BAL) to T 270d (1.060.1 for BAL; 0.760.4 for NER-EIS+BAL) with significant differences according to the Duncan test (p-value,0.05). This observation indicated that a biodegradation process was involved in alkane removal since linear alkanes (such as n-C 18 ) are more sensitive to microbial degradation than alkylated alkanes (such as phytane) [34].

Dynamics of Bacterial Communities
Bacterial community structures were monitored by T-RFLP based on 16S rRNA gene (DNA) and transcript (cDNA) analyses for all incubation conditions throughout the experiment (Fig. 3). Pearson's correlation analysis showed about 70% similarity between the replicates, revealing the spatial heterogeneity within and between microcosms, highlighting the importance of the sampling strategy to obtain representative samples (variability of the BRs is shown by nMDS analysis; Fig. S4A and B in the supplemental material). DNA and cDNA community structures were clearly separated by nMDS analysis (Fig. 3) but similar trends were observed, particularly across time. The bacterial DNA and cDNA community structures significantly shifted as function of time (model I PerMANOVA, Table S3 in the supplemental material), with a pronounced shift of both DNA and cDNA bacterial community structures from day 180 (Fig. 3). Despite the strong effect of time, the BAL addition and the NEREIS addition  Table S3 for a proper estimation of BAL and NEREIS effects). However, the significant three-level interaction (model I, Table S3 in the supplemental material) revealed complex shifts across time and conditions. Due to these complex interactions, the most similar communities across time and conditions were characterized by hierarchical clustering analysis (HCA). According to the similarity observed between biological replicates (BR), thresholds of 69.9% and 70.3% were applied to determine the cluster for 16S rRNA gene and transcript analyses respectively. The analysis of the total bacterial community (DNA analyses) revealed three significantly different clusters ( Fig. 4A; see ANOSIM test values in Table S4 in the supplemental material), with cluster I containing samples from the beginning of the experiment (T 2d ) irrespective of the experiment conditions, and clusters II and III containing oiled samples from T 90d and T 180d respectively. This observation showed that, irrespective of H. diversicolor addition, the presence of oil determined a different bacterial community structure (BCS). Thus, in oiled microcosms, BCS at T 90d and T 180d were not affected by the addition of H. diversicolor, despite the pronounced effect observed in NEREIS treatment (180 days). Nevertheless, in NEREIS and BAL microcosms, BCS modifications were observed from day 7 of incubation, whereas in NEREIS+BAL condition, BCS modification was observed from day 30. At the end of the experiment, specific BCS were observed for each treatment showing a H. diversicolor effect in the presence of oil (T 270d ). Bacterial community dynamics are also presented in nMDS by incubation condition and over time ( Fig. S4 and S5 in the supplemental material). In the BAL and NEREIS+BAL conditions, the most pronounced modification of the community structure was observed at the end of the experiment (T 270d ), testifying to the impact of oil on BCS. nMDS also shows the variability between microcosm replicates, represented by confidence ellipses that take into account the standard deviation of replicates. Although the predictive identification of OTUs has to be taken with caution, the in silico restriction of cloned sequences allowed clear identification at the phylum/class taxonomic level (see Fig. S6 in the supplemental material). At T 180d , BAL and NEREIS+BAL communities were mainly composed of OTUs belonging to Gammaproteobacteria representing about 45% of the profiles, explaining their similarity. The shift between T 180d and T 270d was characterized by a modification of the composition of OTUs related to Gammaproteobacteria, especially in the NEREIS+BAL community.
With the cDNA analyses ( Fig. 4B and Table S4 in the supplemental material), the HCA confirmed the BCS shifts observed by the DNA analyses further highlighting the treatment effects. Nevertheless, in contrast to the BCS shifts observed through DNA analyses, oiled bacterial communities (BAL and NEREIS+BAL) clustered together from day 7 (Fig. 4B, clusters V and VI) revealing a fast response of the bacterial community to oil addition. After 90 days, distinct BCS were observed except at T 270d , with un-oiled samples (CTRL and NEREIS) forming a cluster (cluster VII). Finally, at the end of the experiment (T 270d ), BCS were clearly distinct from T 2d , as also observed in nMDS represented by incubation conditions (Fig. S4B in the supplemental material). Furthermore, in silico digestion of cloned sequences by restriction enzymes (see Fig. S6 in the supplemental material) revealed different taxonomic patterns to those observed with the DNA analyses, particularly the OTUs related to Deltaproteobacteria that were found in high abundance, whereas they were not detected by the DNA analyses. At T 180d, their abundance in the BAL community (21.4%) compared to the NEREIS+BAL community (2.4%) explained the discrepancies observed between these communities. At T 270d , the NEREIS+BAL community exhibited a higher abundance of OTUs related to Gammaproteobacteria (56.3%) while the Alphaproteobacteria were more abundant in the BAL community (19.5%) than in the NEREIS+BAL community (2.3%).
The comparison of T-RFLP analyses of 16S rRNA genes and 16S rRNA transcripts (Fig. 5A) showed that the active population captured at the genomic level (RNA+ DNA+) in BAL and NEREIS+BAL conditions represented approximately 20% at T 270d . However, during the course of the experiment, different trends were observed according to the condition. The OTUs common to DNA and RNA analyses (RNA+ DNA+) increased significantly from day 180 in the BAL condition. In the NEREIS+BAL condition, it increased from day 7, with a concomitant decrease of specific OTUs observed through the DNA analyses (RNA2 DNA+). These observations suggested an early effect when H. diversicolor was added.
The OTUs common to both fractions (RNA+ DNA+) were significantly and negatively correlated with the n-alkanes concentrations (Fig. 5B). The same results were obtained with TPH and PAHs concentrations (data not shown). The fractions of OTUs present only in DNA-based analysis (RNA2 DNA+) and only in RNA-based analysis (RNA+ DNA2) were not significantly related with the hydrocarbon concentrations. These observations indicated that the bacterial active dominant fraction increased concomitantly with hydrocarbon biodegradation.

Bacterial Community Composition by 16S rRNA Transcript Analyses
In order to characterize the bacterial community compositions at T 270d , after hydrocarbon removal was observed, 16S cDNA libraries were constructed for each condition (see Fig. S7 to S9 in the supplemental material). The characteristics of the libraries are shown in Table 1. The coverage of the T 270d libraries from 51.4 to 76.4% indicated that the libraries were representative of the most abundant populations. The analysis of 333 complete sequences, representing 186 OTUs (considering a 97% similarity threshold), showed that 80.9% of OTUs were found specifically in a single condition, which supported the LIBSHUFF library's comparison, The comparison of the libraries' composition at the phylum level highlighted the differences between the bacterial communities (Fig. 6, Table 2). Sequences related to Proteobacteria and Bacteroidetes being the most abundant in all libraries except the NEREIS library (T 270d ), in which only Proteobacteria (98.6%) and Nitrospirae (1.4%) were detected. Among the Proteobacteria, the Gammaproteobacteria was the most abundant class found, especially in the NEREIS library (90.3% of sequences). Deltaproteobacteria, the second most abundant Proteobacteria, were found with similar abundances (20.7% and 22.5%) in the oiled NEREIS+BAL and BAL libraries, while they represented only 6.9% in the NEREIS library. Interestingly, in agreement with T-RFLP analyses, Alphaproteobacteria were found to be most abundant in the CTRL and BAL libraries (7.9% and 13.4% respectively) while they represented less than 4% in libraries from microcosms in which H. diversicolor was present. Betaproteobacteria were found only in the BAL library and one sequence related to Zetaproteobacteria was found in the CTRL library.
Focusing on Gammaproteobacteria and Deltaproteobacteria, the most abundant OTUs, it is important to note that most of the sequences of the four T 270d libraries were associated with microorganisms involved in the sulfur cycle ( Table 2, see Figs. S7 and S8 in the supplemental material). They included sulfur-oxidizing bacteria (SOB) such as purple sulfur bacteria within the Chromatiaceae and Ectothiorhodospiraceae families, sulfur-reducing Desulforomonadaceae and sulfate-reducing bacteria (Desulfobacteraceae and Desulfobulbaceae). The NEREIS library was characterized by 74.9% of sequences related to these families (about 40% in the other libraries). They were mainly associated with SOB (56.9%) as Sulfur-oxidizing bacterium ODIII6, Thiohalophilus thiocyanatoxydans and Olavius algarvensis sulfur-oxidizing endosymbiont. Chromatiaceae were found in lower abundance in the BAL and NEREIS+BAL libraries than in the CTRL and NEREIS libraries, indicating that they were affected by oil. In contrast, we observed that the sulfatereducing Desulfobacteraceae were abundant (15.9% and 13.5%) in oiled libraries. The abundance of Alteromonadacae was higher in microcosms to which H. diversicolor had been added.
The phylogenetic analysis also revealed the presence of sequences associated with bacteria known to be able to degrade hydrocarbons (see Fig. S10 in the supplemental material). They were more abundant in oiled microcosms accounting for 14.5% (6 OTUs) and 18.8% (8 OTUs) of sequences in the BAL and NEREIS+BAL libraries respectively. The sequences in the BAL library were related to Betaproteobacteria and sulfate-reducing strains of Desulfococcus and Desulfosarcina genera while the NEREIS+BAL library was characterized by the presence of sequences related to members of Gammaproteobacteria and sulfate-reducing strains of Desulfobacterium genus, known for their ability to degrade hydrocarbons (Fig. S10).

Discussion
In this study, we found that the structure and diversity of bacterial communities of marine sediment were affected by oil addition and the polychaete. Depending on the addition of the polychaetes Hediste (Nereis) diversicolor different bacterial communi-  ties were obtained without affecting the overall oil removal capacity. However, we did not find evidence to support the hypothesis that burrowing activity affected the bacterial community's oil degradation capacity. We also did not observe any degradation of the heaviest PAHs, suggesting the absence of photo-oxidation. Without excluding other abiotic processes, we assume that microbial biodegradation of hydrocarbon compounds was effective due to the n-C18/phytane ratio. The NEREIS+BAL bacterial community showed higher diversity than the BAL bacterial community reflecting different structures with distinct sequences associated with bacteria known for their ability to degrade hydrocarbons. A few studies have reported that burrowing organisms, through mucus excretion and oxygenation of sediments, influence the bacterial communities in the burrow and its surroundings as well [35][36][37]. Their presence had an impact on the bacterial degradation capacities, stimulating microbial degradation of organic matter [38] and organic contaminants [39] including crude oil compounds [7].
In the present study, we also found evidence that H. diversicolor's burrowing increased as previously reported [40] and that the burrowing was further enhanced in the presence of BAL 110 oil. Similar observations have been reported [5], but this was not the case during a laboratory experiment where H. diversicolor reworking was found to be reduced in the presence of Arabian Light crude oil [41]. In our study, the presence of oil would limit the penetration of oxygen, nutrients, and organic matter that are required for macrobenthic life [42,43]. Thus, the observed enhancement of the reworking activity could be explained by the fact that polychaetes would limit their stay in deeper sediment zones. Regarding the effect of H. diversicolor and its reworking activity on the microbial degradation capacities contradictory observations were reported. It was shown that pyrene degradation was not increased by the presence of H. diversicolor [6] while another report showed that the mineralisation of naphthalene was decreased in the presence of high density meiofauna [44]. These studies, performed with single molecules, did not focus on the composition of active bacterial communities. Thus, to the best of our knowledge, this work represents the first attempt to study the impact of petroleum, a complex mixture, on the composition of active bacterial communities in bioturbated sediments.
T-RFLP (both DNA and cDNA analyses) showed that BAL and NEREIS+BAL bacterial community dynamics were oil-dependent. A first community shift was observed after 90 days of incubation in both contaminated BAL and NEREIS+BAL conditions followed by a second community shift after 180 days of incubation, suggesting that a succession of bacterial communities occurred, which resulted in different bacterial patterns according to the treatment, both showing oil removal capacities. This succession of bacterial communities could be related to the different phases of hydrocarbon removal observed, the earlier degradation of n-alkanes followed by degradation of PAHs. This is supported by other studies showing that bacterial communities are first dominated by n-alkane-degrading bacteria and then by bacterial groups degrading polyaromatic hydrocarbons (PAHs), compounds more difficult to degrade (for a review see [45]).
In the NEREIS+BAL condition, the number of OTUs detected only in the DNA but not in the RNA analyses (RNA2 DNA+; representing inactive dominant fraction of the community) decreased suggesting that some OTUs became much more active to a point where they could be detected on the RNA-based analyses (RNA+ DNA+; representing the active dominant fraction of the community). The OTUs detected only in the RNA analyses (RNA+ DNA2; representing the active minor fraction of the community) were stable during the course of the experiment. The dominant fraction of the NEREIS+BAL community became active earlier than that of the BAL community, showing the role of polychaetes in the stimulation of the bacterial population. Furthermore, the significant correlation observed between the bacterial dominant active fractions and hydrocarbon concentrations may suggest the involvement of dominant active fraction of NEREIS+BAL and BAL communities in hydrocarbon biodegradation.
In oiled microcosms, the comparison of bacterial communities based on T-RFLP and 16S cDNA library analyses revealed that the bacterial communities showed different dynamics in microcosms to which the polychaetes was added, resulting in two distinct bacterial communities. The NEREIS+BAL active community was characterized by OTUs related to Gammaand Deltaproteobacteria while the most active OTUs in the BAL community were related to Alphaand Deltaproteobacteria. These results are consistent with recent Deepwater Horizon oil spill studies, showing that bacterial communities are predominated by members of the Gammaand Alphaproteobacteria in oil-polluted sediments [46] or sand beach [47]. Indeed, sequences retrieved in our study are related to sequences of hydrocarbon-degrading isolates obtained from sediments contaminated by Deepwater Horizon oil spill, especially associated with members of Cycloclasticus and Marinobacter genera (see Fig. S10 in the supplemental material). It is also important to notice that we obtained sequences related to sequences retrieved in oil contaminated sediments, some of them related to sequences retrieved in Prestige oil contaminated sediments [48] (see Fig. S10 in the supplemental material). These observations indicated that hydrocarbon degraders were induced in our sediment. In our study, the phylogenetic affiliation of 16S cDNA sequences indicated that almost all sequences were related to bacterial groups involved in the sulfur cycle; only 14% of total phylotypes were associated with bacteria known to be able to degrade hydrocarbons. Many studies have reported the dominance of bacteria able to oxidize hydrocarbons in oil-contaminated ecosystems [3,4,46,47,[49][50][51]. It is possible that our analysis missed the bacterial bloom related to alkane degraders occurring rapidly after oil contamination [45]. However, since the hydrocarbons were present until the end of the experiment (40% TPH at T 270d ), the low abundance of sequences associated with hydrocarbon-degrading-bacteria in our sediments could be explained by the presence of microorganisms capable of oxidizing hydrocarbons different from those already cultivated, as well as by the fact that the sediments were from a non-polluted environment. Only biogenic hydrocarbons were detected, harboring a bacterial community non-adapted to petroleum compounds. Several reports [52,53] have demonstrated that the impact of crude oil on the bacterial community and the oil degradation capacities are dependent on the environment's history of pollution events [52]. These studies show that oil degradation was faster when bacterial communities were adapted to oil compounds [53]. The slow degradation in our experimental sediments that were begun in the winter season, could be explained by their being maintained at low ambient temperatures. Indeed low temperatures affect the bio-availability of hydrocarbon compounds by increasing their viscosity and decreasing their volatility [54] as well as the microbial activities [55]. In the present study, sequences related to sulfur-oxidizing bacteria (SOB) were found in lower abundance in oiled libraries than in un-oiled libraries. Most of these sequences were related to the cultivable strains Thioalkalivibrio denitrificans and Thiohalophilus thiocyanoxidans [56,57]. The ability of these strains to grow with different electron acceptors such as oxygen and nitrogen oxides has been shown, demonstrating their adaptation capacity to different conditions and their role in the sulfur, nitrogen, and carbon biogeochemical cycles [58].
Sulfate-reducing bacteria (SRB) were found to be dominant in oiled microcosms (BAL and NEREIS+BAL). The simultaneous Unclassified ----- presence of both SOB and SRB communities suggests that the sulfur cycle was active in the microcosms. SRB are known to play an important role in anaerobic oil degradation [59][60][61]. Previous studies have found SRBs to be dominant in oil contaminated sediments [62,63], their activity being stimulated by the presence of hydrocarbons [64]. Interestingly in the present study, their abundance (based on library analyses) was not affected by the addition of H. diversicolor. This observation could be consistent with the increase in sulfate reduction rates observed in bioturbated sediments which other studies explained by the occurrence of reduced microniches through the increased burrowing activity [65]. However, other authors have also reported a decrease in sulfate-reduction rates in bioturbated sediment [12,66]. Clearly, with the results in the present study the debate continues and more research is warranted.
Our results show that the bacterial community structure in mudflat sediments shows evidence of modification after petroleum addition. The modification was dependent on the presence of the added burrowing polychaetes H. diversicolor in the microcosms. Contrary to our initial hypothesis, the overall oil removal capacity was not affected by the addition of polychaetes, although it did show evidence of affecting bacterial community structure. Indeed, the addition of burrowing organisms stimulated the bioturbation activity, resulting in the selection of a particular bacterial community characterized by the dominance of Gammaproteobacteria, as previously described in shrimp-inhabited sediment [36]. In the presence of petroleum, the OTUs related to potential hydrocarbon-degrading bacteria were distinct from those observed in oiled microcosms without added H. diversicolor. Thus, there was a change in the bacterial community to two distinct bacterial communities, both showing a similar overall oil removal capacity. Several hypotheses explaining how environmental perturbations may affect bacterial community composition have been drawn; they can become more resistant, resilient, or functionally redundant [67,68]. The mechanisms underlying functional redundancy are difficult to study because manipulating bacterial communities in the field is far from easy. In our study, the addition of burrowing organisms to sediments maintained near-environmental conditions allowed the manipulation of bacterial community structure and composition, opening the way for such a study. Furthermore, understanding the structuring effect of the reworking processes on bacterial communities will be useful for the management of bacterial metabolisms involved in bioremediation processes of coastal marine sediments.