Impacts of Multiple Stressors on a Benthic Foraminiferal Community: A Long-Term Experiment Assessing Response to Ocean Acidification, Hypoxia and Warming

Ocean chemistry is changing as a result of human activities. Atmospheric carbon dioxide (CO2) concentrations are increasing, causing an increase in oceanic pCO2 that drives a decrease in oceanic pH, a process called ocean acidification (OA). Higher CO2 concentrations are also linked to rising global temperatures that can result in more stratified surface waters, reducing the exchange between surface and deep waters; this stronger stratification, along with nutrient pollution, contributes to an expansion of oxygen-depleted zones (so called hypoxia or deoxygenation). Determining the response of marine organisms to environmental changes is important for assessments of future ecosystem functioning. While many studies have assessed the impact of individual or paired stressors, fewer studies have assessed the combined impact of pCO2, O2, and temperature. A long-term experiment (∼10 months) with different treatments of these three stressors was conducted to determine their sole or combined impact on the abundance and survival of a benthic foraminiferal community collected from a continental-shelf site. Foraminifera are well suited to such study because of their small size, relatively rapid growth, varied mineralogies and physiologies. Inoculation materials were collected from a ∼77-m deep site south of Woods Hole, MA. Very fine sediments (<53 μm) were used as inoculum, to allow the entire community to respond. Thirty-eight morphologically identified taxa grew during the experiment. Multivariate statistical analysis indicates that hypoxia was the major driving factor distinguishing the yields, while warming was secondary. Species responses were not consistent, with different species being most abundant in different treatments. Some taxa grew in all of the triple-stressor samples. Results from the experiment suggest that foraminiferal species’ responses will vary considerably, with some being negatively impacted by predicted environmental changes, while other taxa will tolerate, and perhaps even benefit, from deoxygenation, warming and OA.


INTRODUCTION
The oceans are undergoing significant environmental change. The average sea surface temperature (SST) has increased over the last century (e.g., Hansen et al., 2006), rising atmospheric partial pressure of carbon dioxide (pCO 2 ) is lowering the pH of the oceans (Doney et al., 2009), and the extent and intensity of lowoxygen (i.e., hypoxic) bottom waters is growing, at least in certain regions (e.g., Rabalais et al., 2007;Stramma et al., 2010;Capone and Hutchins, 2013;Altieri and Gedan, 2015). The biological and ecological impacts of these ongoing changes-warming, ocean acidification (OA), and hypoxia-on marine organisms are not yet fully understood.
Great progress has been made understanding the responses of organisms to each of these stressors over the past two decades (e.g., see reviews by Levin et al., 2009;Mathis et al., 2015;Gobler and Baumann, 2016;Smale, 2020). Because many environmental factors typically co-vary in nature, single-stressor studies may not be the best predictor of species' future responses. Also, different stressors can have different impacts on any given organism: they may be additive, synergistic (greater than additive), or antagonistic (offsetting) (reviewed by Breitburg et al., 2015). For instance, one stressor can be beneficial to reproduction of an organism while another can be detrimental to reproduction. In that situation, such antagonistic stressors may have no net impact on reproduction of that organism. To enable better predictions of future organism and ecosystem responses to ongoing global change, analysis of multiple stressors is critical.
Foraminiferan protists are good candidates to study differential responses to OA and hypoxia because the class includes calcareous and non-calcareous species (Sen Gupta, 1999) as well as aerobes and anaerobes (e.g., Risgaard-Petersen et al., 2006). Many marine organisms that secrete carbonate structures have been shown to be negatively impacted by OA (e.g., Courtney et al., 2013;Kroeker et al., 2013;Diner et al., 2015) and the majority of calcareous foraminifera subject to OA manipulations showed a decrease in calcification (e.g., shell normalized weight, growth rate, wall thickness) with increased pCO 2 (reviewed by Keul et al., 2013). From the perspective of oxygen stress, few eukaryotic taxa have both aerobic and anaerobic representatives. While generally thought to be aerobes, foraminifera are known to survive anoxia for extended periods (e.g., Bernhard and Reimers, 1991;Bernhard and Alve, 1996;Moodley et al., 1998;Langlet et al., 2013;Nardelli et al., 2014), with many species able to perform denitrification (Risgaard-Petersen et al., 2006;Koho et al., 2011;Bernhard et al., 2012;Woehle et al., 2018;Orsi et al., 2020).
Many studies have documented foraminiferal response to a single environmental stressor or two stressors (reviewed in Keul et al., 2013;Doo et al., 2014). These studies have addressed single species responses in lab experiments (Nigam et al., 2008;Dissard et al., 2010;Guamán-Guevara et al., 2019), single species in field assessments (e.g., Uthicke and Fabricius, 2012) or multiple species in lab cultures (OA, Haynert et al., 2014), or in field settings (temerature and salinity; Titelboim et al., 2016) (pH, salinity, oxygen;Charrieau et al., 2018a). Because of their conspicuous presence on tropical reefs, perhaps most notable climate change stressor impact studies are those on Larger Benthic Foraminifera (LBF) (Sinutok et al., 2011;Uthicke and Fabricius, 2012;Schmidt et al., 2014;Botté et al., 2020;Oron et al., 2020), which are mm-scale species that typically host algal symbionts. These foraminifera are major contributors to calcification in coral reef and tropical settings (e.g., Langer, 2008). Additional studies of other representatives of the vastly diverse benthic foraminiferal community are warranted to assess climate change impacts to deeper water and/or non-tropical habitats.
Attributing the impacts of climate change on each environmental parameter from field studies are difficult given that many parameters co-vary in the environment. Laboratorybased multiple stressor studies are a way to vary environmental parameters individually to better differentiate impacts. To date, few multiple-stressor laboratory studies on foraminifera have been published. A dual-stressor experiment showed that some calcareous foraminifera are able to survive when carbonate was undersaturated (i.e., < 1) only if salinities approached marine values as opposed to brackish waters (Charrieau et al., 2018b). In this contribution, we focus on so-called smaller benthic foraminifera rather than LBFs or planktonic foraminifera, both of which often have photosynthetic symbionts that may impact their response to environmental changes.
Because there will be ecological "winners and losers" in response to global change (Fabricius et al., 2011), analysis of an entire community or ecosystem might be expected to provide additional insights into future community/ecosystem structuring. We employed the "propagule approach" (Goldstein and Alve, 2011) that exposed an entire "juvenile" foraminiferal community to experiment manipulations. Because these propagules are small (e.g., < 53 µm) and likely transported by currents (Alve and Goldstein, 2010), the propagule bank is composed of both endemic as well as exotic species. If the environment changes, only the propagules that are able to grow in those particular conditions will ultimately compose the living "adult" assemblage. In sum, propagules provide a rapid-response mechanism to environmental perturbations, thereby conferring resilience to this important microeukaryote group.
Applying the propagule approach imposes minimal stress prior to the experiment: most prior foraminiferal environmentalresponse studies used larger ("adult") specimens, individually separated from sediments prior to experiment initiation. Depending on the richness of the sample, such isolations can take days to weeks, potentially introducing laboratory artifacts. Goldstein and Alve (2011) performed a multistressor study assessing impacts of salinity, temperature and disturbance on foraminiferal propagules; their results showed differential response of their intertidal propagule inoculum.
Here, we describe results from a long-term (∼10 month) experiment on New England-shelf benthic foraminifera where oxygen, pCO 2 , and temperature were variables, to mimic a range from pre-Industrial pCO 2 to present day pCO 2 to predicted endof-the-century pCO 2 , aerated to hypoxic conditions, and 2 • C warming. While the life span of bathyal benthic foraminifera is unknown, rapid growth over 3 months in field populations (Corliss and Silva, 1993) implies life spans on the order of 6-12 months. Thus, bathyal benthic foraminifera are good candidate organisms for lab-based incubation studies similar to ours. We show that, while two taxa were ubiquitous across all treatments and samples, resulting foraminiferal assemblages differed among treatments, indicating some species will likely fare better than others as climate change continues.

Field Sampling
All sediment and seawater samples were collected from the R/V Endeavor in May 2013 at one ∼77 m deep site, informally called the "Mud Patch, " on the New England continental shelf [centered approximately at 41 • 30 N, 70 • 30 W (Bothner et al., 1981)]. Bottom waters, which were collected from within 5 m of the seafloor using a CTD-O 2 -equipped Niskin rosette, were analyzed for total alkalinity and Dissolved Inorganic Carbon (DIC) as described and reported by Wit et al. (2016). Bottom-water [O 2 ] was 5.9 mL/L (n = 3; Wit et al., 2016). Additional bottom water was transferred from Niskin bottles into 20-L Cubitainers R and chilled in a temperature-controlled environmental van (7 • C). Water-column conductivity, temperature, and oxygen profiles were obtained on the same casts performed for bottomwater collections.
Sediments were obtained with a Soutar boxcorer. The top ∼2 cm of each box core was siphoned into a stack of serially finer sieves (500, 125, 90, 53 µm) and sieved using chilled bottom water dispensed via a pressurized container (typical gardener's herbicide dispenser). Each sediment-size fraction along with bottom water was kept separately in sealable plastic tubs, which were stored in the environmental van. The same size fractions from 2 to 3 boxcores were pooled into one container. At the end of the cruise, these samples were transported chilled to WHOI where they served as the source of foraminifera for the experiment as described below. Extra bottom water was also transported to the lab in Cubitainers R .

Experiment Overview
The experiment had two stages: incubations were conducted at 7 • C for ∼7.5 months and then at 9 • C for ∼2.5 months. Five different combinations of [O 2 ] and pCO 2 were established in closed, humid chambers, with elevated and lowered concentrations maintained by additions of compressed gases controlled by O 2 and CO 2 sensors. These five treatments included: present day pCO 2 and normoxic (aerated) ("Present Day"; Treatment 1); present-day pCO 2 with hypoxia ("Hypoxia"; Treatment 2); high pCO 2 normoxic ("OA"; Treatment 3); high pCO 2 with hypoxia ("OA + Hypoxia"; Treatment 4); and Pre-industrial pCO 2 normoxic ("Pre-industrial"; Treatment 5) ( Table 1). Very fine (<53 µm) foraminifera-laden sediments were placed in plastic containers (four containers per treatment) open to the treatment atmosphere in each chamber. The experiment began on 5 June 2013 and ended on 4 April 2014. On 14 January, half of each container was processed and designated the 7 • C sample; the remaining half continued through a 2 • C increase. The temperature increase was 0.5 • C for each of four days to allow acclimation. These remaining samples were processed on 4 April and designated as 9 • C samples. Thus, the 9 • C assemblages also experienced 7 • C for ∼7.5 months. Throughout the experiment, weekly feeding, salinity monitoring and adjustment, and bi-weekly sampling for carbonate chemistry were performed, as described in detail below.

Environmental Controllers
The O 2 and pCO 2 conditions for each treatment were set in BioSpherix C-chambers (Parish, New York, United States), which allow installation of CO 2 and/or O 2 sensors (Figure 1; see below). Each C-chamber had an open tray of distilled water installed to maintain high humidity, minimizing evaporation. Additionally, extra plastic containers, which were filled with Mud Patch bottom water, were open to the C-chamber atmosphere. This seawater was allowed to equilibrate for at least 1 week, typically ≥ 2 weeks, prior to being used as fresh seawater for sample containers in that chamber (see below).
Four of the five C-chambers were connected to a compressed gas source(s) via a BioSpherix sensor-controlled valve. To produce hypoxia, a ProO 2 sensor was attached to 100% N 2 and controlled with a BioSpherix ProO 2 controller. To control CO 2 above ambient (i.e., at 2,000 ppmv), a ProCO 2 sensor was attached to 1% CO 2 /99% N 2 and controlled with a BioSpherix ProCO 2 controller. To control CO 2 below ambient (i.e., at 275 ppmv), a ProCO 2 sensor was attached to 100% N 2 and controlled with a BioSpherix ProCO 2 controller. The sensors monitored atmospheric concentrations of the gas(es) and released gases into the C-chamber as necessary. Controllers are rated to ± 2%, after initial calibration by the manufacturer just prior to experiment initiation. It was established previously that the seawater of plastic containers inside C-chambers equilibrated with the chamber atmosphere within about 40 h (McIntyre-Wressnig et al., 2014). Four treatments (Treatments 2-5) were housed in a climate-controlled environmental room. The environmental room remained dark unless laboratory personnel were tending to the experiment or other cultures housed in the facility. Generally, the environmental room remained dark ∼23 h per day.
The fifth C-chamber lacked sensors/controllers because it was to mimic present-day atmospheric conditions (Treatment 1) and due to resource limitations. Thus, we relied on ambient room conditions rather than stringent controls. This Treatment 1 chamber was ventilated periodically, usually daily, and the pCO 2 was measured approximately every 2 weeks with an infrared CO 2 analyzer (Qubit Systems Inc.).

Treatment Values
Values for the variables were selected to invoke responses, not mortality of the entire assemblage. Although the Mud Patch bottom-water temperature at the time of collection was ∼10.5 • C (Wit et al., 2016), the average in situ bottom-water temperatures on the continental shelf encompassing our sampling site range from ∼4 • C in winter, to 6-7 • C in spring, to 7-8 • C in summer, to 10-11 • C in autumn (Richaud et al., 2016). Thus, our cooler temperature (7 • C) was warmer than average winter temperature and approximately equalled average in situ temperatures for half the year. Our warmer phase (9 • C), while lower than the average ), calcite saturation state ( calcite ) and pCO 2 ) were calculated for the stable portion (24 July 2013-19 March 2014) of the whole experiment (both temperatures) and for the stable portion of each temperature (see text). Species richness (S) and Shannon's diversity (H') are also given. autumn temperature, was higher than in situ temperatures over three seasons.
The value selected for acidified treatments (Treatments 3, 4) was 2,000 ppm CO 2 , a value that is known to be tolerated by other benthic foraminifera (e.g., McIntyre-Wressnig et al., 2013 and is known to be within the range of present day seasonal variability of pore waters at our sampling site (Wit et al., 2016) and from shallow-water sites (< 30 m) in a restricted embayment of the Baltic (Haynert et al., 2012). The Pre-industrial pCO 2 selected was 275 ppmv CO 2 , close to inferred concentrations by climate modelers (e.g., Lerman et al., 2011).
The oxygen value selected for hypoxia treatments was 0.7 mL O 2 /L (=∼29 µM), to ensure that a diverse population survived and grew during the experiment. While some foraminifera can live in lower [O 2 ] (e.g., Nonionella stella, Buliminella tenuata, Bolivina argentea, Globobulimina turgida; Bernhard and Reimers, 1991;Bernhard et al., 1997;Koho et al., 2011), we do not know what extended hypoxia exposure would do to most foraminiferal juveniles, so we selected a mid-range hypoxia value (i.e., Levin et al. (2009) defined the threshold of hypoxia as 1.4 mL O 2 /L or 62.5 µM O 2 ). For treatments that did not have oxygen sensors (Treatments 1, 3, 5), dissolved oxygen concentrations were assumed to be in equilibrium with the treatment gas (or, for Treatment 1, air) and calculated using mass balance of the added gases and the known gas concentrations and temperature. Because these three [O 2 ] levels were all well above the hypoxia threshold defined by Levin et al. (2009), all three were considered "aerated."

Experiment Details
Each treatment was represented by four samples (containers), each comprised of fine-grained sediments (<53 µm) in a separate plastic container. In the lab, the fine-grained sediments that passed through the 53-µm screen at sea were re-sieved again over a 53-µm screen with fresh bottom water immediately prior to aliquot dispensation. While still in suspension, 20 mL of this finegrained sediment-seawater slurry was measured via syringe and placed in each 140-mL plastic container and an additional 120-mL of fresh bottom water was added to each container. Twenty containers were filled in this manner. The 20 containers were randomly divided into 5 sets, each of which was placed into one of the 5 treatment chambers. Once settled, the fine-grained sediment was ∼7 mm in thickness. Salinity of the treatments was monitored weekly with a hand-held refractometer and, if necessary, adjusted to 35 by addition of deionized water across all treatments and containers. Because the chambers were humid, salinity was monitored, and seawater was changed weekly, the system configuration was able to maintain near-constant salinity. All containers of each treatment were fed weekly a concentrated mixture (1 mL) of live Dunaliella tertiotecta and Isochrysis galbana.
Because the experimental design had only one chamber per treatment (i.e., combination of [O 2 ] and pCO 2 ), the experiment lacked true replication of treatments. However, having duplicate or more sets (replicates) of sensors, controllers and chambers was logistically and financially unrealistic. We report on the results of the four containers as the containers provide an estimate of the effects of each treatment on the foraminiferal assemblage.

Carbonate System Parameters
At each sampling, 10 mL was removed for alkalinity analyses and 8 mL was removed for DIC analyses. The total volume removed (18 mL) was replaced with pre-equilibrated seawater poised in the same treatment chamber. Using an automated vacuum extraction system on ∼5 mL samples, DIC concentrations were determined manometrically. Using 1-mL samples and certified reference materials obtained from Dr. A. Dickson (Scripps Institution of Oceanography), alkalinity values were determined by automated Gran titration. Using CO2SYS (Lewis and Wallace, 1998) and methods noted in Wit et al. (2016), alkalinity and DIC values were used to calculate carbonate system parameters (i.e., pH, calcite , [CO 3 2− ], pCO 2 ).

Foraminiferal Analyses
On 13 January 2014 (the 7 • C sampling event), each container was gently stirred and half of each container was removed via syringe and placed in another identical vessel containing appropriate pre-equilibrated seawater including the viability indicator CellTracker R Green CMFDA [5-chloromethylfluorescein diacetate; 1 µM final concentration; hereafter referred to as "CellTracker Green, " Invitrogen/Thermo Fisher Scientific ]. All containers were then returned to the appropriate treatment. After ∼20 h, the CellTracker Green-incubated samples were preserved in ∼3% formalin buffered with sodium borate (Borax R ). Benthic foraminifera were then obtained from these samples as described below. To allow acclimation, the temperature of the environmental room housing four of the chambers was increased 0.5 • C per day over four successive days (14-17 January 2014), to a final temperature of 9 • C. Simultaneously, the temperature of the Nuaire incubator housing Treatment 1 (Present Day) was also raised 0.5 • C each day. As noted, the remaining ∼2.5 months of the experiment was conducted at 9 • C. The day prior to experiment termination on 16 April 2014, CellTracker Green was introduced into each remaining sample container as described above and allowed to incubate in the appropriate treatment for ∼20 h, at which time samples were preserved in ∼3% buffered formalin. Just prior to isolation of foraminifera, preserved samples were sieved over a 63-µm screen with artificial seawater. The > 63-µm fraction was examined with a Leica MZ FLIII epifluorescence stereomicroscope equipped for fluorescein detection (∼480nm excitation; ∼520-nm emission). All Rhizarian protists (foraminifera and their close relative gromiids) were isolated, being recorded as either CellTracker-Green labeled (i.e., fluorescent; live) or unlabeled (i.e., non-fluorescent; dead). Thecate and other "soft-bodied" foraminifera were stored in buffered formalin; "hard shelled" calcareous and agglutinated specimens were air dried and sorted on micropaleontology slides. Species identifications are based on previous identifications of the in situ assemblage (Lang, 2013) and on authors' expertise. Identifications of some morphospecies are tentative, especially for most thecate and rare forms.
The most appropriate count data to analyze across the whole experiment is the total assemblage (fluorescent + nonfluorescent specimens) for the 7 • C time point and solely the CellTracker Green-labeled (fluorescent; live) populations for the 9 • C treatment. We used the combination of 7 • C Total and 9 • C Live because (1) all individuals, including dead specimens, in the 7 • C samples must have grown during the 7 • C exposure; (2) dead specimens in the 9 • C samples had an unknown time of death during the experiment-either during 7 or 9 • C exposure, or during the 4 days of warming.

Statistical Analysis
To determine the similarity of the communities and species trends between treatments, experiment-yield data (specimen counts, entered by morphospecies; live + dead for 7 • C treatments; live only for 9 • C treatments) were square-root transformed prior to statistical analysis with Primer version 7.0.11 (Clarke et al., 2014), with the PERMANOVA+ for Primer add on (Anderson et al., 2008). Data was not standardized prior to analyses because we are generally interested in total instead of proportional responses. Also, diversity measures (Shannon's diversity, species richness) were calculated in Primer.

SSU rRNA Sequencing
Specimens of the ubiquitous thecate morphospecies were isolated for sequencing. DNA from a pooled sample of ∼10 specimens was extracted by homogenization in DOC lysis buffer as described in Pawlowski et al. (1994). Small subunit rRNA gene fragments were PCR amplified using the foraminifera-targeting primers S14F3A (Habura et al., 2004) and Rib sB (Pawlowski, 2000). PCR conditions were: 95 • C for 5 min followed by 35 cycles of 30 s at 95 • C, 1 min at 48 • C, and 2 min at 70 • C with a final 10-min extension at 70 • C. PCR products of the expected size were gel purified using the Zymoclean Gel DNA Recovery Kit (Zymo Research) and cloned into the pCR4 vector in the TOPO TA cloning kit (Invitrogen) according to the manufacturer's instructions. Plasmid DNA from a single clone was sequenced bidirectionally with sequencing primers S14F1 and S17 (Pawlowski, 2000) at Beckman Coulter Genomics. Sequences were edited and assembled into a contig using Sequencher (Gene Codes Corporation) and subsequently deposited into the Genbank sequence database under accession number KY129824.

Experimental Design
The experiment was not replicated due to resource limitations: multiple sets of chambers, controllers and sensors were not available. Further, because two environmental rooms were not simultaneously available, the temperature aspect of this investigation was linked to time. The statistical results should be considered with these caveats in mind.

Carbonate Chemistry
The average carbonate chemistry of the treatments differed as expected (Table 1 and Figure 2); the acidified treatments (3, 4) had pCO 2 values > 3 times higher than those of the remaining treatments (1, 2, 5). However, the carbonate chemistry data from the earliest samples revealed that two treatments (1 and 2) required adjustments. The pCO 2 of the Present Day treatment (1) was too high (evidenced by the relatively high DIC for the initial data points; Figure 2B). It became clear that the environmental room had relatively high pCO 2 due to the presence of the other treatments, where chambers with positive pressure resulted in CO 2 leakage into the environmental room atmosphere, and the lower rate of air exchange in environmental rooms in general. Thus, in Week 6, Treatment 1 was moved into a darkened, temperature-controlled incubator (Nuaire) in our main lab, which had better ventilation and no sources of elevated CO 2 . The average calculated pCO 2 for Treatment 1 after it was moved from the environmental room into the Nuaire was 408 (± 70) ppmv (Table 1). Treatment 2 (Hypoxia) had low DIC concentrations for the initial ∼1.5 months of the experiment (Figure 2B). This low pCO 2 concentration was inadvertently caused by the simultaneous loss of CO 2 while lowering [O 2 ] with N 2 , without CO 2 control. A CO 2 sensor and controller were added to Treatment 2 in Week 6.
Because of the aberrant initial data in Treatments 1 and 2, the average carbonate chemistry data reported in Table 1 is based on readings during the stable portion of the experiment (i.e., from 24 July 2013 to 19 March 2014). Also, suspect readings occurred in the last time point for Treatments 2 and 4; these values are also omitted from calculations presented in Table 1, but all of the data are plotted in Figure 2. Calculated pCO 2 values, over the stable period, indicate that Treatments 3 and 4 were considerably below nominal pCO 2 values (21.6 and 23.4% too low, respectively) while Treatments 2 and 5 were generally above nominal pCO 2 values (12.5 and 14.5% too high, respectively) ( Figure 2D).
The alkalinity and DIC of acidified Treatments 3 and 4 (each nominally 2,000 ppm CO 2 ) trended downward over approximately the last half of the experiment (Figures 2A,B). The calcium carbonate content of the inoculum sediments that were used in the experiment was 4.1% dry weight (Wit et al., 2016). Calculations based on the observed alkalinity and DIC elevations in the first months of the acidified treatments and the number and volume of water changes in each treatment suggest the bulk calcium carbonate of the fine-grained sediments supporting the foraminiferal communities was likely dissolved after ∼5 months, approximately half way through the experiment. Thus, the DIC and alkalinity in these high pCO 2 treatments began to decrease after all sedimentary carbonate was dissolved. Additionally, it is likely that calcification of foraminifera growing in these treatments (as well as other treatments) also impacted the carbonate chemistry by removing carbonate ions.
The calcite saturation states ( calcite ) for the Present Day, Hypoxic, and Pre-industrial treatments (1, 2, 5, respectively) exceeded saturation throughout the experiment while the saturation states of the OA and OA + Hypoxia treatments (3, 4, respectively) each approximated 1 over that time frame (Figure 2C).

Community Responses
Due to the fact that we inoculated with propagules (sediments < 53 µm) and harvested only material > 63 µm, we are confident that the reported assemblages grew during the experiment. Hundreds of specimens were recovered from each container (218-1,210 specimens, live and dead combined). CellTracker Green-labeled (i.e., live) specimens comprised approximately half of the yields with 94-628 live specimens per container. A total of 36 foraminiferal and two gromiid morphotypes were present across the full dataset 1 (Bernhard, 2016b). When considering yields from both temperatures, the Hypoxia treatment (Treatment 2) had the highest foraminiferal yield (live + dead combined), with the OA + Hypoxia treatment (Treatment 4) having the second highest total yield. When considering temperature separately, the highest total yield was in the warm Hypoxia treatment. Generally, the non-acidified warm samples had higher abundances than the comparable cool samples. Species richness indicated no clear trend although the lowest richness was in the triple-stress treatment ( Table 1). Shannon's diversity index also showed no clear trend, although the cool OA + Hypoxia treatment had the lowest value (<0.8) and the highest values were in the Present Day and Pre-industrial treatments ( Table 1). Cluster analysis of the 7 • C Total/9 • C Live assemblages (i.e., morphotype (taxon) counts; all containers) revealed two large clusters along with two singletons and a doublet (Figure 3). The larger clusters roughly assembled into aerated vs. hypoxic groupings. Specifically, 21 of the 24 assemblages from the aerated treatments clustered together while there was more distinction among the hypoxic assemblages, with three major groupings: (1) all four cool OA + Hypoxia assemblages (Treatment 4-7) clustered with one warm OA + Hypoxia assemblage (4-9); (2) all four cool Hypoxia assemblages (2-7) clustered with one warm Hypoxia (2-9) assemblage; (3) two warm Hypoxia assemblages (2-9) clustered with two warm OA + Hypoxia (4-9) assemblages. One warm Hypoxia assemblage (2-9) clustered in the larger Hypoxia cluster while one warm OA + Hypoxia assemblage (4-9) plotted individually.
Count data (7 • C Total and 9 • C Live) from each container was averaged within each treatment to ascertain impacts via non-metric Multidimensional Scaling (nMDS) analysis. A twodimensional nMDS plot shows a clear distinction on the primary axis, between hypoxic treatments (Treatments 2 and 4) and aerated treatments (1,3,5) (Figure 4). Additionally, the secondary axis suggests temperature was also an important driving force. Surprisingly, Pre-industrial, Present Day and cool OA assemblages all plotted close together in this projection, indicating that pCO 2 had less impact on these assemblages than did hypoxia and warming. The stress level for the 2-D nMDS was 0.1, which is the threshold between "good ordination with no real prospect of a misleading interpretation" (Clarke et al., 2014). Thus, examination of this data in a 3-D nMDS, which has a stress value that "gives an excellent representation with no prospect of misinterpretation" (Clarke et al., 2014), confirms that hypoxia and temperature were the major drivers while pCO 2 had a lesser impact (Supplementary Figure 1).
Plotting all containers on a two-dimensional nMDS graph provides the opportunity to observe relative positions of each container (Supplementary Figure 2), noting the projection into two dimensions is only marginally representative of reality. Additionally, the stress level of this 2-D nMDS (0.19) is not robust, which is why additional statistical analyses were executed. Regardless and in general, the four containers of each treatment plot in the same vicinity (Supplementary Figure 2).
A comparison of the experiment yields (7 • C Total and 9 • C Live) indicates that approximately 43% of the taxa resulting from the experiment occur at the sampling site. At the species to genus level, over 58% of the in situ foraminiferal community occurred in the experiment yields. FIGURE 3 | Cluster analysis generated from the Bray-Curtis similarity matrix, which was calculated from square-root transformed abundance data of the 7 • C Total and 9 • C Live counts. The red-dashed portions of branches indicate no significant difference between those connected samples (PRIMER SIMPROF analysis, p < 0.05). Treatment 1 (blue) is Present Day (nominally 7.1 mL/L O 2 , 400 ppmv pCO 2 ); Treatment 2 (orange) is Hypoxia (nominally 0.7 mL/L O 2 , 400 ppmv pCO 2 ); Treatment 3 (green) is OA (nominally 5.9 mL/L O 2 , 2,000 ppmv pCO 2 ); Treatment 4 (red) is OA + Hypoxia (nominally 0.7 mL/L O 2 , 2,000 ppmv pCO 2 ); Treatment 5 (purple) is Pre-industrial (nominally 4.9 mL/L O 2 , 275 ppmv pCO 2 ). Filled symbols represent 7 • C samples; open symbols represent 9 • C samples.

Species Responses
Taxa co-occurrences and relative abundances are best illustrated in a shade plot (Figure 5) presented by treatment. Many species showed clear affinity for treatment and/or temperature. For example, Stainforthia fusiformis and Leptohalysis scotti occurred more often in hypoxic samples compared to aerated samples. Additional comments about particular species appear below.
Two morphospecies grew in each container of each treatment. The most abundant ubiquitous morphospecies was a thin cylindrical thecate foraminifer (Figure 6A), which was typically present in high abundances (145-868 total specimens per sample) with a total of 16,554 specimens. This morphotype often, but not always, dominated samples. Basic Local Alignment Search Tool (BLASTN) analysis of a partial SSU rRNA gene fragment (308 nucleotides) recovered from this ubiquitous thecate morphospecies indicated high similarity to the Allogromiidae, with Micrometula sp. being the top score. Thus, we tentatively identify this taxon as Micrometula? sp. This taxon was more abundant in 7 • C samples vs. 9 • C samples (Figure 6A), and proportionally more dead specimens occurred in warm assemblages compared to cool assemblages. Importantly, however, there were Micrometula? sp. living in all four of the triple-stress assemblages (Treatment 4-9).
The other ubiquitous morphospecies, which rarely dominated samples, was the agglutinated Spiroplectammina sp., which had a total of 539 live specimens and 97 dead specimens. Generally, Spiroplectammina sp. abundances were highest in warm samples compared to cool samples from the same treatments ( Figure 6B). Other common agglutinated species were Leptohalysis scotti, Ovammina cf. O. opaca, and Textularia earlandi. L. scotti had higher abundances and more consistent occurrences in both hypoxic treatments compared to aerated treatments ( Figure 7A). Ovammina cf. O. opaca occurred most commonly in samples from both temperatures of Treatment 1 and some OA samples; it rarely occurred in hypoxic samples (with or without OA; Figure 7B). Each of Spiroplectammina sp., L. scotti and T. earlandi had survivors in each of the triple-stress samples (sub-treatment 4-9). Gromia sp., which is a close relative of foraminifera, occurred in all but one sample, including all triple-stress assemblages; it did best in the warm Pre-industrial treatment (5-9; Supplementary Figure 3).
The calcareous Stainforthia fusiformis was the second most common species in the experiment, with 694 live specimens and 202 dead specimens. This taxon was far more abundant in hypoxic assemblages than in aerated assemblages ( Figure 8A). Nonionella iridea, another calcareous species, was similarly more abundant in hypoxic vs. aerated assemblages (Figure 8B). Both of these calcareous species had considerable numbers of survivors in each of the four triple-stress samples (4-9). Although not as common, live Bulimina marginata, another calcareous species, also occurred in all triple-stress assemblages and in all but two samples, generally, with higher abundances in warm assemblages ( Figure 9A).
Elphidium cf. E. excavatum was most abundant in Preindustrial assemblages, as well as in cool acidified assemblages ( Figure 9B). Lesser abundances of this calcareous taxon occurred in Present Day and cool OA + Hypoxia samples. The calcareous FIGURE 6 | Bubble plots projecting abundances in containers (7 • C Total/9 • C Live) by morphospecies on 2-D nMDS plots. Bubbles are color coded to treatment; lighter shade represents warmer assemblages. Circle diameter is keyed to abundance. (A) Micrometula? sp., inset includes reflected light micrograph of partial specimen and SEM image of apertural end; (B) Spiroplectammina sp., inset is SEM image. Treatment 1 (blue) is Present Day; Treatment 2 (orange) is Hypoxia; Treatment 3 (green) is OA; Treatment 4 (red) is OA + Hypoxia; Treatment 5 (purple) is Pre-industrial.
Cornuspira sp. only occurred in Pre-industrial and Present Day samples plus one cool Hypoxia (2-7) sample; it did not occur in any acidified assemblages (Supplementary Figure 4).
Taxa that did not survive in any of the triple-stress assemblages (4-9) included Allogromia spp., Elphidium cf. E. excavatum, Rosalina cf. R. floridana, Quinqueloculina spp., Miliammina fusca, and Buliminella tenuata. While rare, Quinqueloculina sp. did best in the Warm Pre-industrial and cool Present Day treatments, but also occurred in 7 of the 8 Hypoxia assemblages, never appearing in any high pCO 2 samples. A number of taxa had rare occurrences in general (e.g., Epistominella exigua, Paratrochammina challengeri, Textularia spp.) or were absent in all cool Present Day samples (1-7; e.g., Bolivina spp.), so we do not assess species-specific responses for these cases.

DISCUSSION
Our experiment to assess the response of a continental shelf benthic foraminiferal propagule "bank" to warming, deoxygenation and ocean acidification (OA), as single or combined stressors, was successful in that each sample of each treatment produced substantial populations. The densities produced in the experiment containers (up to 60 specimens cm −3 of sediment) approximated the in situ abundances at the Mud Patch collection site, which were ∼10-44 specimens cm −3 of sediment (Bernhard, 2016a) at the time materials were collected for the experiment.

The Propagule Method
The fact that our propagule-inoculum assemblages differed at the end of experimentation depending on environmental conditions harkens the microbial ecology premise by Baas Becking, "Everything is everywhere but the environment selects" (de Wit and Bouvier, 2006). The length of time foraminiferal propagules remain viable in sediments is unknown, but one study showed that propagules grew after 2 years in adverse conditions (Alve and Goldstein, 2010). Additionally, dormancy, which has been documented in foraminifera (reviewed in Ross and Hallock, 2016), can bolster their persistence as a group.
The use of propagules had three important benefits. First, this approach yielded high numbers of individuals, typically in the hundreds per container. Use of propagules also enabled a short time between collection and experiment initiation, an important factor because keeping foraminifera in the laboratory prior to experiments may be stressful to target populations. And finally, as noted above, the propagule method guaranteed that all adult foraminifera found in the final samples had been alive and growing during the course of the experiment. This eliminated the need to rely on methods such as uptake of a fluorescent marker to mark pre-experimental carbonate (Hintz et al., 2006a,b).
There are some limitations to the propagule method, including initiation with an unknown number and composition of propagules and the necessity for relatively long experiments to allow sufficient growth. In addition, some specimens may be omitted from the inoculum due to their large proloculus size, either due to the taxon's overall size (e.g., Globobulimina) or because they are megalospheric individuals (those resulting from asexual reproduction) (Alve and Goldstein, 2003). Depending on the investigator's objective, the advantages of the propagule method may far outnumber its limitations.
Our result that the Pre-industrial assemblages (Treatment 5) clustered with Present Day and OA assemblages (Treatments 1, 3; Figure 3) indicates that pCO 2 is not the major driving factor distinguishing the assemblage response in this experiment. Clustering of Pre-industrial and OA assemblages suggests that (1) species that grew in those samples were highly adaptable, (2) CO 2 did not have a significant impact on survival and growth of our propagule community or subsets thereof, and/or (3) that culture conditions were not ideal. Prior long-term benthic foraminiferal culturing experiments (e.g., Hintz et al., 2004;McCorkle et al., 2008;Barras et al., 2010;Filipsson et al., 2010) indicate that responses and yields can be difficult to predict, with some taxa seemingly more tolerant to varied conditions than others. A prior long-term experiment that varied pCO 2 (430-3247 µatm) showed a shallow-water Baltic foraminiferal assemblage dominated by Ammonia aomoriensis and Elphidium spp. was not greatly affected by relatively high pCO 2 in terms of abundance, growth and reproduction (Haynert et al., 2014). Since our source material was collected from a site that can have high pore-water [CO 2 ] (4,143 ppmv; Wit et al., 2016), we could expect to have foraminiferal survivors, as local variability can invoke adaptation (Vargas et al., 2017).
There was no universal response of calcareous foraminiferal species to higher pCO 2 : one calcareous species was consistently present, albeit in rather low abundances, in the OA treatments (B. marginata; Treatments 3, 4) but other calcareous foraminifera were not (e.g., Quinqueloculina sp.).

Impact of Hypoxia
Of the three variables, oxygen was the most influential factor driving our foraminiferal assemblage response (Figure 3). While some foraminiferal species, such as Cornuspira sp. and Ovammina cf. O. opaca, only rarely occurred in the hypoxic samples, others, such as Stainforthia fusiformis, Leptohalysis scotti and Nonionella iridea, had higher abundances in hypoxic treatments (2, 4) compared to well aerated treatments (1,3,5). This increased abundance in response to deoxygenation is not surprising given many foraminiferal species are known to survive hypoxia and, sometimes, anoxia (e.g., Moodley et al., 1997;Bernhard et al., 2006a;Risgaard-Petersen et al., 2006;Glock et al., 2019). Stainforthia fusiformis has been documented from low oxygen habitats (Høgslund et al., 2008;Leiter and Altenbach, 2010), as have species of nonionellids (Cedhagen, 1991;Bernhard et al., 1997;Høgslund et al., 2008;Leiter and Altenbach, 2010;Fontanier et al., 2014). At least some species of Nonionella are known to store nitrate and denitrify and it is assumed that some Stainforthia denitrify (Risgaard-Petersen et al., 2006;Høgslund et al., 2008). Some species of Gromia are known to denitrify (Høgslund et al., 2017) and, thus, not be significantly challenged by low oxygen. This is consistent with the result that its highest abundances were in the warm Pre-industrial samples, a wellaerated treatment.

Impact of Warming
Yields in our samples subject to warmer conditions must be considered with due caution because our experiment did not allow any propagule aliquots to grow starting at 9 • C. The survivors of the warm treatments had to survive their cooler treatment, so the yields of the warm treatments may be biased, necessitating survival of the experiment's cooler phase. Therefore, the species surviving our warm treatments may not necessarily be actual "winners" in the face of global warming.
In general, many of the 9 • C sub-treatments had higher numbers of conspecifics compared to their respective 7 • C subtreatments (e.g., Spiroplectammina sp., Nonionella iridea). Higher abundances of a given morphotype in the warm sub-treatments compared to colder sub-treatments may be a species-specific response or may merely indicate a longer time to grow beyond the screen size threshold of 63 µm. Of the common taxa, however, only Micrometula? sp. consistently had fewer specimens in samples from the warmer conditions than in samples from the cooler temperature.
A prior long-term temperature-impact experiment indicated that B. marginata grew in treatments between 3.8 and 21.0 • C (Filipsson et al., 2010). Here, B. marginata occurred in all treatments, although in relatively low abundances, with a slight abundance increase in warmer samples. Furthermore, B. marginata's ability to occur in all five treatments, including all triple stressor containers, demonstrates its adaptability over the range of environmental variables tested here.
Considering the full experiment, the highest average abundance occurred in the cooler Hypoxic + OA containers (Treatment 4-7) suggesting that this combination of dual stressors was additive or beneficial to this foraminiferal assemblage. This highest overall yield was mostly due to the plethora of Micrometula? sp. that occurred in Treatment 4-7. Other species that had relatively high yields after exposure to these two stressors were Stainforthia fusiformis and Leptohalysis scotti, although their highest abundance was in other treatments. The lower pH imposed by high pCO 2 appeared to negatively impact S. fusiformis, most likely because it has a calcareous test. The other common species in Treatment 4-7 samples generally lacked calcareous tests.
The other two dual stressor treatments (Treatments 2-9, 3-9) both involved warming. In Treatment 2 (hypoxia), most of the species that grew in hypoxia were not negatively impacted by the 2 • C warming. Species that were more abundant in the hypoxia samples exposed to warming were Spiroplectammina sp. and Nonionella iridea. These results are consistent with natural distributions, where congeners of both these species occur in the oxygen-deficient Santa Barbara Basin, often in high abundance (Bernhard et al., 1997). When combined with hypoxia, warming had a negative impact on Micrometula? sp., which had fewer specimens in Treatment 2-9 compared to Treatment 2-7, suggesting that warming with hypoxia negatively impacted this thecate species.
Few publications exist to date on the impacts of three stressors to foraminifera. For example, Haynert et al. (2014) performed a laboratory-based experiment to assess impacts of different salinity, temperature and OA on the near-coastal foraminifer Ammonia aomoriensis. Results indicated no impact of salinity on growth, but warmer temperatures partially offset decreases in calcification due to OA. Thus, warming and OA were somewhat antagonistic. In a field-based study, Charrieau et al. (2018a) also documented impacts of the same three stressors (i.e., O 2 , salinity, OA) on a shallow-water benthic foraminiferal assemblage in the Skagerrak-Baltic Sea; the highest abundances and species richness was in well aerated, higher pH and higher salinity sites, implying that assemblage will be challenged by future climate change.
Our experiment had one treatment with three stressors, Treatment 4-9. The foraminiferal assemblage of the triplestressed treatment had the lowest number of species as well as lower abundances for most morphospecies. Many taxa did not occur, even rarely, in any of the four containers of the triple-stressor samples. Some of these morphotypes occurred in some of the hypoxic, high pCO 2 samples (Treatment 4-7; e.g., Elphidium cf. E. excavatum, M. fusca, Quinqueloculina sp.) suggesting that those two stressors challenged those taxa and that the added impact of warming proved insurmountable in this experiment.
Conversely, other morphospecies persisted in the triplestressor sub-treatment (4-9), occurring in all four containers, typically in low numbers. The majority of these taxa were agglutinated (L. scotti, E. advena) or thecate (Gromia sp.) forms although some calcareous foraminifers (S. fusiformis, N. iridea, B. marginata) also occurred in all four triple stressed samples. The observation that some calcareous species tolerated and grew in our hypoxic, warmer, and relatively high pCO 2 treatments suggests that carbonate from one foraminiferal species can be analyzed geochemically for proxies of low oxygen (e.g., Mn/Ca, I/Ca), acidification (B isotopes), and for temperature (e.g., Mg/Ca). While many species grew in the triple-stressor sub-treatment, none of the common morphotypes had their highest average abundance in those samples. For Micrometula, the dual stressors of low oxygen and high pCO 2 were favorable in that the taxon produced its largest population, while the addition of the third stressor (warming) was antagonistic, causing the population to decline. Because Spiroplectammina sp. abundances were generally slightly higher in the warmer samples, including their abundances in the triple-stressor samples (4-9) compared to the hypoxic, high pCO 2 cooler samples (4-7), we infer the triple-stressor combination did not severely impact this morphotype over the time period tested. Thus, it is possible that some foraminiferal species may benefit from combined impacts of warming, OA, and deoxygenation.

CONCLUSION
A culture experiment examining three stressors (oxygen, pCO 2 , and temperature) suggests that many species of benthic foraminifera will be challenged in the next few decades, perhaps to the point of extirpation. However, these negative impacts did not affect all of the species cultured, suggesting that some benthic foraminiferal species in continental shelf settings are likely to survive the twenty-first century, unlike the extinction that has been suggested for some benthic foraminiferal assemblages (i.e., shallow-water tropical settings; Uthicke et al., 2013). If the benthic foraminiferal assemblage cultured here is representative of temperate bathyal settings, those community compositions will be altered from currently occurring taxa in response to OA, deoxygenation, and warming, and will likely include higher proportions of agglutinated and, perhaps, thecate taxa vs. those that precipitate carbonate tests. As has been noted in other marine ecosystems, there will be winners and losers in the benthic foraminiferal response to the multiple stressors predicted to result from anthropogenic global change.

AUTHOR CONTRIBUTIONS
JB and DM designed the experiment. JB and JW collected the samples. JW performed the experiment. JW and DM analyzed the water chemistry. JB, VS, and WP performed the statistical analyses. DB did the sequencing. JB, JW, and DM wrote the manuscript. All authors provided manuscript comments and approved the final version.

FUNDING
This work was supported by the US NSF SEES-OA grant OCE-1219948 to JB and the Investment in Science Program at WHOI. DM also received support from the NSF Independent Research and Development Program.