Benthic Prokaryotic Community Response to Polycyclic Aromatic Hydrocarbon Chronic Exposure: Importance of Emission Sources in Mediterranean Ports

The present study aimed to evaluate the impact of polycyclic aromatic hydrocarbons (PAHs) produced by multiple emission sources on prokaryotic communities in sediments chronically affected by anthropogenic pressures. In this context, surface sediments were investigated in three Mediterranean touristic ports over three sampling periods and in different port sectors. The levels of 16 priority PAHs varied over three orders of magnitude (25–49,000 ng g − 1 ) covering the range of concentrations previously reported for Mediterranean harbors. Pyrogenic processes were found to be the dominant emission source of PAHs, with considerable differences among ports. The prokaryotic communities were identiﬁed by using the terminal restriction fragment length polymorphism, targeting the 16S rRNA gene for Bacteria and Archaea as well as the dsrAB gene for sulphate-reducing bacteria (SRB). The structure of the three benthic prokaryotic communities varied consistently among the ports. The structure of Bacteria and Archaea exhibited strong spatiotemporal variations that did not allow us to speciﬁcally link the observed differences in community structures with PAH sources. On the contrary, our study provided, for the ﬁrst time, evidence that the PAH emission sources play a role in structuring benthic communities of SRB. Our ﬁndings indicate that the SRB community can be used as a valuable candidate biotic descriptor for bioremediation monitoring in heavily impacted port sediments.


INTRODUCTION
Pollution by polycyclic aromatic hydrocarbons (PAHs) remains one of the major threats to marine ecosystems and human health due to their toxic, carcinogenic and mutagenic properties (Ghosal et al., 2016). The ubiquitous distribution of PAHs in coastal environments has been mainly attributed to anthropogenic activities (Crain et al., 2009). Because of their hydrophobic nature, PAHs are quickly adsorbed onto suspended particulate material and settle on benthos, thus sediments represent a sink for these pollutants. Thereafter, their persistence in marine sediments is controlled by the degradation activities of benthic prokaryotic communities, which play a key role in PAH removal (Cravo-Laureau and Duran, 2014).
The chemical class of PAHs comprises numerous compounds with two or more fused aromatic rings, characterized by diverse physical, chemical and biological properties. In general, their recalcitrance toward biodegradation as well as their genotoxicity rise with the number of aromatic rings, while volatility and acute toxicity tend to increase with decreasing molecular weight (Ghosal et al., 2016). Moreover, the relative content of PAHs can vary considerably depending on their origin, namely petrogenic and pyrogenic emission sources (Baumard et al., 1998).
Petrogenic PAHs are introduced into the marine environments through accidental spills and leakages of crude oil and its products and are characterized by the predominance of low molecular weight compounds (LPAHs, 2-3 rings) (Neff et al., 2005), while pyrogenic PAHs are formed by the incomplete combustion of organic matter (e.g., fossil fuels, coal, biomass). High temperature processes (e.g., vehicle engines) produce high molecular weight PAHs (HPAHs, >3 rings), which enter the marine environment from different routes, such as atmospheric deposition, riverine inflows or discharge of urban run-offs (Yunker et al., 2002). The process by which PAHs are formed is crucial in determining not only the composition of the contaminant mixtures but also the bioavailability of PAHs in sediments (Akkanen et al., 2012;Duran and Cravo-Laureau, 2016). Petrogenic PAHs are usually more bioavailable than pyrogenic PAHs, which are known to bound more tightly onto sediment particles (Hylland, 2006). Furthermore, the bioavailability of pyrogenic PAHs produced by different processes differs according to the burned materials and the combustion conditions (Hylland, 2006;Akkanen et al., 2012). In highly urbanized and industrialized areas, pyrogenic PAHs are usually produced by multiple emission sources giving rise to complex contaminant mixtures. In this context, harbors and marinas are hot spots of contamination dominated by PAHs produced by multiple combustion processes, with usually higher levels of PAHs than the adjacent coastal zones (Baumard et al., 1998;Yunker et al., 2002;Merhaby et al., 2015;Schintu et al., 2015). However, a coexistence of petroleum and pyrogenic PAHs has been observed in multisectoral harbors where industrial and touristic activities coexist (McCready et al., 2000;Mostafa et al., 2003;De Luca et al., 2004Sprovieri et al., 2007).
The effects of petrogenic PAHs on benthic bacterial communities have been extensively investigated during artificial oil spills (Suárez-Suárez et al., 2011), accidental oil discharges (Acosta-González et al., 2013;Kimes et al., 2013), and natural oil seeps (LaMontagne et al., 2004;Orcutt et al., 2010). The effect of oil spills is very consistent across studies (Kimes et al., 2014;Acosta-González et al., 2015), and generally characterized by negative impacts on benthic bacterial communities (LaMontagne et al., 2004;Orcutt et al., 2010;Suárez-Suárez et al., 2011) and archaeal assemblages (Miralles et al., 2010). Petrogenic contaminations also enhance sulphate reduction and enrich sulphate reducing bacteria (SRB) able to anaerobically degrade hydrocarbons (Orcutt et al., 2010;Suárez-Suárez et al., 2011;Acosta-González et al., 2013;Kimes et al., 2013). In contrast, the effects on microbial communities of chronic exposure to complex mixtures of PAHs, produced by multiple emission sources, is still not well understood. Several studies have dealt with the impacts on prokaryotic communities of PAHs in sediments chronically affected by anthropogenic activities, such as harbors (Slater et al., 2008;Zhang et al., 2008a), anthropized estuaries (Sun et al., 2012(Sun et al., , 2013, industrial sites and densely populated coastal areas (Ben Said et al., 2010;Korlevic et al., 2015;Quero et al., 2015;Jeanbille et al., 2016a,b;Misson et al., 2016). Up to date, only contrasting results have been obtained, which do not allow for any general pattern to be depicted (Acosta-González and Marqués, 2016). In some cases, PAH levels have been related with no changes or minor shifts in bacterial assemblages (Zhang et al., 2008a;Iannelli et al., 2012;Jeanbille et al., 2016b), while other studies have highlighted variations in bacterial communities (Ben Said et al., 2010;Sun et al., 2012Sun et al., , 2013, reduction in microbial network connectivity (Jeanbille et al., 2016a), or an increase in specific populations of SRB (Quero et al., 2015). Usually, these studies have been conducted by comparing severely polluted to more natural (or less contaminated) sediments. However, only few studies have marginally taken into account the emission sources of PAHs entering the sediments (Korlevic et al., 2015;Jeanbille et al., 2016a,b).
The main aim of present study was to evaluate the impact of complex PAH mixtures, produced by multiple emission sources, on the benthic prokaryotic communities in sediments chronically affected by anthropogenic pressures. For this reason, in the framework of the ENPI CBCMED project MAPMED, sediments were collected from three touristic ports located along the Mediterranean Sea and characterized for emission sources and levels of contaminant PAHs. The benthic communities of Bacteria, Archaea and SRB were separately targeted because of their different roles in controlling and determining the fate of PAHs in marine sediments. Figure 1 shows the maps of the touristic harbors selected as case study sites (MAPMED Consortium, 2013). The Port of Cagliari is a large port (2.07 km 2 ) and represents the main harbor of the island of Sardinia (Italy, North-West Mediterranean) and one of the key nodes for the trans-shipment activities in western Mediterranean. The harbor is placed in front of the city center (about 150,000 citizens). Sampling was limited to the Historic Port of Cagliari where tourism-related activities are operated (Figure 1, C). The Port of El Kantaoui (Tunisia, South-West Mediterranean) is a relatively recently built, small, artificial marina (0.04 km 2 ), located on the eastern Tunisian coast within an international touristic center and surrounded by a small permanent population (Figure 1, E). The Port of Heraklion, an intermediate size port (0.87 km 2 ), is the main harbor of the island of Crete (Greece, Eastern Mediterranean) and one of the most important and active ports in the Eastern Mediterranean supporting intense touristic and transport uses. The port is located in close proximity to the city center (about 160,000 citizens) (Figure 1, H). The ports of Cagliari and Heraklion are both characterized by a significant maritime traffic including large passenger ferries, cruise and cargo vessels, while El Kantaoui port offers moorings for smaller fishing boats, luxury yachts and boats for sporting activities (MAPMED Consortium, 2013). The three ports are subject to a variety of potential pollution sources ranging from boats/ships docking and traffic, accidental or operational discharges of oily wastewaters, discharges of sewage effluents, spills from fueling stations, wastes from shipyard activities, and urban run-offs (MAPMED Consortium, 2013). In addition, the Cagliari port also receives freshwater drainage from the surrounding lagoons through an artificial channel. Within each port, three to five sampling stations were selected in order to represent the discrete sectors dominated by different port activities and to achieve good spatial coverage of the whole port area (Table 1 and Figure 1).

Study Sites and Sampling
Three sampling campaigns were carried out in all ports in February, May-June, and September (2012). During the whole sampling year, no periodic maintenance dredging operations were conducted in the three ports. The collected samples were labeled as follows: the first digit of the sample identifier refers to the time period of the sampling campaign (1: February; 2: May-June; 3: September), the second character specifies the port (C: Cagliari; E: El Kantaoui; H: Heraklion), while the third digit identifies the sampling station within each port sector ( Table 1). For analyses of hydrocarbons and microbiological variables, surface sediments (0-1 cm) were collected according to Chatzinikolaou et al. (2018). The surface layer was used to determine the distribution of contaminants since it provides information on the most recently deposited sediment materials (Simpson et al., 2005). The environmental status of the three investigated ports has been recently assessed by Chatzinikolaou et al. (2018), based on physico-chemical parameters measured on samples collected during the coordinated sampling campaigns. In the present work, the following parameters were included (Supplementary Table  S1): salinity measured in surface seawater as well as temperature, redox potential, silt-clay ratio, organic carbon, chlorophyll-a and phaeopigments, levels of aliphatic hydrocarbons (AHs) and unresolved complex mixture (UCM) measured in superficial sediments. The concentrations of 31 individual AHs in superficial sediments (Supplementary Table S2) were determined according to Mandalakis et al. (2014), and the total level of AHs ( AH) was calculated as the sum of individual compounds in the range C10-C40. Total phytopigment concentration was calculated as the sum of chlorophyll-a and phaeopigment concentrations (Danovaro et al., 1999). The ratio of phaeopigments to the sum of chlorophyll-a and phaeopigments (PAP ratio) was determined as an indicator of the "freshness" of the organic material deposited on the sediment (Boon et al., 1998).

Chemical Analysis of PAHs
The surface sediment samples were stored at −20 • C immediately after collection. For determination of PAHs, sediment samples were homogenized and spiked with the multi-standard mixture PAH-Mix 31 (Dr. Ehrenstorfer GmbH, Augsburg, Germany) of five deuterated PAHs (Acenaphthene-D10, Chrysene-D12, Naphthalene-D8, Perylene-D12, Phenanthrene-D10), as surrogate standard. The samples were then subjected to solvent extraction by sonication, sulfur removal by copper addition, and purification of extracts using silica gel column chromatography. The analysis of the 16 priority PAHs issued by the United States Environmental Protection Agency (EPA) was subsequently performed using a gas chromatography-mass spectrometry as described by Mandalakis et al. (2014).

Analysis of Prokaryotic Communities by Terminal Restriction Fragment Length Polymorphism
The sediment samples were stabilized in the field, immediately after collection, by adding an equal volume of RNAlater R reagent (1:1 v/v). The three core samples collected from each station were pooled in order to mitigate the effects of microscale heterogeneity in the sediment habitat. Genomic DNA was extracted from sediment samples (250 mg wet weight) using the NucleoSpin R Soil kit (Macherey-Nagel), following the manufacturer protocol. The DNA concentration was determined on agarose gel using a DNA quantitation standard.
Terminal restriction fragment length polymorphism (T-RFLP) of 16S rRNA genes was used as phylogenetic marker for addressing the bacterial and archaeal communities while the structure of SRB was specifically analyzed by using the gene dsrAB as a functional marker for dissimilatory sulfite reductase (Pérez-Jiménez and Kerkhof, 2005). The bacteriaspecific primer pair 27f-FAM/1492r (Moeseneder et al., 1999) and the archaea-specific primer pair 21f-NED/958r (De Long, 1992) were used for T-RFLP of the 16S rRNA gene of bacterial (16SBact-TRFLP) and archaeal (16SArch-TRFLP) domains, respectively. The amplification reactions were carried out as previously described by De Long (1992). For T-RFLP of the dsrAB gene (dsr-TRFLP), the primer pair DSR1F-FAM/DSR4R was used to amplify about 1.9 kb of the dsrAB gene (Pérez-Jiménez and Kerkhof, 2005).
In order to minimize stochastic PCR bias, the PCR products from three replicate reactions with each primer pair were combined and purified from the agarose gel using the QIAquick Gel extraction kit (Qiagen). The 16S rRNA gene amplicons were digested with the restriction enzyme AluI, which was selected on the basis of an "in silico" digestion using the ERPA function in the MiCA3 software (Shyu et al., 2007). The dsrAB amplicons were digested with MboI (isoschizomer of NdeII) as previously described by Pérez-Jiménez and Kerkhof (2005).
Terminal restriction fragments (T-RFs) were separated by capillary electrophoresis on an ABI 310 sequencer (Applied Biosystems) and sized using GeneScan-500 LIZ as size standard. T-RFLP profiles were analyzed using GeneMarker 2.4.0 software (SoftGenetics). The T-RFs with sizes ranging from 40 to 500 bp and peak heights of more than 35 Relative Fluorescence Unit were included in the analysis. Data were standardized by dividing the area of each peak by the total peak area. The T-RFs comprising less than 0.5% of the total peak area were excluded from the analysis and the percentage values of the remaining T-RFs were recalculated (Rees et al., 2004).

PAH Pollution Descriptors
The following PAH pollution descriptors were investigated and statistically analyzed: (i) the concentrations of the 16-EPA priority PAHs (Supplementary Table S3); (ii) the sum of the concentrations of the 16-EPA priority PAHs ( 16PAH) assigning the values to four pollution categories according to Baumard et al. (1998): low = < 100 ng g −1 ; moderate 100-1,000 ng g −1 ; high 1,000-5,000 ng g −1 ; very high >5,000 ng g −1 ; (iii) the sum of four LPAHs (Phen, Ant, Flu, Pyr, LPAH) and the sum of eight HPAHs (BaA, Chr, BbF, BkF, BaP, Inp, BgP, DBA, HPAH) as defined by Soclo et al. (2000) as descriptors of PAH pollution levels; (iv) the molar ratios of selected PAHs as descriptors of PAH sources (Supplementary Table S4): LPAH/ HPAH, Ant/Ant + Phen, Flu/Flu + Pyr, Frontiers in Marine Science | www.frontiersin.org BaA/BaA + Chr, Inp/Inp + BgP (Yunker et al., 2002). More specifically, the PAH pollution sources were inferred from an increase in the proportion of the less stable, "kinetic" PAH isomers relative to the more stable, "thermodynamic" isomers. We examined four different ratios, since the use of multiple diagnostic ratios provides safer interpretations in case of samples contaminated by PAHs from mixed sources (Yunker et al., 2002). In order for the data obtained to be properly interpreted, the transition points reported by Yunker et al. (2002) were applied: (i) for the Ant/Ant + Phe ratio, a value <0.10 is indicative of petroleum contamination, while a value >0.10 indicates the dominant contribution of combustion sources; (ii) the BaA/BaA + Chr ratio indicates petroleum inputs for values <0.20, mixed sources for values between 0.20 and 0.35, and combustion sources for values >0.35; (iii) a Flu/Flu + Pyr ratio is characteristic of petroleum inputs for values <0.40, combustion of liquid fossil fuels and crude oils for values between 0.40 and 0.50, whereas values >0.50 indicate combustion of grass, wood, or coal; (iv) the Inp/Inp + BgP ratio implies petroleum inputs for values <0.20, combustion of liquid fossil fuels and crude oils for values between 0.20 and 0.50, and combustion of grass, wood and coal for values >0.50.

Statistical Analyses
Four matrices including all the data on the prokaryotic and abiotic variables were constructed: (i) bacterial T-RFs of the 16S rRNA gene (as rows) by sampling stations (as columns); (ii) archaeal T-RFs of the 16S rRNA gene by sampling stations; (iii) T-RFs of the dsrAB gene by sampling stations; (iv) abiotic variables (as row) by sampling stations (as columns). Dissimilarities between sampling stations based on the prokaryotic variables were calculated by means of the Bray-Curtis dissimilarity coefficient, whereas for those based on the abiotic variables the Euclidean distances was used. A Euclidean distance matrix of geographical position was calculated from the geographical coordinates of sampling stations ( Table 1). Dissimilarity matrixes were calculated using the vegdist function in vegan R package (Oksanen et al., 2013).
For the abiotic parameters, the Kruskal-Wallis rank nonparametric ANOVA was used for comparing the median of values measured on all samples collected in each port. The test was performed by applying the kruskal.test function in the stats package of the R software (R Core Team, 2013). In addition, the Wilcoxon rank sum test was performed for the post hoc pairwise comparisons (among ports) by using the wilcox.test function and finally the Bonferroni correction for multiple comparisons was applied by means of the p.adjust function in the stats package.
An ordination analysis of the three targeted prokaryotic communities, that is Bacteria, Archaea and SRB, based on Bray-Curtis dissimilarity index between individual T-RFLP profiles, was performed by means of non-metric multidimensional scaling (nMDS) using the metaMDS function (100 random starts) in the vegan package. The function ordisurf was used to relate the un-constrained community ordination analysis provided by nMDS with the abiotic variables by using a non-linear model (gam model).
For the three targeted prokaryotic communities, the heterogeneity in their structure was estimated using the betadisper function of the vegan package. For each port, the function identifies the group centroid in the ordination space and then calculates the distance of all samples collected within a port from its group centroid (Lee-Cruz et al., 2013). We used the permutest function (999 permutations) to test whether the community structure heterogeneity level significantly differed among ports.
The Bray-Curtis dissimilarity matrices were also tested by permutational multivariate analysis of variance (PERMANOVA) using the adonis function (9,999 permutations) of vegan package in order to estimate the effect on the prokaryotic communities structure of two factors, sampling period and port (Cagliari, El Kantaoui, and Heraklion).
BIOENV and RELATE routines of the PRIMER v.6 software (Clarke and Gorley, 2006) were used to investigate the relation between the structure of the three prokaryotic communities (Bacteria, Archaea, SRB) and the measured abiotic variables. The environmental parameters (silt-clay ratio, salinity, temperature, redox potential, organic carbon, total phytopigment concentration, PAP ratio), the levels of hydrocarbon pollutants ( AH, UCM, Naph, LPAH, and HPAH) and the descriptors of PAH pollution sources ( LPAH/ HPAH, Ant/Ant + Phen, Flu/Flu + Pyr, BaA/BaA + Chr, Inp/Inp + BgP) were included as abiotic variables in the analysis. Highly correlated abiotic variables (Spearman's ρ > 0.9) were mutually excluded from the analysis. For each targeted prokaryotic group, BIOENV routine was used to calculate the Spearman rank correlation coefficient between the Euclidean distance similarity matrix for abiotic variables and the Bray and Curtis similarity matrix computed for T-RFLP profiles. Once the best subsets of variables associated with each biotic community structure were obtained by the BIOENV routine, the RELATE routine was applied on the same subsets to assess the significance of the biotic-abiotic relation.
The Mantel test (Ramette, 2007) was used to assess the overall correlations between the structure of prokaryotic communities and the emission sources of PAHs or the compositions of the pollutant mixtures in superficial sediments. The Mantel test was applied using the mantel function (9,999 permutations) in the vegan package to correlate each Bray-Curtis distance matrix of prokaryotic community structure with the Euclidean distance matrix of: (i) concentrations of individual 16-EPA priority PAHs (Supplementary  Table S4), (iv) geographical coordinates of the sampling stations ( Table 1).
The partial Mantel test was subsequently performed by applying the mantel.partial function in vegan package to evaluate the correlation between each Bray-Curtis distance matrix of community structure and the Euclidean distance matrix of each set of pollutant descriptors (PAHs, AHs, PAH ratios) while controlling the effect of spatial autocorrelation with the Euclidean distance matrix of geographical coordinates of sampling stations (GEO). This allows us to determine the relationship between community structures and pollutant mixtures after the effects of spatial autocorrelation have been removed.

Levels, Composition and Origins of PAH Pollution
The sums of PAH concentrations ( 16PAH, LPAH, and HPAH), as measured in the surface sediments of the different stations and across the three ports along with their statistical analysis, are reported in Table 2 and Supplementary Table S5, respectively. With the only exception of sediments near the shipyard in Heraklion, the concentrations of HPAH were always higher than those of LPAH ( Table 2). More specifically, the compounds with molecular weights 202, 228, 252, and 276 were the most abundant components in all samples. Typically, parent PAHs with molecular weights 202 and 252 presented the highest levels (Supplementary Table S3).
Sediments from Cagliari exhibited significantly higher concentrations and wider variation (CV%) in 16PAH, LPAH and HPAH than those from the other two ports. With regard to Heraklion and El Kantaoui, a significant difference in PAHs was observed only for LPAH, which presented significantly higher levels in Heraklion than in El Kantaoui. In sediments collected from El Kantaoui, 16PAH concentrations ranged from 25 ng g −1 at the port entrance (E3) to 367 ng g −1 near the fuel station (E2) and pollution was classified as "low" (0-100 ng g −1 ) and "moderate" (100-1,000 ng g −1 ), respectively, according to the categories determined by Baumard et al. (1998). In Heraklion port, the 16PAH varied from 89 to 1,987 ng g −1 ranging in the categories from "low" to "high" in the different sectors. Particularly two stations, one located near the docking facility for leisure and fishing crafts (H1; Min 1,418; Max 1,885 ng g −1 ) and one near the shipyard (H5, Min 564; Max 1,987 ng g −1 ), fell in the "high" (1,000-5,000 ng g −1 ) pollution category. The PAH levels in Cagliari varied from 686 ng g −1 at the port entrance (C5) to 48,999 ng g −1 in the sector hosting the military navy vessels (C2), indicating pollution categories of "moderate, " "high, " and "very high" (>5,000 ng g −1 ). Table 2 reports the LPAH/ HPAH ratio as a reliable indicator for discriminating the origin of PAHs. The lower the value of the ratio, the higher the prevalence of pyrolytic on petrogenic PAHs is (De Luca et al., 2004). The values of the LPAH/ HPAH ratio were significantly lower in the sediments from Cagliari and El Kantaoui than in those collected from Heraklion, with the station near the shipyard (H5) exhibiting the highest values. No significant differences were found between sediments of Cagliari and El Kantaoui (Supplementary Table S5).
The scatter plots of the molar ratios Flu/Flu + Pyr versus Ant/Ant + Phen, Inp/Inp + BgP, or BaA/BaA + Chr are showed in Figure 2. Values of the four selected diagnostic ratios at the three ports and their statistical analysis are reported in Supplementary Tables S4, S5, respectively. The Cagliari port showed significantly higher values of Ant/Ant + Phen and BaA/BaA + Chr ratios as compared to the other two port areas. Differences between El Kantaoui and Heraklion were not statistically significant. The station located near the shipyard in Heraklion (H5) was the only one showing an Ant/Ant + Phen ratio typical for unburned petroleum sources. Except for H5 (shipyard) and E3 (port entrance), where the BaA/BaA + Chr ratio suggested a mixed pollution origin, a clear pyrogenic fingerprint was evident in all harbor sediments with the station hosting military navy vessels (C2) characterized by the highest values of both Ant/Ant + Phen and BaA/BaA + Chr ratios.
In order to evaluate whether PAHs in sediments originated from fuel combustion or from combustion of biomasses and coal, the molar ratios of Flu/Flu + Pyr and Inp/Inp + BgP were investigated. With the only exception of the port entrance (E3), sediments from El Kantaoui were always characterized by a predominance of PAHs from petroleum combustion with an excellent consistency between the Flu/Flu + Pyr and Inp/Inp + BgP ratios. In sediments from the Cagliari port, values of both diagnostic ratios were consistent with a PAH origin from combustion of biomasses and coal. In the port of Heraklion, stations hosting passenger (H3) and cargo (H4) ships were characterized by a contribution of petroleum combustion, while a influence from combustion of biomasses and coal was found at the leisure and fishing boatdocking station (H1).

Structure of Prokaryotic Communities
The nMDS ordination plots resulting from the T-RFLP profiles of the three targeted prokaryotic communities are shown in Figure 3. The community structure of Bacteria and Archaea in the three investigated ports revealed to be fairly similar as highlighted by the partial overlap of community profiles (i.e., no evident separation between different colored points), albeit 16SArch-TRFLP showed a greater divergence between port areas. A clear grouping of SRB communities was evident on the basis of the port. Notably, sediments near the shipyard (H5) differentiated from the other Heraklion samples (i.e., positioned distantly in the ordination space) by 16SArch-TRFLP as well as by dsr-TRFLP.
The PERMANOVA analysis (Supplementary Table S6) demonstrated that the "port" was the main factor affecting the three targeted prokaryotic communities (p < 0.01). The "sampling period" was also found to significantly (p < 0.01) affect the T-RFLP profiles for both Bacteria and Archaea but was not a factor significantly separating the SRB communities.
As revealed by TRFLP of bacterial 16S rRNA gene, the port of El Kantaoui showed the highest structure heterogeneity of bacterial communities while the Cagliari port had the lowest one, with the difference between the two ports being statistically significant (p < 0.01) (Figure 4A). On the other hand, the port of Heraklion showed the highest heterogeneity in dsr-TRFLP profiles, but this difference was proven statistically significant (p < 0.05) only in comparison to El Kantaoui ( Figure 4C). No significant differences (p > 0.05) were found among the three ports in terms of structure heterogeneity 2 | Sums of PAH concentrations (ng g −1 dw), pollution levels according to Baumard et al. (1998), and LPAH/ HPAH ratio in surface sediments from the ports of Cagliari, El Kantaoui, and Heraklion. of the archaeal communities as determined by 16SArch-TRFLP ( Figure 4B).

Linking Prokaryotic Communities to Environmental Parameters and Pollution Descriptors
The BIOENV routine, revealed weak (ρ ≤ 0.4) or moderate (0.4 < ρ < 0.6) but significant (p < 0.05) correlations of combined abiotic variables with the community structure, as follows: (i) for Bacteria, two environmental parameters (temperature, redox potential), the total level of AHs ( AH) and the three descriptors of PAH pollution sources, namely Flu/Flu + Pyr, BaA/BaA + Chr, and Inp/Inp + BgP (RELATE ρ = 0.355, p = 0.002); (ii) for Archaea, salinity, temperature, content of organic carbon, the total level of AHs ( AH), as well as the three diagnostic ratios of pollution sources, LPAH/ HPAH, Flu/Flu + Pyr, and BaA/BaA + Chr (RELATE ρ = 0.381, p = 0.001); (iii) for SRB, the total level of AHs ( AH) and three descriptors of PAH pollution sources, FIGURE 2 | Scatter plots between molar ratios of selected PAHs calculated for surface sediments collected at the port of Cagliari (blue), El Kantaoui (green), Heraklion (pink). Flu/Flu + Pyr vs. Ant/Ant + Phe (A), Flu/Flu + Pyr vs. Inp/Inp + BgP (B), Flu/Flu + Pyr vs. BaA/BaA + Chr (C). Dashed lines identify the transition points between different pollution sources according to Yunker et al. (2002). The digit in the symbol identifies the sampling station within each port area.
namely LPAH/ HPAH, Flu/Flu + Pyr and BaA/BaA + Chr (RELATE ρ = 0.479, p = 0.001). Table 3 presents the results of the Mantel test assessing the overall correlations between the structure of prokaryotic communities and the composition profiles of individual PAHs and individual AHs as well as the source-related diagnostic ratios of PAHs in sediments. Significant (p < 0.05) but weak correlations (ρ values ranging from 0.238 to 0.365) were revealed between the structures of the three targeted communities and the molar ratios of PAHs, while no significant correlations were found with the levels of both aliphatic ( AH) and aromatic hydrocarbons ( LPAH, HPAH). However, it is important to stress that the "port" was also a significant factor shaping the structure of the three targeted prokaryotic communities, as demonstrated by the PERMANOVA analysis (Supplementary  Table S6). Therefore, it is reasonable to suppose that other ecologically important variables, not considered in this study, may influence the structure of the prokaryotic assemblages in the investigated sediments leading to the emerging weak correlations between community structure and the measured parameters due to autocorrelation (i.e., purely spatial effect). A second Mantel test confirmed the presence of significant correlations between the geographic locations of the sampling stations and the community structures of Bacteria, Archaea, SRB (spatial autocorrelation). With the only exception of dsr-TRFLP profiles, correlations between community structures and the geographical positions of sampling stations were weak. The effect of sourcerelated PAHs ratios was thus untangled from the one owed to autocorrelation by partial Mantel test (Table 3). Among the three targeted groups, only the structure of SRB was found to be significantly (p < 0.0001) correlated with the sourcerelated diagnostic ratios of PAHs when the effects of spatial autocorrelation were removed (Table 3).
In line with the results obtained by the BIOENV and RELATE routines as well as by the partial Mantel tests, the ordisurf function predicted a significant (p < 0.001) fit between the molar ratios of PAHs and the position of the different samples on the nMDS ordination plot based on the dsr-TRFLP profiles (Figure 5). The curved isolines in Figure 5 roughly illustrated an increase of the molar ratio of Inp/Inp + BgP (R = 0.653) along the diagonal of the nMDS plot, from top left to bottom right. This specific trend is indicative of an increasing contribution of pyrogenic PAHs in the sediment samples along the nMDS plot diagonal. Accordingly, the SRB communities in sediments from Cagliari were located in the lower-right corner of the plots, corresponding to a greater contribution of pyrolytic PAH sources (high values of Inp/Inp + BgP ratio).

DISCUSSION
Polycyclic aromatic hydrocarbons are among the most widespread hazardous organic pollutants in the Mediterranean Sea, which is one of the most severely polluted marine areas worldwide (Abdulla and Linden, 2008;Durrieu de Madron et al., 2011). Our results are in good agreement with the chemical characterization of other harbors and marinas in the Western part (De Luca et al., 2004Mzoughi and Chouba, 2011;Schintu et al., 2015) and in the Eastern part (Mostafa et al., 2003;Barakat et al., 2011;Kapsimalis et al., 2014) of the Mediterranean Sea. Pyrogenic processes are found to be the dominant emission sources of PAHs entering the sediments of the three studied ports. Pollutant mixtures are  dominated by HPAHs having a molecular weight between 202 and 252. This finding has been considered as indicative of anthropogenic contamination since these compounds have been commonly observed in sediments adjacent to or downstream of urbanized and industrialized areas (Yunker et al., 2002). The level and origin of PAHs in superficial sediments vary considerably among the investigated ports. In the port of El Kantaoui, the origin of PAHs is primarily attributed to the incomplete combustion of fuels (i.e., gasoline, diesel, lubricating oils), in agreement with its small dimension, moderate marine traffic, and the narrow spectrum of harbor activities. The pollution levels are comparable to those previously measured in Tunisian coastal sediments (Zaghden et al., 2007) and they are 5-to 12-fold lower than the values recorded by Mzoughi and Chouba (2011) in the Tunisian marina of Sidi Bou Said. The Heraklion port is characterized by the coexistence of different PAH sources, mainly coal combustion and fuel combustion. A contribution of petrogenic sources is also evident in sediments near the shipyard, which could be attributed to the vessels repair and maintenance activities. Contamination levels in Heraklion port are comparable to those previously recorded in other medium-size harbors of the Mediterranean basin, such as Porto Torres (De Luca et al., 2004). In the port of Cagliari, the primary PAH source is the burning of coal and biomass, which could be explained by atmospheric deposition and street run-offs emitted from the adjacent city. Emission sources and pollution levels in Cagliari are comparable to those frequently recorded in other relevant touristic and commercial Mediterranean harbors (Barakat et al., 2011;Merhaby et al., 2015). The only exception was the sector hosting the military navy vessels, where we found one order of magnitude higher PAH levels. Pollution concentrations higher than those measured in this sector of the Cagliari port have been previously recorded in a limited number of highly urbanized areas along the Mediterranean coasts, including the harbors of Alexandria (Mostafa et al., 2003;Barakat et al., 2011), Naples (Sprovieri et al., 2007), and Genova Punta Vagno (Bertolotto et al., 2003). Overall, the investigated ports enable us to evaluate the link between benthic prokaryotic communities and complex mixtures of PAHs produced by multiple pyrolytic sources and over a wide concentration range (25-49,000 ng g −1 dw). At temporal scale, temperature, driven by season, seems to play an important role in shaping the bacterial and archaeal communities in superficial sediments of the investigated ports. In contrast, changes in SRB assemblages observed in the present study were not ascribed to temporal variations. In line with this finding, previous works have demonstrated that psychrotolerants and mesophiles predominate in SRB communities from temperate sediments (Robador et al., 2009(Robador et al., , 2016, resulting in a seasonally stable community structure (Mussmann et al., 2005;Robador et al., 2009). At spatial scale, the structure of the three benthic prokaryotic communities varies among ports. Moreover, the geographical location of ports, as delimited by their coordinates, is associated with variations in the community structure of Bacteria, Archaea and SRB. Our results are in good agreement with previous studies where the benthic SRB communities have been found to be significantly different among sites at global scale (Pérez-Jiménez and Kerkhof, 2005;Robador et al., 2016). On the other hand, Chatzinikolaou et al. (2018) have recently demonstrated that the three studied ports are markedly different in terms of their physico-chemical parameters, as a result of different geographical positions and local factors. On one side, the Mediterranean Sea is characterized by wellknown longitudinal gradients, such as a west to east increase in salinity and temperature, as well as a decrease in productivity (e.g., Coll et al., 2010). On the other hand, ports are typically characterized by varying hydromorphological properties due to physical alterations (i.e., dredging, confinement) caused by human activities to support specific uses (e.g., navigation). These alterations induce changes in the local environmental conditions (e.g., low hydrodynamism, low oxygen, high organic content) and account for port designation as Heavily Modified Water Body (WFD 2000/60/EC).
Despite the different statistical approaches applied, no direct link was found between any of the three targeted prokaryotic communities and the concentrations of individual PAHs or the total levels of LPAHs and HPAHs. These findings can be well understood in the framework of previous studies demonstrating the combustion process by which PAHs are formed determines not only the composition of the contaminant mixture but also its bioavailability in sediments (Shrestha et al., 2010;Akkanen et al., 2012;Duran and Cravo-Laureau, 2016). Similarly, Jeanbille et al. (2016b) have recently given the first evidence that the structure of bacterial communities in sediments of five highly urbanized Mediterranean coastal areas contaminated by pyrolytic and mixed petrogenic/pyrogenic PAHs was significantly correlated with diagnostic PAH ratios, rather than the actual content of PAHs. However, these correlations were weak since most of the variance observed by those authors remained unexplained, implying that other variables or processes may also play an important role in shaping prokaryotic communities. Accordingly, we found a weak relation between the emission sources of PAHs entering the surface sediments of the three studied ports and the structure of bacterial and archaeal communities. This relation faded away when the strong effect of geographic position was taken into account suggesting that other ecologically important variables (e.g., the mineral composition of sediments, the presence of other organic pollutants) or processes (e.g., geochemical, sediment reworking by benthic meio-and macrofauna, etc.) may shape the structure of bacterial and archaeal assemblages in the sense of what resources are readily available and under which environmental conditions those resources can be exploited by the prokaryotic communities in Mediterranean port sediments. Differences in the emission sources of PAHs among the three investigated ports might be associated with the observed variations in the structure of SRB communities despite the pronounced effect of geographical location. Even though other important environmental variables could contribute to the spatial variation in structure among ports, our results indicate that PAH contamination may be a driver of SRB communities in superficial sediments chronically impacted by multiple anthropogenic pressures. Indeed, we found the highest heterogeneity of SRB communities in the Heraklion port, characterized by the broadest spectrum of PAH emission sources (i.e., petroleum products, fuel combustion and coal burning) among the studied Mediterranean ports. Moreover, the total level of AHs can also play a role in shaping the observed structure of benthic SRB community. Nevertheless, it's crucial to stress that the response of SRB communities cannot be directly attributed to any field measured variable until confirmed by laboratory tests under controlled experimental conditions (Cravo-Laureau and Duran, 2014).

CONCLUSION
The present study support the notion that SRB communities significantly differ between sites with diverse anthropogenic impacts and geographical locations, while temporal variation seems to play only a marginal role as a driver of SRB community change in temperate sediments. According to the published literature, this study provides for the first-time evidence that the emission sources of PAHs entering superficial sediments could influence the specific response of the benthic SRB communities in the context of a chronic contamination dominated by pyrolytic PAHs. The influence of PAH emission sources on SRB assemblages seems to contribute to the local-specific differences in community structure among ports. Currently, PAH source identification is a prerequisite for the evidencebased decision-making in environmental management of port systems, such as the designation of optimal strategies for source control and remediation. The findings of this study indicate the SRB community is a valuable candidate target descriptor for the monitoring of bioremediation processes in port sediments.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
ET and CA developed the concept for this study. CA, TD, ECh, and ET carried out the field work. ET, CA, ECh, GM, and ECa designed the experiments. FV, GS, and SS performed the microbiological analyses. MM performed analysis of hydrocarbons. ET, CA, ECa, and GM contributed reagents, materials, and analysis tools. FV, CA, ET, and RL analyzed the data. ET, FV, and MM drafted the manuscript. All authors edited and reviewed the final version of the manuscript and approved before submission.

FUNDING
This work has been supported with the financial assistance of the European Union under the ENPI CBC Mediterranean Sea Basin Programme in the framework of the project "MAnagement of Port areas in the MEDiterranean Sea Basin" (MAPMED). The contents of this document are the sole responsibility of UNICA, HCMR, and UNIFI and can under no circumstances be regarded as reflecting the position of the European Union or of the Programme's management structures.

ACKNOWLEDGMENTS
The authors would like to thank the representatives of the Port Authority of Cagliari, El Kantaoui, and Heraklion for their constant support during the implementation of the sampling campaigns.