Spatial scales of bacterial community diversity at cold seeps (Eastern Mediterranean Sea).

Cold seeps are highly productive, fragmented marine ecosystems that form at the seafloor around hydrocarbon emission pathways. The products of microbial utilization of methane and other hydrocarbons fuel rich chemosynthetic communities at these sites, with much higher respiration rates compared with the surrounding deep-sea floor. Yet little is known as to the richness, composition and spatial scaling of bacterial communities of cold seeps compared with non-seep communities. Here we assessed the bacterial diversity across nine different cold seeps in the Eastern Mediterranean deep-sea and surrounding seafloor areas. Community similarity analyses were carried out based on automated ribosomal intergenic spacer analysis (ARISA) fingerprinting and high-throughput 454 tag sequencing and were combined with in situ and ex situ geochemical analyses across spatial scales of a few tens of meters to hundreds of kilometers. Seep communities were dominated by Deltaproteobacteria, Epsilonproteobacteria and Gammaproteobacteria and shared, on average, 36% of bacterial types (ARISA OTUs (operational taxonomic units)) with communities from nearby non-seep deep-sea sediments. Bacterial communities of seeps were significantly different from those of non-seep sediments. Within cold seep regions on spatial scales of only tens to hundreds of meters, the bacterial communities differed considerably, sharing <50% of types at the ARISA OTU level. Their variations reflected differences in porewater sulfide concentrations from anaerobic degradation of hydrocarbons. This study shows that cold seep ecosystems contribute substantially to the microbial diversity of the deep-sea.

Introduction

Cold seeps are seafloor ecosystems fueled by chemical energy originating from microbial transformation of hydrocarbons and sulfide. They are characterized by the largest biomass and highest productivity of all deep-sea ecosystems (Jørgensen and Boetius, 2007; Levin and Sibuet 2012; Boetius and Wenzhöfer 2013). Because of the localized, spatially and temporarily dynamic supply of hydrocarbons from deep subsurface reservoirs, cold seeps have a highly fragmented distribution along continental margins and represent isolated habitats on the vast deep-sea floor (Sibuet and Olu, 1998; Levin, 2005). Main questions about biodiversity patterns, dispersal capabilities and life history of fauna and microbes restricted to these chemosynthetic habitats need to be answered in order to understand the interconnectivity and resilience of these dynamic ecosystems (Tyler and Young, 1999; Pradillon , 2007). Moreover, as biodiversity changes over multiple spatial scales (Green and Bohannan, 2006; Ramette and Tiedje, 2007; Martiny ; Hanson ), a better understanding of the interconnectivity between these fragmented ecosystems needs to be based on quantitative assessments of spatial patterns of species diversity.

Members of the domains Bacteria and Archaea contribute the main functions, that is, hydrocarbon degradation, sulfide production and consumption and chemosynthetic CO2 fixation (Jørgensen and Boetius, 2007) at cold seep and at other types of chemosynthetic ecosystems. The prevailing energy sources at seeps select for specific types of bacteria and their metabolic functions (Knittel and Boetius, 2009; Goffredi and Orphan, 2010), which differ from those of the surrounding oligotrophic, particle-flux dependent deep-sea ecosystem (Orcutt ; Quaiser ). But still little quantitative information is available about the similarity and spatial variation of bacterial communities at cold seeps (Pop Ristova ) and at other types of chemosynthetic ecosystems (Huber ; Bienhold ), as compared with the deep-sea floor (Jacob and references therein). At the global scale, it was shown that chemosynthetic ecosystems displayed the highest heterogeneity in bacterial community structure across all oceanic ecosystems (Zinger ). Spatial turnover, that is, the relative change in community composition along a spatial gradient (Anderson ), of bacterial communities in the deep-sea was previously found to occur not only at large spatial scales in pelagic sediments (>1000s of km; Schauer ; Zinger ) but also to lesser extent at small scales (<1 km; Pop Ristova ; Meyer ; Ruff ) among reduced habitats of individual systems, mainly due to variations in energy availability and environmental heterogeneity, as well as spatial separation.

Here we investigated bacterial diversity patterns of seep ecosystems of the oligotrophic Eastern Mediterranean deep-sea, with the specific aim to understand how bacterial communities of fragmented, highly reduced sediment habitats vary with spatial scales. Bacterial communities of nine locations distributed over three seep regions were studied using a combined approach, including bacterial community fingerprinting and 454 massively parallel tag sequencing (MPTS), porewater geochemistry, quantification of element fluxes and microbial consumption rates. This allowed for a comparison of microbial diversity at geographical scales ranging from tens of centimeters to hundreds of kilometers. The main aims were to better understand (i) how bacterial communities of reduced methane-seeping habitats compare with those inhabiting surrounding oxidized deep-sea sediments, (ii) which biogeochemical factors shape the bacterial diversity of cold seeps and (iii) if and how bacterial communities of cold seep ecosystems vary with spatial scales.

Material and methods

Description of sampling sites and of analytical procedures

During two consecutive expeditions MSM13/3 and 4 in 2009, on board the RV Maria S Merian, similar types of reduced sediments were sampled from three cold seep regions in the Eastern Mediterranean deep-sea, that is, the Amon Mud Volcano (AMV), the Central Pockmark (Pock) area of the Nile Deep-sea Fan and the Amsterdam Mud Volcano (AmsMV), at depths of 1120–2030 m, located at 128–355 km distance from each other (Table 1, Supplementary Table S1 and Figure 1). In each of the three regions, 2–4 sites were sampled upon visual identification as typical patches of greyish to blackish color, some of which were covered by thin, white thiotrophic bacterial mats (Figure 2b). Geochemically, all seep sites were characterized by methane-seepage and high sulfide concentrations (Table 2). Of these, the AMV_seep_2 samples were obtained from a recently gas-vented seafloor area in the active center of the Amon MV, with not yet fully established active methane-oxidizing community (Felden ). Additional non-seep sediment samples were obtained from the immediate vicinity of the cold seep sites (Table 1) and also from distant sites outside of the seepage regions sampled during previous missions (Table 3, also see Supplementary Information).

Table 1

Overview of the main sampling sites investigated in this study

Cold seep	Location	Habitat	Sampling site ID	Longitude; latitude	Sampling time
Amon MV (AMV; water depth 1120 m)	Center	Reduced sediment with bacterial mat	AMV_seep_1	32°22.132272'N	October–November 2009
				31°42.654951'E
		Reduced sediment	AMV_seep_2	32°22.174029'N	October–November 2009
				31°42.627174'E
	Sulfur band	Reduced sediment with bacterial mat (M14 marker)	AMV_seep_3	32°22.045089'N	October–November 2009
				31°42.276642'E
		Reduced sediment with bacterial mat (M16 marker)	AMV_seep_4	32°22.045956'N	October–November 2009
				31°42.265083'E
		Beige non-seep sediment	AMV_non-seep_1	32°22.3002'N	October–November 2009
				31°41.9598'E
Amsterdam MV (AmsMV; water depth 2030 m)		Reduced sediment (M5 marker)	AmsMV_seep_1	35°20.034018'N	November–December 2009
				30°16.167342'E
		Reduced sediment (M6 marker)	AmsMV_seep_2	35°20.079390'N	November–December 2009
				30°16.129803'E
	Center	Reduced sediment (M10 marker)	AmsMV_seep_3	35°19.945893'N	November–December 2009
				30°16.098681'E
	E Rim	Beige non-seep sediment	AmsMV_non-seep_1	35°19.9602'N	November–December 2009
				30°16.8198'E
	NE-E Rim	Beige non-seep sediment	AmsMV_non-seep_2	35°20.2398'N	November–December 2009
				30°17.2902'E
	Center	Beige seep-influenced non-seep sediment	AmsMV_non-seep_3	35°20.0598'N	November–December 2009
				30°16.1202'E
Pockmark (Pock; water depth 1690 m)		Bacterial mat at marker M10	Pock_seep_1	32°31.99'N	November–December 2009
				30°21.12'E
		Bacterial mat at marker M9	Pock_seep_2	32°31.98'N	November–December 2009
				30°21.15'E
		Beige seep-influenced non-seep sediment	Pock_non-seep_1	32°32.07'N	November–December 2009
				30°21.39'E

Detailed list of all samples and measurements performed at each sampling sites are shown in Supplementary Table S1. Physical markers (M) were deployed at almost all the investigated habitats for easier visual recognition.

Figure 1

Overview maps of all three investigated cold seep regions, including seep and non-seep sites in the Eastern Mediterranean Sea (a) the Amsterdam MV (b), the Pockmark area (c), and Amon MV (d). The main sampling and measurement sites are depicted with symbols, that is, push core and multicore sampling (green triangles), benthic chamber incubations (yellow squares) and microprofiler measurements (black circles).

Figure 2

Underwater photographs of reduced seep habitats investigated in this study: Upper panel (a, b, c) photographs were taken at the Amon MV and lower panel (d) photographs at the Amsterdam MV and (e, f) at a Pockmark seep. (a) Center of the Amon MV characterized by gassy sediments and patchy bacterial mats, (b, c) a lateral mud flow with extensive bacterial mats at the outer rim of the volcano. (d) Typical methane-seeping habitat marked by black sediment patches encountered at the Amsterdam MV. Overview (e) and a more detailed photo (f) of a bacterial mat habitat at the Pockmark cold seep.

Table 2

Biogeochemical characterization of seep and nearby non-seep sites, including average depth integrated (0–10 cm) sulfate reduction (SR) rates, total oxygen uptake (TOU), methane efflux, diffusive oxygen uptake (DOU), oxygen penetration depth (OPD), minimum detected concentration of sulfate and sulfate penetration depth

Cold seep	Sampling site	SR	TOU	CH4 efflux	DOU	OPD	Min. SO4	SO4 depth
		mmol m−2 d−1	mmol m−2 d−1	mmol m−2 d−1	mmol m−2 d−1	mm	mmol l−1	cm
Amon MV	AMV_seep_1	16	115	83; 85; 49	47	2	<1	10
	AMV_seep_2	2	5	1; 70; 24	11	3	29	>19
	AMV_seep_3	28	119	<1	3; 6	1; 3.9	<1	5
	AMV_seep_4	1	13	ND	5	45	29	>18
Amsterdam MV	AmsMV_seep_1	35	101	1175	ND	ND	<1	15
	AmsMV_seep_2	8	23	196	ND	ND	4	>15
	AmsMV_seep_3	13	89	1169	ND	ND	<1	8
Pockmark	Pock_seep_1	20	153	688	ND	ND	<1	12
	Pock_seep_2	50	228	36	ND	ND	1	20
Amon MV	AMV_non-seep_1	0a	5a	ND	12a	>40a	29b	>13b
Amsterdam MV	AmsMV_non-seep_1	0	ND	ND	ND	ND	27	>30
	AmsMV_non-seep_2	1	ND	ND	ND	ND	28	>25
	AmsMV_non-seep_3	56	ND	ND	ND	ND	27	>17
Pockmark	Pock_non-seep_1	0	1a	ND	0.5a	>10a	27	>29

Abbreviation: ND, not detected. For further details on the Amon MV biogeochemistry, see Felden .

Intermediate sites with subsurface seepage include AMV_seep_4 and AmsMV_non-seep_3.

Grünke .

Data derived from M70/2a cruise.

Table 3

454 MPTS—derived sample characteristics and species richness estimates

Sample ID	Number of sequences	Number of OTU0.03	Number of OTU0.03 singletons	Relative number of OTU0.03 singletons (%)	Chao1	ACE	Shannon index
AMV_seep_1	10988	1478	700	47	3592	6794	5.6
AMV_seep_3	9062	1892	1073	57	6164	11943	5.5
AMV_seep_4	11661	2250	1198	53	5041	9526	6.1
AmsMV_seep_1	9929	1216	451	37	2863	4780	5.1
AmsMV_seep_2	8598	1192	460	39	3047	6060	4.9
AmsMV_seep_3	13186	1487	623	42	2793	4833	5.0
Pock_seep_1	15774	1640	740	45	3015	5778	5.1
Pock_seep_2	7925	1023	418	41	3120	5919	5.0
Emed_non-seep_region_1a	19369	1535	527	34	ND	ND	5.3
Emed_non-seep_region_2a	27606	2280	864	38	ND	ND	5.2
Emed_non-seep_region_3a	15422	1522	538	35	ND	ND	6.0

Abbreviations: ND, not determined; OTU, operational taxonomic unit.

Three regions without cold seep systems (Emed_non-seep_region_1–3) were included to test for differences in biodiversity at larger distances.

Sediment samples from non-seep regions were taken at 1375–2968 m water depth in the deep Eastern Mediterranean Sea (Boetius and Lochte, 1996).

Total oxygen fluxes and methane effluxes were determined in situ with a remotely operated vehicle (ROV QUEST, Marum, Bremen, Germany) benthic chamber module, in which changes in solute concentrations were monitored with an oxygen optode (Aanderaa, Bergen, Norway) and a preprogrammed syringe system sampling bottom waters for methane concentrations (Felden and references therein). A deep-sea modular microprofiler was used to carry out high-resolution (200 μm) microsensor measurements of oxygen concentrations and to calculate diffusive oxygen uptakes, as described elsewhere (Glud ). Total numbers of single cells (1-cm resolution, about 20-cm sediment depth) were determined by applying the acridine orange direct count method according to Boetius and Lochte (1996). Porewater was extracted at a resolution of 1 cm using the Rhizon moisture samplers (Seeberg-Elverfeldt ; pore size 0.1 μm), and subsamples were immediately fixed for different types of analyses, that is, NH4+, total H2S, SO42−, PO43−, dissolved inorganic carbon (DIC) and alkalinity, in the home laboratory. Sulfate reduction (SR) rates were determined ex situ by the whole-core injection method with a 35SO4 radiotracer (Kallmeyer ; Felden ). Detailed description of all analytical procedures can be found in the Supplementary Information. All biogeochemical data are available open access via the PANGAEA database (Pop Ristova ; http://doi.pangaea.de/10.1594/PANGAEA.830241).

Bacterial community structure at seep and nearby non-seep sites was assessed with automated ribosomal intergenic spacer analysis (ARISA; 1-cm layers, 5-cm sediment depth). In short, DNA was amplified in triplicates with the ITSR and HEX-labeled ITSF primers (Cardinale ), targeting the internal transcribed spacer region as described elsewhere (Pop Ristova , also see Supplementary Information). In addition, bacterial communities of seep sites were analyzed by 454 MPTS of the V4–V6 region of the 16 S rRNA gene (10-cm sediment depth) from the same DNA extracts, using the 518 F and 1064 R primers, and following the protocols published on http://vamps.mbl.edu. Data were processed with mothur (v.1.29) following the standard operating procedure (Schloss ; also see Supplementary Information for detailed information). SFF files of all sequence data have been submitted to the GenBank Sequence Read Archives (http://www.ncbi.nlm.nih.gov) under BioProject ID: PRJNA244496. For the comparison of the community composition of seep regions (V4–V6-based sequences) to non-seep regions (V6-based sequences), only the V6 hypervariable region was compared, and data sets were processed together using the mothur (see Supplementary Information for detailed information). Comparability between the derived V6-excised and the original data set was high and significant at the level of beta-diversity (change in community structure between sites), Shannon diversity index and the taxonomical level of class (Supplementary Table S2 and Supplementary Figures S1–S3).

Statistical analyses

ARISA peaks were binned (2-bp bin size; Interactive Binner function, http://www.ecology-research.com; Ramette, 2009), and ARISA PCR triplicates were merged into a consensus profile by keeping only operational taxonomic units (OTU) that appeared in at least two of the PCR replicates and by averaging the corresponding relative peak intensities. 454 MPTS sequences were clustered at 97% sequence's similarity level (OTU0.03). Singletons were counted as OTU0.03 represented by only one sequence in the whole data set. Shannon diversity index was used to compare the diversity of seep and non-seep regions, as this index was shown to be reliable and stable across studies and data sets based on different primer sets (He ). Additionally, variations in the diversity of seep sites were compared by calculating the Chao1 and ACE richness estimates (same primer set). For both analyses, sample diversity was calculated on normalized 454 MPTS subset based on the sample with the least number of sequences (Pock_seep_2, 7925), using mothur v.1.29.1 (Schloss ). Percentage of shared OTUs was calculated by pairwise comparison to assess similarity between sites at small (<10 cm; within individual sites), intermediate (50–630 m, between sites within one cold seep region) and large spatial scales (128–355 km, between sites of different cold seep regions). All statistical analyses were performed in the R statistical computing language (v.2.15), using vegan (Oksanen ) and labdsv (Roberts, 2013). Indicator taxa characterizing each cold seep structure were identified via indicator species analysis according to Dufrêne and Legendre, (1997). Dissimilarity matrices based on community (pooled ARISA replicates or 454 MPTS) and on porewater data (DIC, H2S and SO4) were calculated using Bray–Curtis and Euclidean distances, respectively. For all analyses, ARISA profiles of individual depth sediment layers were pooled in silico, except for the analyses of correlation between diversity and porewater geochemistry, as well as investigation of bacterial community structure of seep and non-seeps. Community differences were visualized with non-metric multidimensional scaling analyses and tested for significance using the analysis of similarity test. Mantel test with 999 Monte-Carlo permutations was used to test for significance of Spearman and Pearson correlations between community and environmental (for example, differences in porewater concentrations) distance matrices. Mantel P-values were corrected for multiple testing using the Bonferroni's correction (Ramette, 2007).

Results

Geochemical characterization of seep versus non-seeps habitats

At all three cold seep systems, the Amon Mud Volcano (AMV), Amsterdam Mud Volcano (AmsMV) and the Pockmark area (Pock), we sampled visually and biogeochemically distinct patches of highly reduced, methane-seeping sulfidic sediments of dark grey to blackish color, which were separated by non-reduced oxygenated seafloor areas (Figure 1). The areas directly surrounding the reduced patches consisted of beige sediments, showed little to no influence of methane seepage and were sampled as ‘non-seep sites'. Microsensor measurements confirmed that all non-seep sites were fully oxygenated until the penetration depth of the sensors (>1–4 cm sediment depth) (Figure 3f). Throughout the top 20 cm, these had sulfate and DIC porewater concentrations of 28 and 2.3 mmol l−1, respectively, reflecting bottom water concentrations (Table 2, Supplementary Figures S4 and S5). AMV_seep_4 and AmsMV_non-seep_3 were intermediate sites with some signs of subsurface seepage. Total oxygen uptake (TOU) and diffusive oxygen uptake were low (1–12 mmol m−2 d−1), and SR rates barely detectable (Table 2, Supplementary Figure S6). Ammonium and phosphate concentrations in the porewaters were also low, with, on average, 2 μmol l−1 and 1 μmol l−1, respectively (Supplementary Figure S7). Elevated NH4+ and subsurface SR rates measured at AmsMV_non-seep_3 indicate horizontal methane migration, placing this sample biogeochemically between a seep and a non-seep site. At all non-seep sites, cell numbers of 0.2–1.7 × 109 cm−3 were detected in the top 5-cm surface layers, as well as a decline with sediment depth (Supplementary Figure S8).

Figure 3

Plots of porewater concentration profiles at methane-seeping, reduced sites (a–c; Pock_seep_1) and at nearby non-seep sites from the Pockmark region (d–f; Pock_non-seep_1). (a, d) Depth profiles of total sulfide (closed circles) and sulfate concentrations (open circles); (b, e) DIC concentration (open circles) and alkalinity (closed circles); (c, f) in situ oxygen concentration determined with microsensors, the insert in panel (c) shows the same profile at higher depth resolution. Further detailed geochemistry plots of all regions can be found in the Supplementary Figures S1–S5.

Seepage of dissolved methane to the bottom water was detected at all seep sites, with maximum effluxes of 1169 and 1175 mmol m−2 d−1 measured at two of the Amsterdam MV sites (Table 2). Methane efflux measured at the AMV_seep_2 site varied between 1 and 70 mmol m−2 d−1 within days of sampling, reflecting high temporal and/or spatial variations in the seepage of methane at this location. At most of the seep sites, oxygen penetration was limited to the topmost 0.4 cm of sediment (Figure 3c), except at the AMV_seep_4 where oxygen was detected down to 4 cm depth (Table 2). Diffusive oxygen flux was elevated at all seep sites compared with the non-seep sites and reached a maximum value of 47 mmol m−2 d−1 at the AMV_seep_1 site. Both Pockmark sites had the highest benthic TOU (153 and 228 mmol m−2 d−1) detected in this study. All other investigated seep sites had more variable TOU (5–119 mmol m−2 d−1; Table 2).

Integrated (0–10 cm) SR rates at seep sites showed similar patterns as TOU, with the highest rates detected in the Pockmark region (20 and 50 mmol m−2 d−1), while the other sites had more variable sulfate consumption (1–35 mmol m−2 d−1; Table 2). Concentrations of total sulfide were elevated (max. 24 mmol l−1) at all seeps and increased with sediment depth, except at two of the Amon MV sites (AMV_seep_2 and AMV_seep_4) where no free sulfide was detected throughout the complete investigated sediment depth (Figure 3a, Supplementary Figure S4). At all other seep sites, sulfate concentration profiles steeply decreased with depth, and sulfate was often completely consumed (Table 2). Alkalinity and DIC concentrations were substantially elevated (max. 33–37 mmol l−1) and increased with depth, except at the AMV_seep_2 and AMV_seep_4 where background values (2.4 mmol l−1) were measured through the complete core depth (Figure 3b, Supplementary Figure S5). Ammonium concentrations varied substantially between different sites and mainly increased with sediment depth, to reach a maximum value of 4.4 mmol l−1 at the AMV_seep_3 (Supplementary Figure S7). Phosphate concentrations (average 3 μmol l−1) were similar at all seeps, exceeding those of non-seep sites. All of the seeps, except AMV_seep_2, had higher bacterial cell numbers compared with the non-seep sites, with maximum values of 0.3–6.7 × 109 cm−3 in the topmost 5-cm sediment depth. Differences in the cell numbers between non-seeps and seeps (excluding the newly formed AMV_seep_2 site) were statistically significant (Mann–Whitney W=30, P=0.01). Cell numbers at all seeps declined with sediment depth to reach 0.1 × 109 cm−3 at 15 cm (Supplementary Figure S8).

Comparison of bacterial communities between seep and non-seep regions

The Shannon diversity index did not differ between seep regions (5.3±0.4 s.d., n=8) and non-seep regions of the Eastern Mediterranean Sea (5.5±0.5 s.d., n=3; Mann–Whitney W=7, P-value=0.4; Table 3). Within the seep regions, observed ARISA OTU numbers were also not different between methane-seeping sites (241±31 s.d., n=9) and non-seeping sites (212±55, n=5; Mann–Whitney W=28, P=0.5).

Pairwise comparison of methane-seeping sites and the corresponding nearby non-seeping sites (0.01–1.8 km apart) revealed that, on average, between all pairs only 36% of the ARISA OTUs were shared (Supplementary Table 5a). Overall, seep and non-seep sites shared only 18% of a total of 158 ARISA OTUs and had communities with significantly different structures (analysis of similarity R=0.8, Bonferroni's P=0.001; Figure 4).

Figure 4

Non-metric multidimensional scaling analysis depicting differences in the bacterial community structure between seep (black symbols) and non-seep sites (red symbols) of three cold seep regions, as derived from ARISA-based data. Convex hulls depict significant differences between the groups, as determined by analysis of similarity (R=0.8, Bonferroni's P=0.001). Pooled (0–5 cm) samples per site were used to calculate Bray–Curtis dissimilarity matrix.

Pooled 454 sequences of all seep sites were dominated by Deltaproteobacteria (26%), Gammaproteobacteria (17%), Anaerolineae (11%), Epsilonproteobacteria (7%) and Caldilineae (7%) (Table 4). In contrast, most abundant taxa of Eastern Mediterranean deep-sea sediments from non-seep regions were dominated by Betaproteobacteria (17%), Gammaproteobacteria (14%), Synergistia (12%), Alphaproteobacteria (8%) and Flavobacteria (8%). Furthermore, the relative sequence abundance of the Epsilonproteobacteria, encompassing many types of microbes that gain energy via oxidation of sulfide, differed between seep (7% of all sequences) and non-seep regions (1%) (Table 4). Among the three most abundant OTU0.03 in the whole seep data set, two belonged to Epsilonproteobacteria and were closely related to sulfide-oxidizing Sulfurovum. Putative sulfur reducers and oxidizers of the Desulfobacteraceae, Desulfobulbaceae, Thiotrichales and Acidithiobacillales, along with aerobic methane oxidizers of the Methylococcus, were the most abundant microorganisms within the Deltaproteobacteria and Gammaproteobacteria at all seep sites (Supplementary Table S3). At non-seep regions, facultative anaerobic heterotrophs of the Alteromonadales and Oceanospirillales, as well as Desulfobacteraceae and Desulfobulbaceae, were the most abundant taxa within the Gammaproteobacteria and Deltaproteobacteria, respectively.

Table 4

Ten most sequence-abundant bacterial classes

AMV_seep_1	AMV_seep_3	AMV_seep_4	AmsMV_seep_1	AmsMV_seep_2	AmsMV_seep_3
Deltaproteobacteria (36)	Deltaproteobacteria (36)	Deltaproteobacteria (22)	Deltaproteobacteria (30)	Deltaproteobacteria (27)	Deltaproteobacteria (30)
Epsilonproteobacteria (16)	Anaerolineae (17)	Gammaproteobacteria (18)	Gammaproteobacteria (16)	Gammaproteobacteria (18)	Gammaproteobacteria (16)
Anaerolineae (11)	Caldilineae (12)	Anaerolineae (16)	Anaerolineae (15)	Anaerolineae (8)	Epsilonproteobacteria (10)
Gammaproteobacteria (11)	Gammaproteobacteria (6)	Caldilineae (10)	Sphingobacteria (5)	Epsilonproteobacteria (7)	Anaerolineae (7)
Sphingobacteria (5)	Epsilonproteobacteria (4)	Epsilonproteobacteria (5)	Caldilineae (4)	Milano-WF1B-44 (7)	Flavobacteria (5)
Holophagae (2)	Phycisphaerae (2)	Holophagae (3)	Milano-WF1B-44 (4)	Sphingobacteria (5)	Sphingobacteria (4)
Clostridia (2)	Sphingobacteria (2)	Milano-WF1B-44 (2)	Holophagae (2)	Flavobacteria (5)	Milano-WF1B-44 (4)
Caldilineae (2)	Clostridia (1)	Actinobacteria (2)	Actinobacteria (2)	Caldilineae (4)	Clostridia (4)
Actinobacteria (1)	Actinobacteria (1)	Sphingobacteria (2)	Epsilonproteobacteria (2)	Holophagae (3)	Caldilineae (2)
Flavobacteria (1)	Holophagae (1)	Flavobacteria (1)	Clostridia (2)	Actinobacteria (2)	Holophagae (2)

The relative percentage of sequence abundance for each seep site is given in parentheses. For comparison, three samples from non-seep regions of the deep Eastern Mediterranean sea were included.

Bacterial communities of cold seeps and factors shaping their structure

Among the seep sites, the richness estimates varied twofold (2793–6164 Chao1 index). Highest variation in the estimated richness existed between seep sites of the Amon MV (Table 3). Samples of this MV had significantly higher estimated richness than the Pockmark and Amsterdam MV samples (Mann–Whitney W=15, P=0.04. Similar patterns were obtained when using ARISA OTUs instead (Supplementary Figure S9).

The majority of OTU0.03 (32%) occurred at only 1–2 sites, and only 92 of the total 8117 OTU0.03 (that is, 4% of all types, excluding singletons) were common to all seep sites and comprised the core seep community with 59% of all sequences (Supplementary Figure S10). Most of the core OTU0.03 were affiliated to Deltaproteobacteria (23%), Caldilineae (14%), Gammaproteobacteria (11%), Anaerolineae (11%), Sphingobacteria (8%) and Epsilonproteobacteria (5%) (Supplementary Table S4). The three different cold seep regions had significantly different community structures (analysis of similarity R=0.3–0.6, Bonferroni's P⩽0.003, Figure 5c), and only the AMV_seep_1 showed higher similarity to the Amsterdam MV samples (Figures 5a and b). Indicator species analysis revealed that most indicator OTU0.03 of Amon and Amsterdam MV belonged to the class of Deltaproteobacteria and Gammaproteobacteria. The Pockmark sites had no unique OTU0.03 in the Gammaproteobacteria but in the Anaerolineae, Epsilon- and Deltaproteobacteria.

Figure 5

Non-metric multidimensional scaling analysis plots (based on Bray–Curtis dissimilarity) depicting differences in the bacterial community structure between different seep regions, as calculated from depth-pooled ARISA (a) and 454 MPTS (b), as well as depth-unpooled ARISA data (c). Seep sites are colored according to the seep region: Amsterdam MV sites in red, Amon MV sites in black, and Pockmark sites in green. Different sites within the same seep region are depicted with different symbols. Convex hulls in (c) depict significant differences between the groups, as determined by analysis of similarity (R=0.3–0.6, Bonferroni's P⩽0.003).

Sample diversity of seep habitats was inversely related to seepage activity, as revealed by the significant negative correlation between CH4 effluxes and Chao1 and ACE estimated richness (Spearman's coefficient R=−0.9, P⩽0.01). No significant correlation between beta-diversity (full 454 data set or only Deltaproteobacteria sequences) with seep activity, defined as SR, TOU or CH4 efflux, could be observed (Mantel test R<0.05, P>0.4). On the other hand, beta-diversity based on ARISA (data set, including all sites and sediment depth layers) was significantly and positively correlated with differences in porewater concentrations of sulfide, DIC and sulfate between seep sites (R=0.5–0.3, Bonferroni's P=0.01–0.04; Figure 6).

Figure 6

Plot depicting positive relation between beta-diversity (derived from ARISA) and differences in sulfate (a), sulfide (b) and DIC concentrations (c). Solid line represents a LOWESS curve (locally weighted scatterplot smoothing).

Spatial variations of seep-influenced bacterial communities

Pairwise analyses based on shared ARISA or 454 OTUs revealed decreasing similarity between samples with increasing distance for seep and non-seep sites within the scale of 0–350 km investigated here. On the smallest investigated spatial scale (<10 cm; within site comparison), a pair of samples shared, on average, 76% OTUs for seep (±4% s.d., n=9) and non-seep sites (±8% s.d., n=5). The sharpest decrease in the percentage of shared ARISA OTUs (54±8% s.d., n=3) was evident at intermediate spatial scales (20–630 m; between sites within one region), after which on the largest spatial scales (100 s of km; between cold seep regions) the trend leveled off at 51% shared OTUs for non-seep (n=2) and 52% (± 3% s.d.; n=3) for seep samples (Supplementary Table S5a). Pairwise comparison based on the 454 MPTS data set revealed a similar pattern but with lower values (Supplementary Table S5b). Accordingly, species accumulation analyses showed a linear increase, without reaching a plateau, in OTU0.03 diversity with increasing number of seep samples considered (Figure 7).

Figure 7

OTU accumulation curves, based on the 454 MPTS (a) and ARISA (b) data. Arrow indicates the number of seep sites needed to recover 95% of the observed OTUs.

Discussion

Cold seeps as islands with distinct microbial communities at the deep-sea floor

Previous investigations have indicated that cold seeps host higher biomass and lower diversities of animals than non-seep sediments (reviewed in Levin, 2005). As to microbial diversity, it was suggested that extreme environments in general host less diverse microbial communities than more stable, homogenized environments (Frontier, 1985). Results of this study reveal that cold seep habitats have similar levels of bacterial diversity (determined as Shannon index; Table 3) and host significantly higher cell numbers compared with typical oxygenated sediments in the oligotrophic deep Eastern Mediterranean. In line with this finding, but using other methods, it was previously shown that other cold seeps in the Eastern Mediterranean harbor similar bacterial diversity to background sediments with no sulfide production (Heijs ).

Investigation of the whole community structure using high-resolution fingerprinting techniques revealed that communities of seep sediments in the Eastern Mediterranean deep-sea differed also in their composition and relative abundance compared with adjacent and distant non-seep sediments. This expands the findings of previous studies, which have identified differences in the presence and proportion of key functional groups of the Deltaproteobacteria, Gammaproteobacteria and Epsilonproteobacteria involved in aerobic and anaerobic hydrocarbon degradation, as well as sulfur cycling (Orphan ; Knittel ; Heijs ; Orcutt ; Ruff ). Differences were also found at the level of the whole community between seep and non-seep regions, for example, by the high sequence abundances of the chemo-organotrophic Anaerolineae and Caldilineae of the Chloroflexi at seeps and their virtual absence in non-seep sediments (Table 4; Polymenakou ). In line with our results, it was previously found in the Pacific that deep-biosphere hydrate-bearing sediments harbor distinct bacterial communities compared with hydrate-free sites (Inagaki ). Hence, at the level of beta-diversity, bacterial community patterns match those of seep animals (Sahling ; Levin , 2010; Olu ; van Gaever ; Menot ; Ritt ), in that their community structures differ significantly from background communities colonizing deep-sea sediments.

Seep sites harbor distinct bacterial communities when compared with nearby non-seep sites (<1.8 km), sharing 22–45% of their OTUs with neighboring habitats, and 18% over the complete ARISA data set (Supplementary Table S5). In comparison, at other seep sites at the northeastern Pacific margin macrofauna communities were found to share a much higher proportion of species between seep and non-seep sites (Levin ). High adaptation of communities to their respective environments and different types of energy sources, as well as the variable tolerance to toxic conditions caused by sulfide production at cold seeps, are regarded as main factors selecting for different macrofaunal communities at seep and non-seep sites (Sibuet and Olu, 1998; Levin ). The same proposed factors and mechanisms could be also responsible for the differences in the bacterial communities of seeps and non-seeps observed here, but further studies are necessary to determine this relationship.

Factors shaping the bacterial diversity at cold seep ecosystems

Cold seeps are highly heterogeneous ecosystems (Cordes ), with patchy emission pathways of the main sources of energy—methane, other hydrocarbons and sulfide (Barry ; de Beer ; Omoregie ; Treude ; Pop Ristova ). Accordingly, only few bacterial groups were found in common to all the investigated cold seeps, of which sulfate reducers belonging to the Desulfobacteraceae and Desulfobulbaceae, as well as sulfide oxidizers Helicobacteraceae, were the most abundant and diverse taxa (Supplementary Table S4). Hence we propose that microorganisms affiliated to these families can be regarded as bacterial indicators of seep ecosystems, with overlaps to whale and wood falls (Bienhold ). Previous studies have identified differences in the distribution of key functional groups of the Deltaproteobacteria, Gammaproteobacteriaand Epsilonproteobacteria involved in hydrocarbon degradation and sulfur cycling in relation to environmental features (Orphan ; Knittel ; Heijs ; Orcutt ; Ruff ). At the level of the bacterial community structure, patterns could be related to contrasting local geochemical conditions within one seep structure (Omoregie ; Grünke ; Pop Ristova ). In this study, bacterial communities were compared across a range of different biogeochemical activities related to seepage (Table 2). Bacterial communities exhibited a patchy distribution from local to regional scales and increasing dissimilarity with increasing difference in seep sediment geochemistry (Figures 5 and 6). Similarity among seep sites was high at the level of bacterial classes, but seep habitats differed substantially regarding their relative abundance and composition at the OTU0.03 level (Table 4). The highly dynamic Amon MV system (Felden ) with very different types of seep settings—that is, gas versus brine driven fluid flow—showed significantly higher bacterial richness and variability compared with the Pockmark and Amsterdam MV sites (Table 3, Figures 5a–c). Underlying seep geochemistry was responsible for shaping both the alpha- and beta-diversity of cold seeps. Positive relationships between richness and energy availability have been revealed for other marine ecosystems (Pommier ; Fuhrman ; Bienhold ). A negative trend revealed here between methane efflux and microbial richness in seeps may be explained by the finding that the highest methane efflux rates came from relatively recent seep sites from within a disturbed mud volcano center (Felden ). Non-metric multidimensional scaling analyses based on a subset of the community, including only classes such as Deltaproteobacteria, Gammaproteobacteria, Epsilonproteobacteria, Caldilineae and Anaerolineae, revealed highly similar patterns correlated to those based on the whole bacterial community (Supplementary Table S6). This implies that, in terms of diversification, not only key functional groups involved in hydrocarbon degradation such as methane-oxidizing, sulfur-reducing and -oxidizing microorganisms but also other bacterial types such as the Chloroflexi may respond similarly to geochemical variations at cold seep ecosystems. Members of the Chloroflexi have been found at cold seeps, deep hydrocarbon-bearing sediments and, in general, organic-rich sediments (Kormas ; Inagaki ; Pachiadaki ), where they comprise high portion of the total bacterial community (30%, this study). The role of this bacterial group at cold seep ecosystems is largely unknown and should be addressed in future studies.

Effect of spatial scales on bacterial diversity of cold seep ecosystems

Knowledge on the biodiversity patterns of deep-sea organisms across space and time is important to anticipate responses of deep-sea ecosystems to current and future environmental and anthropogenic changes (Danovaro , 2010; Levin and Sibuet, 2012). Results of this study show that spatial patterns of bacterial communities of cold seeps comply with the taxa–area relationship, that is, an increase in species richness with increasing size of sampled area (Figure 7). This pattern has been reported for all domains of life and across different habitats, including meiofauna and macrofauna seep organisms (Cordes ), and hence is thought to be one of the few universal laws in ecology (Rosenzweig, 1995; Lawton, 1999; Zinger ). A general trend of increasing dissimilarities in the bacterial communities with geographic distance has been identified across various ecosystems on large spatial scales (Papke ; Whitaker ; Martiny ; Ramette and Tiedje, 2007; Schauer ; Zinger ). In our study, the percentage of shared ARISA bacterial types was the highest within sites (76%), substantially declined between sites (54%) and remained at the same level between different cold seep regions (52%) (Supplementary Table S5). Hence, highest spatial turnover occurred on intermediate spatial scales (<1 km), within individual seep regions, among locations spatially separated by few meters to hundreds of meters. With approximately one-third of the microbial diversity being unique to the individual seep habitats, that is, patches of reduced sediments of only a few meters in diameter, these localities significantly contribute to the biodiversity of deep-sea, explaining why cold seep systems can be regarded as biodiversity hotspots. The small, reduced habitats may have a critical role for the connectivity and diversity of margin communities and might provide resilience to perturbations caused by increased anthropogenic impact. This is a relevant finding with regard to spatial management of deep-water ecosystems, such as for deep-sea oil and gas resources; however, this also suggests that cold seep communities could be vulnerable to habitat losses at the scale of the individual seep ecosystem (Ramirez-Llodra ; van Dover ).