Exposure to CO₂ influences metabolism, calcification and gene expression of the thecosome pteropod Limacina retroversa

Thecosome pteropods, a group of marine planktonic shelled gastropods, contribute substantially to food webs and both carbon and carbonate flux in many pelagic ecosystems (Berner and Honjo, 1981; Fabry and Deuser, 1991; Noji et al., 1997; Hunt et al., 2008; Manno et al., 2010). Serving as food for zooplankton, fish and whales, they also contribute an estimated ∼12% of the global calcium carbonate export flux (Fabry and Deuser, 1991; Bednaršek et al., 2012a). Thecosomes have recently received increased research interest as pelagic calcifiers as they are sensitive to anthropogenic increases in CO₂ concentration that result in declining ocean pH and calcium carbonate saturation state (ocean acidification, OA). Because of these changes, acidification has the potential to influence both physiological and biomineralization processes (Fabry et al., 2008; Doney et al., 2009).

OA has been shown to cause complex responses in various groups of marine calcifiers, but generally results in trends towards impaired biomineralization, survival and development (reviewed in Kroeker et al., 2013; Wittmann and Pörtner, 2013). Calcifier sensitivity seems to be partially attributable to increased energetic expenditures associated with biomineralization under less favorable saturation states (Waldbusser et al., 2015b). The particular sensitivity of pteropods is assumed to be a consequence of their thin aragonitic shells, which are highly responsive to under-saturation (Comeau et al., 2009; Lischka et al., 2011; Bednaršek et al., 2012b, 2014; Bednaršek and Ohman, 2015), and due to their prevalence in polar and upwelling regions, which are known to respond rapidly to anthropogenic influences (Hunt et al., 2008; Bednaršek et al., 2012a). Although capable of shell repair and shell accretion in under-saturated conditions (Comeau et al., 2010; Lischka et al., 2011; Peck et al., 2016), it can be assumed that increased energetic expenditures limit the success of pteropods over prolonged periods of exposure or during repair, as has been suggested for other calcifiers (Melzner et al., 2011). No explicit studies of energetics have been conducted in pteropods; consequently, as with many other organisms the broader physiological and ecological implications of OA exposure on thecosomes are unclear. For example, a few thecosome species respond to CO₂ exposure with changes in aerobic respiration, but respiration in others is minimally affected, unless a co-stressor is involved (Comeau et al., 2010; Maas et al., 2012; Seibel et al., 2012; Lischka and Riebesell, 2016).

To provide an understanding of the processes sensitive to acidification, full transcriptomic responses of various biomineralizing organisms have recently been investigated by high throughput RNA sequencing (RNA-seq; i.e. Moya et al., 2012; Li et al., 2016). These techniques have been applied to thecosomes (Koh et al., 2015; Maas et al., 2015; Moya et al., 2016), but due to differences in regional chemistry or treatment severity, they have explored very different types of CO₂ stress – some under-saturated with respect to aragonite, others not. Neural function, ion transport and energetic metabolism seem to be frequently influenced, but comparisons across studies remain difficult. These studies consistently show that exposure to elevated CO₂ results in changes to gene expression that are probably related to biomineralization, but the results are sometimes contradictory in direction of response.

Using paired assessments of respiration rate and gene expression, two previous pteropod studies suggested that despite a lack of measurable organismal level response, metabolism may be influenced by CO₂ exposure, indicated by altered expression of genes associated with the mitochondrial inner membrane complex (Maas et al., 2015; Moya et al., 2016). Gene expression profiling can thus reveal responses to OA prior to the emergence of gross physiological changes, increasing our sensitivity in detecting subtle effects or revealing resource reallocation. This is important as keeping healthy pteropods in captivity is difficult (Howes et al., 2014), making long-duration behavioral and physiological studies experimentally challenging. It is important, however, to attempt longer exposures as we seek to predict responses to the chronic stress that OA represents. Previous transcriptomic (Moya et al., 2012, 2015; Li et al., 2016) and physiological (i.e. Hochachka and Somero, 2002; Langenbuch and Pörtner, 2004) studies have consistently indicated that the duration and severity of CO₂ exposure can influence the mechanism of acidosis response. Although the application of paired molecular and organismal metrics in pteropods has begun to address the response of these organisms to CO₂ stress, several important unknown factors remain regarding the influence of severity and duration of exposure.

This study aims to address these specific questions using the species Limacina retroversa (Fleming 1823), a temperate to sub-polar pteropod species found year-round in the Gulf of Maine. This region has recently been identified as the most sensitive to acidification stress along the US East Coast (Wang et al., 2013). During certain seasons the deep basins already reach under-saturation with respect to aragonite, although the typical vertical distribution of the pteropods presently never intersects these corrosive waters (Wang et al., 2017). These surface waters are, however, predicted to reach under-saturation within 100 years (IPCC, 2013; Wang et al., 2017). In response to acidification L. retroversa veligers from the Gulf of Maine exhibit increased mortality and delayed development (Thabet et al., 2015), and adults show changes in shell condition and sinking speed (Bergan et al., 2017); the influence of OA on adult physiology has, however, not been thoroughly studied in this region. Limacina retroversa shells from the Arctic and Northern Norwegian Seas have been shown to be sensitive to CO₂ exposure (Lischka and Riebesell, 2012, 2016; Manno et al., 2012), but mortality, respiration rate and swimming were unaffected without co-exposure to other stressors, such as warming or low salinity.

In the present study, response metrics were interrogated following both short-term (1–3 days) and medium-term (7–21 days) exposure to enhanced CO₂. Three levels of CO₂ exposure, intended to achieve over-saturated, near-saturated and under-saturated conditions, were used. In addition to examining how duration and severity of exposure influence physiological processes, analyses of transcriptome-wide responses were paired with respiration and shell analyses to investigate how any changes at an organismal level corresponded to changes in gene expression.

Adult pteropods (parapodia fully developed) and water for animal maintenance were collected from the research vessel Tioga in the western Gulf of Maine off Provincetown, MA, USA (at ∼42°N and 70°W, ∼50–250 m water depth). Water was collected on 25 April 2014 from ∼30 m depth using a submersible pump and filtered with 63 µm mesh. For later animal maintenance some of this water was transferred into 1 litre jars, which were pre-chilled in a refrigerator at 8°C. The rest of the water was stored for transport in large (∼150 l) covered plastic cans. A large batch of this water (>400 l) was transferred on the evening of 25 April to a walk-in cold-room at the Woods Hole Oceanographic Institution. The cold-room was maintained at a temperature of 8°C, to approximate the temperatures at which animals are found during this season. This water was deposited into a holding tank where it recirculated through a 1 µm filter for the remainder of the experiment. The temperature of the cold-room was monitored continuously with a HOBO logger (Onset, Bourne, MA, USA).

After ∼12 h of filtration, batches of this water transported from the field collection site (hereafter in situ water) were transferred into three ∼100 litre covered pre-equilibration tanks and bubbled with three different gas concentrations. Local air (nominally ∼400 ppm) was compressed and controlled with in-line pressure valves to provide the ambient treatment. Gas concentrations of 800 or 1200 ppm CO₂ were generated with mass flow controllers (Aalborg, Orangeburg, NY, USA), combining compressed ambient air and CO₂ using the system described in White et al. (2013). These are referred to as the medium and high treatments, respectively. These concentrations were intended to yield over-saturated, near-saturated and under-saturated conditions with respect to aragonite.

On 27 April 2014, pteropods were collected via slow tows (ca 5–10 m min⁻¹ vertical and vessel speeds of 1–2 knots) of a Reeve net with 333 µm mesh and a large (∼20 l) cod end designed for gentle treatment of live organisms. Upon retrieval, L. retroversa were isolated into 1 litre glass jars at densities of 20–40 individuals l⁻¹ in chilled pre-filtered local seawater (see above). These jars were initially kept in an 8°C refrigerator and then transported in coolers with ice packs back to the laboratory cold-room facility on the evening of 27 April.

On the same evening that animals were returned to the laboratory, the pre-equilibrated in situ water was transferred into three, 12 litre glass experimental carboys per treatment (a total of nine carboys placed in a randomized pattern in the cold-room), where gas equilibration was continued using micro-bubbling stones. All collected animals were pooled, and then individually counted into the pre-bubbled experimental carboys in a randomized manner reaching a density of ∼20 individuals l⁻¹. After this transfer, pteropods were fed with a mixture of Rhodomonas lens and Heterocapsa triquetra. Animals were fed the same mix of algae once every 2 or 3 days and immediately after each water change. Water was changed at weekly intervals using in situ seawater from the holding tank following the same protocol of pre-bubbling each treatment, as above. Mortality was not explicitly measured, but appeared similar among treatments; survival was initially high, but decreased substantially in the second week.

To monitor the carbonate chemistry during the experiments, sub-samples of water (50 ml) from each carboy were taken every 2–3 days. The pH was measured spectrophotometrically as described in White et al. (2013) using 2 mmol l⁻¹ m-Cresol Purple indicator dye and an Ocean Optics USB4000 spectrometer with an LS-1 light source and a FIA-Z-SMA-PEEK 100 mm flow cell (Ocean Optics, Dunedin, FL, USA). pH is reported on the seawater scale (pH_sw). Salinity and seawater density were measured using a seawater refractometer (model 96822, Hanna Instruments, Woonsocket, RI, USA) at each pH reading. Temperature in the cold-room was 7.9±0.27°C (mean±s.d.) for the duration of the experiment (HOBO logger; Onset), while salinity of all treatments was consistently 34.

To more fully characterize the carbonate chemistry of the experimental treatments, discrete samples for analysis of dissolved inorganic carbon (DIC) and total alkalinity (TA) were also collected from the pre-equilibrated water at weekly intervals, prior to its transfer to the carboys (days 0 and 7) and from the water in each carboy just before it was replaced during water transfers (days 7 and 14). This provided both the starting and ending DIC and TA of each batch of water for weeks 1 and 2, and the starting conditions for week 3. These were measured following the methods of Wang and Cai (2004).

For DIC and TA measurements, water was sampled from carboys using a plastic syringe soaked with deionized (DI) water and a hose rinsed thoroughly with sample water three times. Sampling from the pre-equilibration tanks was done by submerging the glass bottle in the plastic can, with three rinses prior to the final fill. The sample jars were 250 ml Pyrex borosilicate glass bottles, and air headspace of about 1% of the bottle volume was left in each sample bottle to allow room for expansion. Each sample was poisoned with 100 µl of saturated mercuric chloride, capped with an Apiezon-L greased stopper, thoroughly mixed, and then tied with a rubber band over the glass stopper.

DIC samples were analysed on a DIC auto-analyser (AS-C3, Apollo SciTech, Bogart, GA, USA) via acidification, followed by non-dispersive infrared CO₂ detection (LiCOR 7000; Wang and Cai, 2004). Three aliquots were measured from each sample to obtain replicate measurements. The instrument was calibrated with certified reference material (CRM) from Dr A. G. Dickson at the Scripps Institution of Oceanography, and has a precision and accuracy of ±2.0 µmol kg⁻¹. Total alkalinity was measured using an Apollo SciTech alkalinity auto-titrator, a Ross combination pH electrode, and a pH meter (ORION 3 Star) based on a modified Gran titration method, which has a precision and accuracy of ±2.0 µmol kg⁻¹ (Wang and Cai, 2004). Additional carbonate chemistry parameters were estimated using concurrently collected DIC and TA pairs for ingoing or outgoing water at the weekly water changes, or using the more frequent pH_sw measurements with the TA pair from the closest water change. The saturation state of aragonite (Ω_A) and P_CO2 were calculated with CO2SYS software (Pierrot et al. 2006), using constants K₁ and K₂ by Dickson and Millero (1987), and the KHSO₄ dissociation constant from Dickson (1990).

At each sampling time point individuals were randomly selected from the ambient, medium and high treatments using a long plastic pipette. Respiration experiments were nominally conducted on days 1, 3, 7 and 14, but as measurements were conducted over 24 h periods, these experiments began ∼12 h and 2.5, 6.5 and 13.5 days after initial exposure. To allow for gut clearance, pteropods were first transferred to a treatment-specific 1 litre beaker of water containing filtered (0.2 µm) treatment-specific CO₂-equilibrated in situ water for 8–12 h. After this period, the nine healthiest looking individuals (actively swimming or with parapodia extended) were transferred to custom small volume (2–3 ml) glass respiration vials with fresh treatment-specific 0.2 µm filtered CO₂-equilibrated in situ water and an optically sensitive oxygen sensor (OXFOIL; PyroScience, Aachen, Germany). As a control for bacterial respiration, every fourth chamber was filled with water but was left without an animal. The oxygen concentration in each chamber was then measured non-invasively at the start of the experiment using a FireSting optical oxygen meter (PyroScience) following similar methods to Maas et al. (2016a). At the conclusion of each 24 h respiration experiment, the O₂ concentration was again measured in all chambers. The chambers were then weighed wet, emptied and weighed dry to more accurately estimate the water volume with a correction for density of 1.024 (measured with the seawater refractometer). Each organism was visually inspected to confirm survival and then was briefly rinsed with DI water, placed in a pre-weighed aluminium dish, and weighed on a Cahn microbalance (C-33; 1 µg precision). Final oxygen consumption rates were calculated (µmol O₂ wet h⁻¹) based on the change in oxygen between the final and initial oxygen measurements, corrected for the bacterial respiration from the control chambers (which proved to be low, ∼0.32 nmol O₂ h⁻¹). Statistical analyses used this corrected metabolic rate as the independent variable with wet mass as a covariate to account for size differences in respiration rate.

Solutions of calcein (50 mg l⁻¹) in in situ 0.2 µm filtered seawater were mixed in quart jars and then bubbled with CO₂ gas at the same concentration as the three treatments. Ten live L. retroversa were removed from each treatment on days 7, 14 and 21 and placed into the jars with the calcein solution for 1 h. During the time that a live pteropod is in calcein solution, the fluorescent pigment binds to calcium ions in cells and can be used to visualize the calcium carbonate that is being actively precipitated. The animals were then rinsed with DI water and imaged under a fluorescence and differential interference contrast microscope. The images were analysed in MATLAB using a custom code following similar methods to Bergan et al. (2017). Firstly, the shell region was identified from the background by making the image binary. Secondly, the average intensity of pixels comprising the shell in a grayscale version was calculated normalized to a maximum intensity of 255 to result in a scale of 0 (no fluorescence) to 1 (maximum brightness in all pixels). This calculation provides an estimate of the calcein-derived fluorescence and is used as an indicator of active calcification. It is important to note, however, that in L. retroversa, if the period of photography is less than 24 h after staining (as was the case in this experiment) other regions with high calcium ion concentration, such as the gut and muscle, also fluoresce (Thabet et al., 2015). Thus some of the changes in intensity relative to duration of exposure could also be linked to changes in movement or feeding and associated differences in calcium use.

After 1, 7 and 14 days of exposure, pteropods were preserved in RNAlater for analysis of gene expression. Total RNA was extracted from samples of five to nine pteropods pooled from each treatment using the E.Z.N.A. Mollusc RNA Kit (Omega Biotek, Norcross, GA, USA). Three RNA samples were selected from each day and treatment combination (27 samples in total) based on spectral profile and quantity, and these were submitted to the University of Rochester Genomics Research Center for sequencing. Libraries were constructed using TruSeq reagents and then sequenced on three lanes of an Illumina HiSeq 2500 (nine samples per lane) as a High Output version 4 125 bp PE project. The sequencing facility used Trimmomatic (version 0.32; Bolger et al., 2014) to eliminate adapter sequences (2:30:10) and removed low-quality scores using a sliding window (4:20), trimming both trailing and leading sequences (13) and leaving only sequences with a minimum length of 25 for downstream use.

One replicate per day and treatment was randomly chosen for transcriptome assembly (nine RNA samples consisting of 59 individuals), and de novo assembly was performed with Trinity (r2014-07-17; Grabherr et al., 2011) using default parameters (k-mer=25, cut-off of 200 bp). The de novo transcriptome was annotated through blastx searches against both nr (e-value <1.0 e⁻⁵) and Swiss-Prot (e-value <1.0 e⁻⁶) to allow for direct comparison with previous pteropod transcriptome assemblies (Koh et al., 2015; Maas et al., 2015; Johnson and Hofmann, 2016; Moya et al., 2016; Thabet et al., 2017). Annotations identified as unknown, hypothetical or uncharacterized proteins were not included in calculations of annotation success but were retained for downstream differential expression analyses.

Reads from individual RNA samples were then aligned to the transcriptome using Bowtie2 (version 2.2.3; Langmead and Salzberg, 2012), estimates of abundances were made with RSEM (Li and Dewey, 2011), and edgeR analysis of differential expression (DE) was performed with R version 3.0.1 (Robinson et al., 2010), using the pipeline packaged with Trinity (Haas et al., 2013). Genes were defined as DE if the false discovery rate and the P-value were both <0.05, and the log₂-fold change was >2 (corresponding to a fourfold change in expression). To explore gross patterns of gene expression among days and treatments, gene counts in each sample were compared via a Bray Curtis Similarity Matrix using the FATHOM toolbox (http://www.marine.usf.edu/user/djones/matlab/matlab.html) in MATLAB.

Gene ontology (GO) annotation and an Interpro scan were performed in Blast2GO (Conesa et al., 2005). The GO terms from these annotations were then merged for GO enrichment analyses, applying a Fisher's exact test with false discovery rate (FDR) (Benjamini and Hochberg) using a P-value and FDR cut-off of <0.05 within Blast2GO. One analysis was done for each of the DE gene lists associated with captivity (generated by edgeR, see above) and for each associated with CO₂. Results were reported using the most specific GO term available. DE genes were additionally annotated beginning with a Blastx search against nr with an e-value <1.0e⁻¹⁰.

The effect of captivity alone was assessed by identifying changes of gene expression over time (days 1, 7 and 14) in the ambient treatment. The effect of CO₂ on gene expression was tested within each day and then the differentially expressed transcripts that were shared among analyses were queried to determine whether there were consistent patterns in transcriptomic response among treatments and across days.

Quantitative RT-PCR (qPCR) was used to measure transcript expression of four genes that were differentially expressed in the Illumina dataset and two reference genes that showed stable expression in the Illumina dataset. All primer sequences are given in Table S1. Residual genomic DNA was removed from total RNA using the RNA Clean & Concentrator Kit (Zymo Research, Irvine, CA, USA). Sufficient RNA was available for reverse transcription from a total of 25 samples; some of these samples had also been used in RNA-seq analyses and others were distinct samples collected within the same experiment. Complementary DNA (cDNA) was synthesized of 450 ng total RNA in 20 µl reactions using the Iscript cDNA Synthesis Kit (Bio-Rad, Hercules, CA, USA).

qPCR reactions were run in duplicate in 20 µl volume using 10 µl iTaq Supermix (Bio-Rad), 8 µl molecular biology grade water, 1 μl cDNA and 1 μl of 10 μM primers in an iCycler iQ real-time PCR detection system (Bio-Rad). Cycling conditions were: 95°C for 1 min followed by 40 cycles of 95°C for 15 s and 60°C for 25 s. After cycling, the products were subjected to melt curve analysis to ensure that only a single specific product was amplified. Amplification of control samples (lacking either template or reverse transcriptase) was negligible. Amplification efficiency was calculated within each sample using the LinReg program, and produced mean values ranging from 97 to 104%. Transcript expression was calculated in arbitrary fluorescence units by LinReg. Expression of genes of interest was divided by the geometric mean of the two reference genes. Pearson's coefficient was calculated to assess the correlation between mean Illumina and qPCR expression using seven points (day–treatment combinations) where both analyses were conducted.

There was a clear effect of captivity evident in all the response metrics within the ambient treatment, with a substantial increase in the number of DE genes over time, as well as decreases in respiration rates and calcification. Separate from this captivity effect, however, CO₂ exposure significantly influenced respiration, calcification and gene expression. Pteropods in the medium and high CO₂ treatments had significantly higher respiration rates after one week and calcification after two weeks when compared with the ambient treatment. Respiration showed an interactive effect, with increasing divergence in the respiration rates when comparing ambient and medium/high CO₂ treatments as the exposure duration increased. These organismal responses were reflected in the gene expression patterns, with increasing differences between the ambient and CO₂ treatments as the exposure duration exposure increased.

Carbonate chemistry measurements averaged over the course of the exposures suggest that target values were achieved throughout the experiment (Table 1). Ambient P_CO2 was somewhat higher than expected based on the global average level (i.e. ∼460 ppm rather than 400 ppm), probably owing to local respiratory processes near the compressor in-line, but within the normal annual P_CO2 that the pteropods experience in surface waters (Wang et al., 2017). Although calculations of P_CO2 based on measurements of water chemistry (DIC/TA) at the initiation of the experiment suggest insufficient pre-equilibration of the medium and high treatments, subsequent pH measurements indicate that target levels were reached within a day of bubbling (Table S2, Fig. S2). Based on an average of the DIC/TA sampling pairs from the two weekly water changes, and excluding these low initial measurements, the ambient treatment was always over-saturated with respect to aragonite (1.56), while on average the medium treatment was slightly under-saturated (0.96) and the high treatment was strongly under-saturated (0.71; Table 1).

Both duration and severity of CO₂ exposure significantly influenced respiration (Fig. 1). Over the course of the experiment, the average mass declined from 0.45 mg on day 1 to 0.28 mg on day 14, with a more substantial decrease in the high treatment (data not shown). Pteropods ranged in size from 0.05 to 0.7 mg (0.3–0.8 mm in length), and as respiration rate correlated significantly with mass, mass was included as a covariate in subsequent analyses. A two-factor (duration and treatment) generalized linear model with a Bonferroni post hoc test revealed a significant decrease in respiration with increasing exposure time (F_3,89=10.783, P<0.001). CO₂ treatment also had a statistically significant effect (F_2,89=11.355, P<0.001), with the medium and high treatments having a higher respiration rate than ambient treatment after 7 and 14 days of exposure. These effects of duration and treatment interacted significantly, with an increasing difference between ambient and the other treatments as the duration of exposure increased (F_6,89=4.805, P<0.001). Data met the assumptions of normality and homogeneity of variance.

Imaging of calcein-stained shells indicated higher calcification in pteropods from the medium and high treatments, and also showed an overall decreasing trend relative to duration of exposure (Fig. 2). Although there were regions of enhanced fluorescence intensity, calcein binding occurred over the entire shell and not just at the leading edge of the last whorl where growth would be expected to occur (Fig. 2A). Fluorescence intensity differed significantly among the CO₂ treatments (Fig. 2B) after exposure durations of 2 weeks (one-way ANOVA: P=0.001) and 3 weeks (one-way ANOVA: P<0.001), but not during week 1 (one-way ANOVA: P=0.053), although all three time points showed similar trends of increased fluorescence in the medium and high treatments. For weeks 2 and 3, pairwise analysis revealed the ambient treatment to have shells with significantly lower fluorescence than the medium and high treatments (Holm–Šidák method: P<0.05), but that the medium and high treatments were not different from each other. Both treatment and exposure duration affected fluorescence intensity (two-way ANOVA: P<0.001), and there was no interaction between factors (P=0.753). Calcein fluorescence after 3 weeks of exposure was significantly lower than after 1 or 2 weeks (Holm–Šidák method: P<0.05), while weeks 1 and 2 were not significantly different from each other. Data met the assumptions of normality and homogeneity of variance.

The assembly metrics were comparable to recent pteropod de novo transcriptomes and sufficient to allow for DE analysis (Table 2). Trimming produced ∼624.5 million paired-end reads with 23±1.8 million reads per library (mean±s.d.). A transcriptome was assembled from 199.3 million reads (9-library subset), resulting in 116,910 transcripts constituting 95,145 ‘genes’ with an N50 of 737 bp. Of the transcripts, 34,941 (29.9%) were annotated by a Blastx search against Swissprot (e-value <1.0 e⁻⁶) and 48,870 (41.8%) against nr (e-value <1.0 e⁻⁵). On average, 71.73% of all sequences used in the DE analysis mapped back onto this assembly, and of these 100% mapped uniquely to one gene.

Comparing the transcriptomic profiles across all CO₂ treatments, there was greater similarity among samples from the same day (15–20%) than among samples from the same CO₂ treatment across days (23–26%), with decreasing similarity within a day as the duration of exposure progressed, irrespective of CO₂ treatment (Fig. 3). The least similarity was between the day 1 and day 14 samples (35%). Analysis of DE in the ambient treatment based on the day of sampling (hereafter captivity analysis) shows an increasing number of DE genes over time (Table 3). Overall, captivity influenced between ∼2 and 10% of the genes in the transcriptome.

GO enrichment analyses indicated that within the first 7 days of captivity there were significant changes in physiology, with increased DE in transcripts associated with translation, energetic metabolism in the mitochondria, the ribonucleoprotein complex and stress responses (immune/inflammatory, and apoptosis signalling; Table S3A). After two weeks (i.e. day 1 versus day 14), GO enrichment analysis showed similar but more pronounced responses. There was also an increase in the number of GO terms associated with neuronal signalling, oxidative stress and protein metabolism/degradation. During this period there was also a decrease in the number of GO terms associated with DNA replication, and signalling associated with mitochondrial energetic metabolism (related to GTPases and the cAMP pathway). There were no significant differences between the day 7 and day 14 captivity DE based on GO enrichment analyses. To further understand these results, we examined the BLAST annotations of transcripts within the enriched GO categories (Table S3B). The stress response transcripts were particularly represented by upregulated heat shock proteins. An increasing number of mitochondrial genes (ATP synthase, cytochrome oxidases, NADH dehydrogenases/ubiquinone and GTP-binding proteins) were upregulated as the duration of captivity increased (13 in day 1 versus day 7, 23 in day 7 versus day 14, and 62 in day 1 versus day 14).

Despite being an analysis of only the ambient treatments over time, multiple classes of transcripts that have previously been associated with CO₂ responses were upregulated after 7 and 14 days of captivity, including genes associated with ion transport, neural signalling and biomineralization (chitin binding/metabolism, cytoskeletal organization). In contrast, carbonic anhydrase (which has many physiological functions, but is also associated with biomineralization) showed a U-type response curve with multiple forms showing greatest expression on day 7, followed by day 14.

Compared with the effect of captivity, fewer genes were DE in response to CO₂ exposure (<1.5%). This was, however, similar to transcriptomic responses to CO₂ that have been seen previously in ecologically similar pteropod species (1.2%, Koh et al., 2015; 2.6%, Moya et al., 2015), and far greater than has been observed in Clio pyramidata, a vertically migratory pteropod species that probably experiences a wider range of CO₂ conditions naturally (0.001%, Maas et al., 2015) or other pelagic organisms such as copepods (0.003%, Bailey et al., 2017). Although this is low in comparison with OA response in benthic, longer-lived organisms such as coral (∼19%, Moya et al., 2012), and oysters (∼10%, Li et al., 2016), these studies on larger organisms were conducted with specific tissues [e.g. Li et al. (2016) used pearl oyster mantle], or a discrete population [e.g. Moya et al. (2012) used larvae generated from crosses of 15 coral colonies] which decreases variability and increases statistical power in gene expression analyses. Furthermore we feel that the population size and life history of pelagic organisms contribute to greater variability in both genotype and phenotype within captured planktonic populations, and subsequently result in fewer significantly DE genes.

There were no significantly enriched GO terms among any of the CO₂ comparisons. Increasing numbers of DE transcripts were detected in response to CO₂ as both the duration and level of exposure increased (Table 4). The greatest number of DE transcripts was consistently associated with the comparisons between the ambient and high treatments. While no transcripts were consistently DE across all days, there were some consistent patterns within gene families, particularly when comparing the ambient and high treatments. At all three time points there were genes associated with collagen, fibrocystin, zinc metalloproteinase, the extracellular region, calcium ion binding, and serine-type endopeptidase activity were upregulated in the higher CO₂ treatments. Beyond those transcripts that have a link to respiration and biomineralization, many DE transcripts were associated with inflammation pathways (with an overall downregulation) including modification of the eicosanoid signalling, prostaglandin synthesis, phospholipase, annexin and plasminogen transcripts.

On day 1, CO₂ treatment primarily resulted in upregulation of transcripts within the high relative to the ambient treatment (129 of 133 DE transcripts). GO annotation identified transcripts associated with the extracellular matrix (plasminogens, collagens, serine protease), lipid metabolic processes (allene oxide synthase-lipoxygenase, lipases), carbohydrate, metal ion and/or calcium ion binding properties (zinc metalloproteinase-like; Table S4). Many of these genes were also DE in the comparison between the high and medium treatments (32 of 133 transcripts).

On day 7, an increasing proportion of transcripts were downregulated in response to enhanced CO₂ (76; 37%), and the overall number of DE transcripts increased (205). Those transcripts that were downregulated in response to CO₂ were often shared (31) in the contrasts of ambient to both the medium and high treatments and included transcripts with potential biomineralization function (alkaline phosphatase, carbonic anhydrase, ceramide kinase, a sodium bicarbonate cotransporter and exchanger, and a transcript associated with chitin binding/metabolic process). The medium and high treatments shared a number of upregulated transcripts in contrast to ambient (31), including fibrocystins, collagens, cytochrome oxidase subunits and a zinc metalloproteinase.

On day 14 the number of DE transcripts in the medium and high treatments in comparison with ambient had substantially increased (566 and 1304, respectively), while the number of transcripts differing between the medium and high treatments remained similar (20). The number of genes downregulated in response to CO₂ remained similar in proportion of the DE from day 7 (731 of 1504; 24% in the ambient versus medium, and 47% in the ambient versus high treatments). Once again the downregulated and upregulated transcripts were substantially shared in the contrasts of ambient to both the medium and high treatments (25 and 344, respectively). A few of the downregulated biomineralization-associated transcripts that were DE in the day 7 contrast (alkaline phosphatase, carbonic anhydrase) remained downregulated, while a number of additional transcripts became downregulated, including those with GO terms associated with microtubule processes, oxidoreductase activity/processes, metabolic process, and metal or calcium ion binding. In comparison with the ambient treatment, fibrocystins, cilia/flagella associated and biomineralization-associated (i.e. calmodulin, calponin, mucin) transcripts were consistently upregulated in the medium and/or high treatments. There were also a number of sequences that have previously been implicated as being important in biomineralization (ferric-chelate reductase, chitin-binding, chitinase) that were upregulated on day 14 in the high CO₂ treatment. Finally, a set of lipid metabolic processes became substantially influenced by day 14 in the high CO₂ treatment, with upregulation of lipases, downregulation of phospholipases associated with inflammatory pathways, and upregulation of lipoxygenases and other signalling transcripts possibly associated with prostaglandin synthesis.

It is important to note that some genes that have previously been implicated in CO₂ stress, such as those associated with the energetic metabolism in the mitochondria (i.e. cytochrome oxidases, and NADH dehydrogenase), belonged to multi-gene families from which various transcripts showed contrasting responses to CO₂ exposure. Genes in these families were also substantially upregulated as the duration of captivity increased, thus there may be interactive effects of duration and exposure level. Overall, many genes associated with mitochondrial metabolism were generally upregulated with CO₂ by day 7 and strongly downregulated on day 14.

The qPCR measurements broadly corroborated results from Illumina-based sequencing (Fig. S1). For example, both methods showed downregulation of carbonic anhydrase in response to CO₂ exposure (Fig. 4A,B). The qPCR measurements also supplemented the Illumina analyses, showing that downregulation of carbonic anhydrase begins within 3 days of exposure. For exposure durations and CO₂ levels where both measurements were made, the two measurements were strongly correlated (Fig. 4C).

Finally, when comparing the DE genes from the CO₂ analysis with those from the captivity analysis, there were a number of transcripts that co-occurred (8–17% of the DE genes in the captivity analysis were also DE in at least one comparison between ambient and a CO₂-exposed group). GO enrichment analyses based on the co-occurring transcripts identify these as being associated with cytochrome c oxidase activity, oxygen binding, microtubules and the structural constituents of the cytoskeleton (Table S5).

Together our results suggest that L. retroversa collected from the Gulf of Maine during spring responded to CO₂ exposure with complex physiological changes that included an increased respiration rate, increased calcification activity and patterns of gene expression that implicate changes in biomineralization, metabolism, neural function and ion binding (summarized in Fig. 5). Captivity, however, produced a similar or greater effect on physiology, resulting in decreased respiration, reduced calcification and a substantial change in gene expression, requiring careful interpretation of the combined findings. Despite some overlap in gene expression with the captivity response, a number of the patterns of DE appear to be solely attributable to the exposure to CO₂ as they were consistent among days and treatment levels.

The number of transcriptome elements (genes or isoforms) generated in molluscan studies varies substantially (∼20,000–700,000), and understanding of function remains limited by low annotation success (14–36%; reviewed in Harney et al., 2016). Our study falls mid-range for pteropod studies (Table 2), with 30% of assembled sequences identified using Swissprot and e-value <1.0 e⁻⁶, and 42% using the nr database and an e-value <1.0 e⁻⁵. The relatively high N50 (∼700), the ∼70% mapping success of all the treatments, and reasonable annotation success suggest that this is an adequate assembly for the purposes of DE analysis.

Captivity resulted in the induction of stress response genes, including heat shock proteins, and transcripts associated with apoptosis, immune and inflammatory responses, and energetic metabolism in the mitochondria. The upregulation of these genes, the lack of response of lipid metabolism genes, as well as the downregulation of those in the cAMP and G-protein signalling pathways in the captivity response is in direct contrast to patterns observed in studies of other starving invertebrates (Zinke et al., 2002; Suo et al., 2006) and fish (Salem et al., 2007). These patterns lend support to the conclusion that nutrient deprivation did not dominate the observed changes in physiology. In the ambient treatment, the respiration rate declined over time, with the latest time point (14 days) significantly reduced compared with day 1. This is in apparent direct contrast with the upregulation of transcripts associated with energetic metabolism (oxidative phosphorylation) in the mitochondria. The set of genes that was influenced, including NADH dehydrogenase, ubiquinol and cytochrome c subunits, however, also play important roles in the oxidative stress response and apoptosis (Kadenbach, 2003; Jiang and Wang, 2004). As these functions were also substantially upregulated in response to captivity, it is likely that these transcripts represent a malfunction of the mitochondria in relation to stress, which has been shown to result in a decrease in organismal respiration rate (Kadenbach, 2003; Murphy, 2009; Salin et al., 2015). Upregulation of transcripts associated with the energetic metabolism of the mitochondria may then be compensation for the inefficiency of the malfunctioning mitochondria and replacement of cytochrome c lost to apoptosis signalling.

Captivity also affected calcification, with a decrease in intensity of calcein fluorescence across all treatments over time, which may be connected to the observed reduction in respiration rate and resource allocation away from active biomineralization. In contrast to the calcein findings, however, transcripts potentially associated with biomineralization, such as ion transporters, chitin binding proteins and cytoskeletal organizational elements became upregulated over the duration of captivity, suggesting more biomineralization effort. These transcripts are not directly related to calcium carbonate deposition, so they may represent changes in the organic matrix, or other functions that are not visualized by the fluorescent stain.

Several DE genes were shared between the captivity and CO₂ exposure analyses. These include transcripts related to the ribosome (site of protein synthesis), mitochondrial energetic metabolism, stress response, biomineralization and the cytoskeleton. These transcripts may represent a general stress response, be the spurious result of random chance due to multiple statistical analyses, or occur in the CO₂ analyses only due to an interactive effect of captivity. It is likely that some shared transcripts are attributable to each of these factors. If the protocol of capturing and culturing increases the sensitivity of the organisms to CO₂ exposure, this could also increase the number of shared transcripts. Shell repair (due to mechanical damage) is one possible mechanism for interaction between the stressors. It has been shown that thecosome pteropods can repair their shell after physical damage (Lischka et al., 2011; Thabet et al., 2015; Peck et al., 2016), and some of the gene expression patterns associated with this sort of damage are probably similar to those in response to dissolution.

These results emphasize that the culturing protocol was not optimal. Although this is not particularly surprising, as keeping the negatively buoyant planktonic pteropods in captivity has consistently been shown to be a difficult problem (reviewed in Howes et al., 2014), it must be kept in mind when interpreting the results of the experiments. Concurrent with the experiments reported here, a number of culturing tests were conducted that successfully maintained the same species from eggs for up to 6 months, during which time they grew into large adults and spawned, in identical containers and food conditions, but in seawater drawn from shallow waters local to the laboratory and without the presence of bubbling (reported in Thabet et al., 2015). As such, it is likely the individuals had sufficient and appropriate food under these laboratory conditions. Thus it seems likely that the ‘captivity’ patterns reported in this study are driven by effects of the bubbling itself, which was used to maintain the acidified conditions, or the use of stored filtered water, which provided consistent oceanic alkalinity levels. If bubbling reduced feeding capabilities, then nutritional state, which has been implicated in pteropod OA responses (Seibel et al., 2012), may have contributed to the observed physiological changes. The gene expression patterns associated with captivity in this study are not consistent with starvation, but atypical feeding conditions may still have played a role.

This captivity analysis is important as it suggests that during longer term experiments with pteropods (>7 days), some of the observed reductions in metabolic rate may be the consequence of culturing conditions (Lischka and Riebesell, 2016). It also suggests that longer duration exposure to CO₂ may be necessary to detect changes in respiration, possibly resolving some of the variability observed in previous experiments which ranged from ∼12 h to 9 days of exposure (Comeau et al., 2010; Maas et al., 2012; Seibel et al., 2012; Lischka and Riebesell, 2016). All gene expression studies to date have been conducted on individuals exposed for less than a week and, for these, patterns of gene expression associated with cytochrome c oxidase activity, oxygen binding, microtubules and the structural constituents of the cytoskeleton may be influenced by captivity. Importantly, there appears to be no interactive effect of captivity on shell quality analyses, emphasizing the utility of biomineralization of pteropods as a useful bio-indicator of ocean acidification response (Bednaršek et al., 2017).

Despite the apparent declining health of the organisms, there was a clear influence of CO₂ on respiration, calcification and gene expression. From day 3 onwards the medium, and to a greater extent the high, CO₂ treatments consistently showed a higher mass-specific respiration rate in comparison with ambient treatment. This finding is in contrast with previous work done with this species in the Arctic (Lischka and Riebesell, 2016), which found no significant effect of enhanced CO₂ on the respiration rate for L. retroversa after 7 or 9 days of exposure, but a significant suppression in overwintering juveniles after 18 h of exposure at under-saturated conditions (Ω_Ar<0.6).

These two populations are quite separate, and based on differences in water chemistry, those in the Arctic are much more likely to experience conditions of high CO₂, low pH and aragonite under-saturation with regard to aragonite in their natural vertical distribution. Furthermore, the two studies were done during different phases in the life cycle, with the Arctic pteropods exposed during their overwintering period when they would be experiencing low food availability, in situ saturation states as low as ∼0.75, and possibly a period of ceased growth (Lischka and Riebesell, 2012). In contrast, April in the Gulf of Maine coincides with the spring bloom when surface saturation state exceeds 1.8 in the upper 50 m where pteropods are found (Wang et al., 2017), and the difference between environmental and laboratory conditions may have resulted in the more pronounced physiological response in our study region. Our findings suggest that exposure to enhanced CO₂ resulted in physiological compensation in this population, at least during this season.

In our analysis of L. retroversa gene expression under enhanced CO₂, changes in the respiration rate again appear to be in direct contrast to changes in mitochondria-associated genes, including DE of cytochrome oxidase and NADH dehydrogenase transcripts. There was a strong and consistent pattern of upregulation of these genes in association with the captivity analysis, and an associated clear pattern of reduction in respiration rate. In the CO₂ analysis the same genes show a more erratic pattern, often with multiple transcripts with similar annotation occurring concurrently with contrasting expression patterns. It was, however, more common for these transcripts to be upregulated on day 7 and downregulated by day 14 in the higher CO₂ treatments. For respiration there was a significant interaction between treatment level and duration, with a statistically significant increase in respiration only by day 7. Thus the change from upregulation to downregulation in these mitochondrial genes could represent changes in energetic metabolism of the thecosomes between days 7 and 14 that are a reflection of this interaction between duration and level of CO₂ stress.

A possible explanation for this pattern is that the pH changes in the environment are prompting a change in mitochondrial membrane proton permeability, which is associated with decreased efficiency of oxidative phosphorylation and increased oxidative stress (Capitanio et al., 1996; Kadenbach, 2003; Murphy, 2009). The observed upregulation by day 7 may be an attempt to restore metabolic rate in these unfavorable conditions and to cope with the oxidative stress, while by day 14 the organisms may be so compromised in the high CO₂ condition that they are shutting down production of mitochondrial components to reduce oxidative stress. A number of previous studies have reported oxidative stress associated with CO₂ exposure (Tomanek et al., 2011; Hu et al., 2015; Klanian and Preciat, 2017), and one measured a concomitant increase in respiration (Klanian and Preciat, 2017). These studies, however, showed increases in anti-oxidant protein level or enzyme activity (superoxide dismutase, catalase, glutathione peroxidase, etc.) rather than implicating transcripts associated with mitochondrial oxidant production as seen in this study. Genes associated with the mitochondrial inner membrane complex have also shown patterns of downregulation (as in our day 14 comparison) in association with exposure to CO₂ in C. pyramidata (Maas et al., 2015), Heliconoides inflatus (Moya et al., 2016), urchins (Todgham and Hofmann, 2009; O'Donnell et al., 2010) and in coral (Moya et al., 2012), and are often interpreted as being associated with a reduction (rather than increase) in oxidative metabolism. Of the previous studies that reported a change in mitochondrial inner membrane complex genes, those that directly measured respiration rate did not find a substantial change (Maas et al., 2015; Moya et al., 2016). The significant influence of CO₂ on the respiration rate of L. retroversa in this study may be a result of the larger sample size, or phylogenetic differences. A study of the adaptive influence of OA on populations of copepods showed selection on transcripts associated with mitochondrial genes (including cytochrome oxidase) and oxidative phosphorylation after OA exposure (De Wit et al., 2015), and similar selection processes may contribute to the observed patterns in pteropods. Overall the linkage between this mitochondrial energetic transcript expression pattern, respiration rate and oxidative stress requires more directed study.

Enhanced binding of calcein relative to the ambient control was consistently observed in the medium and high treatments throughout the experiment. Calcein binding was seen over the entirety of the shells in all treatments, although it was not uniformly distributed (Fig. 2). Calcein incorporation over much of the shell was seen previously in juvenile Limacina helicina (Lischka et al., 2011) exposed to ∼1100 µatm (Ω_Ar 0.71) for 29 days. That study examined cross-sections of the shell near the aperture and found calcein fluorescing along the internal edge of the shell. These results are in contrast to a study by Comeau et al. (2009), which found that there was a distinct growth band of calcein dye incorporated only along the growing edge of L. helicina exposed to either ambient (350 µatm; Ω_Ar 1.90) or high (765 µatm; Ω_Ar 1.00) CO₂ conditions for 5 days. Calcein staining of normally developing early-stage L. retroversa under ambient CO₂ showed calcification occurring over much of the shell, with greater dye incorporation along the leading edge in late-stage veligers and juveniles (Thabet et al., 2015). As calcification was occurring throughout the shell of the present live animals, and was more pronounced in the higher CO₂ treatments, the fluorescent regions would seem to be sites of shell repair.

Unlike respiration, calcification patterns were clearly and directly linked with gene expression in response to CO₂ exposure. Upregulation of transcripts associated with biomineralization, including calcium-binding proteins, metalloproteinases, mucins, chitin and collagen in our enhanced CO₂ exposures, was consistent with the CO₂-associated increase in calcein staining, beginning on day 1 and continuing throughout the experiment. These DE results were similar to the findings in H. inflatus, which showed upregulation of a similar set of calcification-associated genes in response to high CO₂ (Moya et al., 2016), and to observations in C. pyramidata, which also responded to high CO₂ with increases in transcripts associated with collagen, calcium ion binding and the extracellular region (Maas et al., 2015). In contrast, Koh et al. (2015) found downregulation of mucins, chitinase and collagen under the highest CO₂ treatments, although some chitin and collagen-associated transcripts were both downregulated and upregulated in the two CO₂ treatments. As aragonite saturation state was not reported by Koh et al. (2015), it is difficult to assess if these differences among treatments and studies are due to dissimilarities in carbonate chemistry or species sensitivity.

In contrast with many other biomineralization-associated transcripts, carbonic anhydrase and alkaline phosphatase, which have both been posited to be important in the biomineralization of molluscs (Gaume et al., 2011; Xiang et al., 2014), were downregulated in the longer duration enhanced CO₂ treatments in our study. Alkaline phosphatase is implicated in the initiation of biomineralization and/or the secretion of structural organic matrix, while some forms of carbonic anhydrase accelerate bicarbonate formation. These processes have been shown to be influenced by OA in pearl oysters (Li et al., 2016), where lower pH (7.5; Ω_Ar ∼1.0, similar to our medium treatment) resulted in decreased alkaline phosphatase activity and downregulation of carbonic anhydrase precursors. In our study, the similar downregulation of both of these transcripts in the longer high-CO₂ exposures may reflect changes in the structural organic matrix and the initiation of calcification, which are not visualized by the organism-level calcein staining. This is supported by a companion study that found changes in transparency of the L. retroversa shell in the presence of CO₂ (Bergan et al., 2017). The downregulation of these two transcripts, in addition to our calcein staining results and upregulation of extracellular matrix-associated genes, are in direct contrast to the response documented by Moya et al. (2016) in H. inflatus who found increased carbonic anhydrase and alkaline phosphatase gene expression, decreased biomineralization (⁴⁵Ca uptake), and under-representation of GO terms associated with the extracellular matrix with their higher CO₂ treatment. It is important to note that the carbonate chemistry treatments of Moya's experiment were substantially different from the conditions experienced by L. retroversa in our experiments or by the pearl oyster in the experiment by Li et al. (2016). The high treatment in Moya's experiment had a pH of 7.9 and Ω_Ar ∼2, which is more similar to our ambient condition than to our medium or high treatments, which both reached under-saturation. As a consequence it seems likely that the reversal in biomineralization pattern between the studies is a direct consequence of differences in saturation state. More generally, differences in the basic water chemistry between studies (total alkalinity, salinity, etc.) still remain a hurdle in comparisons of DE responses as it remains unclear which carbonate chemistry parameters are most important for biomineralization stress (reviewed in Cyronak et al., 2015a,b; Waldbusser et al., 2015a)

Our results support the notion that pteropods can undergo rapid gross changes in biomineralization-related processes (Comeau et al., 2009; Lischka et al., 2011; Bednaršek and Ohman, 2015; Bergan et al., 2017), probably of both the aragonitic and proteinaceous components of the shell, as a consequence of CO₂ exposure. Both borderline and strongly under-saturated conditions had a similar influence, particularly as duration of exposure increased. The increased calcification, indicated by calcein staining, in pteropods within the medium and high CO₂ treatments provides organism-level validation of the observed transcriptional changes and is particularly consistent with the increased calcium-ion binding activity. A more detailed understanding of the particular mechanisms would be facilitated by RNA-seq of mantle-specific expression, as has been done for other molluscs (Artigaud et al., 2014; Hüning et al., 2016; Li et al., 2016), but the large number of individuals this would require and the technical difficulty of the associated dissection prevented any attempt during this study.

Several other transcripts that showed DE patterns in relation to CO₂ exposure were consistent with previous studies of thecosome pteropods (Koh et al., 2015; Maas et al., 2015; Moya et al., 2016), including an upregulation of transcripts associated with serine peptidase activity, a transcript with putative neuropeptide signalling and calcium ion binding function (annotated as a polycystic kidney disease protein), fibrocystins, zinc metalloproteinases and coagulation factors. These gene families are thus good candidates for greater mechanistic exploration as they appear to be broadly influenced within pteropods, although their specific physiological function in the group remains unclear.

Overall, our results suggest that both borderline and strong aragonite under-saturation substantially affect gene expression in L. retroversa and result in increased respiration rate and calcification effort. Of particular note, we saw similar changes in calcification under both CO₂ treatments over time, but complex changes in respiration and the DE transcripts as the severity and duration of the treatment increased. There were transcripts with DE in the mid-level CO₂ treatment that were not DE in the high-level CO₂ treatment, and, in some gene sets (particularly mitochondrial associated transcripts), increasing DE followed by decreasing DE as the duration of exposure increased, suggesting that some responses are either non-linear or that there are isoforms that respond specifically to the different carbonate chemistry conditions. As this is the first transcriptomic study we know of that explicitly tests for the effect of duration of CO₂ exposure, these findings, and their corresponding organismal-level metrics, may be useful to compare among earlier studies of differing exposure length.

When viewed as a whole, our results emphasize the differences between short-term laboratory experiments and the longer-term reality of acidification exposure. In marine environments organisms may periodically experience under-saturation, such as during diel vertical migrations, upwelling events or via riverine discharges. Under these sorts of conditions (brief but acute), adult and juvenile pteropods typically show changes in shell quality (i.e. Comeau et al., 2010; Bednaršek et al., 2012c; Lischka and Riebesell, 2012), but no change in respiration rate, and a limited transcriptional response (Maas et al., 2012, 2015, 2016b). From our results, biomineralization remains a clear early indicator of OA stress in the group and there are a number of gene patterns that are immediately influenced, re-emphasizing the idea that thecosome shell quality and associated transcripts could be a potential bio-indicator of environmental acidification (Bednaršek et al., 2017). These short-term responses are important to understand, but it is clear from this study that longer-term exposure to high CO₂, even if it is close to saturation, has a larger and more significant influence on pteropods than shorter duration exposure. Thus the insights afforded by the present analysis may suggest a better assessment of how this species will respond to real-world OA.

We would like to thank R. Galat, S. Gallager, D. McCorkle, M. White and C. Zakroff for assisting in the set-up of the culturing and CO₂ exposure facilities. We greatly appreciate the insight of D. McCorkle and the collaboration of Z. A. Wang and K. Hoering on the carbonate chemistry measurements. We much appreciate the hard work and dedication of Captain K. Houtler and Mate I. Hanley and would like to thank them for their excellent support aboard the research vessel Tioga. At-sea sampling was supported by S. Chu, T. Crockford, M. Edenius, A. Thabet and P. Wiebe. Special thanks is owed to P. Alatalo for critical assistance in maintaining cultures of pteropods and phytoplankton, providing insight and advice, and for consistent hard work during experiments.