Past, Present, and Future: Performance of Two Bivalve Species Under Changing Environmental Conditions

Steeves, Laura E.; Filgueira, Ramón; Guyondet, Thomas; Chassé, Joël; Comeau, Luc

doi:10.3389/fmars.2018.00184

ORIGINAL RESEARCH article

Front. Mar. Sci., 25 May 2018

Sec. Marine Fisheries, Aquaculture and Living Resources

Volume 5 - 2018 | https://doi.org/10.3389/fmars.2018.00184

This article is part of the Research Topic Coastal and Marine Aquaculture in Climate Change View all 3 articles

Past, Present, and Future: Performance of Two Bivalve Species Under Changing Environmental Conditions

$\r\nLaura E. Steeves*$ Laura E. Steeves¹^*

Ramón Filgueira²

Thomas Guyondet³

Joël Chassé³

Luc Comeau³

¹Biology Department, Dalhousie University, Halifax, NS, Canada
²Marine Affairs, Dalhousie University, Halifax, NS, Canada
³Science Branch, Fisheries and Oceans Canada, Gulf Fisheries Centre, Moncton, NB, Canada

Globally, the production of marine bivalves has been steadily increasing over the past several decades. As the effects of human population growth are magnified, bivalves help provide food security as a source of inexpensive protein. However, as climate change alters sea surface temperatures (SST), the physiology, and thus the survival, growth, and distribution of bivalves are being altered. Challenges with managing bivalves may become more pronounced, as the uncertainty associated with climate change makes it difficult to predict future production levels. Modeling techniques, applied to both climate change and bivalve bioenergetics, can be used to predict and explore the impacts of changing ocean temperatures on bivalve physiology, and concomitantly on aquaculture production. This study coupled a previously established high resolution climate model and two dynamic energy budget models to explore the future growth and distribution of two economically and ecologically important species, the eastern oyster (Crassotrea virginica), and the blue mussel (Mytilus edulis) along the Atlantic coast of Canada. SST was extracted from the climate model and used as a forcing variable in the bioenergetic models. This approach was applied across three discreet time periods: the past (1986–1990), the present (2016–2020), and the future (2046–2050), thus permitting a comparison of bivalve performance under different temporal scenarios. Results show that the future growth is variable both spatially and interspecifically. Modeling outcomes suggest that warming ocean temperatures will cause an increase in growth rates of both species as a result of their ectothermic nature. However, as the thermal tolerance of C. virginica is higher than M. edulis, oysters will generally outperform mussels. The predicted effects of temperature on bivalve physiology also provided insight into vulnerabilities (e.g., mortality) under future SST scenarios. Such information is useful for adapting future management strategies for both farmed and wild shellfish. Although this study focused on a geographically specific area, the approach of coupling bioenergetic and climate models is valid for species and environments across the globe.

Introduction

With a growing human population expected to surpass 9 billion by 2050, the cultivation of bivalves helps to provide food security as a source of inexpensive protein (Shumway et al., 2003; Godfray et al., 2010). Globally, the production of bivalves has been steadily increasing over the past several decades and is expected to continue to do so (FAO, 2016). Considered ecosystem engineers as well as keystone species (Gutiérrez et al., 2003; Zippay and Helmuth, 2012; Han et al., 2017; Sorte et al., 2017), bivalves interact with their environment through both top-down and bottom-up processes (e.g., Coen et al., 2007; Rice, 2008). Top-down control via filter-feeding may significantly curtail phytoplankton populations (Cranford et al., 2003; Newell, 2004; Forsberg et al., 2017) potentially affecting bivalve performance itself (Dame and Prins, 1998; Bacher et al., 2003; Strohmeier et al., 2005), but also impacting other filter-feeders and grazers (Kluger et al., 2017). Filtration activity can also play an important role in regulating water quality and depth of light penetration (Gallardi, 2014; Guyondet et al., 2015; Petersen et al., 2016). Bivalves can also impose bottom-up control on plankton communities by altering fluxes of nutrients (Menge, 1992; Newell, 2004). Bivalves change geographical patterns in energy transfers from the pelagic environment to the benthos, since feces and pseudofaeces contribute to augmenting organic matter directly below farms (Newell, 2004; Cranford et al., 2007, 2009) and filtration acts in reducing organic matter deposition elsewhere in the system (Guyondet et al., 2015).

Interactions between bivalves and their environment are bidirectional, meaning the abiotic environment also imposes effects on bivalves. Temperature plays an important role in physiology, gene expression, distribution, and fitness of bivalves (Zippay and Helmuth, 2012; Shelmerdine et al., 2017). The internal body temperature of bivalves usually matches external water temperatures, except when intertidal species are subject to aerial exposure (Zippay and Helmuth, 2012). Although the effect of temperature on physiology is species-specific, generally as temperature increases, physiological rates will increase until a threshold is met, at which point performance will decline (Kooijman, 2010). Due to this relationship between physiological functions (e.g., filtration rates and oxygen consumption), short-term changes (days/weeks) in temperature will predictably impact their survival (e.g., Malham et al., 2009; Rinde et al., 2016) and long-term changes (years) will impact their reproductive timing and effort, and consequently their spatial distribution (Bayne et al., 1976; Philippart et al., 2003; Kittner and Riisgård, 2005; Toupoint et al., 2012; Zippay and Helmuth, 2012; Filgueira et al., 2014). Accordingly, the relationship between temperature and physiology of bivalves is particularly relevant within the context of climate change.

The Intergovernmental Panel on Climate Change (IPCC) has reported that since 1971 global ocean surface temperatures have increased on average 0.11°C (±0.02 SD) per decade (Pachauri and Mayer, 2015). Although temperature will not be the only change to impact bivalves, it has been shown to be the most deterministic factor influencing shellfish growth and distribution, and is one of the most widely studied abiotic factors related to climate change (Zippay and Helmuth, 2012; Rodrigues et al., 2015; Buckley and Huey, 2016; Filgueira et al., 2016). Warming oceans are modifying current natural ranges of many marine species (e.g., Diederich et al., 2005; Jones et al., 2009, 2010; Rinde et al., 2016; Shelmerdine et al., 2017). Despite inherent uncertainty, modeling techniques which incorporate both climate data and organismal bioenergetics are the only tools available to explore the effects of future climate change scenarios on animal populations.

Climate models are quantitative representations of natural processes that make up Earth's conditions, and are often used to predict the effects of climate change (Pachauri et al., 2014). Currently, global emissions of CO₂ are the best predictor of Earth's surface warming, and are directly related to both human population and economic growth (Pachauri et al., 2014). Outputs from these models include estimated surface, air, and water warming, ice cover, and change in circulation patterns (Pachauri et al., 2014). Refining the scale of global climate models promotes understanding how climate processes and conditions will change on local scales at highly detailed spatial resolutions. This has been done for the Northwest Atlantic shelf region of Canada (Brickman and Drozdowski, 2012), integrating atmospheric and oceanic information to estimate future sea surface temperature (SST) and salinity of the Scotian Shelf and Gulf of Saint Lawrence (Long et al., 2016).

Regarding bioenergetics, Dynamic Energy Budgets (DEB) provide a mathematical method for modeling energetic flows through individual organisms (Kooijman, 2010). DEB models breach interdisciplinary boundaries by merging the principles of thermodynamics, physiology, and theoretical biology. The mechanistic nature of DEB models permits its application to a wide range of environmental conditions. DEB has been parameterized for several bivalve species (e.g., Pouvreau et al., 2006; Crassostrea gigas; van der Veer et al., 2006; Macoma balthica, Mya arenaria, Cerastoderma edule, Mytilus edulis, Crassostrea gigas; Filgueira et al., 2014; Crassostrea virginica) and used to predict their growth (e.g., Lavaud et al., 2017; Crassostrea virginica), and reproductive effort (Montalto et al., 2016; Brachidontes pharaonis, Mytilaster minimus, Mytilus galloprovincialis). The coupling of climate and growth models is being used under the context of climate change to explore the effect of predicted temperatures on the performance and distribution of organisms (e.g., Sarà et al., 2011, 2013; Brachidontes pharaonis; Thomas et al., 2011, 2015; Mytilus edulis, Crassostrea gigas; Klinger et al., 2017).

The development of the bivalve aquaculture industry is contextualized by climate change and its concomitant effects on the oceans, given its reliance on natural environmental conditions. Climate change is generating uncertainty around future production levels of bivalves, which makes creating resilient government policies increasingly difficult (e.g., Rodríguez-Rodríguez and Roberto, 2017). On the Atlantic coast of Canada, it has been recognized that SST is increasing at a rate higher than the global average (IPCC, 2013; Saba et al., 2016). As these changes are predictable, impacts to cultivated bivalve species can be anticipated, and their effects could be mitigated with management plans. In the present study, the future growth of two widely cultivated bivalve species in Atlantic Canada (Mytilus edulis and Crassostrea virginica) are estimated by coupling their bioenergetics to a high resolution climate model. In this way, bivalve growth and performance can be predicted for the coming decades, to prepare for the impacts of climate change on the industry. The findings of this study are relevant for planning bivalve farming, in terms of both site- and species-selection.

Methods

Description of Study Area

The study area for this research is the Scotian Shelf and Gulf of Saint Lawrence, the bodies of water surrounding Nova Scotia (NS) and Prince Edward Island (PEI) (Latitude: 42.7130–49.0416, Longitude: −67.1065 to −59.0403). The study area contains widely varying temperatures both seasonally and spatially. Colder waters in the Bay of Fundy reach an average summer high of 13°C, however water temperatures in the sheltered Northumberland Strait (e.g., the body of water between NS and PEI) may exceed 20°C (Feindel et al., 2013). Furthermore, this area is extensively used for bivalve aquaculture, primarily the eastern oyster (Crassostrea virginica), and the blue mussel (Mytilus edulis). Culture methods consist primarily of long-lines, floating bags, and oyster bottom culture (DFO, 2003a,b).

Climate Change Model

The climate change model used in this study was produced by the Department of Fisheries and Oceans Canada, as a part of the Canadian Government's Aquatic Climate Change Adaptation Services Program. Unlike other Global Climate Models (GCMs), this is a high resolution (1/12°) model of regional climate dynamical downscaling system of the Gulf of Saint Lawrence, Scotia Shelf, and Gulf of Maine (Long et al., 2016). The model domain, for the use in this research was restricted to the waters surrounding NS and PEI (Latitude: 42.7130–49.0416, Longitude: −67.1065 to −59.0403). The model is constructed of the atmospheric Canadian Regional Climate Model (CRCM), and the oceanic model the Canadian Océan PArallélisé (CANOPA) model. The CANOPA model was developed at the Bedford Institute of Oceanography (Brickman and Drozdowski, 2012), based on the Océan PArallélisé model, version 9 (OPA 9.0; Madec et al., 1998), and the Louvain-la-Neuve ice model, version 2 (LIM2; Fichefet and Morales Maqueda, 1997; Bouillon et al., 2009). The model covers the time period from 1970 to 2100 under the A1B and Representation Concentration Pathway (RCP) 8.5 scenario (Long et al., 2016). The output is produced in grid cells with dimensions of 5–6 km horizontally. River inputs are included in the model; however tidal forcing is not used. A full description of the model can be found in Long et al. (2016)

Bioenergetic Model: Dynamic Energy Budget (DEB)

The bivalve models used in this study are based on DEB theory (Kooijman, 1986, 2010). DEB models define how energy is moved through an individual by describing them with three standard state variables: reserves, structure, and maturity/reproduction (depending on the life stage). Energy moves as organisms assimilate food, which is first stored as reserves. Subsequently, a fraction of reserves (κ) is directed toward growth and maintenance of the structure of the individual, and the rest (1 − κ) is allocated to reproductive maturity in juveniles, and reproduction in adults (specifically, gametes). The notations and symbols used in this text are consistent with Kooijman (2010) where square brackets [] denote quantities per unit structural volume, braces {} denote quantities per unit surface-area of the structural volume, and rates are defined by dots above their symbol. The model equations are described briefly in Table 1, and further details of the model can be found in Pouvreau et al. (2006). The models were run for a full year, beginning on January 1st, and were initialized with the same length for both species (2 cm) and dry flesh mass (DFM), including structural weight only, of 0.013 g and 0.02 g for C. virginica and M. edulis respectively.

TABLE 1

Table 1. Equations used in the Dynamic Energy Budget (DEB) model.

The parameterization of the DEB models (Table 2) followed existing studies, Lavaud et al. (2017) for C. virginicia and Rosland et al. (2009) and Saraiva et al. (2011) for M. edulis. The models were calibrated using the scaled functional response (f) as a simplified proxy for food availability. This parameter originates from the Holling Type II response used in DEB theory (Kooijman, 2010):

\begin{array}{l} f = \frac{X}{X + X_{k}} \end{array}

Where X is food density (or availability) and X_k is the half-saturation coefficient, which represents the time and energy an individual allocates to searching for food. As the value of f moves from 0 to 1, increasingly less time and energy is spent looking for food. Using this relationship, f was simplified as a proxy ranging between 0 and 1 for food availability where low values of f (i.e., f = 0) reflect low food availability, and high values of f (i.e., f = 1) reflect food availability at saturation (Kooijman, 2010). This proxy does not incorporate potential negative effects of high food availability. This proxy was used to calibrate the growth rates of both C. virginica and M. edulis. Calibration was completed to ensure that growth rates were within biologically reasonable ranges. This was done by comparing C. virginica and M. edulis growth from the literature with model outputs (Table 3). Growth rates were calculated following Clausen and Riisgård (1996):

\begin{array}{l} μ = (\frac{W_{t}}{W_{0}}) x t^{- 1} \end{array}

Where W_t and W₀ are the average dry weight or shell length on Day 0 and Day t respectively.

TABLE 2

Table 2. Standard DEB parameters for Crassostrea virginica (from Lavaud et al., 2017) and Mytilus edulis (from Rosland et al., 2009¹; Saraiva et al., 2011²).

TABLE 3

Table 3. Review of observed growth rates of Mytilus edulis and Crassostrea virginica used to validate the DEB model output.

Coupled Model

The DEB models for C. virginica and M. edulis were coupled to the high resolution ocean climate model (Long et al., 2016) following an off-line scheme. The forcing variable used from the climate model was seawater temperature. Temperature data was extracted from the model between the years of 1980 and 2050 for each grid square in the region of interest, and then the top two depth measurements (0 and 12 m) were daily-averaged. Model outputs were available twice a month, and linear interpolation was used to estimate data points between observations. Three periods were studied by averaging 5 year periods, 1985–1990, 2016–2020, 2046–2050, representing the past, present, and future scenarios respectively. Averages were used to minimize the impact of potential outliers.

Data Analysis and Numerical Experiments

To determine how temperature changed both spatially and temporally, past SST data (1986–1990) was subtracted from future SST data (2046–2050), for the entire region being analyzed. This was done for two periods during the year, by averaging the temperature for January and August, the coldest and warmest months, respectively. To determine how thermal stress may change over time, the number of consecutive days exceeding physiologically relevant thermal thresholds (°C) were calculated for comparison between the three times periods. This was done for warm temperatures relevant for the upper thermal threshold of M. edulis and cold temperatures relevant for the lower thermal threshold for C. virginica. The upper thermal threshold used for M. edulis in this DEB model is 23°C (Rosland et al., 2009), however increased mortality rates have been observed in laboratory conditions at 22°C (Clements et al., 2018). Similarly, the lower thermal threshold of C. virginica used for this model was 2°C (Lavaud et al., 2017), however behavioral changes (e.g., filtration rates decrease by 50%) have been observed at 9°C (Comeau et al., 2008). This was done for the warmest and coldest areas relevant for bivalve aquaculture contained in the model domain. For this, only coastal areas were considered. To explore the species-specific effects of SST changes, growing degree days (GDD) were calculated for both C. virginica and M. edulis across the three time periods, in the warmest area within the model domain relevant for bivalve aquaculture. GDD were calculated as follows:

\begin{array}{l} GDD = \sum_{i = s t a r t d a t e}^{e n d d a t e} (T_{i} \geq T_{T L}, T_{i} \leq T_{T H}) * Δ d \end{array}

Where T_i is the SST (°C) on day i, and T_TL and T_TH are the predetermined lower and upper threshold temperatures respectively (°C), and was a determine time-step (1 day). GDD is in units of °C × day, and was calculated for 1 year (January 1–December 31). Lower and upper thermal tolerances for C. virginica were 2°C and 35°C (Lavaud et al., 2017), respectively, and −1.8°C and 23°C (Saraiva et al., 2011) were used for M. edulis.

To estimate the physiological performance of both species, shell length (SL), dry weight (DW), and gonadosomatic index (GSI), were estimated for the present time period after 1 year of growth, for the warmest and coldest areas relevant for bivalve aquaculture within the model domain. GSI was calculated as a ratio of reproductive tissue to total dry weight. To determine how these indices changed spatially and temporally, percent change was calculated for both SL and GSI in the entire model domain as follows:

\begin{array}{l} \frac{F u t u r e V a l u e - P a s t V a l u e}{P a s t V a l u e} x 100 % \end{array}

To estimate changes in phenology across the model domain, the earliest spawning dates of both species were estimated using a combination of temperature and GSI thresholds, wherein both conditions must be met for a spawning event to occur. The temperature thresholds for C. virginica and M. edulis were 17°C (DFO, 2003a), and 14°C (Newell and Moran, 1989), respectively. Both of these temperature thresholds reflect relatively conservative estimates. Additionally, GSI thresholds of 0.2 (Choi, 1992) and 0.28 (Troost et al., 2010) for C. virginica and M. edulis, respectively. On the first day that both thresholds were met, spawning was triggered in the model by emptying reproductive reserves.

All statistical analyses were performed in Rstudio version 3.1.2. For all parametric analyses, tests for normality and homogeneity of variance were performed with Shapiro-wilk and Barlett's tests respectively. No data transformations were required; all parametric assumptions were met. For GDD comparisons, 2-way ANOVAs were run. When factors yielded significant effects, post-hoc testing was done using a Tukey test. All α levels were 0.05. No significant interactions between factors were found.

Results

Climate Model

The climate model indicates that SST warming will be spatially dependent within the model domain. Differential warming rates were observed between the past (1986–1990) and the future (2046–2050) scenarios (Figure 1). For example, some areas experienced average August temperatures of 2.5°C higher in the future, than in the past (Figure 1). Smaller changes were observed in January, with the highest absolute change being 1.5°C. SST warming may be more impactful to the upper thermal limits of M. edulis than the lower thermal limits of C. virginica. In the warmest coastal area within the model domain, the number of consecutive days exceeding the upper thermal threshold of M. edulis (23°C; Rosland et al., 2009) increased over the three time periods (between past and future for all temperatures p < 0.05) (Figure 2A). In the past (1986–1990), 6 days in a row exceeding 22°C were expected, compared to the future, where 39 consecutive days could exceed this temperature (Figure 2A). For temperatures closet to the lower thermal threshold of C. virginica (2°C; Lavaud et al., 2017), no significant changes were observed in the number of consecutive days falling below these temperatures (between past and future for all temperatures, p > 0.05; Figure 2B). These were calculated for the coldest coastal area relevant for bivalve farming within the model domain. Further linking temperature and physiology, in the warmest coastal area GDD increased between the past and the future for both species (p < 0.05, Figure 3). In the coldest coastal location, no changes in GDD were observed for either species (Figure 3).

FIGURE 1

Figure 1. Difference in average temperature (°C) between the past (1986–1990) and the future (2046–2050) for the Gulf of Saint Lawrence and Scotian Shelf areas within the model domain. Temperatures were calculated by averaging the SST values for the month of January (A) and August (B), the coldest and warmest months, respectively.

FIGURE 2

Figure 2. (A) Number of consecutive days that exceeded M. edulis' upper thermal threshold (21°C) (Rosland et al., 2009). Temperature data was taken from the warmest coastal pixel within the model domain. Bars not sharing the same letter indicate significantly different groups (p < 0.05). (B) Number of consecutive days falling below C. virginica's lower thermal threshold (8°C) (Lavaud et al., 2017). Temperature data was taken from the coldest coastal pixel in the model domain. Within each temperature group, no significant effect of time was found (p > 0.05).

FIGURE 3

Figure 3. Growing-degree days for C. virginica and M. edulis for each time period analyzed. Upper thermal tolerances were 35 and 23°C, and lower thermal tolerances were 2 and −1.8°C, respectively. This was done for the warmest coastal pixel within the area analyzed. Bars not sharing the same letters indicate significantly different groups (p < 0.05).

DEB Model Validation

The performance of the DEB models for both M. edulis and C. virginica were assessed by comparing growth rates (calculated using shell length) to those observed in the literature for similar geographic regions as those used in this study (Table 3). Using the present output (2016–2020), M. edulis growth rates (0.00233 ± 3.65 × 10⁻⁵ cm/day; mean ± SD) matched those observed in the literature (0.00265 ± 0.00137 cm/day) (Table 3). For C. virginica, growth rates produced for the present (0.00198 ± 5.1 × 10⁻⁵ cm/day) were below those observed in the literature (0.00340 ± 9.2 × 10⁻⁴ cm/day) (Table 3).

Coupled DEB-Climate Model

The performance, as measured by DW, SL, and GSI of both species varied by both speices and location (Figures 4A–H). Generally, the performance of both species was greater (i.e., higher growth rates) in the area with warmer temperatures. Performance as measured by SL and GSI of both species was also visualized for the entire geographic region (Figure 5). The predicted temperature increases are likely to have more observable effects on C. virginica than M. edulis in terms of both length and GSI, over 1 year of growth (Figure 5). In the first year of growth, it was found that C. virginica did not meet the GSI thresholds needed to spawn (Figures 4D,H) either in the present or the future. For M. edulis, the spawning date occurred earlier in the future than in the past, with the largest differential being observed in the Bay of Fundy region (Figure 6).

FIGURE 4

Figure 4. Plots of dynamic energy budget (DEB) output indicating C. virginica and M. edulis growth in dry weight (B,F), length (C,G), and gonadosomatic index (D,H). Over 1 year (beginning on January 1st). This was replicated for areas with two differing thermal regimes, where (A–D) are model outputs from an area with higher average water temperatures than (E–H).

FIGURE 5

Figure 5. Difference (%) in length (SL) (A,B) and gonadosomatic index (GSI) (C,D) of C. virginica and M. edulis between the past (1986–1990) and the future (2046–2050), calculated after 1 year of growth.

FIGURE 6

Figure 6. The spawn date of M. edulis indicated in Julian day for the past (A) (1986–1990) and the future (B) (2046–2050). Spawn date was calculated as combination of temperature (14°C) (Newell and Moran, 1989) and GSI (0.28) (Troost et al., 2010) thresholds. Spawn dates above day 365 indicate that spawning did not occur in the first year of growth.

Discussion

Water temperature impacts the physiology of bivalves in species specific ways for C. virginica and M. edulis. Downstream effects of temperature changes on these species have the potential to impact distribution, performance, and phenology (Zippay and Helmuth, 2012; Shelmerdine et al., 2017). In this study, a high resolution climate model was coupled to bioenergetic models of C. virginica and M. edulis to explore the differential effects of predicted SST changes on these commercial species in the coming decades. By simulating the growth of both of these species over a year, throughout different time periods, the effect of temperature on performance in terms of dry weight, shell length, and reproductive effort was examined for both species.

Climate Model

Results from the climate model indicated variable rates of warming within the model domain. Temporally, between the past (1986–1990) and the future (2046–2050), SST warmed for all areas. For seasonal extremes (January and August) in terms of absolute warming, larger changes are predicted for the summer compared to the winter (Figure 1). Spatially, SST increased differentially; in January, the southern Scotian Shelf region experienced the most warming, and in August higher rates of absolute warming were observed in coastal areas around northern NS and the northern shore of PEI (Figure 1). These results match predictions for Canada's mid-high latitude waters in the Scotian Shelf and Gulf of Saint Lawrence (Feindel et al., 2013; Bush et al., 2014). The climate model used in this research has a high spatial resolution, however it does not capture inner coastal bays where aquaculture is carried out. Although it is difficult to create climate change models at high enough resolution to capture the temperatures in these bays (e.g., Stobart et al., 2016), adding 1°C to coastal temperatures in PEI may provide a closer estimation to current bay/inlet temperatures (Filgueira et al., 2015). This variability in SST warming both spatially and temporally indicates that climate change, in terms of ocean warming, will have variable impacts on bivalve bioenergetic processes.

DEB Model

DEB models, as mechanistic tools, currently require calibration to the local environmental conditions where they are applied (e.g., Bernard et al., 2011; Picoche et al., 2014). The application of DEB to future scenarios is then limited to the availability of projected climate data, namely temperature and food availability. Inherently, these estimations are impacted by uncertainties as outlined by Skogen et al. (2014): scenario uncertainty (the unknown future socioeconomic climate), model uncertainty (flaws within model estimations), and internal uncertainty (inability to detect change until variability of a signal flattens out). Although these uncertainties have been shown to decrease as models are applied to more local scales (Hawkins and Sutton, 2009), they still limit the ability to integrate environmental data into locally calibrated DEB models. This may explain the lower growth rates produced by the C. virginica model in this study, compared to those observed in the literature (Table 3). The DEB model used in this study was originally calibrated using field data from the Gulf of Mexico (Lavaud et al., 2017). The scaled functional response was the only parameter calibrated in this study. Potential ecophysiological (Dittman et al., 1998) or genetic (Murray and Hare, 2006) variability from oyster populations from different latitudes could be missed in the current parameterization of the model, and consequently explain these lower growth rates. This lower growth may also be exacerbated by the conservative temperature estimate of coastal bays, as discussed above.

Although methods for estimating some environmental variables (e.g., temperature) are clearly defined, others (e.g., food availability) are less clear. Defining food availability for bivalves is frequently done using proxies such as phytoplankton concentration (e.g., Riisgård et al., 2012), often represented by chlorophyll (e.g., Lesser et al., 2010). These methods for defining food availability are associated with inherent model uncertainty for current conditions (Smith, 1980), and become even more uncertain when estimating future chlorophyll a concentrations (Elliott et al., 2005). Filgueira et al. (2016), using a spatially explicit model combining the physical environment, aquaculture practices, and climate change drivers, have suggested that climate change (increased SST) may cause a decrease in chlorophyll a concentration in coastal embayments, where bivalve aquaculture is present. There is also potential for the abundance of specific size classes of phytoplankton to shift temporally (Agirbas et al., 2015), indicating the possibility for both the quantity and quality of food availability to change over time. Although temperature has been suggested to be the most deterministic factor for bivalve performance under climate change conditions (Filgueira et al., 2016), forcing predicted food availability onto the bioenergetic model would provide a more comprehensive understanding of how bivalve performance will change over time. As predicting chlorophyll a in climate change modeling is inherently difficult, food availability was held ad libitum for the entire year in this model, following the methods of Lavaud et al. (2017). By limiting the forcing variables in the model to temperature, the results of this study are restricted to the impact of temperature on bivalve physiology.

Coupled Climate-DEB Model

Predicted changes in SST coupled to bivalve physiology have shown to benefit to C. virginica over M. edulis in terms of growth, due to their differing thermal physiologies. For M. edulis a maximum SL increase of 8.7% was predicted for the future (2046–2050), compared to 16.0% for C. virginica (Figure 5). Growth rates of these ectothermic species are highly temperature dependent (Zippay and Helmuth, 2012; Feindel et al., 2013), and although growth rates are predicted to increase, thermal stress associated with ocean warming must be considered. The upper thermal limit of M. edulis used in this model was 23°C (Saraiva et al., 2011), but temperatures above 20°C are associated with increases rates of mortality (Newell and Moran, 1989; Mallet and Myrand, 1995), and significantly reduced growth rates (Gonzalez and Yevich, 1976). Specifically, in the southern Gulf of St Lawrence, mortality associated with summer heat stress has recently been observed (Myrand et al., 2000). Although mortality is not predicted by the model, an analysis of predicted temperature could be used as an indicator of mortality risk. For mussels, the number of consecutive days per year at temperatures above 22°C would increase from ~7 in the past (1986–1990) up to ~50 in the future (2046–2050) (Figure 2A). The potential for increased mortality rates at this temperature is highly relevant; in a recent laboratory experiment, mussels held at 22°C for a 30 day period experienced significantly higher mortality rates than a control group held at 16°C (Clements et al., 2018). This suggests significantly increased mortality rates for M. edulis in the future, compared to past conditions (Figure 2A). Note that, as discussed above, mortality risk would become even higher in aquaculture areas, that is, sheltered bays which cannot be captured with the spatial resolution of the climate model. In this way, these predictions are a conservative estimate of mortality risk in highly coastal, sheltered bays.

Contrastingly, C. virginica has the potential to grow to market size at an expedited rate under future scenarios (Figure 5), with no adverse effects from high temperatures due to their high thermal tolerance. Although there was no significant decrease observed in the number of consecutive days below temperatures relevant to the lower thermal threshold of C. virginica (Figure 2B), increased growth rates are still likely due to warming summer temperatures. C. virginica growth becomes observable in water temperatures around 9°C (Shumway, 1996), a temperature threshold which is likely to be met earlier in the year in the future. Food availability (recorded as timing of phytoplankton blooms) is currently not temporally synced to the seasonal initiation of C. virginica growth. In waters contained within the model domain, spring phytoplankton blooms occur in water temperatures between 4 and 9°C, creating a mismatch between the timing of peak food availability and the seasonal initiation of C. virginica growth (Pernet et al., 2007; Comeau et al., 2008; Feindel et al., 2013). As the temporal match between food availability (phytoplankton blooms) and temperatures suitable for growth could become more closely linked under climate change scenarios, the production potential for C. virginica would be further benefitted.

In addition to the effects on growth and mortality, SST can also affect phenology, in terms of both reproductive timing and distribution (Thackeray et al., 2010; Feindel et al., 2013). Warming ocean temperatures and the expanded northward distribution of the oyster Crassostrea gigas has been extensively studied in European waters (Laugen et al., 2015; Shelmerdine et al., 2017). Spawn dates were estimated in this study using a combination threshold of both temperature and GSI, a common method to trigger spawning in DEB modeling (e.g., Bourlès et al., 2009). The application of this threshold suggested that in the future M. edulis will spawn earlier in the year (Figure 6B), and there is potential for oysters to demonstrate a similar trend, however they did not spawn in the first year (Figure 5). For all areas in which spawning occurs, the dates fall within the currently observed spawning times for the Gulf region (between May and August) (DFO, 2003b). GSI calculations (Figure 5) indicate that for both species, proportionally more energy will be allocated to reproduction in the future, although a greater increase was observed in C. virginica.

Limitations

The predictions of these DEB models are forced by temperature, the most significant variable for bivalve ecophysiology (Filgueira et al., 2016); however, the effects of other drivers such as ocean acidification, food availability and ecosystem dynamics should not be ignored (Feindel et al., 2013). Ocean acidification, although difficult to predict, has been incorporated into DEB models (e.g., Klok et al., 2014). The waters included within this model have experienced changes in their pH, showing an average pH decline of 0.1–0.2 units since the 1930s (Stewart and White, 2001). Gledhill et al. (2015) have also suggested that the coastal waters of NS have a reduced buffering capacity due to significant freshwater inputs. Acidification has the potential to negatively affect fertilization, larval settlements, and spat shell formation (Curren and Azetsu-Scott, 2013; Gurney-Smith, 2015). Negative impacts have been reported on M. edulis larvae in terms of slow growth rates and shell deformities (Gazeau et al., 2010, 2013), and also for recruitment (Brown et al., 2016). For C. virginica and M. edulis in North America, negative impacts to shell calcification are expected under marine pH conditions predicted for 2050 and 2100, respectively (Gazeau et al., 2007; Whitman-Miller et al., 2009). Despite this, most studies show few impacts of elevated CO₂ on growth or mortality rates of adult bivalves (Keppel et al., 2015; Clements et al., 2018). The impacts of acidification on larval mortality may be particularly important for the aquaculture industry in NS and PEI, where spat are primarily harvested from the wild and therefore rely on unbuffered water (Feindel et al., 2013). However, rearing spat in hatcheries may help to mitigate the impact of acidification on bivalve aquaculture (Clements and Chopin, 2016).

Culture conditions, where bivalves are grown at high densities are also not represented in these models, which are individual based. Local effects caused by high culture density, such as competition, can have significant impacts on bivalve performance (Cubillo et al., 2012). The primary culture method of bivalves in NS and PEI is in-shore, which increases the risk of thermal stress, and requires cultivation at high densities due to limiting space. Off-shore aquaculture avoids high summer temperatures experienced in coastal bays, potentially reducing thermal stress on M. edulis (Myrand et al., 2000; DFO, 2017). Additionally, less competition for space could permit cultivating bivalves at lower stocking densities. Although off-shore aquaculture lacks the physical protection of sheltered bays, by lowering long-lines several meters into the water column, abrasive wave action can be avoided, and more stable temperatures can be achieved (DFO, 2017; Klinger et al., 2017). Preliminary results from M. edulis grown off-shore in Newfoundland indicated that off-shore growth rates were comparable to in-shore growth, and that spawning occurred less often, but was more predictable (DFO, 2017). Deeper cooler waters have the potential to avoid thermal stress for M. edulis under climate change conditions, compared to sheltered inlets (Figure 2A). Compounding stressors such as low food availability, reproduction, and summer heatwaves can act to simultaneously increase mortality rates. For example, previous summer mortalities of M. edulis has been explained by the combined effects of high water temperatures, low food availability, and high reproductive output (Tremblay et al., 1998; Myrand et al., 2000). For C. virginica, as mentioned, it should also be considered that this model was calibrated in the Gulf of Mexico (Lavaud et al., 2017), with a differing thermal regime to the one in the model domain. Not accounting for possible differences in the ecophysiology of C. virginica latitudinally may introduce error into the model.

Conclusions

Results of this study, although subject to the uncertainties previously discussed, indicate that SST in coastal waters will undergo differential rates of warming both spatially and seasonally. Due to the different thermal physiologies of C. virginica and M. edulis, these predicted changes in SST will create species specific risks and opportunities in terms of growth and phenology. Summer heat stress may pose an increased threat to M. edulis, due to their lower thermal tolerance. This can impact performance in terms of growth, but also potentially increase the risk of mortality. C. virginica, with a higher thermal tolerance, is unlikely to experience negative impacts to growth rates, but instead is likely to show an increase in performance. Additionally, as SST warms there may be an increased temporal match between food availability and the growing period of C. virginica, which may in turn impact growth and phenology. When planning for the future of cultivating bivalves in NS and PEI, these results provide information on which areas will become stressful for M. edulis, and concomitantly where opportunities will arise to cultivate C. virginica more effectively. The continued success of this industry is dependent upon the careful selection of species and farming sites. To avoid increased temperature related mortality of M. edulis from thermal stress, off-shore aquaculture could be considered as a cultivation method. Furthermore, C. virginica could be opportunistically grown in areas where warming SST could promote increased growth rates compared to colder areas. To build upon these results, and reduce uncertainties, future studies should incorporate environmental variables such as seasonal food availability and ocean acidification. Factors such as these can be explored under climate change scenarios, and further integrated into bioenergetic models. bioenergetic models.

Author Contributions

Climate model data was provided by JC, and DEB modeling expertise provided by RF. All authors contributed to the improvement of the manuscript.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Funding

This work was prepared as part of a graduate research project at Dalhousie University, Canada. Funding was provided through a graduate NSERC scholarship to LS, and an NSERC discovery grant to RF.

References

Agirbas, E., Martinez-Vicente, V., Brewin, R. J. W., Racault, M.-F., Airs, R. L., and Llewellyn, C. A. (2015). Temporal changes in total and size-fractioned chlorophyll-a in surface waters of three provinces in the Atlantic Ocean (September to November) between 2003 and 2010. J. Mar. Syst. 150, 56–65. doi: 10.1016/j.jmarsys.2015.05.008