Quantification of the effects of ocean acidification on sediment microbial communities in the environment: the importance of ecosystem approaches.

To understand how ocean acidification (OA) influences sediment microbial communities, naturally CO2-rich sites are increasingly being used as OA analogues. However, the characterization of these naturally CO2-rich sites is often limited to OA-related variables, neglecting additional environmental variables that may confound OA effects. Here, we used an extensive array of sediment and bottom water parameters to evaluate pH effects on sediment microbial communities at hydrothermal CO2 seeps in Papua New Guinea. The geochemical composition of the sediment pore water showed variations in the hydrothermal signature at seep sites with comparable pH, allowing the identification of sites that may better represent future OA scenarios. At these sites, we detected a 60% shift in the microbial community composition compared with reference sites, mostly related to increases in Chloroflexi sequences. pH was among the factors significantly, yet not mainly, explaining changes in microbial community composition. pH variation may therefore often not be the primary cause of microbial changes when sampling is done along complex environmental gradients. Thus, we recommend an ecosystem approach when assessing OA effects on sediment microbial communities under natural conditions. This will enable a more reliable quantification of OA effects via a reduction of potential confounding effects. © FEMS 2016.

INTRODUCTION

The increase in atmospheric CO2 concentrations has led to a continuous increase in the partial pressure of CO2 and decrease in pH in the oceans since pre-industrial times, a process called ocean acidification (OA). Since then, the ocean water pH has declined from approximately 8.2 to 8.1, and will reach pH 7.8 by the year 2100 if we continue on the current and predicted CO2 emissions trajectories (IPCC 2013). Such changes in the carbonate chemistry can have dramatic effects on marine organisms, e.g. reducing biogenic calcification (Hofmann et al. 2010), modifying fish sensory perception (Munday et al. 2014), and enhancing seagrass and macroalgal growth (Koch et al. 2013). Marine microbial populations are also expected to respond to OA. As key players in nutrient cycling and remineralization of organic matter, changes in marine microbial communities and in the services they provide have the potential for far-reaching consequences (Liu et al. 2010; Joint, Doney and Karl 2011).

Previous work on the effects of OA on marine microbes focused mostly on planktonic microbes, and has so far reported rather variable and inconsistent outcomes. Several incubation experiments documented a change in bacterial community composition, involving e.g. Gammaproteobacteria and Flavobacteria, as well as function under decreased pH conditions, e.g. increased nitrogen fixation rates and decreased nitrification (Beman et al. 2011; Kitidis et al. 2011; Krause et al. 2012; Lomas et al. 2012). There is also evidence that OA may increase microbial carbon degradation rates by increasing enzyme activity (Piontek et al. 2010). On the other hand, in many mesocosm studies the composition of the planktonic microbial community was mostly stable over varying pH levels without changes in biogeochemical functions (Newbold et al. 2012; Lindh et al. 2013; Roy et al. 2013; Oliver et al. 2014). Similarly divergent results were obtained from incubation experiments with microbial biofilms or marine sediments, which predicted either rapid responses of the microbial community to OA (Witt et al. 2011; Laverock et al. 2013; Tait, Laverock and Widdicombe 2013; Braeckman et al. 2014), or only minor impacts of OA on community composition and function (Tait and Laverock 2013; Gazeau et al. 2014).

To venture beyond the limited scope of incubation experiments and to assess long-term OA effects on microbial community composition and functions in their natural environment, the usage of naturally CO2-rich sites as OA analogues has become increasingly popular. Previous observations on microbial communities in sediments at naturally CO2-rich sites focused mostly on sites in the Mediterranean Sea and in Papua New Guinea (Kerfahi et al. 2014; Taylor et al. 2014; Raulf et al. 2015). Pronounced changes in microbial richness and community composition were observed along acidification gradients. However, studies disagreed on the direction of the change in microbial richness: Kerfahi et al. (2014) and Raulf et al. (2015) detected an increase, whereas Taylor et al. (2014) observed a decrease of microbial richness with decreasing pH. Furthermore, reports on the prevalence of dominant microbial taxa in marine sediments (i.e. Gammaproteobacteria, Alphaproteobacteria, Bacteroidetes) under acidified conditions were inconsistent, with evidence pointing towards both increases and decreases as well as no changes in relative abundances (Kerfahi et al. 2014; Taylor et al. 2014; Raulf et al. 2015). So far these contrasting results have not been reconciled.

Many naturally CO2-rich sites are associated with hydrothermal activity, which is the primary source of CO2 at these sites (Hall-Spencer et al. 2008). However, the hydrothermal character often introduces confounding factors that make it difficult to specifically assess the impact of high pCO2 and reduced pH (Vizzini et al. 2013). For instance, most CO2-rich sites exhibit increased levels of methane, sulfide, temperature and various trace elements (Wenzhöfer et al. 2000; Meyer-Dombard et al. 2012; Vizzini et al. 2013; Burrell et al. 2015). Especially the sediments at naturally CO2-rich sites may be strongly affected as they are influenced from above by the acidified water as well as from below by the hydrothermal fluids (Wenzhöfer et al. 2000). The increased temperatures and altered chemical composition of hydrothermal fluids compared with seawater can change pore water geochemistry and affect geochemical gradients in the sediment (Wenzhöfer et al. 2000; German and von Damm 2003). As such, sediments at naturally CO2-rich sites constitute a highly complex system with a multitude of environmental parameters that shape microbial habitats and influence microbial community composition and associated functions.

In general, research on microbial communities at shallow-water hydrothermal seeps, including naturally CO2-rich sites that were used as OA analogues, is very limited (Giovannelli et al. 2013) and often environmental measurements are not paired spatially and temporally with microbiological characterizations. Apart from recent work by Molari et al. (unpublished), previous conclusions were based on a limited number of samples and recorded environmental parameters (Kerfahi et al. 2014; Taylor et al. 2014; Raulf et al. 2015). In several instances, water column parameters were used to describe microbial communities without comprehensively assessing the environmental conditions in the sediment, which the microbes were subjected to (Kitidis et al. 2011; Kerfahi et al. 2014; Raulf et al. 2015). Insufficient environmental characterization and subsequent poor selection of naturally CO2-rich sites as OA analogues may constitute a potential cause for the diverging results regarding the impact of OA on sediment microbial communities.

Here, we quantified the impact of OA on sediment microbial communities at hydrothermal CO2 seeps, while taking the complex environment associated with the CO2 seepage into account. For this, we investigated natural shallow-water CO2 seeps in Papua New Guinea, which have been extensively used in OA research (Fabricius et al. 2011; Morrow et al. 2014; Raulf et al. 2015). We used molecular community fingerprinting (Automated Ribosomal Intergenic Spacer Analysis) of more than 100 samples as well as next generation amplicon sequencing of the 16S rRNA gene to assess bacterial and archaeal community composition. We accompanied the molecular analyses with a comprehensive characterization of environmental conditions, including bottom and pore water chemistry, sediment permeability as well as carbon and nitrogen content. Through these analyses we further identified microbial taxa that may constitute ‘losers’ and ‘winners’ under the conditions at the CO2 seeps.

MATERIALS AND METHODS

Sampling

Sediment samples were collected at two CO2 seeps at Upa Upasina, Normanby Island (Reef 1: S 9.82, W 150.82) and Dobu Island (Reef 2: S 9.74, W 150.86), Papua New Guinea (Supplementary Fig. S1). Both reefs were characterized by a pH gradient created by hydrothermal CO2 seepage and have been previously used as OA analogues to study sediment microbial communities (Raulf et al. 2015). Along this pH gradient, 14 sites were sampled at Reef 1 and six sites at Reef 2. At each site three independent replicate cores (diameter: 2.5 cm) were taken and the sediment of the first 2 cm (0–2 cm) and the next lower 2 cm (2–4 cm) was preserved in RNAlater Solution (Ambion) for molecular analysis.

To characterize the environmental conditions in the sediment along the pH gradient at Reef 1 the following parameters were measured: in situ oxygen concentrations, temperature, redox potential and pH using microsensor profiles (temperature sensor: UST Umweltsensortechnik GmbH, Gschwenda, Germany; pH sensor: Microelectrodes Inc., Bedford, NH, USA; http://doi.pangaea.de/10.1594/PANGAEA.858091); total organic and inorganic carbon and total nitrogen content of the sediment; grain size, porosity and permeability (http://doi.pangaea.de/10.1594/PANGAEA.858091); and pore water geochemistry analysed from 6 ml pore water (http://doi.pangaea.de/10.1594/PANGAEA.858033). These parameters will be referred to as ‘sediment parameters’. Additionally, bottom water samples were collected approximately 5 cm above each of the sediment cores at Reef 1 and 2 to measure carbonate chemistry (pH, dissolved inorganic carbon (DIC), total alkalinity (TA)) and nutrient concentrations (SiO4, PO4, NOx, NH4). Data on bottom water carbonate chemistry and nutrient concentrations will further be referred to as ‘bottom water parameters’ and are available in the Pangaea database (http://doi.pangaea.de/10.1594/PANGAEA.854018). An overview of the sampling design and the measured parameters is provided in Supplementary Table S1.

Automated Ribosomal Intergenic Spacer Analysis

From each sample DNA was extracted from 1 g of sediment using the UltraClean Soil DNA extraction kit according to the manufacturer's instructions (MoBio Laboratories Inc., Carlsbad, CA, USA). The DNA was eluted in TE buffer (100 mmol L−1 Tris-HCl, 10 m mmol L−1 EDTA) and quantified photometrically with a NanoQuant (Tecan, Crailsheim, Germany). To screen for changes in microbial community structure we used the community fingerprinting technique Automated Ribosomal Intergenic Spacer Analysis (ARISA) with universal bacterial primers (Fisher and Triplett 1999). As previously described in Ramette (2009), triplicate reactions of the ARISA-PCR were run for each sample in an Eppendorf MasterCycler using the PeqLab PCR kit (PeqLab Biotechnology GmbH, Erlangen, Germany). Each PCR reaction contained 10–15 ng DNA, 1× buffer S, 0.25 mmol L−1 dNTPs, 0.1 g L−1 bovine serum albumin (BSA), 0.4 μmol L−1 fluorescently labeled forward primer (ITSF: 5′-GTCGTAACAAGGTAGCCGTA-3′) and reverse primer (ITSReub: 5′-GCCAAGGCATCCACC-3′), an additional 1 mmol L−1 MgCl2 and 0.05 U μL−1 Taq polymerase in a total reaction volume of 25 μL. PCR conditions were 3 min at 94°C followed by 30 cycles of 45 s denaturation at 94°C, 45 s annealing at 55°C, and 90 s elongation at 72°C, with a final elongation step for 5 min at 72°C. Fragment sizes were determined on a capillary sequencer. Fragments between 100 base pairs (bp) and 1000 bp were binned into operational taxonomic units (OTUs) using custom R scripts available at http://www.mpi-bremen.de/en/Software_4.html.

Amplicon sequencing and sequence processing

To taxonomically classify the microbial community, we sequenced the hypervariable regions V3–V4 and V4–V6 of the bacterial and archaeal 16S rRNA gene, respectively, using the bacterial primers S-D-Bact-0341-b-S-17 (5′-CCTACGGGNGGCWGCAG-3′) and S-D-Bact-0785-a-A-21 (5′-GACTACHVGGGTATCTAATCC-3′), and the archaeal primers Arch349F (5′-GYGCASCAGKCGMGAAW-3′) and Arch915R (5′-GTGCTCCCCCGCCAATTCCT-3′; Amann et al. 1990; Klindworth et al. 2013). The community of the upper sediment layer of one replicate core was sequenced from 13 sampling sites at Reef 1. Sequences were generated on the Illumina MiSeq platform (CeBiTec Bielefeld, Germany), in a 2 × 300 bp paired-end run. Primer sequences were removed from the raw paired-end reads with cutadapt (Martin 2011). The primer-trimmed sequences are available on ENA (PRJEB11384). Sequences were quality trimmed with a sliding window of 4 bases and a minimum average quality of 15 with trimmomatic v0.32 (Bolger, Lohse and Usadel 2014), merged with PEAR v0.9.5 (Zhang et al. 2014), clustered into OTUs using swarm v2.0 (Mahé et al. 2014) and taxonomically classified with SINA (SILVA Incremental Aligner) v1.2.10 using the SILVA rRNA project reference database (release 119) at a minimum alignment similarity of 0.9 and a last common ancestor consensus of 0.7 (Pruesse, Peplies and Glöckner 2012). OTUs that were unclassified on domain level as well as those matching chloroplast and mitochondrial sequences were excluded from the analysis. The final OTU tables are accessible at http://doi.pangaea.de/10.1594/PANGAEA.854018. Throughout this article taxon names are used to designate sequence affiliation to the respective taxon and ‘abundance’ of bacterial and archaeal taxa refers to the sequence abundance of Illumina amplicons.

Statistical analysis

Principal component analysis (PCA) was used to classify the sampling sites at Reef 1 into categories of hydrothermal influence based on observed sediment parameters. To assess the covariation among environmental parameters, pairwise Pearson correlation coefficients were calculated.

α-Diversity indices were calculated based on repeated random subsampling of the amplicon data sets to assess richness and evenness of the microbial community, namely OTU number, Chao1 and abundance-based coverage estimator, inverse Simpson, percentage of absolute (occurring only once in the complete data set) and relative singletons (occurring only once in the sample) as well as absolute doubletons (occurring twice in only one sample in the complete data set). Significant differences in α-diversity indices between hydrothermal influence categories were determined by analysis of variance (ANOVA) using permutation tests. P-values of subsequent pairwise tests were adjusted using the Benjamini–Hochberg (BH; Benjamini and Hochberg 1995) correction procedure at the significance level of P = 0.05.

The change in community structure (β-diversity) between samples was quantified by calculating Bray–Curtis and Jaccard dissimilarity from relative OTU abundances. The former was used to produce non-metric multidimensional scaling (NMDS) plots and to test for community similarity between hydrothermal influence categories (ANOSIM tests); the latter was used to calculate the number of shared OTUs between samples. For the amplicon data sets, we excluded the rare biosphere by retaining only those OTUs that were present with more than two sequences in more than 10% of the samples. This reduction of the data sets did not change β-diversity patterns (Mantel test, r > 0.9, P < 0.001).

The contribution of environmental parameters to explaining the variation in community structure was calculated using redundancy analysis (RDA) and variation partitioning of centered log-transformed relative OTU abundances. Prior to significance testing, parameters were excluded using forward model selection until the minimum Akaike Information Criterion (AIC) value was reached. Of highly correlated parameters only one parameter was kept in the final models. Differentially abundant OTUs were detected using the R package ALDEx2 (Fernandes et al. 2014) at a significance threshold of 0.05 for BH-adjusted parametric and non-adjusted non-parametric P-values. The significance threshold for planned parametric post-hoc tests was 0.1 (BH-adjusted).

All statistical analyses were conducted in R using the core distribution with the additional packages vegan (Oksanen et al. 2015), compositions (Van den Boogaart, Tolosana and Bren 2014), ALDEx2 (Fernandes et al. 2014) and FactoMineR (Husson et al. 2015).

RESULTS

Environmental conditions at sampling sites

The bottom water pH gradients created by the CO2 seeps at Normanby and Dobu Island ranged from approximately 6.8 to 8.3 at both Reefs (Table 1). In addition to the changes in pH, bottom water silicate concentrations increased with decreasing distance to the main CO2 seepage area from approximately 4 μmol L−1 at reference sites to more than 10–50 μmol L−1 at the seep sites. Reef 2 was further characterized by the presence of microbial mats in close proximity to the CO2 seeps and by a strong sulfidic smell emanating from the CO2 seeps. The presence of sulfide at the CO2 seeps at Reef 2 convinced us to focus our sampling effort on Reef 1 as a more likely OA analogue.

Table 1.

Environmental conditions at the sampling sites on Normanby and Dobu Island, Papua New Guinea. Values are given as mean ± standard error. For pH, oxygen concentration and temperature microprofiles mean values per hydrothermal influence category were calculated based on median values per sediment layer (0–2 cm, 2–4 cm).

		Reef 1			Reef 2
	Sediment layer	Reference	Medium	High	Reference	Medium	High
Water depth (m)		3.87 ± 0.08	4.30 ± 0.48	3.74 ± 0.56	3.50 ± 0.00	2.50 ± 0.08	2.30 ± 0.23
pH	0–2 cm	7.86 ± 0.23	7.15 ± 0.24	6.67 ± 0.21
	2–4 cm	7.75 ± 0.20	6.64 ± 0.40	6.53 ± 0.31
O2 (μmol L−1)	0–2 cm	16.28 ± 2.25	5.78 ± 4.70	0.41 ± 0.33
	2–4 cm	0.17 ± 0.42	0.12 ± 0.08	0
T (°C)	0–2 cm	29.80 ± 0.35	30.26 ± 0.14	30.49 ± 1.55
	2–4 cm	29.86 ± 0.38	30.69 ± 0.20	31.56 ± 2.57
Porosity	0–2 cm	0.49 ± 0.01	0.46 ± 0.02	0.47 ± 0.01
	2–4 cm	0.45 ± 0.01	0.44 ± 0.02	0.45 ±0.03
Permeability (m2)	0–4 cm	(8.24 ± 0.81) × 10–11	(3.57 ± 0.87) × 10–11	(2.96 ± 0.02) × 10–11
pH	Bottom water	8.24 ± 0.02	7.83 ± 0.08	7.53 ± 0.13	8.33 ± 0.00	7.56 ± 0.05	6.78 ± 0.02
TA (mmol L−1)	Bottom water	2.51 ± 0.02	2.50 ± 0.02	2.60 ± 0.03	2.42 ± 0.02	2.47 ± 0.02	2.52 ± 0.03
SiO4 (μmol L−1)	Bottom water	3.06 ± 0.18	5.17 ± 0.97	14.92 ± 3.17	4.26 ± 0.04	20.06 ± 3.48	50.95 ± 4.90
PO4 (μmol L−1)	Bottom water	0.04 ± 0.01	0.07 ± 0.01	0.08 ± 0.01	1.16 ± 0.65	0.10 ± 0.02	1.69 ± 0.72
NO3 (μmol L−1)	Bottom water	0.28 ± 0.02	0.46 ± 0.08	0.48 ± 0.07	3.37 ± 1.59	0.50 ± 0.17	0.96 ± 0.29
NO2 (μmol L−1)	Bottom water	0.02 ± 0.01	0.02 ± 0.00	0.03 ± 0.01	0.03 ± 0.01	0.02 ± 0.01	0.05 ± 0.01
NH4 (μmol L−1)	Bottom water	3.15 ± 2.55	2.47 ± 0.73	1.27 ± 0.79	0.76 ± 0.25	2.56 ± 1.56	1.33 ± 0.20

At Reef 1, the median pore water pH in the first 2 cm of the sediment and the following 2 cm calculated from in situ microprofiles ranged from 8.1 at reference sites to 5.9 at seep sites. Generally, median sediment pH was lower than bottom water pH and declined more rapidly when moving towards the main CO2 seepage area. However, at a seagrass patch influenced by CO2 seepage this relationship was reversed with a bottom water pH (7.1) that was lower than the median pore water pH (7.8). With decreasing sediment pH, median permeability dropped from approximately 8 × 10–11 to 3 × 10–11 m2. Temperature profiles further showed increased temperatures at some of the seep sites, especially in the sediment layer from 2 to 4 cm, where maximum temperatures of more than 34°C were reached compared with approximately 30°C at the majority of the sites (Table 1). Pore water geochemistry of sediments with active fluid seepage showed increases in cations that are typically present in hydrothermal fluids, such as lithium, silicon and manganese, and decreases in cations that are typically depleted in hydrothermal fluids such as magnesium.

Sediment parameters were used to conduct a PCA to identify sampling sites at Reef 1 with similar characteristics based on all of the following parameters: pH, TA, oxygen concentration (O2), redox potential (redox), temperature, porosity, grain size, permeability, total nitrogen (TN), total organic carbon (TOC), total inorganic carbon (TIC), nutrients (NH4, NOx), anions (Cl, Br, SO4) and cations (Na, Li, Si, B, Mg, K, Fe, Mn, Sr, Rb, Ba, Ca, Cs). We identified three main clusters of sampling sites: reference and seep sites were separated along principal component 2 (PC2). The seep sites were further separated along PC1. The parameters contributing to PC1, which accounted for 40% of the variation in the data, were mostly concentrations of signature elements for hydrothermal activity such as lithium, silicon and manganese. Permeability, grain size, pH, nitrogen and organic carbon content of the sediment were among the parameters that contributed to PC2 (16% of the total variance in the data; Fig. 1 and Supplementary Fig. S2). Based on the pattern in the PCA the sampling sites were classified into three categories: (i) reference sites characterized by ambient pH, (ii) sites characterized by low pH and low hydrothermal influence, (iii) and sites characterized by low pH and high hydrothermal influence, including pronounced temperature increases. Sampling sites falling within these three hydrothermal influence (HI) categories will be referred to as ‘reference’, ‘medium HI’ and ‘high HI’ sites thereafter.

Figure 1.

Principal component analysis (PCA) of the environmental conditions in the sediment at the sampling site at Reef 1 on Normanby Island to classify the sampling sites according to hydrothermal influence. The data are available at Pangaea (http://doi.pangaea.de/10.1594/PANGAEA.858033, http://doi.pangaea.de/10.1594/PANGAEA.858091). For oxygen, pH, temperature and redox potential microprofiles, the PCA was calculated with median values for each sediment layer. Missing values were replaced by the mean of the respective parameter for the calculation of the PCA. The arrows show the loadings of the environmental parameters scaled to 4 times their value for better visualization. DIC, dissolved inorganic carbon; TA, total alkalinity; TIC, total inorganic carbon; TN, total nitrogen; TOC, total organic carbon.

Only sediment parameters where data were available for at least eight of the 14 sampling sites at Reef 1 were included in the RDA models to explain the variation in microbial community structure. Together with bottom water parameters and spatial information these were reef, position along the reef, water depth, bottom water silicate, phosphate, nitrate, nitrite, ammonium concentrations, TA and pH, sediment porosity, permeability, temperature, pH, redox potential, and oxygen concentrations. Pairwise linear correlations among environmental parameters that were used for the RDA models revealed a complex pattern of covariation. Several parameters were highly correlated, e.g. bottom water silicate concentrations and bottom water pH or permeability and the position of the sampling sites along the reef where absolute correlation coefficients were approximately 0.8. On the other hand, bottom water and sediment pH did not show any linear correlation (Supplementary Fig. S3).

Variation in microbial community structure

A total of 117 sediment samples were analysed with ARISA to identify overall changes in microbial community structure (Fig. 2). Microbial community structure showed distinct clusters associated with hydrothermal influence categories. Those clusters were confirmed to be significantly different from each other by ANOSIM (Supplementary Table S2). Whereas the microbial community structure of reference samples from both reefs was quite similar, there was a divergent trend at each of the two reefs from reference to medium to high HI sites. The microbial communities at the seep sites at Reef 1 and 2 were significantly different from each other (ANOSIM, R > 0.8, BH-adjusted P < 0.05). Furthermore, there was no apparent influence of sediment layer on the microbial community composition (Fig. 2). The average number of shared OTUs between samples from different hydrothermal influence categories at each reef was in most cases less than 40%. Generally more OTUs were shared between reference and medium HI sites or between medium and high HI sites compared with reference and high HI sites (Supplementary Table S2). Even among samples within hydrothermal influence categories, the microbial community was very heterogeneous with about 50% shared OTUs between any two samples. The patterns of change in bacterial and archaeal community composition in the first 2 cm of the sediment based on 16S sequences were very consistent with ARISA results (Supplementary Table S2). The number of shared 16S OTUs was slightly lower with on average 12–28% shared OTUs between hydrothermal influence categories and about 30–50% within categories.

Figure 2.

Non-metric multidimensional scaling (NDMS) plot based on the Bray–Curtis dissimilarity matrix of the microbial community based on automated ribosomal intergenic spacer analysis (ARISA). Larger symbols mark the subset of samples that was used in the redundancy analysis (RDA) models with sediment parameters (Table 2).

Table 2.

Contribution of observed environmental parameters to explaining the variation in microbial community structure based on redundancy analysis (RDA)-based variation partitioning. To compare the explanatory power of different sets of environmental parameters, several RDA models were tested: the complete ARISA data set was analysed using bottom water parameters and hydrothermal influence categories; a subset of the ARISA data was analysed using bottom water, sediment parameters and hydrothermal influence categories. The bacterial and archaeal amplicon data sets were analysed using hydrothermal influence categories.

Data set	Modela	Source of variation	Total R2 adjusted	Pure R2 adjusted	F	df	P-valueb	Covariationc
ARISA	Bottom water	All	0.277		5.126	9,88	<0.001
complete	(AIC = 574.1)	Reef	0.067	0.037	5.555	1,88	<0.001	0.030
		Position along reef	0.051	0.034	5.163	1,88	<0.001	0.017
		Water depth	0.049	0.034	5.149	1,88	<0.001	0.015
		SiO4	0.055	0.029	4.540	1,88	<0.001	0.026
		PO4	0.028	0.025	4.068	1,88	<0.001	0.003
		NO3	0.017	0.015	2.827	1,88	0.028	0.002
		NO2	0.018	0.008	1.962	1,88	0.272	0.010
		TA	0.033	0.015	2.900	1,88	0.025	0.017
		pH	0.041	0.029	4.547	1,88	<0.001	0.012
	Categories	All	0.198		9.789	3,104	<0.001
	(AIC = 638.5)	Reef	0.059	0.061	8.999	1,104	<0.001	−0.002
		Hydrothermal influence	0.137	0.139	10.185	2,104	<0.001	−0.002
ARISA	Bottom water	All	0.352		5.665	5,38	<0.001
subset	(AIC = 248.6)	Position along reef	0.080	0.056	4.378	1,38	0.005	0.024
		Water depth	0.068	0.059	4.531	1,38	<0.001	0.009
		SiO4	0.097	0.066	4.994	1,38	<0.001	0.031
		NO3	0.066	0.041	3.442	1,38	0.024	0.025
		pH	0.045	0.043	3.579	1,38	0.025	0.003
	Sediment	All	0.436		5.156	8,35	<0.001
	(AIC = 244.84)	Water depth	0.068	0.057	4.632	1,35	<0.001	0.011
		SiO4d	0.097	0.117	8.486	1,35	<0.001	−0.020
		Porosity	0.029	0.043	3.728	1,35	<0.001	−0.014
		Permeability	0.069	0.043	3.721	1,35	<0.001	0.026
		O2	0.016	0.012	1.763	1,35	0.002	0.004
		Redox potential	0.038	0.011	1.693	1,35	0.037	0.027
		Temperature	0.051	0.052	4.314	1,35	<0.001	−0.001
		pH	0.088	0.038	3.457	1,35	<0.001	0.050
	Categories (AIC = 251.8)	Hydrothermal influence	0.259		8.500	2,41	<0.001
16S Bacteria	Categories (AIC = 120.5)	Hydrothermal influence	0.375		4.606	2,10	<0.001
16S Archaea	Categories (AIC = 100.2)	Hydrothermal influence	0.267		3.188	2,10	<0.001

The Akaike Information Criterion (AIC) is given for each model as goodness-of-fit statistic.

Only the significance of the whole model (all) and the pure effects of the respective parameters (accounting for the effects of all other factors in the model) were tested. P-values were calculated based on restricted permutations.

Covariation constitutes the amount of variation that can be explained by more than the parameter of interest.

Bottom water silicate concentrations were used as proxy for pore water concentrations, because of a high correlation based on point measurements of selected samples.

To estimate and compare the contribution of different sets of environmental parameters to explaining the patterns in microbial community composition, the ARISA data were analysed with several RDA models: (i) the complete ARISA data set using bottom water parameters, (ii) a subset of the ARISA data set from Reef 1 where sediment parameters were available that was also analysed with (iii) bottom water parameters for comparability, (iv) the complete data set, and (v) the ARISA subset using solely hydrothermal influence categories as explanatory variable (Table 2).

All ARISA RDA models as well as their individual factors significantly explained variation in microbial community structure, with the exception of bottom water nitrite concentration in the complete ARISA model, which was not significant but retained in the model selection procedure. The total amount of explained variation (R2) varied from 20% to 44% between the RDA models based on the ARISA dataset. The models based only on hydrothermal influence categories explained approximately 10% less variation than the models based on numeric environmental parameters. The model with the ARISA subset from Reef 1 and bottom water parameters explained more variation than the model with the complete ARISA data set that also included samples from Reef 2 with an R2 of 35% compared with 28%. Including sediment parameters further increased the total amount of explained variation resulting in by far the best model with a R2 of 44% and the lowest AIC of the models with the ARISA subset. pH alone was able to explain 4.5% (total R2) in the bottom water models and almost 9% when pore water pH was used. Accounting for the variation explained by the other parameters in the models (pure R2), pH explained much less of the changes in community structure due to covariation with the other parameters in the model. Despite the high degree of covariation, the contribution of pH to explaining the microbial community composition was significant in all models (Table 2). The variation explained by the other factors in the models was similar to that of pH ranging from an R2 of 1.1 to 11.7%. Water depth, silicate concentrations, temperature and permeability were among the factors with the highest effects on community changes. For the amplicon data sets, only the significance of hydrothermal influence categories was tested due to the reduced number of samples. For both bacterial and archaeal communities, hydrothermal influence significantly explained the variation in community structure with an R2 of 38% and 27%, respectively (Table 2).

Taxonomic composition of microbial communities

The sediment at the CO2 seeps hosted a very diverse community: estimated species richness based on 16S OTUs ranged from 13 000 to 48 000 for bacteria and 475 to 4100 for archaea (Supplementary Fig. S4). Median Chao1 richness was highest in samples from reference sites with 42 000 for bacteria and 2700 for archaea followed by 35 000 and 2000 at medium HI sites and 14 000 and 960 at high HI sites, respectively. However, within each hydrothermal influence category, Chao1 estimates varied considerably over a range of more than 20 000 bacterial and 2000 archaeal species, without significant difference between groups. Bacterial and archaeal inverse Simpson index showed a weak decreasing trend from reference to medium and high HI sites. Both bacterial and archaeal communities consisted of a large proportion of rare OTUs, which made up approximately 40% (absolute singletons), 25% (relative singletons) and 10% (absolute doubletons) of all OTUs in each sample (Supplementary Fig. S4).

The bacterial community was dominated by Gamma- and Deltaprotobacteria with a total relative abundance of 18% and 17%, respectively (Fig. 3A). The next most abundant bacterial phyla were Bacteroidetes (9%), Actinobacteria (9%) and Chloroflexi (7%). Even at low taxonomic resolution levels, several bacterial taxa were identified as being differentially abundant between hydrothermal influence categories. The number of differentially abundant bacterial taxa increased as higher taxonomic resolution levels were considered (Supplementary Table S3). Generally between 20 and 70% of the sequences per sample belonged to differentially abundant bacterial taxa. In the next sections, we will focus on dominant bacterial taxa. The full list of the differentially abundant bacterial taxa is provided in Supplementary Table S4.

Figure 3.

Taxonomic composition of the 10 most abundant bacterial phyla (A) and archaeal genera (B) per sample at reference, medium and high HI (hydrothermal influence) sites. For Proteobacteria class-level resolution is shown; for taxa that were unclassified on the respective level of resolution, the next higher level classified taxonomic rank is shown. Asterisks mark significantly differentially abundant taxa between hydrothermal influence categories. SCG: Soil Crenarchaeotic Group, MHVG: Marine Hydrothermal Vent Group, DHVEG-6: Deep Sea Hydrothermal Vent Group 6.

Among the dominant differentially abundant phyla, Chloroflexi, Chlorobi and Deferribacteres increased at medium and high HI sites compared with reference sites. The changes in relative abundance within the Chlorobi and Chloroflexi were caused by several taxa within these phyla, i.e. of the classes Ignavibacteria for Chlorobi and Anaerolineae, Ardenticatenia and Caldilineae for Chloroflexi. The genus Caldithrix was mostly responsible for the increase of Deferribacteres towards the seep sites. The difference in the relative abundance of Chloroflexi was most apparent towards high HI sites and was not statistically significant between reference and medium HI sites. Chlorobi and Caldithrix also showed significant increases between reference and medium HI sites.

Although not detected as differentially abundant at the phylum level, Cyanobacteria, predominantly of the Subsections I and II were significantly less abundant at the seep sites than at reference sites. In all cases, this decrease was significant between reference and medium HI sites. The genus Pleurocapsa of Subsection II showed the strongest decrease among the Cyanobacteria from on average 2.3% at reference sites to disappearing almost completely at the seep sites. Alphaproteobacteria decreased towards the seep sites. This decrease was not uniform in all differentially abundant taxa of the Alphaproteobacteria, e.g. Rhodobacteriaceae increased from on average 1.7 to 3.4% at medium and high HI sites, whereas Rhodospirillaceae decreased from 8 to 3% and 1.5%, respectively. Flavobacteria were most abundant at reference and medium HI sites with on average 6% and 4%, respectively, and decreased toward high HI sites (1.8%). Similar to the case of Alphaproteobacteria this trend was not uniform, e.g. Zeaxanthinibacter already decreased significantly from reference sites (1.7%) to about 0.3% at medium and high HI sites. Within the Deltaproteobacteria several subgroups were differentially abundant, exhibiting divergent trends: whereas the order Desulfuromonadales (reference: 0.6%, medium HI: 1.1%, high HI: 2.9%) increased towards the seep sites with the strongest increase towards high HI sites, the family Desulfobacteraceae decreased severely towards high HI sites with an average relative abundance of 1.4% as compared with 6% at both reference and medium HI sites.

Differentially abundant OTUs largely fell into the previously mentioned differentially abundant bacterial taxa. Additionally, a large number of differentially abundant OTUs belonged to the actinobacterial OM1 clade (15 OTUs) and the gammaproteobacterial family JTB255 (25 OTUs) with a total relative abundance of 4% and 1.4%, respectively. Yet, these OTUs did not follow a consistent trend with highest relative abundance at sites of any of the three hydrothermal influence categories.

The archaeal community was dominated by the Marine Group I archaeon Candidatus Nitrosopumilus, which accounted for 36% of the total number of sequences in the archaeal 16S data set (Fig. 3B). The next most abundant groups were Deep Sea Hydrothermal Vent Group 6 (18%), other Marine Group I archaea (15%), Marine Benthic Group E (10%) and the genus Cenarchaeum (9%) also of the Marine Group I. The candidate genus Nitrosopumilus was generally more abundant at high HI sites where they reached maximum relative abundances of more than 60%. At reference and medium HI sites, they were much less abundant with on average 21% and 13%, respectively. A full list of the differentially abundant archaeal taxa is proved in Supplementary Table S5.

DISCUSSION

Usage of naturally CO2-rich sites as OA analogues

We used sampling sites at the CO2 seeps in Papua New Guinea (PNG) as analogues for OA research. The detailed characterization of the sediment based on 29 different environmental parameters showed that carbonate chemistry alone, specifically the pH gradient, was insufficient to describe the environmental conditions at these sites. Whereas some of the observed environmental parameters, e.g. sediment permeability, grain size, and nitrogen, organic and inorganic carbon content, supported the classification of sites along the pH gradient, mainly pore water element concentrations varied among sampling sites with comparable pH. The enrichment or depletion of certain elements in the pore water can be used as an indicator of hydrothermal activity (German and von Damm 2003). Here, the increased concentrations of lithium, manganese and silicon pointed towards a stronger hydrothermal influence in several of the sampling sites associated with the CO2 seeps, which were therefore termed ‘high HI’ sites.

The hydrothermal character of natural CO2 seeps has been recognized before as a potentially confounding influence for OA research (Vizzini et al. 2013). Vizzini et al. (2013) suggest referring to natural CO2 seeps as low pH environments rather than analogues for OA. Here, the distinction between medium and high HI sites allowed us to define a level of hydrothermal influence that may still be considered a reasonable OA scenario. Because of their strong hydrothermal character, we propose that high HI sites exhibit conditions that are less likely to occur under future OA than those at medium HI sites. Consequently, changes in the microbial communities that were only specific to high HI sites should be interpreted with caution because they may bias the assessment of OA impacts. Especially in cases where the community analysis is based on categories and not numerical environmental variables, accounting for the hydrothermal influence at naturally CO2-rich sites is of paramount importance.

Because of logistic constraints, environmental data in OA studies on sediment microbial communities is often collected from the water column and rarely from the pore water and sediment directly (Kerfahi et al. 2014; Raulf et al. 2015). Here, particularly at the CO2 seep sites, we observed strong deviations of the same parameters, such as pH, between bottom and pore water measurements. Additionally, we detected pronounced differences in pore water element concentration between sampling sites, which were not recorded in previous observations on water column chemistry (Fabricius et al. 2011). Such deviations question the usage of bottom water measurements as a proxy for conditions in the sediment at such sites, and will strongly affect our ability to accurately assess the driving environmental forces behind patterns of microbial richness and community composition in marine sediments. In this study, the divergence between bottom water and pore water parameters is evident in their explanatory power regarding the variation in microbial community composition. Noteworthy, basing the description of the microbial community on sediment rather than bottom water parameters increased the model fit and the percentage of explained variation of the microbial community composition by approximately 25%.

In both bottom water and sediment models, pH was among the factors significantly explaining changes in microbial community composition, which in the case of bottom water pH was also observed previously in PNG (Raulf et al. 2015). Yet, many of the other observed parameters were statistically equally important, e.g. temperature and permeability. Furthermore, as commonly found in natural systems, it was difficult to account completely for confounding factors and to isolate the effects of interest due to a high degree of covariation among environmental parameters (Sunagawa et al. 2015). Therefore, it is possible that changes in the microbial community that were related to pH, as observed by RDA, could also be related, directly or indirectly, to silicate concentrations or organic carbon content, since these factors were strongly correlated. Only a small, although significant, percentage of explained variation was attributed to pH alone when accounting for covariation. These circumstances make the interpretation of the data in a biological context more difficult because changes in microbial community structure can often be assigned to several factors (Hanson et al. 2012).

Furthermore, there may be causal relationships between environmental parameters. Long-term acidification may reduce the accumulation of coarse carbonate sediments, resulting in finer sediments composed of silicate sands (Artur Fink, personal communication). The lower permeability of finer sediments may limit water circulation and the supply of oxygen and organic matter into the sediment and may thus change the microbial habitat immensely, e.g. in terms of energy availability (Schöttner et al. 2011). Therefore, changes in the microbial community attributed to permeability may constitute indirect pH effects. The interdependencies of environmental factors may also help to explain the contrasting results regarding the effects of reduced pH on sediment microbial communities from different natural CO2-rich sites: before the onset of pH decrease, the initial conditions in the sediment may have varied among naturally CO2-rich sites even if present pH values are similar. To improve the comparability of naturally CO2-rich sites and to reconcile divergent microbiological results from different naturally CO2-rich sites, knowledge of sediment parameters such as permeability is therefore required.

As an alternative to the detailed analysis based on measured environmental parameters, the analysis of sediment microbial communities can be based on acidification categories, which summarize the conditions at the sampling sites (Kerfahi et al. 2014; Taylor et al. 2014). This approach was implemented here by using hydrothermal influence categories to describe the sediment microbial community. Although hydrothermal influence categories had generally less explanatory power than numeric environmental parameters, we were still able to explain about 20–38% of the total variation in microbial community composition. Values within this range of explained variation are not uncommon in microbial ecology (Hanson et al. 2012).

Changes in microbial community composition

The present study constitutes the second time the same sampling area in PNG was investigated to assess OA impacts on sediment microbial communities, with the initial exploration taking place in 2010–11 (Raulf et al. 2015). By using a larger sample size and a more comprehensive environmental characterization, we were able to confirm the general observations made in the previous study, thus underlying the reproducibility of the results obtained under natural conditions. Indeed, we found a similar overall taxonomic composition of the bacterial and archaeal communities as well as consistent shifts in microbial community composition and similar turnover rates of ARISA OTUs from reference to seep sites (Raulf et al. 2015). However, we showed that the composition of the microbial community at the seep sites was not comparable between Reef 1 and Reef 2, which was not detected previously (Raulf et al. 2015). Furthermore, trends in α-diversity of the microbial community here and from Raulf et al. (2015) diverged markedly: in Raulf et al. (2015) estimated bacterial richness and the percentage of rare bacterial OTUs increased as pH decreased, whereas here we detected a decreasing trend in richness that was, however, not statistically significant. A decrease in microbial richness would be consistent with similar observations in the Mediterranean (Taylor et al. 2014). Whether those contrasting results from PNG were caused by technical biases or by natural temporal variation still needs to be further investigated.

We used differential sequence abundance of bacterial and archaeal taxa to identify potential ‘winners’ and ‘losers’ under low pH conditions. In many samples, the relative abundance of differentially abundant taxa accounted for the majority of sequences in the sample, suggesting a large impact of hydrothermal influence on the microbial community. Among potential winners were the bacterial phyla Chloroflexi, Chlorobi and Deferribacteres. The mostly thermophilic Chloroflexi only showed a pronounced increase at high HI sites. Isolates of the three dominant differentially abundant classes Anaerolineae, Ardenticatenia and Caldilineae have been shown to prefer pH conditions around 7 or below (Yamada et al. 2006; Kawaichi et al. 2013), which were observed at the seep sites. Unlike other Chloroflexi, these three classes most likely do not photosynthesize, and grow chemoheterotrophically (Yamada et al. 2006; Kawaichi et al. 2013). Chloroflexi have also previously been found to be more abundant at decreased pH in sediments and associated with corals and sponges in PNG (Morrow et al. 2014; Raulf et al. 2015).

Chlorobi, which are also known as green sulfur bacteria, have previously been detected at shallow-water hydrothermal seeps exhibiting decreased pH (Maugeri et al. 2009; Giovannelli et al. 2013). Unlike other Chlorobi, representatives of the class Ignavibacteria, which contained the most dominant differentially abundant OTU in our data set, were characterized as chemoheterotrophic, strictly anaerobic bacteria lacking genes for photosynthesis and sulfur oxidation (Iino et al. 2010; Liu et al. 2012). Bacteria of the Deferribacteres, specifically the genus Caldithrix, are most likely anaerobic chemoorganotrophs known to be associated with hydrothermal seeps and decreased pH (Miroshnichenko et al. 2003, 2010). Chloroflexi, Chlorobi and Deferribacteres may constitute important carbon degraders at PNG seep sites.

The alphaproteobacterial family Rhodobacteraceae was another potential winner at lower pH. Rhodobacteraceae, also called purple sulfur bacteria, have previously been reported to increase in abundance at the seep sites in PNG as well as other naturally CO2-rich sites (Taylor et al. 2014; Raulf et al. 2015). The most abundant OTU of the Rhodobacteraceae was closely related to Rhodovulum species, which are able to oxidize iron and reduced sulfur compounds for anoxygenic photosynthesis preferably around pH 7 or lower (Straub, Rainey and Widdel 1999). Both dissolved iron and sulfide were present at the seep sites in PNG, although elevated sulfide concentrations were only detected in sediment layers deeper than 3 cm (Artur Fink, personal communication). These conditions might have facilitated the increased abundance of Rhodobacteraceae in general and Rhodovulum in particular. In coral reef ecosystems Rhodobacteraceae have been implicated in coral and seaweed diseases and their abundance in these associations was even higher at lower pH, supporting the hypothesis of a decline in reef health under OA (Meron et al. 2011; Bourne et al. 2013; Egan et al. 2013).

Rhodospirillaceae, an alphaproteobacterial family known as purple non-sulfur bacteria, were among the potential losers at decreased pH, with some OTUs disappearing completely already at medium HI sites. Rhodospirillaceae constitute a physiologically very diverse bacterial family with many uncultured types, complicating the functional classification of the Rhodospririllaceae in general. Among the OTUs that could be classified at the genus level, Defluviicoccus is a chemoheterotroph that grows at temperatures below 30°C and at a variable pH range (Maszenan et al. 2005), which suggests that temperature might be the limiting factor for this genus at the seep sites.

Other potential losers were the cyanobacterial genus Pleurocapsa and the deltaproteobacterial family Desulfobacteriaceae. Although OA is expected to boost photosynthesis, Cyanobacteria in general are expected to benefit less from the increased availability of CO2 when co-occurring with eukaryotic photosynthetic organisms (Koch et al. 2013). Preliminary data indicated that bulk net photosynthesis did not vary significantly between sampling sites in PNG (Artur Fink, personal communication), suggesting that there might be a replacement of cyanobacterial photosynthetic taxa at the seep sites by other photosynthetic organisms. Given the increased silicate concentration, diatoms are likely to dominate primary production at the seep sites. The genus Pleurocapsa is further capable of calcification (Krumbein and Giele 1979; Rippka et al. 1979), which may increase its sensitivity to decreased pH.

Desulfobacteraceae reduce sulfate to degrade organic carbon anoxically (Kuever, Rainey and Widdel 2005). Although Desulfobacteraceae are expected to be able to adapt to decreased pH (Koschorreck 2008), they have been shown to decrease at CO2 seeps (Raulf et al. 2015). This decrease was also observed here in the abundance of Desulfobacteraceae with the strongest decrease towards high HI sites. Measurements of sulfate reduction rates further documented a strong decline in sulfate reduction from reference to high HI sites (Artur Fink, personal communication). The decrease in Desulfobacteraceae is accompanied by an increase in Desulfuromonadales, which can reduce iron or manganese as alternative electron acceptors to sulfate (Vandieken and Thamdrup 2013). Since iron and manganese were enriched at the seep sites, Desulfuromonadales may outcompete Desulfobacteraceae as carbon degraders.

Other bacterial taxa did not show a consistent trend, either in the monotony of the trend from reference to medium to high HI sites, or when considering the trends at higher taxonomic resolution levels. Flavobacteria, with OTUs mostly related to aerobic carbon degraders (Bernardet, Nakagawa and Holmes 2002), were most abundant at reference and medium HI sites and decreased severely towards high HI sites. At higher taxonomic resolution levels, divergent trends emerged between reference and medium HI sites with several flavobacterial taxa either increasing or decreasing. The strong decrease at high HI sites was common to most differentially abundant flavobacterial taxa and may be attributed to the lower oxygen availability at these sites resulting from a shallower oxygen penetration depth as compared with reference sites. The flavobacterial genus Zeaxanthinibacter was among the taxa that displayed a continuous decrease from reference to medium and high HI sites. The decrease in Zeaxanthinibacter may also be related to decreased zeaxanthin pigment concentrations at the seep sites (Artur Fink, personal communication; Asker, Beppu and Ueda 2007). Members of this genus are characterized as strictly aerobic bacteria, which grow at a range of pH values and can degrade long-chain carbohydrates (Asker, Beppu and Ueda 2007). Their aerobic lifestyle supports the hypothesis that oxygen may be the limiting factor restricting the distribution of at least some Flavobacteria. At high HI sites, Flavobacteria might be replaced by Chlorobi, Chloroflexi and Deferribacteres as major carbon degraders under anoxic conditions. Our results were only partially consistent with previous studies at naturally CO2-rich sites, which have so far unanimously supported an increase in Flavobacteria at decreased pH potentially linked to an increased availability of organic matter (Kerfahi et al. 2014; Taylor et al. 2014). However, in our study the organic carbon content of the sediment was lower at the seep sites (Artur Fink, personal communication), which may contribute to explaining these divergent findings.

The dominant archaeal taxon, Candidatus Nitrosopumilus, has previously been described in PNG to increase in abundance as pH decreases (Raulf et al. 2015). Its role as ammonia oxidizer may contribute to maintaining nitrification in reef sediment under reduced bacterial ammonia oxidation (Laverock et al. 2013; Raulf et al. 2015). We also detected an increase in Candidatus Nitrosopumilus at decreased pH, which was, however, restricted to high HI sites and may therefore be less likely to occur under future OA.

In summary, we were able to confirm the majority of the trends in differentially abundant taxa, which were previously observed either in PNG or at other naturally CO2-rich sites. However, in specific instances, we detected taxa that have not been identified in OA research before, or which did not display the same trends as reported previously. These inconsistencies with previous studies may have several reasons. (i) The distinction between medium and high HI sites allowed us to identify taxa only changing in abundance towards sites exhibiting conditions that may be less likely under future OA scenarios. The most drastic changes in microbial community composition were observed at the high HI sites and thus are likely correlated with parameters other than pH. (ii) In many instances, pH may not be the primary factor limiting the distribution of a microbial taxon even at medium HI sites. Other factors such as oxygen availability, temperature, organic carbon content and the presence of alternative electron acceptors may be equally or more likely to explain the differential abundance of microbial taxa between hydrothermal influence categories. (iii) It is important to consider the taxonomic resolution when analysing microbial communities. Especially in physiologically highly diverse taxa, a comparison at low taxonomic resolution levels might be insufficient to recover more detailed patterns in microbial community composition. (iv) As already suggested by Raulf et al. (2015) the sampling area in PNG may be subject to a large temporal variation in microbial community composition. To better account for this phenomenon and to resolve inconsistencies with previous results, repeated monitoring studies, as conducted here, would be necessary.

CONCLUSION

We conclude that naturally CO2-rich sites such as the CO2 seeps in PNG may continue to be used for OA research on sediment microbial communities under the premise that the environmental conditions in the sediment are well documented. Furthermore, we recommend caution in attributing changes in microbial communities to acidification without a careful consideration of other environmental parameters. Thus, for future research on the impact of OA on sediment microbial communities, we strongly recommend performing a detailed assessment of the environmental conditions in the sediment and of other parameters than those of the carbonate system. This will enable a more reliable selection of naturally CO2-rich sites as analogues for OA scenarios and will also improve the comparability between studies from different naturally CO2-rich sites.

Click here for additional data file.