Forty-two sediment samples were collected in different sites along the Mediterranean and the Atlantic French coasts (Figure 1, Supplementary Table S1). An a priori contamination criterion was determined according to the vicinity of potential input of hydrocarbons, like urbanized areas, industrial activities, industrial or recreational shipping activities. According to this criterion, contaminated and non-contaminated surface sediment samples (1–2 cm) were collected within each site. Non-contaminated samples were collected at least 2 km away from contaminated sites. Sites from Le Havre (samples MA1 and MA2), Lorient (GC1 and GC2), La Rochelle (LR1), Port-Vendres (PV1 to PV3), Thau lagoon (M5 and M6), Toulon (M3 and M4), Beirut (L1 and L2), Bizerte (BI1) and Ajaccio (M1 and M2) were subjected to industrial activities and industrial shipping (Figure 1). Sites from Banyuls-sur-Mer (BA1 et BA4) and Quiberon (GC1) were subjected to recreational shipping activities. Finally, GC6 and GC7 samples were collected in the same location before (August 1999) and after (July 2000) the Erika oil spill, which occurred in December 1999.

Samples were frozen immediately and transported at -18°C to the laboratory, where they were conserved at -80°C until further analysis. Salinity, temperature, particle size distribution (PSD, % < 63 μm), total organic carbon content (%TOC) and PAH concentrations were different among sites (Supplementary Table S1). PSD was measured using laser diffraction, according to Azoury et al. (2013) and %TOC was determined according to Tronczynski et al. (2005), on the decarbonated sediment samples by flash combustion (900–1000°C) with thermal conductivity detection (TCD). PSD was not determined for Bizerte samples, two samples from Lebanon and samples from La Rochelle. %TOC was not determined for two samples from Lebanon, one from the Gulf of Lion, three from Corsica and two samples from La Rochelle (Supplementary Table S1). PAH analyses were processed according to Stauffert et al. (2013) using Gas-Chromatography coupled to Mass-Spectrometry (GC-MS). As PAH concentrations were significantly correlated (rho > 0.75, p < 0.05), we were able to calculate PAH sums and compare them with several Sediment Quality Guidelines (SQG, Buchman, 2008): PEL (Probable Effect Low), ERL (Effect Range Low) and T50 (50% toxicity concentration). Two PAH diagnosis ratios [phenanthrene/anthracene (P/A) and fluoranthene/pyrene (F/P)] were considered in order to evaluate the source of contamination across sites. Fluoranthene and phenanthrene are found in fossil fuels and are therefore petrogenic PAHs, while anthracene and pyrene are produced during combustions and are pyrogenic PAHs. Values below 10 and above 1 for P/A and F/P, respectively, were shown to indicate a pyrogenic source of PAHs (Budzinski et al., 1997).

Total genomic DNA was extracted from 250 mg of sediment samples using the PowerSoil DNA isolation kit (MO-BIO Laboratories Inc.). PCR amplifications were performed using primers 27f (5′-AGRGTTTGATCMTGGCTCAG-3′) and 519r (5′-GTNTTACNGCGGCKGCTG-3′) targeting the V1–V3 region of the bacterial 16S rRNA gene (Weisburg et al., 1991), primers 344f (3′-ACGGGGYGCAGCAGGCGCGA-5′) and 915r (3′-GTGCTCCCCCGCCAATTCCT-5′) targeting the V3–V4 region of the archaeal 16S rRNA gene (Raskin et al., 1994) and primers 1560f (3′-TGGTGCATGGCCGTTCTTAGT-5′) and 2035r (3′-CATCTAAGGGCATCACAGACC-5′) targeting the VR4 region of the eukaryotic 18S rRNA gene (Hardy et al., 2010). Pyrosequencing was performed at Mr. DNA – Molecular Research (Shallowater, TX, USA), using a Roche 454 GS-FLX Titanium instrument (Roche, NJ, USA). Pyrosequencing resulted in 693687, 290816, and 110044 raw reads for bacterial, archaeal and eukaryotic datasets, respectively, (Supplementary Table S2). Raw reads from the three datasets were independently quality-based, trimmed, and aligned on the May 2013 Greengenes reference alignment for prokaryotes (McDonald et al., 2012) and SILVA 119 database for eukaryotes (Quast et al., 2013). Chimera were subsequently detected and eliminated using UCHIME (Edgar et al., 2011). All statistical analyses were performed on a random subsample of 3238, 829 and 804 sequences for bacterial, archaeal and eukaryotic datasets, respectively, corresponding to the smaller number of sequences per sample in the datasets, after trimming and quality processing. OTUs were clustered at a 90% identity threshold following the Mothur 454 Standard Operating Procedure (Schloss et al., 2009) with few modifications. The 90% threshold for OTUs is close to the taxonomic level of the family (Konstantinidis and Tiedje, 2007). This threshold was used to get sufficiently abundant OTUs for correlation analysis (i.e., enough occurrence and abundance to get statistically robust Spearman correlations). Another advantage of using this threshold is to circumvent potential taxonomic misclassifications due to sequencing anomalies (Barberán et al., 2012; Fillol et al., 2016). This resulted in 6034 bacterial OTUs, 680 archaeal OTUs and 774 eukaryotic OTUs. Taxonomic classification of each OTU was performed using the May 2013 Greengenes reference alignment for prokaryotes (McDonald et al., 2012) and the SILVA 119 database for eukaryotes (Quast et al., 2013). A representative sequence (i.e., the closest sequence of all other sequences) for each OTU was selected and a phylogenetic tree was reconstructed for each domain of life using the ARB software (Ludwig et al., 2004). Bacterial, archaeal and eukaryotic base frequency ARB filters were applied to exclude highly variable positions before sequences were added using the ARB parsimony insertion tool to the May 2013 Greengenes and SILVA 119 backbone trees, for prokaryotes and eukaryotes, respectively. Resulting phylogenetic trees were exported and used to run UniFrac weighted matrix calculation for each domain of life (Lozupone et al., 2006). Weighted UniFrac distance matrices were further used in beta-diversity analysis. We used PERMANOVA based on 1000 permutations (McArdle and Anderson, 2001) with function Adonis of the vegan package in R (Oksanen et al., 2013) to assess the source of variation in the weighted UniFrac matrices. Each dataset encompass the 42 samples of this study. The complete dataset was deposited in the NCBI Sequence Read Archive (SRA) database under study Accession no SRP063723.

Taking independently the Atlantic and Mediterranean regions (see Results), contaminated and non-contaminated site groups were defined by calculating and plotting Euclidean distances based on 11 PAHs as a dendrogram based on hierarchical ascendant classification (HAC, Supplementary Figure S1). Mean total PAH concentrations between contaminated and non-contaminated sites from each geographic region (i.e., Atlantic and Mediterranean regions) were compared using the Wilcoxon–Mann–Whitney U test.

Collinearities in independent environmental variables were tested before running permutational Manova (PERMANOVA), species abundance distribution (SAD) and network analysis. Variables with collinearity up to 0.75 were grouped together, and proxies were used as explanatory variables. Fluroranthene was the representative predictor for 10 other PAH (i.e., Fluorene, Pyrene, Anthracene, Benz[a]anthracene, Benzo[a]pyrene, Benzo[g,h,i]Perylene, Chrysene, Dibenz[a,h]Anthracene, Indeno[1,2,3,c-d]Pyrene and Phenanthrene) and salinity for latitude and temperature (Supplementary Table S3). Salinity and fluoranthene concentration were log transformed, and F/P ratio was squared transformed to improve the linearity and homoscedasticity of residuals.

The OTU abundance table was used to examine the SAD patterns of each OTU and partition the SAD into core and satellite microbial OTUs. For this purpose, the index of dispersion for each OTU was calculated as the ratio of the variance to the mean abundance (VMR) multiplied by the occurrence. This index was used to model whether lineages follow a Poisson distribution (i.e., stochastic distribution), falling between the 2.5 and 97.5% confidence interval of the χ² distribution (Krebs, 1999).

Associations between bacterial, archaeal and eukaryotic OTUs were inferred from an undirected co-occurrence network. Pairwise score between 90% core OTUs were computed using Spearman’s rank correlations. Only co-occurrences corresponding to correlations with a coefficient (rho) > 0.6 and a statistical significance (P-value) < 0.001 were considered for further analysis. Non-random co-occurrence patterns were tested with the checkerboard score (C-score) under a null model preserving site frequencies (Stone and Roberts, 1990; Gotelli and McCabe, 2002). A C-score calculated for each pair of microbial OTUs was compared with the C-score computed for 5000 randomly assembled null matrices. In order to avoid biases affecting raw C-score values (i.e., OTU number, abundance…), the standardized effect size (SES) was calculated (Sridhar et al., 2012). Because the C-score is an inverse indicator of the frequency of co-occurrence, positive SES values indicate less co-occurrence than expected by chance (i.e., predominance of segregation within communities) and vice versa for negative values (i.e., predominance of facilitation). If co-occurrences were not different from what was expected by chance, values of SES should fall between -2 and 2.

The network was visualized with the Gephi software (Bastian et al., 2009). Nodes represented prokaryotic and eukaryotic OTUs at 90% identity and edges represented the significant correlations between them. Network characterization for comparisons between contaminated and uncontaminated sediments was performed using a set of overall network topological indices (i.e., average node connectivity, average path length, average clustering coefficient and modularity; Newman, 2003). After module detection, each module was represented by its singular value decomposition (SVD) of abundance profile by using module eigengene analysis (Langfelder and Horvath, 2007). Module response to environmental and contaminant parameters was assessed by correlations between module eigengenes and environmental factors.

All analyses were run using the R packages vegan (Oksanen et al., 2013), igraph (Csardi and Nepusz, 2006) and WGCNA packages (Langfelder and Horvath, 2008).

Particle size distribution, %TOC, salinity, and temperature averages of the water column were significantly greater in the Mediterranean sites than in the Atlantic sites (Supplementary Figure S2). The salinity, temperature, and latitude gradients were highly correlated (Spearman rho = -0.82, p-value = 0.001, Supplementary Table S3). Because of the important environmental heterogeneity among the sediments from the two regions, we analyzed the PAHs contamination gradient in both regions individually. PAHs concentrations were strongly correlated to each other (rho > 0.75, p < 0.05, Supplementary Table S3). In order to compare the PAH concentrations to the Sediment Quality Guidelines thresholds (SQG; Buchman, 2008), which are available for 9 PAHs among the 11 PAHs of our dataset, the sums of these 9 PAHs concentrations (Σ PAH₉) and of the corresponding SQGs were calculated (i.e., for fluoranthene, pyrene, anthracene, benz[a]anthracene, benzo[a]pyrene, chrysene, dibenz[a,h]anthracene, fluorene, phenanthrene; Figure 2). The SQG used are the following: PEL, ERL, which are low trigger concentrations and T50 (50% toxicity concentration) values.

The HAC analysis confirmed the a priori classification (i.e., pre-sampling classification) for most of the samples, and particularly for Mediterranean samples (Supplementary Figure S1). Figure 2 shows the differences of mean Σ PAH₉ concentrations between Mediterranean and Atlantic regions. Contaminated samples determined using the HAC analysis exhibited significantly higher mean PAH concentrations in both Mediterranean and Atlantic regions (Figure 2). However, the difference between contaminated and non-contaminated samples was much greater in the Mediterranean sites. When comparing mean Σ PAH₉ concentrations, PEL, ERL, and T50 separated Mediterranean contaminated sites from non-contaminated sites, while only PEL could be a threshold between contaminated and non-contaminated Atlantic sites (Figure 2). SQGs were originally defined regarding macrofauna ecotoxicological assays (Buchman, 2008), and little is known about their reliability for the determination of PAHs toxicity on microorganisms (Sun et al., 2012). However, bacterial community structure from chronically contaminated estuarine sediment was previously found to respond to low trigger concentrations (Sun et al., 2012). These structure modifications were mainly related to change in core OTU abundances (Sun et al., 2013).

We used diagnosis ratios of PAH concentrations (i.e., fluoranthene to pyrene (F/P) and phenanthrene to anthracene (P/A)) to get an indication of the dominant source of contamination by PAHs (i.e., pyrogenic versus petrogenic; Budzinski et al., 1997; Neff et al., 2005). Most of the samples, from both the Mediterranean and the Atlantic regions, exhibited P/A values lower than 10 and F/P values above 1, indicating a pyrogenic contamination (Supplementary Figure S3, Supplementary Table S1). In contrast, none of the sites exhibited clear petrogenic contamination according to these ratios. PAHs derived from combustion in coastal sediments occurs through atmospheric deposition (Lipiatou et al., 1997; Golomb et al., 2001), which constitutes a chronic PAH input since the Industrial Revolution (Lima et al., 2003; Azoury et al., 2013). Considering the diagnosis ratios, we could therefore infer that the samples analyzed here were thus facing mainly chronic contamination by pyrolytic PAHs.

Traditionally in macro-ecology, the SAD in spatially or temporally separated sites can be divided into two groups: the core species/taxa, which are widely distributed and generally abundant, and the satellite species/taxa, which mostly occur in low abundance and at a reduced number of sites (Hanski, 1982; Magurran and Henderson, 2003; Ulrich and Zalewski, 2006). Previous studies showed that this division of the SAD was not artifactual and demonstrated that environmental factors shaped the relative abundance of the core species/taxa, whereas random dispersal mostly explained satellite species/taxa community assembly (Magurran and Henderson, 2003; Ulrich and Zalewski, 2006; Magurran, 2007). The concept has been reinforced in subsequent studies, which suggested that relevant aspects of the SAD could be overlooked without such a partition (Ulrich and Ollik, 2004; Gray et al., 2005; Dolan et al., 2009). Despite a growing number of surveys interested by rare or core communities in microbial ecology (Galand et al., 2009; Pedrós-Alió, 2012; Hugoni et al., 2013; Cariveau et al., 2014), the mechanistic framework developed by Magurran and Henderson (2003) to partition the SAD has been very rarely used (Fillol et al., 2016).

For this purpose, 7488 microbial OTUs (6034 bacteria, 680 archaea, and 774 eukaryotes OTUs at 90% cutoff) from the 42 sediments sampled in this study were analyzed in an abundance-vs.-occurrence plot (Figure 3A). In agreement with one of the most robust trend in macro-ecology (Gaston et al., 2000), we found a significant and positive relationship between the mean abundance and the occurrence. However, due to a high number of OTUs, the typical discontinuity that separates core from satellite OTUs was not appreciable (Figure 3A). In order to partition the SAD, we compared the index of dispersion for each OTU to a model assuming a stochastic distribution (Poisson model) falling between the 2.5% and 97.5% confidence limit of the χ² distribution (Krebs, 1999). Overall, 859 OTUs (11% of the initial number of OTUs) were found to form the core microbial community of coastal sediments (Figure 3B). They represented 80% of the reads. In contrast, the remaining 89% of the OTUs fell below the 2.5% confidence limit line, indicating random distribution through sediment samples (Figure 3B). The ranking abundance plot confirmed that this method represent a tractable and reproducible approach to partition the SAD of microbial communities (Figure 3B). It also highlights the fact that the definition of core or satellite OTUs should not be only based on abundance or occurrence/permanence but on both, as suggested by Hubbell (2001).

Core and satellite community compositions were found to differ significantly for the three domains of life (Figure 4). A high percentage of satellite OTUs remained unclassified at the phylum level, particularly for eukaryotes with more than 50% of the OTU that could not be affiliated to any taxa. For bacteria, the core community was dominated by OTUs belonging to the Gammaproteobacteria (16% of the total abundance) and Deltaproteobacteria (9%), which are known to be dominant clades in marine sediments (both contaminated and pristine) and to encompass a large number of hydrocarbon-degrader strains (Paissé et al., 2008; Zhang et al., 2008; Rojo, 2009; Kostka et al., 2011; Acosta-González et al., 2013; Sun et al., 2013). The dominance of Deltaproteobacteria OTUs within the core community is in agreement with results found along a pollution gradient in estuarine sediments (Sun et al., 2013). This class includes most of the sulfate-reducing genera, which play a crucial role in the anaerobic degradation of organic matter (Leloup et al., 2009), but also in hydrocarbon degradation (Suárez-Suárez et al., 2011; Stauffert et al., 2014a). Archaeal core community was the most divergent from the satellite community, with domination by OTUs belonging to the Thaumarchaeota, Bathyarchaeota and Thermoplasmata (Figure 4). The latter groups, which represented 22 and 14% of total abundance, respectively, were found to be part of the core microbiota in a meta-analysis regrouping more than 200 sediment studies (Fillol et al., 2016). Both archaeal groups seemed to form specific consortia in sediment ecosystems, suggesting a potentially relevant trophic connection related to the degradation of aromatic compounds and detrital proteins (Lloyd et al., 2013; Meng et al., 2014; Fillol et al., 2016). In contrast, distribution and potential functions of the Parvarchaeota (previously known as Haloarchaea), which represented 64% of the total abundance in the archaeal satellite community, are currently unknown.

For eukaryote a large proportion of the OTUs of the core and satellite communities remained unclassified (i.e., 25% and 57%, respectively), highlighting our ignorance about eukaryotic communities inhabiting coastal sediments. With a growing number of environmental surveys investigating the structure and diversity of eukaryotic communities using the 18S rDNA gene and NGS (next generation sequencing), this result underpinned the urge for the development of eukaryotic taxonomy databases similar to those for prokaryotes. The remaining OTUs were dominated both in the core and satellite communities by members of the Holozoa kingdom and the Alveolata superphylum with 45 and 18% of the total abundance, respectively.

We conducted PERMANOVA analyses in order to identify whether the core or the satellite community would show the highest sensitivity to changes in environmental and contaminant variables (Table 1). As expected, we observed a strong filtering effect on the most abundant and frequent OTUs forming the core community, with salinity (associated with latitude and temperature) representing the main environmental parameters driving the community composition. %TOC and particle size distribution (% > 63 μm) together explained an additional 10% of the variance in community composition, while fluoranthene explained only 2% of this variance. The response to the measured predictors was weak for the OTUs forming the satellite community, with only 2% of the overall variance explained by salinity. This result may indicate a minimal influence of satellite OTUs on biogeochemical cycles as previously postulated (Gobet et al., 2012), and contrary to core OTUs that are thought to contribute actively to the biomass and the nutrient cycling within the community (Campbell et al., 2011; Gaidos et al., 2011; Pedrós-Alió, 2012). However, Gammaproteobacteria and Deltaproteobacteria, which encompass many species able to degrade hydrocarbons (Prince et al., 2010), were found in similar proportions in the rare and core communities (Figure 4), suggesting the potential role of the rare benthic biosphere in hydrocarbon degradation. Bacterial hydrocarbon-degraders, such as Cycloclasticus, were depicted as members of the rare community in non-polluted environment, and became dominant in the presence of hydrocarbon contamination (Teira et al., 2007; Sauret et al., 2014). This fluctuating state of specialized taxa from the rare community could facilitate the adaptation of the whole community to chronic hydrocarbon contamination. Overall, despite the predominant response of the core community to environmental changes, we cannot exclude the potentially central role of members of the rare community in chronically hydrocarbon-polluted sediments.

As illustrated by recent studies, microbial co-occurrence patterns can help to unveil ecologically meaningful interactions between taxa (Horner-Devine et al., 2007; Steele et al., 2011; Fillol et al., 2016) and network properties can be used to study community stability in altered conditions (Zhou et al., 2011; Deng et al., 2012; Faust and Raes, 2012; Sun et al., 2013). To get a first insight in the potential microbial consortia affected by PAH contaminations and the potential effect of these contaminants on the core community integrity, we constructed a co-occurrence network based on strong and significant Spearman correlations. The dataset containing core OTUs was filtered for optimizing network sensitivity and specificity. The final network included 131 nodes and 566 edges (Figure 5). We used the SES of the C-score metric to assess the existence of non-random co-occurrence patterns in the network (Gotelli and McCabe, 2002; Horner-Devine et al., 2007; López et al., 2013). We observed a SES value of 2.56 (P < 0.001, C-score_norm = 0.88) indicating non-random network structure. The network was largely dominated by bacterial OTUs, which represented 76% of the nodes (Figure 5). Eukarya and Archaea represented 15 and 9% of the nodes, respectively. Correlations between domains were common within the network although intra-domain connections tended to be more numerous. These co-occurrences might reflect OTUs performing similarly or harboring complementary functions, but they may also be caused by shared preferred environmental conditions (Steele et al., 2011; Deng et al., 2012). Network approaches showed that co-occurring OTUs are often organized into groups, or modules, of functional significance (Chaffron et al., 2010; Barberán et al., 2012; Vick-Majors et al., 2014). Modularity analysis revealed a complex structure composed of 11 modules (Figure 6A). Modules are groups of highly connected OTUs within the group but with very few connections outside the group. Figure 6B shows the response of modules to environment changes, as assessed by means of correlations between modules and environmental parameters. Five of the 11 modules detected were significantly correlated with the environmental parameters measured in this study (Figure 6B). Modules 1, 2, and 6 were correlated to salinity (and associated latitude and temperature), confirming the strong effect of these variables on microbial consortia. In agreement with previous works (Chaffron et al., 2010; Freilich et al., 2010; Faust and Raes, 2012), modules could be considered as ecological and/or functional niches as suggested by modules 1, 2 and 6, which represented low salinity/high latitude and high salinity/low latitude sub-networks, respectively. Module 5 was the only microbial consortium responding (negatively) to fluoranthene and associated PAH concentrations, indicating that PAH contaminations might have a negative impact on the members of this module. This module was composed of four eukaryotic and two bacterial OTUs connected by an unusual number of negative correlations between OTUs belonging to Actinobacteria and OTUs belonging to the Alveolata taxa (Figure 7), suggesting predator – prey relationships (Faust and Raes, 2012). This assumption is in agreement with the fact that Alveolates, and more particularly Ciliates, feed on a variety of prey but that many ciliate groups are bacterivores (Meng et al., 2012).

Comparing network properties and overall topology can be used to study community stability/vulnerability in altered conditions. In order to assess the influence of hydrocarbon contaminants on network properties, we split the core OTU abundance table based on the classification provided by the HAC (Supplementary Figure S2) and constructed separated microbial networks for non-contaminated and contaminated sediment communities (Table 2). Non-contaminated and contaminated networks were generated for the whole dataset and for each region (i.e., Mediterranean and Atlantic) in order to limit the influence of environmental variables (i.e., salinity, temperature, latitude, PSD and %TOC, which were significantly different between Mediterranean and Atlantic regions). Network comparisons are commonly achieved by calculating the average for each node of four network indices (i) connectivity or degree distribution, which represent the number of edges of a node toward other nodes; (ii) path length, which is the shortest path between two nodes; (iii) the clustering coefficient, which describes how well a node is connected to its neighbors; (iv) modularity, a measure of the structuration of the network into modules (Deng et al., 2012). Table 2 showed that the three pairs of networks calculated here presented the typical topology for microbial networks (Chaffron et al., 2010; Steele et al., 2011; Barberán et al., 2012; Deng et al., 2012): scale-free (i.e., node connectivity distribution not different from a power law model, data not shown), small-world (clustering coefficient ranging from 0.56 to 0.89 and average path from 2.22 to 2.75) and modular (modularity ranging from 0.72 to 0.86). Overall and except for the Atlantic samples, significant differences between contaminated and non-contaminated networks were detected in terms of average connectivity, average shortest path and average clustering coefficient (Table 2). These results indicated that hydrocarbon contaminations might have a significant and negative effect on the interactions among microbial taxa in coastal sediments. Indeed, the higher values for average shortest path in contaminated networks may imply lower efficiency and speed in the transmission of information, energy or material in the system, potentially altering the response of the microbial community to environmental perturbations (Zhou et al., 2010; Deng et al., 2012). In addition, the lower connectivity and lower average clustering coefficient in contaminated networks may make the networks less resistant to change (Albert et al., 2000; Montoya et al., 2006). This may also indicate a potential perturbation of the network and the loss of important community members. A lower clustering coefficient is indeed an indication of the loss of highly connected nodes (i.e., hubs), which have been related to the concept of keystone species (Paine, 1969; Steele et al., 2011). Such keystone species loss could cause a greater fragility of the network, which ultimately could lead to its fragmentation and to secondary extinctions of species relying on the hubs (Solé and Montoya, 2001).

Altogether, these results suggested that hydrocarbon contaminations might have altered the microbial network in coastal sediments, rendering the microbial communities inhabiting these habitats more vulnerable to environmental perturbations.