Ocean acidification effects on mesozooplankton community development: Results from a long-term mesocosm experiment.

Ocean acidification may affect zooplankton directly by decreasing in pH, as well as indirectly via trophic pathways, where changes in carbon availability or pH effects on primary producers may cascade up the food web thereby altering ecosystem functioning and community composition. Here, we present results from a mesocosm experiment carried out during 113 days in the Gullmar Fjord, Skagerrak coast of Sweden, studying plankton responses to predicted end-of-century pCO2 levels. We did not observe any pCO2 effect on the diversity of the mesozooplankton community, but a positive pCO2 effect on the total mesozooplankton abundance. Furthermore, we observed species-specific sensitivities to pCO2 in the two major groups in this experiment, copepods and hydromedusae. Also stage-specific pCO2 sensitivities were detected in copepods, with copepodites being the most responsive stage. Focusing on the most abundant species, Pseudocalanus acuspes, we observed that copepodites were significantly more abundant in the high-pCO2 treatment during most of the experiment, probably fuelled by phytoplankton community responses to high-pCO2 conditions. Physiological and reproductive output was analysed on P. acuspes females through two additional laboratory experiments, showing no pCO2 effect on females' condition nor on egg hatching. Overall, our results suggest that the Gullmar Fjord mesozooplankton community structure is not expected to change much under realistic end-of-century OA scenarios as used here. However, the positive pCO2 effect detected on mesozooplankton abundance could potentially affect biomass transfer to higher trophic levels in the future.

1 Introduction

Continuous burning of fossils fuels is causing an increase of atmospheric carbon dioxide (CO2), and current atmospheric pCO2 values (ca. 400 μatm) are projected to reach levels of up to 1000 μatm in less than 100 years [1]. Approximately one-third of the anthropogenic CO2 has been taken up by the oceans [2] leading to a reduction in pH (hence the term “ocean acidification” [3, 4]) and shifts in seawater carbonate chemistry [5]. Coastal marine ecosystems may be less sensitive to increased CO2 than open ocean regions, as the natural CO2 fluctuation in these areas is already substantial [1, 6]. However, ocean acidification (OA) can interact with other natural and anthropogenic environmental processes such as warming [7], eutrophication [8], and deoxygenation [9], making it a potential threat in conjunction with other stressors. Furthermore, OA may affect zooplankton not only directly by decreases in pH, but also indirectly via trophic pathways [10-12]. Consequently, both direct pH as well as pCO2 effects on primary production [13] may travel up the food web [10] therefore altering ecosystem functioning and community composition (e. g. [14]).

Elevated pCO2 in seawater may have positive effects on primary production, but at the same time impact marine organisms both via changes in calcification rates [15, 16], and via disturbance to acid–base (metabolic) physiology [17]. Calcified secretions in marine fauna and flora are not only limited to skeletal CaCO3 (thus, calcifiers sensu stricto) but there are other calcium-based structures that might be a target for low pH effects, such as, for example, the equilibrium organs (statoliths) in gelatinous zooplankton [17]. These organs are calcium magnesium phosphate crystals which may be affected by lowering pH [18], as reported for statoliths of scyphomedusae [19].

Copepods are the most abundant marine planktonic metazoans and, together with microzooplankton, are the major primary consumers in most marine food webs, sustaining secondary consumers such as fish and jellyfish [20, 21]. Copepods typically prefer larger and moving prey, i.e. they feed primarily on ciliates and dinoflagellates than on diatoms [22, 23], with preferred sizes between 20 and 200 μm ([24] and the references therein). As a result, they often switch from phytoplankton to microzooplankton over the course of a phytoplankton bloom [22] as larger prey items typically only become available later in the phytoplankton bloom, and even predate their offspring when resources are scarce [25].

Previously, copepods were considered to be relatively tolerant to OA [26, 27], but several processes in copepods may in fact be affected by low pH, including metabolism [28], pH balance [29], reproduction [30], development [31], growth [32] and survival [33]. Furthermore, diverse sensitivities to OA exist between different species and even between life stages within species [34]. Early life stages are most sensitive, resulting in a potential negative effect on survival and/or development (e. g. [29, 30, 35]). Different sensitivities to OA might also be related to copepod habitats, thus those copepod species more exposed to natural pH fluctuations (as vertical migrators or coastal species) might be more tolerant to OA than others [33, 36].

During the last decade, numerous studies dealing with the potential effects of high CO2 on single species were published (e. g. [35, 37]), while ecosystem-level impacts have attracted less attention. In order to assess future OA effects on natural communities, studies focused on ecological interactions (e. g. [38-41]), as well as long-term multigenerational experiments [42-44] are of paramount importance. To investigate the effects of end-of-century pCO2 levels on coastal pelagic ecosystems, we conducted a long-term mesocosm experiment in a boreal fjord. The present paper is part of the BIOACID II long-term mesocosm study PLoS Collection [45]. Here we focus on the natural mesozooplankton community, in particular on copepods and hydromedusae as the most abundant taxa. Testing the null hypothesis of no-effect, we assessed (1) mesozooplankton community development along the winter-to-summer plankton succession and the OA effects on the community interactions as well as (2) temporal trends and high-CO2 effects on species abundances, supported by two onshore experiments in the case of the most abundant copepod species, Pseudocalanus acuspes.

2 Materials & methods

2.1 Mesocosms setup and experimental design

Within the framework of the BIOACID II project (Biological Impacts of Ocean ACIDification), this study was part of the”BIOACID II long-term mesocosm study”, which was conducted from January to July 2013 in the Gullmar Fjord (58°15’ N, 11°28’ E), on the Swedish Skagerrak coast [45]. We deployed ten mesocosms (KOSMOS, M1-M10: “Kiel Off-Shore Mesocosms for future Ocean Simulation”, [46, 47]) in the fjord to study the effect of changing carbonate chemistry conditions on mesozooplankton community development. The experimental units consisted of large enclosed water volumes (~50 m3), five of them used as controls (ambient pCO2 levels = ca. 380 μatm), and the other five were CO2-enriched in levels adjusted to realistic end-of-century scenarios (RCP 6.0 [1]). Mesocosms were sealed by sediment traps, installed at the bottom of each mesocosm bag. Target pCO2 was reached at the beginning of the experiment by adding CO2 saturated seawater to the mesocosms. Subsequent additions were made on a regular basis in the course of the experiment (day 17, 46, 48, 68 and 88) to compensate for CO2 loss through outgassing. We established realistic end-of-century pCO2 levels (average = ca. 760 μatm) over the study period (see Fig 1a, [45]). Regular sampling every 2nd day included CTD casts, water column sampling, and sediment sampling. Water column samples were collected with integrating water samplers (IWS, Hydrobios), which collect a total volume of 5 L from 0–17 m depth evenly through the water column. This water was used for nutrient analyses, pigment analysis, and microzooplankton microscopy. All analyses are described in detail in [45] within this PLoS Collection. Briefly, nutrient (NO3-+ NO2-) concentrations (Fig 1b, [45]) were measured with a SEAL Analytical QuAAtro AutoAnalyzer and a SEAL Analytical XY2 autosampler. Pigment extracts were used for analysis by means of reverse phase high performance liquid chromatography (HPLC) (Fig 1c, [45]). Every eight days, microzooplankton samples were taken from the IWS carboys, immediately fixed with acidic Lugol’s solution and stored dark until identification (Fig 1d, [48]). Results presented here correspond to t1 (10th March) up to t103 (20th June) of the 113 days that the mesocosms experiment lasted [45].

Fig 1

Abiotic and biotic factors potentially affecting mesozooplankton community along the experiment.

a) in situ pCO2 levels, b) nutrients (NO3-+ NO2-), c) chlorophyll a, and d) microzooplankton abundances (ciliates and heterotrophic dinoflagellates). Colour code: red = treatment (~760 μatm pCO2), blue = control (ambient conditions). Solid lines = mean values; striped area = standard error of the mean.

2.2 Mesozooplankton sampling

The mesozooplankton community was sampled in the mesocosms and the fjord by vertical net hauls with an Apstein net (55μm mesh size, 17 cm diameter) equipped with a closed cod end, sampling a total volume of 385 L. Sampling depth was restricted to the upper 17m to avoid resuspension of the material accumulated in the sediment traps, at 20m depth. One net haul per mesocosm was taken once every eight days, within a narrow time-window (1 to 3 p.m.) to avoid differences in the community composition caused by diel vertical migration. Note that sampling frequency was lower than for other water column samples to avoid overharvesting of the plankton community. Samples were rinsed on board with filtered sea-water, collected in containers and brought to the laboratory, where samples were preserved in 4% formaldehyde buffered with sodium tetraborate. For transportation during summer time, the samples were placed in cooling boxes until fixation of the organisms.

During analysis, organisms were sorted using a stereomicroscope (Olympus SZX16) and classified to the lowest possible taxonomical level, including gender in the case of adult copepods. Copepodites and adults were classified to species level whereas nauplii from different species were pooled together. Taxonomical analyses were carried out focusing on copepods [49-52] and hydromedusae [53-55] as the most abundant groups. Every sample was sieved through 50 μm mesh, rinsed with tap water and poured into a calibrated beaker, where organisms were well mixed before taking a 5% aliquot with a Hensen Stempel pipette [56]. Counting was restricted to 5% (one aliquote) or 10% (two aliquots) of the total sample for the most abundant groups (nauplii, P. acuspes adults and P. acuspes copepodites) when more than 200 individuals were counted in the first aliquot. Otherwise the subsampling procedure was repeated, counting a maximum of a 15% of the total sample for all species.

Since some organisms characteristic to a winter-to-summer succession might not have been included when the experiment started, the community within the mesocosms was enriched by the addition of 22 L of fjord water every fourth day [45]. Likewise Atlantic herring (Clupea harengus) eggs and green sea urchin (Strongylocentrotus droebachiensis) gastrulae were artificially added to each mesocosms on t48 and t56 respectively [45] according to the time of the year that these groups would have been part of the natural fjord community. Densities of herring eggs introduced in the mesocosms were ~70–108 eggs per m3 and peak egg-hatching was estimated to occur around t63, with a final number of 1608 ± 237 hatched larvae per mesocosms, i. e. ~27–37 larvae per m3 [57]. These larval densities are within the natural range for the North Sea [58]. Sea urchin gastrulae were obtained in the onshore laboratory, introduced in the mesocosms (~110 sea urchin gastrulae per m3) and subsequently monitored from the mesozooplankton net tows on a weekly basis. An in depth analyses of Atlantic herring and green sea urchin larvae development are provided by Sswat et al. [57] within the framework of this PLoS Collection and Dupont et al. (unpubl. data).

2.3 P. acuspes condition experiments

Copepods were the most abundant group within the mesozooplankton community during the whole experiment, and the calanoid copepod P. acuspes was the most abundant species. To gain insights in P. acuspes’ physiological response to simulated OA we conducted two additional incubation experiments during the pre-bloom (March, t19) and senescence phase (May, t59) of the phytoplankton community (Fig 1). Every mesocosms was sampled by an extra net haul (see 2.2), and P. acuspes females were sorted immediately and subsequently incubated in a cold room adjusted to the average in situ temperature (t19: 3°C and t59: 5°C [45]) for offspring viability monitoring (n = 12) and respiration measurements (n = 5), or preserved for carbon content analyses (n = 20). Normally swimming females with undamaged eggs (60 females per treatment) were selected and initial clutch sizes were noted prior incubation to assess hatching rates. We aimed to incubate 12 females per mesocosms (i. e., 60 females per treatment), but this was not achieved in all cases due to the scarcity of egg carrying females within some samples or due to mortality of the females after 24h. Considering that incubation in small volumes does not affect egg production [59], females were incubated for 48h in 6-well plates, one female per well, in starvation and simulated field temperature. No additional pCO2 treatment was necessary because the aim of this side experiment was to analyse the memory effects of increased pCO2 on females in the mesocosm rather than effects on the eggs themselves. Clutch size and survival of the females were recorded each day during the condition experiments. Prosome length of all incubated females was measured upon termination of the experiment.

Respiration rates of five egg-carrying females per mesocosm (i. e. 25 animals per treatment) were measured in the cold room. Females were transferred to 1.6 mL vials equipped with fluorescent O2 foil discs (PSt3 spots, PreSens Precision Sensing, Germany) and filled with seawater adjusted to the pCO2 levels from corresponding mesocosms, based on the immediately preceding carbonate chemistry measurements in the mesocosms [45]. Vials were then sealed with Teflon caps and O2 concentrations were measured at 0, 3, and 6 hours using a Fibox 3 optode system. Respiration rates were calculated by subtracting the average oxygen depletion rate measured in five controls from the oxygen depletion rate in the vials holding copepods, multiplying by vial volume and dividing by number of individuals in each vial. Prior testing of the optode system at 5°C showed a 2 min 95% reaction time, i.e. the period of time taken before the output reached within 5% of the final oxygen concentration value (as estimated by exponential regression). Therefore, at every sampling, oxygen concentrations were read for three minutes, and an average of values read during the last minute was used for calculations.

To analyse carbon content, 20 non-ovigerous P. acuspes females were sorted from each mesocosm sample (i. e. 100 animals per treatment). The females were briefly rinsed in Milli-Q water to remove the excess of salt, and preserved in pre-weighted tin cups, which were in time dried (60°C) and preserved in an desiccator until analysed. Weights were obtained with a microbalance (Sartorius SC2). A Vario MICRO cube CHN analyser (Elementar) was used to measure carbon content.

2.4 Statistical analysis

To study Gullmar Fjord’s mesozooplankton community we firstly calculated species diversity for every mesocosm, which were compared using general linear models (GLMs) to detect any differences among treatments (high-pCO2, ambient). Subsequently, we analysed total abundances and abundances from the most frequent mesozooplankton species using general additive mixed models (GAMMs) to analyse the effect of the treatments as well as temporal trends. We compared the development of the community between treatments by a non-metric multidimensional analysis (NMDS) followed by a similarity analysis (ANOSIM). Finally, focusing on the most abundant species in the mesocosms (P. acuspes), we compared productivity and females’ condition between treatments by using GLMs.

Mesozooplankton diversity in mesocosms was calculated by using the Simpson’s Diversity Index (D) for finite communities. This index ranges from 0 to 1, and it is adapted to the form 1-D for a more intuitive interpretation of the results, thus higher values indicate higher sample diversity. Males, females and copepodites of the same copepod species were pooled together. Nauplii were assumed to be P. acuspes since this species accounted for > 90% of the copepod abundance during the whole experiment. General linear models (GLMs) were fitted to the Simpson’s indices to determine the dependence of diversity 1-D on time and pCO2. Calculations of D were performed in the vegan package of the R environment [60].

A multivariate analysis (NMDS) was used to describe the changes in the mesozooplankton community throughout the mesocosm experiment. NMDS is an ordination technique which represents, in an n-dimensional space, the dissimilarities obtained from an abundance data matrix [61]. NMDS takes a rank based approach, being more robust to datasets like the one used here, but as a consequence all the information about the magnitude of distances is lost. NMDS is therefore useful to represent the dissimilarities, and assess the influence of the treatment in the evolution of the community. However, due to the lack of magnitude, this technique is not ideal to evaluate the influence of environmental gradients on community changes [62]. The treatment effect was assessed by using permutation tests on the community position in the NMDS space, by checking if the area of clusters formed by the treatment in the NMDS were smaller than randomized samples of the same size [62]. In a complementary approach, we applied an ANalysis Of SIMilarity (ANOSIM) test [63] as a post-analysis to compare the mean of ranked dissimilarities between treatments (high-pCO2, ambient) to the mean of ranked dissimilarities within treatments. This analysis tests the assumption of ranges of (ranked) dissimilarities within groups are equal, or at least very similar [64].

Only those species that were present in at least one of the mesocosms for more than nine sampling days (2/3 of the number of days sampled) were used for temporal trends and multivariate analyses. By this criterion, the species selected for the analyses were: the hydromedusae Aglantha digitale and Hybocodon prolifer, and the females, males and copepodites of the copepod species Oithona similis, Temora longicornis, and P. acuspes. The aggregated copepod nauplii (pooled in one group and not identified to species level) were also included in these analyses.

To describe the temporal trends of each species during the mesocosm experiment we used GAMMs [61, 65] with a Poisson distribution and with a logarithmic transformation. Four different kinds of models were fitted to each abundance group (Table 1). Each of these models allowed the temporal trends to vary differently between treatments, representing (a) no difference between treatments (α + f), (b) differences in temporal trends but not in abundance (α + f) (c) difference in absolute abundance but not in temporal trends (α + f) and (d) difference both in absolute abundance and temporal trends (α + f). In this way potential differences between pCO2 and ambient mesocosms could be detected as either increase/decrease of overall abundance or changes in phenology. All models were fitted with an autocorrelation structure of first order to account for temporal autocorrelation in the data, and the specific mesocosm was used as a random intercept as the focus of the analyses was not the differences between mesocosms, but between treatments [61]. The models were compared by means of the Akaike Information Criterion (AIC). AIC takes into account both the goodness of fit of the model and model complexity, with lower AIC values indicating models with a better ratio between the explained variance and the number of variables [65]. For each species, the model with the lowest AIC was considered to better represent the temporal trends during the experiment, while avoiding overfitting the data.

Table 1

Generalized Additive Mixed Model (GAMM) structures.

α + f	Temporal trend and absolute abundances are treatment-independent (Model Trtmt_indep)
α + fT	Temporal trends depend on the treatment, but absolute abundances are treatment independent (Model Trtmt_trend)
αT + f	Absolute abundances depend on the treatment, temporal trends are treatment independent (Model Trtmt_absAb)
αT + fT	Both absolute abundances and temporal trends are affected by the treatment (Model Trtmt_absAb_trend)

In the case of copepods, we analysed the effects of the end-of-century pCO2 treatment on P. acuspes productivity by estimating a nauplii-to-adult ratio. Afterwards, GLMs were fitted to these ratios. The differences in the physiological and reproductive condition of P. acuspes females were analysed by GLMs comparing the potential effect of treatment and month in respiration rates, carbon content, prosome length, clutch size and hatching success. The effect of the time of the year (March and May), treatment and their interaction was considered in the models.

We used R (version 3.0.2, [66]) to fit abundances data with the GAMMs and GLMs. The significance level for all statistical analysis was set to p < 0.05.

3 Results

3.1 Mesozooplankton community: Composition, diversity and development

The mesozooplankton community comprised 27 different species and taxonomic groups (for a complete taxon list, see Table 2). The morphological classification of the most abundant groups (copepods and hydromedusae) was consistent with the genetic analyses conducted during the experiment (see [55] for more details). Copepods were the most abundant group throughout the experiment, representing 93–97% of the total abundances. P. acuspes was the dominant species in terms of abundance; based on the sum of adults and copepodites, P. acuspes represented 99.9% of the total copepod population at the beginning of the experiment and 33.6% at the end. Together with P. acuspes, only two other copepod species (T. longicornis, O. similis) and two hydromedusae (A. digitale, H. prolifer) were regularly recorded in our quantitative analyses. Other copepods and hydromedusae, polychaetae, chaetognatha, and appendicularians, as well as echinodermata, pteropoda, fish (larvae, eggs), bivalvia, cirripedia, and cladocera were rare (counted in less than 2/3 of the number of days sampled) or very rare (recorded in less than 3 sampling days during the experiment) in the studied community.

Table 2

Complete list of species and taxa present in the mesocosms registered throughout the study period.

Based on our records, species were classified as common (recorded on at least 9 sampling days, hence used for the GAMM analyses), rare (counted on 3 to 9 sampling days) or very rare (on less than 3 sampling days). C = common, R = rare, VR = very rare.

	Taxonomic groups	Records
1	Aglantha digitale	C
2	Hybocodon prolifer	C
3	Sarsia tubulosa	VR
4	Rathkea octopunctata	VR
5	Obelia sp.	VR
6	Phialella quadrata	VR
7	Bivalvia	VR
8	Pteropoda	R
9	Polychaeta	R
10	Evadne sp.	R
11	Podon sp.	R
12	Copepod nauplii	C
13	Pseudocalanus acuspes	C
14	Temora longicornis	C
15	Oithona similis	C
16	Acartia clausi	R
17	Tisbe sp.	R
18	Centropages cf. hamatus	R
19	Calanus sp.	VR
20	Monstrilla sp.	VR
21	Ectinosoma sp	R
22	Parasagitta elegans	R
23	Cirripedia	R
24	Ophiopluteus larvae	VR
25	Sea urchin larvae and juveniles	R
26	Oikopleura dioica	R
27	Teleostei (fish larvae)	VR

Mesozooplankton abundances (Fig 2A) increased after the first phytoplankton built-up (t17), and decreased during the phytoplankton post-bloom phase (t41-t77) and before microzooplankton increase (t81) (Fig 1C and 1D). GAMM analysis showed a treatment effect in total mesozooplankton abundances, which were higher under acidification scenarios (Trtmt_abdAb, Table 3). Averaged total catch (M1-M10) at the beginning of the experiment (t1) was 14571 ± 2857 individuals per m3, reached maximum in t49 (136342 ± 24451 individuals per m3), to decrease until minimum levels at t103 (9497 ± 3111 individuals per m3). Mesozooplankton biodiversity (1-D) was low during the experiment (Fig 2B), with average values of 0.094 ± 0.018 in ambient conditions and 0.098 ± 0.043 in the high-pCO2 mesocosms. No differences between ambient conditions and high-pCO2 treatment were observed (non-significant effect of treatment in a GLM). Independently from the pCO2 treatment, Simpson’s index (1-D) stayed below 0.1 in both treatments until t81. Then the index increased, with maxima on t103 (0.552 ± 0.045 in ambient and 0.535 ± 0.126 in high-pCO2, respectively).

Fig 2

Mesozooplankton community.

A) Mesozooplankton abundances. Solid lines = prediction from Generalized Additive Mixed Models (GAMMs) (smoother trends p-value < 0.05) with ambient and high-pCO2 mesocosms separately; striped area = confidence interval. B) Simpson’s Diversity Index (1-D) in relation to pCO2 levels within the mesocosms along the study period. Symbols and colours (blue = ambient; red = high-pCO2 treatment) identify each mesocosm.

Table 3

Mesozooplankton community models selection.

Generalized Additive Mixed Models (GAMMs) for the mesozooplankton community: a) α + f, no difference between treatments (Model Trtmt_indep), b) α + f, pCO2 treatment effect on temporal trends but not in abundance (Model Trtmt_trend), c) α + f, pCO2 treatment effect on absolute abundance but not on temporal trends (Model Trtmt_absAb) and d) α + f, treatment causes differences both in absolute abundance and seasonal trends (Model Trtmt_absAb_trend). Only those species that were present in at least one of the mesocosms more than 9 days (2/3 of the number of days sampled) and only convergent models were used for this analyses. The smoother of all selected models had a p-value < 0.05. For each species, the model with the lowest AIC (boldface) was considered to better represent the temporal trend during the experiment. Hyphens (-) indicate non-convergent models.

Taxa	Model type	R2	AIC	Taxa	Model type	R2	AIC
nauplii	Trtmt_indep	0.855	257.797	T. longicornis copepodites	Trtmt_indep	0.123	544.681
	Trtmt_trend	0.855	278.645		Trtmt_trend	0.127	540.113
	Trtmt_absAb	0.859	258.568		Trtmt_absAb	0.169	544.147
	Trtmt_absAb_trend	0.854	279.925		Trtmt_absAb_trend	0.122	536.422
P. acuspes ♀	Trtmt_indep	0.441	189.89	O. similis ♀	Trtmt_indep	0.558	463.501
	Trtmt_trend	0.491	195.135		Trtmt_trend	0.583	445.861
	Trtmt_absAb	0.443	191.887		Trtmt_absAb	0.552	465.903
	Trtmt_absAb_trend	0.5	197.739		Trtmt_absAb_trend	0.582	448.497
P. acuspes ♂	Trtmt_indep	0.564	282.254	O. similis ♂	Trtmt_indep	0.605	484.982
	Trtmt_trend	0.586	307.326		Trtmt_trend	0.635	482.307
	Trtmt_absAb	0.573	283.754		Trtmt_absAb	0.599	482.24
	Trtmt_absAb_trend	0.586	310.298		Trtmt_absAb_trend	0.633	479.176
P. acuspes copepodites	Trtmt_indep	0.727	210.277	O. similis copepodites	Trtmt_indep	0.767	447.67
	Trtmt_trend	0.752	232.495		Trtmt_trend	0.759	469.749
	Trtmt_absAb	0.76	209.844		Trtmt_absAb	0.766	449.509
	Trtmt_absAb_trend	0.75	234.226		Trtmt_absAb_trend	0.758	471.615
T. longicornis ♀	Trtmt_indep	-	-	A. digitale	Trtmt_indep	0.118	735.989
	Trtmt_trend	-	-		Trtmt_trend	0.114	734.663
	Trtmt_absAb	0.044	635.237		Trtmt_absAb	0.11	736.248
	Trtmt_absAb_trend	0.197	668.866		Trtmt_absAb_trend	0.11	739.801
T. longicornis ♂	Trtmt_indep	0.157	614.175	H. prolifer	Trtmt_indep	0.083	811.073
	Trtmt_trend	-	-		Trtmt_trend	0.151	764.543
	Trtmt_absAb	0.148	615.588		Trtmt_absAb	0.19	812.093
	Trtmt_absAb_trend	0.069	614.303		Trtmt_absAb_trend	0.173	764.455
Total catch	Trtmt_indep	0.852	92.57
	Trtmt_trend	0.867	104.36
	Trtmt_absAb	0.868	91.95
	Trtmt_absAb_trend	0.866	106.35

The 2-dimensional representation of the community did not show different patterns between treatments (Fig 3). Permutation tests (with 999 permutations) did not show the areas (i. e. clusters of samples) representing the treatment to be significantly smaller than randomized areas, indicating no treatment effect in the ordination. On the contrary, areas representing the sampling day (Fig 3) were significantly smaller than randomized areas using the same test. This result indicates clear community differences throughout the study period. Results from the ANOSIM test (p-value = 0.322) matched with the NMDS, suggesting that there was no significant difference between the community development under the high-pCO2 treatment and the ambient conditions.

Fig 3

Non-metric Multidimensional Scaling analysis (NMDS) of the mesozooplankton community (stress value = 0.17).

Colour code: red = treatment (~760 μatm pCO2), blue = control (ambient conditions). Sampling days represented as t-day; lines represent patterns. The underlying data implemented in the analysis are shown in Fig 1.

3.2 Species abundances

Temporal trends of the selected species were analysed by using GAMMs (Figs 4 and 5; Table 3). The model selection procedure discerned whether there was a difference in the temporal trends and abundances in between the two different treatments (i.e. high or ambient pCO2).

Fig 4

Copepod abundances along the study period.

A) P. acuspes females, B) P. acuspes males, C) P. acuspes copepodites, D) T. longicornis females, E) T. longicornis males, F) T. longicornis copepodites, G) O. similis females, H) O. similis males, I) O. similis copepodites, J) nauplii. Colour code: red = treatment (~760 μatm pCO2), blue = control (ambient conditions). M = mesocosms. Solid lines = prediction from Generalized Additive Mixed Models (GAMMs) (smoother trends p-value < 0.05) with the ambient and high-pCO2 mesocosms shown separately; striped area = confidence interval. Black lines indicate that the prediction of the model for high-pCO2 treatment and ambient conditions are the same.

Fig 5

Hydromedusae abundances along the study period.

A) A. digitale, B) H. prolifer. Colour code: red = treatment (~760 μatm pCO2), blue = control (ambient conditions). M = mesocosms. Solid lines = prediction from Generalized Additive Mixed Models (GAMMs) (smoother trends p-value < 0.05), with the ambient and high-pCO2 mesocosms shown separately; striped area = confidence interval.

There was no pCO2 effect on the abundance of adult P. acuspes and T. longicornis but copepodite stages of both species responded to increased pCO2. P. acuspes adults did not show differences in abundances nor in temporal trends between treatments (Table 3 Trtmt_indep for both males and females; Fig 4A and 4B). However, the absolute abundance of P. acuspes copepodites differed between treatments, being higher under the high-pCO2 treatment (Table 3 Trtmt_absAb; Fig 4C). Abundance of T. longicornis adults did not show a difference between treatments (Fig 4D and 4E); even though the selected model showed slightly higher abundances of T. longicornis females in the high-pCO2 mesocosms (Table 3 Trtmt_absAb; Fig 4D), the confidence intervals of the modelled abundances were overlapping throughout the study period. This indicates that the difference were small, and probably caused by extreme values at the end of the experiment. Only T. longicornis copepodites (Table 3 Trtmt_absAb_trend; Fig 4F) showed different absolute abundances and a different temporal trend between treatments, being more abundant in the ambient pCO2 mesocosms, particularly during the last 20 days of the study. O. similis adults negatively responded to the elevated pCO2 conditions with an earlier abundance decrease towards the end of the experiment (Fig 4G and 4H). In case of O. similis males the absolute abundance and the temporal trend were negatively affected by the high-pCO2 treatment (Table 3 Trtmt_absAb_trend). However, this effect was not detected on O. similis copepodites (Table 3 Trtmt_indep; Fig 4I), which showed no significant difference between both treatments. Copepod nauplii, the most abundant group in the mesozooplankton (Fig 4J), did not show a difference in temporal trends nor abundance between treatments (Table 3 Trtmt_indep). When analysing abundances in certain time-points, we could detect different pCO2 effects that were not detected by the GAMMs. In the case of P. acuspes, adult copepods were significantly more abundant on t81 (t-test, p-value = 0.010), but the effect disappeared afterwards. Different responses were also observed on nauplii abundances, which were significantly higher under high-pCO2 conditions between t49 and t65 (t-test, p-value = 0.03), whilst we did not detect differences in abundances between treatments when analysing abundances from t65 until the end of the experiment (t-test, p-value = 0.622).

In the case of both hydromedusa species, we also detected species-specific pCO2 effects (Fig 5, Table 3). Under the high-pCO2 treatment, H. prolifer abundance was lower; the GAMM detected an effect not only on the temporal trend, but also on the abundances of this species (Table 3 Trtmt_absAb_trend). The model representing A. digitale also showed a different temporal trend between treatments (Table 3 Trtmt_trend) despite of the confidence intervals overlapping of both patterns.

To sum up, after analysing the abundance of each species under high-pCO2 conditions during the whole study period we observed positive (P. acuspes copepodites, A. digitale), negative (T. longicornis copepodites, H. prolifer, O. similis adults) and no effects of elevated pCO2 (nauplii, P. acuspes and T. longicornis adults, O. similis copepodites). It is worth mentioning that the predictive power (R2) of these models was low in some cases (see Table 3) due to the complete absence of some species in some mesocosms. However, the models represented well the overall trend differences between treatments (Figs 4 and 5). Differences between treatments were at times significant for specific time periods.

3.3 P. acuspes: Productivity and females’ condition

Copepod productivity was assessed by computing the ratio between nauplii and adults for the most abundant species, P. acuspes. We calculated the nauplii-to-adult ratio from t17 until the end of the experiment, since the fraction < 200 μm was preserved only from t17 on. At a significance level of 0.05, no differences in this ratio between the ambient and high-pCO2 treatment (GLM, p-value = 0.576), but a significant effect of time (GLM, p-value < 0.001) was detected. Productivity increased from the beginning of the experiment until t65 or t73 independently of the pCO2 treatment (see Fig 6), and rapidly decreased afterwards. A second increase in the productivity was detected from t97, with the highest ratios in some of the high-pCO2 mesocosms.

Fig 6

P. acuspes productivity in relation to pCO2 levels along the study period.

Symbols and colours (blue = ambient; red = high-pCO2 treatment) identify each mesocosm. Production estimated as the ratio between nauplii and adults. P. acuspes nauplii abundances were estimated from the relative abundances of P. acuspes in relation to total copepod abundances per sampling day and mesocosm.

Regarding the P. acuspes females’ condition, none of the physiological and reproductive parameters investigated (respiration, carbon content, prosome length, clutch size, hatching success) showed a significant difference between treatments, nor in the interaction between month and treatment (p-value > 0.05; Fig 7, Table 4). However, significant differences between the first (March, t19: first phytoplankton bloom) and the second experiment (May, t59: second phytoplankton bloom) were observed. Respiration rate (Fig 7A) was lower during May compared to March (p-value = 0.001). Females’ carbon content and prosome length, as well as the hatching success after 48h incubation (Fig 7B, 7C and 7E) were not different between months, nor between pCO2 conditions. Yet, at the beginning of the incubations (0h), clutch size (Fig 7D) was significantly higher in May (p-value = 0.021). None of the interactions between pCO2 treatment and month rendered in a significant effect on the studied variables.

Fig 7

P. acuspes females’ condition.

General Linear Models (GLMs) comparing the potential pCO2 effect on P. acuspes females: A) respiration rate, B) carbon content, C) prosome length, D) clutch size at the beginning of the incubation (0h), E) hatching success after 48h incubation. Error bars represent standard deviation. Colour code: red = treatment (~760 μatm pCO2), blue = control (ambient conditions). March = t19 (first phytoplankton bloom), May = t59 (decline phase of the second phytoplankton bloom).

Table 4

Results from P. acuspes females’ condition experiment.

Generalized Linear Models (GLMs) based on two laboratory experiments (March, May), n = 120 females per experiment. Boldface represent p-values < 0.05.

Respiration	Estimate	Std.Error	t-value	p-value
(Intercept)	5.035	0.786	6.406	0
pCO2 treatment	0.553	0.37	1.492	0.154
month	-0.786	0.185	-4.246	0.001
Carbon content
(Intercept)	5.586	0.958	5.829	0
pCO2 treatment	0.541	0.452	1.198	0.247
month	-0.056	0.226	-0.246	0.808
Prosome length
(Intercept)	0.728	0.039	18.875	0
pCO2 treatment	-0.016	0.018	-0.895	0.383
month	-0.005	0.009	-0.536	0.599
Clutch size (0h)
(Intercept)	2.394	1.103	2.17	0.044
pCO2 treatment	-0.563	0.52	-1.082	0.294
month	0.661	0.26	2.542	0.021
Hatching success
(Intercept)	11.465	9.875	1.161	0.262
pCO2 treatment	0.275	4.655	0.059	0.954
month	-0.823	2.328	-0.354	0.728

4 Discussion

During this winter-to-summer experiment on the effect of ocean acidification on plankton communities, we did not detect an effect of pCO2 on either the diversity of the mesozooplankton community, nor on its development as a whole. At first sight, this may seem surprising as some taxa showed a response to OA, where others did not. The most parsimonious explanation for this apparent contradiction is the strong dominance of the copepod P. acuspes. As a result, changes in the relative composition of the community were small and were not be picked up by relatively coarse indicators such as Simpson’s Diversity or rank-based methods such as NMDS. Only on the last two sampling days, when P. acuspes abundances declined strongly, a trend towards a higher diversity under high-pCO2 conditions became visible (Figs 2B and 3), and the communities under the two treatments diverged (observed also for microzooplankton [48]). Potentially this indicates a long-term effect of high pCO2 on the communities, but this is impossible to say as, at that time the mesocosm set-up started to deteriorate and the experiment was terminated.

Unlike previous mesocosms studies focusing on the effect of OA on natural coastal plankton communities in the Arctic [67] and the Baltic [39], we detected a positive pCO2 effect on the total mesozooplankton abundance from Gullmar Fjord. This effect was mostly caused by the CO2-driven increase in the abundances of P. acuspes copepodites. This was somewhat unexpected, as previously work on the same species from the same location [42, 68] found significant negative pCO2 effects on egg production and metabolism. The two studies cited above were highly controlled laboratory experiments, where the copepods were cultured under uniform environmental conditions (except for the pCO2 treatments) and offered identical prey in all treatments. Thus, the effects observed were directly caused by changes in carbonate chemistry of the water as all other environmental factors were identical. In semi-natural experiments such as the one described here, these effects are easily masked, either through bottom-up effects (changes in the availability or quality of the food), or as a result of top-down effects (changes in predation rates). In our two condition experiments we excluded the latter effects, and focused on the effects of the overall growing conditions in the mesocosms. In contrast to the laboratory experiments cited above, we did not find significant differences in the physiological condition of P. acuspes females between ambient and high-pCO2 treatments (Fig 7). Secondary production in P. acuspes followed a temporal trend, with higher clutch sizes and nauplii abundances on t59 (May), responding to higher phytoplankton concentration (chla) and microzooplankton biomass. However, this increase in food quantity might not have been coupled with food quality to maintain the copepod population in the mesocosms, which increased from ~260 ± 5 copepods L-1 (t19) to ~1245 ± 32 copepods L-1 (t59). This could explain lower respiration rates in May than in March [69, 70]. Potential food items for copepods on t19 (March) consisted mainly of phytoplankton between 5 and 40 μm and microzooplankton biomass below 2 μg C L-1 before the first phytoplankton bloom in the mesocosms [48, 71]. On t59 the entire mesocosms system was dominated by Coscinodiscus concinnus (representing 47% of the biomass) and the nanophytoplankton fraction (accounting for 21%) [71], both largely outside the food spectrum of P. acuspes. Microzooplankton biomass was ~12 μg C L-1 on t59 [48], but might not have been enough to supply the whole P. acuspes population, so copepods might have searched for alternative food sources such as sinking material. In fact, the decrease in adults from t97 in all mesocosms matched high resolution images taken from sediment trap material, where high abundances of adult P. acuspes were found (Tim Boxhammer, pers. comm.). This observation suggest that, towards the end of the experiment, copepods might have migrated downward searching for food and stayed close to the sediment traps, as previously observed in a mesocosms experiment in a Norwegian fjord [72].

In view of the result of the two laboratory experiments, where we observed no effects of pCO2 on egg production, the most plausible explanation for the higher P. acuspes abundances under the high-pCO2 treatment is a community CO2-driven bottom-up effect [10, 12, 73]. This is not a contradiction, as in the laboratory experiments we specifically looked at the memory pCO2 effect on the clutch, which was not expected to be affected by the 48h food deprivation regime [74]. Thus, the higher abundance of P. acuspes copepodites was probably fuelled by phytoplankton community responses to high-pCO2 conditions during our mesocosms experiment. Higher primary production [75] and higher chla levels under high-pCO2 [45] resulted in higher copepodite abundances. Interestingly, this CO2-driven increase in copepodite abundances did not result in higher abundances of adults later in the season except on t81, when adult P. acuspes were significantly more abundant under high-pCO2 conditions. The most plausible explanation for this trend in adult P. acuspes abundance after t81 is, apart from the potential downward migration as indicated above, that the level of top-down control through herring larvae was different, with higher predation pressure in high-pCO2 mesocosms. As detailed in Sswat et al. [57], after hatching on ~t63, herring larvae would have gradually switched from endogenous to exogenous feeding, preying then firstly on nauplii and ciliates, afterwards increasing the size of their prey gradually with their own body size until they reached copepodites (~t65-t81) and finally adults (~t81-t105)[76-78]. From t77 (14th day post-hatching, DPH) survival of herring larvae was significantly higher in the high-pCO2 mesocosms [57], which would imply higher grazing pressures on P. acuspes. Since consumption rates of smaller larvae are much lower than those of larger ones, we would have only detected a top-down effect of the herring larvae on adult abundance at the end of the experiment. This, together with a more intensive feeding activity by herring larvae because of the higher larvae survival rates under the acidic treatment [57], could have caused lower abundances of adult P. acuspes relative to the opposite pattern in the copepodites.

In the case of T. longicornis, no effects of pCO2 were observed on the adults but copepodites were more abundant under ambient conditions, especially during the last 20 days of the experiment (Table 3, Fig 4D–4F). This finding fits to the last two sampling days divergence between treatments in the NMDS analysis (Fig 3), which points to a different development of the community under ambient and high-pCO2 conditions. The particular tolerance in T. longicornis female reproductive fitness to end-of-century pCO2 scenarios had already been described by McConville et al. [27]. However, the higher abundances of T. longicornis copepodites observed in ambient conditions suggest that this tolerance might be diminished in early life stages, as previously observed in other calanoid copepods [29, 79].

Our results suggest a negative effect of pCO2 on adult O. similis, which were more abundant under ambient conditions when considering the whole experimental period. The explanation for O. similis’ sensitivity to OA observed in adults might be in the life history of this copepod. According to Lewis et al. [33] there is a correlation between sensitivity to OA and vertical migration behaviour. Species that do not exhibit diel vertical migration behaviour (as O. similis) are typically less exposed to variation in pCO2 levels compared to other copepods and more prone to be sensitive to OA [33, 80]. For O. similis, these researchers detected reduced adult and naupliar survival under 700 and 1000 μatm pCO2. Our study would support this observation by lower O. similis adult abundances under high-pCO2 conditions. Towards the end of the experiment, however, we observed an increase in O. similis abundance, likely reacting to the increase in ciliates and dinoflagellates biomass [48]. Adults showed a significant reaction to OA with firstly higher and subsequently lower abundances in the high-pCO2 treatment. As also observed on adult P. acuspes, the differential decrease in adult O. similis within treatments from t81 might respond to herring larvae abundance and the size-dependent feeding activity [57, 77]. Thus considering that during the last two sampling days adults would probably be in the preferred size range for the herring larvae, the release in preying pressure on copepodites and the built-up of protozooplankton [48] might explain the final increase in copepodite abundance in both treatments.

Whilst the connection between jellyfish blooms (scyphomedusae, hydromedusae, siphonophores and ctenophores) and anthropogenic climate change remains unclear (e. g. [81, 82]), the effects of changing seawater carbonate chemistry on planktonic gelatinous species have been rarely tested. However, all results on different gelatinous zooplankton groups (schyphomedusa ephyrae [19, 83, 84], coelenterate records [85]) point to the tolerance of jellyfish to future changes in pCO2. In this study we showed for the first time the species-specific sensitivity of hydromedusae to OA. Thus H. prolifer (Anthomedusa) reacted negatively to high pCO2 by lower abundances, while A. digitale (Trachymedusa) was more abundant in the high-pCO2 treatment (Table 3, Fig 5). This result was unexpected, given the fact that A. digitale has statoliths, which could be a target for lower pH (as Richardson and Gibbons [85] also noted). Our findings suggest that hydromedusae with statoliths are not necessarily more sensitive than those without these calcium-based structures, and consequently hydromedusa statoliths might not be sensitive to OA, at least in realistic end-of-century scenarios. Further ecophysiological analyses, however, are still required for these and other hydromedusae species to confirm this hypothesis.

Conclusion

During this study, we observed species-specific sensitivities to pCO2 in copepods and hydromedusae abundance. In the case of copepods, responses to elevated pCO2 depended also on the life-stage of the individuals, copepodites generally being the most sensitive stage. Our results point that OA could positively affect the calanoid P. acuspes by a bottom-up effect in pCO2-fuelled food webs. Nonetheless, the effect of OA on single species was not detectable in the structure or diversity of this community, probably due to the overwhelmingly dominance of P. acuspes in the studied community. Hence, under a realistic end-of-century OA scenario, the Gullmar Fjord mesozooplankton community structure is not expected to change much, although it could well be that the OA effect on copepodites would potentially affect biomass transfer to higher trophic levels in the future.

Ethic statement

No specific permission was required for activities related to field sampling. The field location was not privately owned or protected, and neither endangered nor protected species were involved. Fish larvae experiment [57] was conducted under the ethical permission (number 332–2012 issued by the Swedish Board of Agriculture "Jordbruksverket"). Animal welfare was assured by minimization of stress from handling and treatment. Specimens were therefore anaesthetized before handling using Tricaine methanesulfonate MS-222. The CO2 concentrations used in this study are far below the lethal level.