Ervilia castanea (Mollusca, Bivalvia) populations adversely affected at CO 2 seeps in the North Atlantic Science of the Total Environment

castanea CO 2 seeps and and net-calci CO 2 levels. at high CO 2 sites. this at the CO 2 seeps. size of E. castanea were positively correlated with Chl-a in sediment. SiteswithnaturallyhighCO 2 conditionsprovideuniqueopportunitiestoforecastthevulnerabilityofcoastaleco- systems to ocean acidi ﬁ cation, by studying the biological responses and potential adaptations to this increased environmental variability. In this study, we investigated the bivalve Ervilia castanea in coastal sandy sediments at reference sites and at volcanic CO 2 seeps off the Azores, where the pH of bottom waters ranged from average oceanic levels of 8.2, along gradients, down to 6.81, in carbonated seawater at the seeps. The bivalve population structurechangedmarkedlyattheseeps.LargeindividualsbecamelessabundantasseawaterCO 2 levelsroseand were completely absent from the most acidi ﬁ ed sites. In contrast, small bivalves were most abundant at the CO 2 seeps. We propose that larvae can settle and initially live in high abundances under elevated CO 2 levels, but that high rates of post-settlement dispersal and/or mortality occur. Ervilia castanea were susceptible to elevated CO 2 levels and these effects were consistently associated with lower food supplies. This raises concerns about the effects of ocean acidi ﬁ cation on the brood stock of this species and other bivalve molluscs with similar life history traits.


H I G H L I G H T S
• The bivalve Ervilia castanea was studied at volcanic CO 2 seeps and reference sites. • Abundance, size and net-calcification were inversely related to CO 2 levels. • Large individuals were scarce or absent at high CO 2 sites. • Recruitment of this bivalve was highest at the CO 2 seeps. • Abundance and size of E. castanea were positively correlated with Chl-a in sediment.

G R A P H I C A L A B S T R A C T
a b s t r a c t a r t i c l e i n f o

Introduction
Increasing atmospheric CO 2 concentrations have caused the pH of the surface ocean to fall by 0.1 units since the pre-industrial age (Bindoff et al., 2019;Dupont and Pörtner, 2013;Feely et al., 2009;Keeling et al., 2017) and further decreases, between 0.1 (RCP2.6) or 0.3 (RCP8.5) pH units, are expected by 2100 (Bindoff et al., 2019). This process of ocean acidification affects calcification of marine organisms, leading to permanent damage and reduced viability (Shirayama, 2005) and causing major changes to marine ecosystems (Doney et al., 2012). Shallow-water CO 2 seeps provide analogues for the effects of chronic ocean acidification, at spatial and temporal scales that are difficult to mimic in aquaria or mesocosms (Hall-Spencer and Harvey, 2019). Care is needed to assess potential anomalies in temperature, alkalinity, metals and H 2 S which are frequently found at these volcanic sites and can also influence biological responses, but are not associated with ocean acidification (Vizzini et al., 2013;Mishra et al., 2020).
This type of in situ approach has shown that ocean acidification can lead to ecosystem simplification and reorganization of marine communities (Fabricius et al., 2015;Kroeker et al., 2013aKroeker et al., , 2011Agostini et al., 2018). The structure of benthic marine communities changes greatly when they are exposed to carbonated seawater (Widdicombe et al., 2009), with decreases in biodiversity, biomass and trophic complexity (Hall-Spencer and Harvey, 2019). Communities tend to change from being dominated by calcareous organisms to non-calcareous ones (Baggini et al., 2014;Christen et al., 2013;Hall-Spencer et al., 2008), granting the latter a perceived advantage in a future acidified ocean (Kroeker et al., 2013b;Vihtakari et al., 2013).
Organisms with calcareous shells and skeletons, such as coralline algae, molluscs and echinoderms, are expected to be adversely affected by increasing CO 2 levels (e.g., reductions in abundance or local extinctions; Kroeker et al., 2011;Ricevuto et al., 2012) due to susceptibility of the calcification process to reduced carbonate availability (Bindoff et al., 2019). Garilli et al. (2015) found two gastropod species dwarfed as an adaptation to carbonated seawater, which allowed them to maintain calcification and partially repair shell dissolution. Other studies show that pH (or pCO 2 ) on their own (apart from the carbonate saturation state) can cause physiological effects that result in less energy available for mineralization (Cyronak et al., 2016;Palmer, 1992;Roleda et al., 2012;Waldbusser et al., 2015).
These results indicate that the responses of the early life stages of bivalves to increased CO 2 go well beyond the detrimental effects on shell calcification.
Responses to ocean acidification are species-specific and even closely related species can react differently, depending on various factors, such as life history traits, habitat, sex, nutritional status, adaptation potential and methodological conditions Doney et al., 2012;Dupont and Thorndyke, 2009;Garilli et al., 2015;Kroeker et al., 2013aKroeker et al., , 2011Marčeta et al., 2020;Pörtner, 2008;Vihtakari et al., 2013). Individual responses to ocean acidification are of key importance for any given species, but they also have cascading effects along trophic relations and ecological interactions, making them important drivers in shaping community structures and phase shifts (Kroeker et al., 2011;Range et al., 2014b). Consequently, the direct effects of CO 2 on competitors, prey, predators and the quantity and quality of food supply can be translated to density control mechanisms (Metaxas, 2015;Range et al., 2012) and mediated by adaptation or acclimation (Kroeker et al., 2011).
The species targeted in this study, Ervilia castanea (Montagu, 1803), is a semelid bivalve (Morton, 1990) found from south of the British Isles into the Mediterranean (Babío and Bonnin, 1987). It has a brownish, oval and elongated shell, which reaches a maximum length of 15 mm (Macedo et al., 1999;Morton, 1990). In the Azores, where the majority of the bivalve species tend to be~50% smaller than their continental conspecifics (Morton et al., 2013), the maximum length reported is 6.7 mm (Morton, 1990). This generalized dwarfism of bivalves is commonly found around oceanic islands and has been related to the low availability of food in these oligotrophic waters (Morton et al., 2013). Individuals of this species have a high potential reproductive output, with maturation occurring at a length of 3.5 mm for females and 5.5 mm for males (Morton, 1990). Ervilia castanea occupies subtidal coarse-grained, well-sorted sediments overlain by clean, nutrient poor, oceanic waters, to depths exceeding 40 m (Morton, 1990). In the Azores it is numerically dominant on sandy substrata and mainly predated upon by naticid gastropods (Morton, 1990;Morton et al., 2013).
In this study, we collected in situ evidence about the effects of elevated seawater CO 2 on populations of the bivalve Ervilia castanea, along submarine volcanic CO 2 gradients in the Azores. We recorded abundance, size-structure, biomass and tissue composition. Three hypotheses were tested: we expected the abundance (Hyp. 1) and netcalcification (i.e., accretion minus dissolution, Hyp. 2) of Ervilia castanea to decrease with rising CO 2 levels and for the size structure of the population (Hyp. 3) to be modified along the CO 2 gradients.

Selection of the study sites
The Azores are nine volcanic islands located on the triple junction of the Eurasian, North-American and Nubian tectonic plates. Shallowwater volcanic CO 2 seeps occur around several of the Azorean islands and seamounts (Couto et al., 2015;Rajasabapathy et al., 2014;Tribollet et al., 2018;Viveiros et al., 2016;Wallenstein et al., 2013Wallenstein et al., , 2009. Two of these coastal sites, located off the islands of São Miguel and Faial, were chosen for this study, based on their shallow depth and accessibility, as well as previous knowledge about the intensity and composition of the gas emissions. Ponta da Lobeira (São Miguel, 37°43′31.8″N 25°19′ 01.6″W, Supplementary Video 1) is a site with strong degassing at 7 m depth, located <2 km from the village of Ribeira Quente, on the south flank of Furnas Volcano. Ponta da Espalamaca (Faial, 38°32′36.4″N 28°35′48.4″W, Supplementary Video 2) is a site with more diffuse degassing, at 37 m depth, located 3 km from Horta port (Fig. 1). Supplementary videos one and two show the type of benthic habitat and the degassing fluxes at Ponta da Lobeira and Ponta da Espalamaca, respectively.

Experimental design
Off each island, three sampling sites, with similar depth, bottom typology exposure to currents and wave action, were selected: a significant degassing site (High CO 2 ); a transitional site (Intermediate CO 2 ), where the seawater carbonate chemistry was still affected, but with no visible CO 2 release; and a reference site (Reference), unaffected by the CO 2 emissions. Test dives with a pH multiprobe (YSI6600) were used to define these sites. In São Miguel, sampling was done between 3 and 7 July 2014, at depths ranging between 6 and 10 m. All sites in São Miguel were relatively sheltered from the dominant wave action. Sampling in Faial was done between 10 and 15 July 2014, at depths between 28 and 37 m. Given the proximity to the Faial-Pico channel, all the sites in Faial were strongly exposed to tidal currents.

Free gas composition and flux
Free gas samples were collected at the High CO 2 sites using Giggenbach bottles (Supplementary Video 3) and gas composition was analysed by titration and gas chromatography (for additional details, Ferreira et al., 2005). Inverted plastic funnels, connected to volumetric flasks, were used to measure gas emission rates. In São Miguel (Ponta da Lobeira), this sampling was stratified into low, medium and high flux areas, defined by visual inspection of bubbling, with three replicate samples collected at each level. In Faial (Ponta da Espalamaca), intensity of bubbling was more homogeneous, so the different levels of emission could not be established and two replicates were collected in each of three haphazardly selected points of emission (A, B and C).

Physicochemical properties of seawater
Seawater temperature, salinity, dissolved oxygen and total pH (pH t ) were measured at least once per day at each sampling site, over the course of 3 or 4 days. A YSI6000 multiprobe was used for continuous (5 second intervals) in situ measurements during vertical profiles (3 days per site) and diving operations (1 or 2 days per site). A Niskin bottle was also used during the vertical profiles, for collecting discrete samples of sub-surface (0-2 m) and near-bottom (6-10 m in São Miguel and 15-38 m in Faial) seawater, for the determinations of dissolved inorganic nutrients, total pH (pH t ) and total alkalinity (TA). pH t was measured from samples without headspace using a glass electrode (WTW 340i pH meter) and calibrated with a Tris seawater buffer provided by the A. Dickson laboratory. TA samples were filtered through a 0.2 μm pore size, stored at 4°C and measured within 48 h by potentiometric titration according to Dickson et al. (2003) using a Metrohm 848 Titrino Plus equipped with Metrohm 869 Compact Changer Sampler. TA measurements were corrected with certified reference material (Dickson, 2010), at 20 μmol·kg −1 accuracy and 2 μmol·kg −1 precision. Samples for the determination of dissolved inorganic silicate and phosphate were filtered using a 0.2 μm polyethersulfone syringe filter and stored at −20°C until analysis. Concentrations of silicate and phosphate were measured following Hansen and Koroleff (1999), by means of a spectrophotometer (Cary 60, Varian). Seawater carbonate chemistry speciation was calculated from temperature, salinity, silicate, phosphate, multiprobe pH t values and TA using CO 2 SYS (Lewis and Wallace, 1998) with the equilibrium constants determined by Mehrbach et al. (1973) as refitted by Dickson and Millero (1987), by interpolating discrete variables from Niskin sampling (silicate, phosphate and TA) with continuous multiprobe records (pH t , temperature and salinity). Accordingly, the calculations were done using the average values of TA, silicate and phosphate for each vertical profile (i.e., average of bottom and surface samples per sampling site and per day).

Macrofauna
Sediment macrofauna samples were collected by scientific divers, from an area of 0.25 m 2 excavated to 15 cm depth, using an airlift fitted with a 500 μm mesh bag (Supplementary Video 4). Replicate samples were spaced by 2 to 6 m and haphazardly positioned within each site, except in the High CO 2 sites, where the points of gas emission were explicitly avoided. In São Miguel, nine replicates were collected at each site. In Faial, time restrictions and operational limitations while diving at greater depths and in stronger currents, enabled the collection of only six samples of macrofauna in the Reference and Intermediate sites and seven in the degassing site.
All organisms retained in the bag were preserved in 5% buffered formaldehyde, hand-sorted and later preserved in 70% ethanol until taxonomic determination and counting. All the Ervilia castanea were individually measured and subsequently pooled per replicate and shell length (1 mm size-classes, 0 to 8 mm). These pools of bivalves were dried for a minimum of 36 h at 60°C (dry weight -DW) and subsequently burnt in a muffle furnace for 4 h at 550°C (ash weight -AW). Ash-free dry weight (AFDW) and relative ash weight (%AW) were estimated, per replicate and size-class using Eqs. (1) and (2), respectively:

Sediments
Hand-held PVC cores (3.5 cm diameter per 15 cm length) were used to collect two sediment samples, adjacent to each macrofauna sample. The first core was used for determination of granulometry and organic matter/carbonate content. A set of six sieves (63 μm, 125 μm, 250 μm, 500 μm, 1 mm and 2 mm) was used for sediment pooling per grain size. After sieving, each fraction was weighed and related to the total weight of the sample, as a percentage. The organic matter content and carbonate content in sediments was determined by the loss-ofignition (LOI) method (Heiri et al., 2001). Sediment was dried in an oven for 72 h at 100°C (DW 100 ), from which, 10 g of each sample were burnt for 8 h at 550°C (DW 550 ) and subsequently for 2 h at 900°C (DW 900 ). Organic matter content (LOI 550 ) and carbonate content (LOI 900 ) were then determined using Eqs. (3) and (4), respectively.
The second core was used for carbon/nitrogen ratio and photosynthetic pigment analyses. Sediment samples were frozen and kept in the dark for preservation of chloroplastic pigments. Approximately 12 g of frozen sediment were lyophilized in a Telstar LyoAlfa 15 for 48 h, to get the 5 g of material needed. For total organic carbon (TOC) and total nitrogen (TN) samples were homogenized and acidified with dilute hydrochloric acid until complete decarbonization. TOC was measured with a Thermo Scientific Flash 2000 elemental analyser. For analyses of the chloroplastic pigments, samples were homogenized and pigments extracted in 90% acetone. Pigment separation was done using reverse phase high performance liquid chromatography (HPLC) and measured with a Gilson fluorescence detector according to Wright and Jeffrey (1997).

Statistical analysis
Two-way analyses of variance (ANOVA) were used to test for differences in abundance and biomass of Ervilia castanea within each island (São Miguel and Faial). CO 2 level (High, Intermediate and Reference) and size-class were considered as fixed orthogonal factors. The data was untransformed and Cochran's test was used to test for heterogeneity of variance. When appropriate, multiple comparisons for means were done on significant main effects using Student-Newman-Keuls (SNK) tests, with adjusted p-values for each pairwise comparison. Univariate analysis and graphs were done using SigmaPlot software (Version 13.0, Systat Software, Inc.) and R software (Version 3.1.1).
Distance-based linear modelling (DistLM) was used to investigate the relationship between sediment variables and the size structure of Ervilia castanea (Bray-Curtis similarity matrix). The significance of the relationship was first evaluated for individual variables with marginal tests (999 permutations). Significant variables were included in the model selection using BEST routine (i.e. all possible variable combinations). Distance-based redundancy analysis (dbRDA) was used for the ordination and visualization of the best overall DistLM solution, according to Akaike's Information Criteria (AIC). All the multivariate analyses were done using PRIMER 6 (Version 6.1.13) statistical package with the PERMANOVA add-on (Version 1.0.3).

Free gas composition and flux
Free gas composition at both vent sites was mainly carbon dioxide. In São Miguel, CO 2 accounted for >99 mol% of gas composition, while in Faial CO 2 accounted for between 98 and 99 mol%. Nitrogen, argon and oxygen followed with values ranging between 0.02 mol% and 1.27 mol%. The remaining gases (helium, hydrogen sulfide, hydrogen and methane) had very low concentrations, with some, mainly in Faial, being absent or below detectable limits (Table 1).

Physicochemical properties of seawater
Seawater temperature and salinity varied between the two islands and with depth, but not among sampling sites within each island. Average salinity was higher in Faial (37.69 ± 0.05), relative to São Miguel (36.83 ± 0.05). Average sea surface temperature did not vary between the two Islands at the time of sampling (20.61 ± 0.14), but the bottom layer was slightly warmer in São Miguel (20.08 ± 0.20), relative to Faial (18.90 ± 0.19). The carbonate chemistry variables showed a similar trend in both islands, consistent with the predefined gradients. In São Miguel, pH t ranged from 6.92 (High CO 2 site) to 8.16 (Reference site). In Faial pH t ranged between 6.81 (High CO 2 site) and 8.2 (Reference site). The lowest pH t values were always found at the degassing sites, while the highest values were similar across all the sampling sites (Fig. 2). pCO 2 and HCO 3 2− increased along the gradients, while CO 3 , OH − and CaCO 3 saturation (calcite and aragonite) showed the inverse pattern. Undersaturation (Ω < 1) of aragonite and calcite was observed for the degassing sites in both islands and for the intermediate site in Faial. These values were, however, restricted to the extremes of the frequency distribution ( Fig. 2 and Supp. Figs. 1 and 2) and did not affect the central tendency metrics (Table 2). Undersaturated conditions were, therefore, relatively infrequent, even at the high CO 2 sites.

Population structure of Ervilia castanea
A total of 3134 Ervilia castanea individuals were collected during this study, representing 58% of the overall abundance of macrofauna. E. castanea was the dominant organism in all sampling sites, except for the Intermediate site in Faial, where Nematodes of the genus Metoncholaimus were numerically dominant. Densities of E. castanea up to 1140 ind·m −2 were recorded in São Miguel, which generally had larger abundances, relative to Faial (Fig. 3). In Faial, the sizefrequency distribution was clearly skewed towards smaller individuals at the High CO 2 site, while in São Miguel the sites tended to have a more similar distribution. Average abundances (across all size classes) consistently peaked at reference sites and this effect was significant in both Islands (Table 3). In contrast, the smallest size class (<1 mm), composed of recently settled postlarvae (i.e., plantigrades, sensu Carriker, 1961) and juveniles (spat), showed consistently larger abundances in sites affected by CO 2 than in reference sites, which explains the significant interaction term for São Miguel. The dominant size class in Faial (1-2 mm) also showed a similar pattern, reaching peak abundance in the high CO 2 site (Fig. 3).

Biomass of E. castanea
The smallest (<2 mm) and largest (>5 mm) size classes were excluded from the biomass analyses, due to lack of precision in weighing and absence in the high CO 2 sites, respectively. The average biomass per individual (DW, AW and AFDW) increased with shell length for the three intermediate size classes (between 2 and 5 mm shell length, Fig. 4). This pattern was visible in all sampling sites, except for the intermediate site in Faial, where individuals larger than 2 mm were absent. Relative ash weight differed significantly among sites (%AW , Table 3), with the dominant size classes showing consistently larger values in reference sites (92% in both islands), than in sites affected by CO 2 emissions (83% in São Miguel and 75% in Faial). Although no quantification of shell dissolution or strength was done in this study, alterations were clearly apparent on the few larger individuals found at Intermediate and High CO 2 areas, particularly around the umbo, where it was usually associated with some degree of erosion to the periostracum ( Fig. 5B and D).

Relationship between sediment characteristics and population structure of E. castanea
DistLM marginal tests indicated that 13 out of the 15 sediment variables considered were significantly related to the distribution and size structure of E. castanea (Table 4). The only exceptions were silt (sediment finer than 63 μm) and total nitrogen (TN). These two variables were, therefore, excluded from the model selection. The best overall DistLM solution included 4 sediment variables: 2 mm grain size, carbonate content, chlorophyll a and other pigments (diadinoxanthin, diatoxanthin, zeaxanthin and lutein). The coarser sediments were strongly associated with the intermediate and high CO 2 sites in Faial (Fig. 6), which also had the smallest abundances of E. castanea above 2 mm shell length (Fig. 3). Chlorophyll a was the dominant photosynthetic pigment in all sampling sites, Table 1 Gas composition (mean values) of samples collected at each of the two degassing systems (Ponta da Lobeira -São Miguel, n = 3; Ponta da Espalamaca, Faial, n = 2). Carbon dioxide (CO 2 ), hydrogen sulfide (H 2 S), methane (CH 4 ), hydrogen (H 2 ), helium (He), nitrogen (N 2 ), oxygen (O 2 ), argon (Ar). B.d.l. (below detection limit).

Discussion
Gas released from both seep systems was mainly CO 2 (>98.6%), with N 2 , O 2 , Ar and He also present in the dry phase. Residual fractions of other typical hydrothermal gases (H 2 S, H 2 and CH 4 ) were also detected at Ponta da Lobeira (Table 1). The sulfuric component was, however, only present when the gas flux was medium to high (≥10 L/min/m 2 ) and in relatively small concentrations (<10 ppm), when compared to the sub-aerial fumaroles of the Furnas Volcano (Caliro et al., 2015) or to other submarine vents (Tarasov et al., 2005). Boatta et al. (2013) recorded much higher concentrations of H 2 S (400 ppm) in the vent bubbles at Vulcano Island, Italy, concluding that only a small proportion of sulfide enters or remains in the aquatic phase and that the dissolved levels dropped below the detection limit at 5 m distance from the emission point. Accordingly, although we cannot dismiss the possibility of Fig. 2. Violin plots for pH (total scale), CO 2 partial pressure (pCO 2 ), saturation of calcite (Ω Ca ) and aragonite (Ω Ar ); the frequency distribution (grey area), median (white dot) and quartiles (black bars) are presented for each site in each island, across all depths and dates of sampling. toxic gases, other than CO 2 , affecting the bivalves at Ponta da Lobeira, this potential confounding factor was clearly less important here than in most previously studied shallow vent sites. Other types of contaminants usually present in vent emissions, such as trace metals, were not assessed in our study. Their adverse biological effects are, however, usually localized to a few hundred meters from the vents (Vizzini et al., 2013), so it is unlikely that they could have affected the Intermediate CO 2 sites.
Seawater temperature, salinity and total alkalinity varied within narrow ranges among the different sampling sites and these variations were independent of the CO 2 emissions. Seawater pH, on the other hand, varied substantially among sites. This was mainly due to the increased frequency of extremely low pH values (down to 6.8) at the sites affected by CO 2 emissions, rather than a consistent decrease of the average pH at those sites. Despite this, the gradient of CO 2 partial pressure in São Miguel was clearly within the range predicted by the IPCC for the end of this century, which varies between 421 (RCP2.6) and 936 (RCP8.5) ppm (Hoegh-Guldberg et al., 2014;Pörtner et al., 2014). The gradients in seawater carbonate chemistry observed in São Miguel can, therefore, be unambiguously applied to infer on the response of benthic organisms to present conditions (Reference site) and future scenarios of ocean acidification (Intermediate and High CO 2 sites). The gradients in carbonate chemistry were less clear for Faial, due to the greater depth, stronger currents and the fact that pH data was collected along the entire water column (37 m). These conditions complicate the interpretation, as variations exceed the average scenarios for ocean acidification in the foreseeable future. Nevertheless, the regulation of carbonate chemistry in coastal waters involves a multitude of drivers (i.e., watershed processes, nutrient inputs, changes in ecosystem structure and metabolism), apart from anthropogenic CO 2 emissions, so variations of this magnitude are not unusual for this type of habitat . In fact, the potential adaptations of coastal ecosystems to this increased environmental variability, induced by interactions between global and local drivers, are critically important to Table 2 Seawater carbonate chemistry variables (mean ± SE) at each sampling site (R -Reference, I -Intermediate CO 2 , H -High CO 2 ) in São Miguel (SM) and Faial (FA). The variables represented are: temperature (Temp), salinity (Sal), pH (total scale), total alkalinity (TA), CO 2 (total scale), CO 2 partial pressure, bicarbonate (HCO 3 − ), carbonate (CO 3 2− ), ion hydroxide (OH − ), saturation of calcite (Ω Ca ) and aragonite (Ω Ar ). Each entry represents the average of bottom and surface measurements per sampling site and per day (n = 3).   Fig. 3. Abundance (mean ± SE) of the bivalve Ervilia castanea per size-class and CO 2 level in each of the islands; bars with different letters denote significant differences among CO 2 levels within each size class (SNK tests, p < 0.05); values in parenthesis represent mean abundances, across all size classes, at each CO 2 level in each of the islands. forecast their sensitivity and vulnerability to future ocean acidification and cannot be ignored. Hyp. 1 and Hyp. 2 were corroborated: Ervilia castanea significantly decreased in abundance and net-calcification, measured as relative ash weight, with increasing CO 2 levels. Individuals of all sizes were represented in reference sites, while larger individuals (above 3 mm shell length) were scarce or absent in sites affected by CO 2 emissions. This steep decline in the abundance of larger individuals was probably caused by high levels of post-settlement mortality, as a consequence of delayed or unsuccessful metamorphosis of recently settled pediveligers into postlarvae and juveniles (Bayne, 1965;Gobler, 2009, 2010), reduced shell calcification (Clements and Hunt, 2017;Green et al., 2013Green et al., , 2009Green et al., , 2004, growth impairment (Parker et al., 2010) or a combination of these factors (Waldbusser et al., 2010).
The average biomass per individual (DW, AW and AFDW) varied among size-classes, but not among sampling sites. Given the relatively small contribution of the soft tissues (AFDW) to the overall biomass (DW), the pattern observed for relative ash weight (%AW) seems to have been mainly driven by variations in the mineral fraction of the shells (i.e., calcium carbonate content). Negative net-calcification rates have frequently been observed for infaunal and epifaunal marine bivalves, in response to elevated CO 2 (see reviews by Gazeau et al., 2013 andHunt, 2017). Although the relative contributions of accretion and dissolution to this net-effect cannot usually be disentangled, several authors suggest increased post-deposition dissolution as the main cause for the observed losses of shell mass (Findlay et al., 2011;Nienhuis et al., 2010;Range et al., 2012;Rodolfo-Metalpa et al., 2011;Tunnicliffe et al., 2009 4. Individual body weight of the dominant size classes of Ervilia castanea per CO 2 level in each Island; dry weight (outline, means ± SE) and ash weight (filled area, mean); percentage values on top of the bars denote relative ash weight (mean). Fig. 5. Microphotographs of specimens of Ervilia castanea from different size classes and CO 2 levels: A (1-2 mm, reference), B (1-2 mm, high CO 2 ), C (2-3 mm, reference), D (2-3 mm, high CO 2 ). Scale bar is 1 mm.
carbonate saturation of the surrounding water and the physiological status of the organism, which is largely determined by food availability, whereas dissolution relates solely to the external saturation state (Marshall et al., 2019). Dissolution of the shell is recognized as an important source of post-settlement mortality in bivalves, having been associated with the exponential loss of individuals after the transition from the larval to the benthic stage (Green et al., 2004). The type of shell damage observed at sites affected by CO 2 emissions has been previously reported in several laboratory (Bressan et al., 2014;Gazeau et al., 2014;Range et al., 2012) and in situ studies (Rodolfo-Metalpa et al., 2011;Thomsen et al., 2010). Other studies have also found that increased CO 2 levels are often associated with decreases in the mechanical properties of the shell, such as, thickness (Bressan et al., 2014;Gaylord et al., 2011;Talmage and Gobler, 2010), density (Klok et al., 2014;Rühl et al., 2017), hardness and fracture resistance (Beniash et al., 2010;Dickinson et al., 2013;Gazeau et al., 2013;Welladsen et al., 2010). Thinner or weaker shells have also been shown to increase bivalve mortality through several pathways, mediated by interactions with other environmental stressors (Dickinson et al., 2012), predators (Grey et al., 2007;Sanford et al., 2014), competitors (Coleman et al., 2014) or parasites (MacLeod and Poulin, 2015).
A number of recent studies have shown that the detrimental effects of elevated CO 2 on the calcification and growth of bivalves can be mitigated or avoided when the food supply is not limiting (Melzner et al., 2011;Pansch et al., 2014;Ramajo et al., 2016bRamajo et al., , 2016aThomsen et al., 2013). We did not assess if food limitation occurred during our study, but sediments at reference sites had greater amounts of chlorophyll a and inorganic carbon, relative to sites affected by CO 2 emissions (Fig. 6, Table 4). This implies that reference sites were better feeding grounds for infaunal bivalves. Concentrations of chlorophyll-a in seawater showed a similar pattern, with the smallest values consistently being found at the High CO 2 sites (unpublished data). Chlorophyll-a was, therefore, negatively correlated with CO 2 level and positively correlated with abundance and size of E. castanea. This strongly suggests that the indirect effects of CO 2 , mediated by trophic interactions (i.e., food availability), may be equally or even more important than its direct effects on growth and survival of these clams.
Exposure to elevated CO 2 during development can also lead to reduced fitness and increased post-settlement mortality, due only to reductions in growth rates and energetic reserves (Hettinger et al., 2012;Waldbusser et al., 2015). Some of these mechanisms probably contributed to the observed effects on the population densities of E. castanea. Inferring on the relative importance of direct and indirect effects would, however, require manipulative experiments, which was clearly beyond the scope of the present study.
Hyp. 3 was also corroborated. In contrast to the overall pattern, abundances of the smaller size classes of E. castanea (<1 mm in São Miguel and <2 mm in Faial) generally increased with rising CO 2 levels. This pattern is clear, despite the abundance of spat (<1 mm) probably being underestimated, given that size at settlement is smaller than the Table 4 Sediment characteristics (means ± SE) by island and location. Granulometry (<63 μm, 63 μm, 125 μm, 250 μm, 500 μm, 1 mm, 2 mm and 4 mm), organic matter (OM), carbonate content (CaCO 3 ), total organic carbon (TOC) and total nitrogen (TN) are expressed as percentages. Other pigments include diadinoxanthin, diatoxanthin, zeaxanthin and lutein. C:N ratio was not independently determined, but calculated, so it was not considered in the DistLM models.

São Miguel Faial
Reference (n = 9) Intermediate (n = 9) High (n = 9) Reference (n = 6) Intermediate (n = 6) High (n = 10) 500 μm mesh used for sampling (Bayne, 1965;Chıćharo and Chıćharo, 2001). Although larval settlement was not directly measured, this suggests that settlement rates were higher in sites affected by CO 2 emissions than in reference sites. The fact that newly settled postlarvae were found only in high CO 2 sites (Fig. 5B) also seems to corroborate this hypothesis. Other populations of the same genus are known to have extreme temporal variations in abundance (i.e., boom-and-bust dynamics - Albano et al., 2016), which implies highly variable settlement rates. Potential spatial or temporal differences in settlement among experimental sites could, therefore, also have influenced the abundance and size-frequency distribution of the populations. The preponderance of small bivalves in naturally acidified environments has previously been observed and interpreted as evidence of continuous recruitment (Metaxas, 2015). As far as we are aware, however, there is no experimental evidence for positive effects of CO 2 on larval settlement of calcifying benthic marine invertebrates (Gazeau et al., 2013;Parker et al., 2013), while neutral (Bechmann et al., 2011;Crim et al., 2011) or negative effects (Cigliano et al., 2010;Ricevuto et al., 2012) have been observed. Accordingly, a direct enhancement of bivalve settlement by elevated levels of CO 2 seems extremely unlikely. Any positive effect was probably indirect, regulated by density-dependent control mechanisms.
Intraspecific interactions, like competition for resources (i.e., food or space) or predation (i.e., oophagy and larviphagy), have previously been recognized as important regulating mechanisms for population size in filter-feeding bivalves, by imposing a negative relationship between recruitment success and densities of adult conspecifics (Andre et al., 1993;Andre and Rosenberg, 1991;Comtet and Desbruyères, 1998;Vânia et al., 2014). It seems likely, therefore, that the high densities of adults in reference sites could lead to reduced rates of recruitment.
Chemical cues released by adult conspecifics are known to induce larval settlement and metamorphosis in several species of bivalves (Gosling, 2015;Porri et al., 2007). While the terms settlement and metamorphosis are often used interchangeably (Hadfield and Paul, 2001), they actually designate separate processes: metamorphosis is a definitive morphogenetic event; settlement, on the other hand, is reversible, as pediveligers and plantigrades of some species can swim up from the benthos and settle again in a new location (Bayne, 1964;Hadfield and Paul, 2001). Given the absence of E. castanea adults, chemical cues may be lacking from sites with elevated CO 2 , which could have contributed to delayed metamorphosis and/or increased postsettlement dispersal, while having no effect on larval settlement itself. Although these potential mechanisms were not investigated in our study, the simultaneous occurrence of high levels of intraspecific predation (in reference sites) and the lack of chemical cues from adult conspecifics (in sites affected by CO 2 ) could explain the observed trends in the abundance of postlarvae and juveniles.
According to Morton (1990), Ervilia castanea in the Azores reach sexual maturity at a shell length of 3.5 mm for females and 5.5 mm for males. The survival threshold in sites with elevated levels of CO 2 was clearly under that size range, suggesting that individuals at those sites are not producing offspring. Being an open system, however, larvae may originate from non-acidified sites, and only after settlement become subjected to the effects of elevated CO 2 (Dupont and Pörtner, 2013;Kroeker et al., 2013b). Overall, and even though reproduction was not investigated in our study, these results suggest that populations of E. castanea in high CO 2 sites are probably sink populations.

Conclusions
One of the advantages of conducting ocean acidification research using natural CO 2 gradients is that the responses observed are not restricted to the direct effect of carbonate chemistry, but also include the indirect consequences of trophic interactions and changes in habitat complexity (Fabricius et al., 2015;Hale et al., 2011;Range et al., 2012). In the future, transplant experiments including measurements of mineralization and dissolution of the shell, physiological, metabolic and genetic endpoints could enlighten and complement the responses observed in this study. There is also a need to consider other stressors (e.g., temperature, nutritional conditions, pollution), as future ecosystem changes will be driven by the synergistic effects of multiple stressors, instead of their isolated action (Baggini et al., 2014;Hale et al., 2011;Kroeker et al., 2013b). Future work should also include environmental benthic variables, such as sediment pH and temperature, as these can significantly vary from water column values (Clements et al., 2016;Dashfield et al., 2008).
We have shown that shallow-water volcanic CO 2 emissions in the Azores can be used to assess the effects of ocean acidification in central Atlantic benthic ecosystems, complementing similar studies in the region (Viotti et al., 2019) and in other parts of the world (Fabricius et al., 2015;Hall-Spencer et al., 2008;Harvey et al., 2018). Our results also emphasized the importance of food availability (i.e., chlorophyll a) in mediating the effects of CO 2 on filter feeders. Although Ervilia castanea is not commercially exploited, these findings raise concerns about the effects of ocean acidification on coastal bivalve populations with similar life history traits and on coastal economies that depend on them for food provision or commercial exploitation.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.  [FWO, grant number: 1242114N]. Sampling was done under the permit CCPJ015/ 2014, issued by the Regional Government of the Azores. We are grateful to Ana Navarro-Campoy, Joana Cruz, Lucia Moreno, Ana Amaral, Suzana Nicolau, Nélson Baião, Domingas Quiatuhanga, "Bombeiros Voluntários de Vila Franca do Campo" and the boats crews in Ponta Delgada ("Toninha Pintada" and "Cachalote I") and Horta ("Pintado" and "Águas Vivas"), for their support during this study. Comments by Dr. Jeff Clements and two anonymous referees substantially improved the original manuscript.