Cusk (Brosme brosme) and climate change: assessing the threat to a candidate marine fish species under the US Endangered Species Act

Abstract

Hare, J.A., Manderson, J.P., Nye, J.A., Alexander, M.A., Auster, P.J., Borggaard, D.L., Capotondi, A.M., Damon-Randall, K.B., Heupel, E., Mateo, I., O'Brien, L., Richardson, D.E., Stock, C.A., and Biege, S.T. 2012. Cusk (Brosme brosme) and climate change: assessing the threat to a candidate marine fish species under the US Endangered Species Act. – ICES Journal of Marine Science, 69: 1753–1768.

In the Northwest Atlantic Ocean cusk (Brosme brosme) has declined dramatically, primarily as a result of fishing activities. These declines have led to concern about its status, which has prompted reviews under the US Endangered Species Act and the Canadian Species at Risk Act. Changes in distribution and abundance of a number of marine fish in the Northwest Atlantic have been linked to climate variability and change, suggesting that both fishing and climate may affect the status of cusk. Our goal was to evaluate potential effects of climate change on Northwest Atlantic cusk distribution. Coupling a species niche model with the output from an ensemble of climate models, we projected cusk distribution in the future. Our results indicate cusk habitat in the region will shrink and fragment, which is a result of a spatial mismatch between high complexity seafloor habitat and suitable temperature. The importance of habitat patch connectivity for cusk is poorly understood, so the population-level consequences of climate-related habitat fragmentation are uncertain. More broadly, climate change may reduce appropriate thermal habitat and increase habitat fragmentation for other cold-water species in the region; thereby, increasing the potential for regional overexploitation and extirpation.

Introduction

To protect against declines in global biodiversity, there is a growing recognition of the need to include climate information in determinations of whether a species is threatened or endangered (Ruhl, 2008; Sommer, 2009). There are a number of international agreements and national legislations aimed at protecting endangered species, each with different definitions and methods of assessment. Using the US Endangered Species Act as an example, an endangered species is one that is in danger of extinction throughout all or a significant portion of its range, and a threatened species is one that is likely to become endangered within the foreseeable future throughout all or a significant portion of its range (Endangered Species Act, 1973). A species is determined to be threatened or endangered based on one or more of the following five factors: (i) present or threatened destruction, modification, or curtailment of a species' habitat or range; (ii) overutilization for commercial, recreational, scientific, or educational purposes; (iii) disease or predation; (iv) inadequacy of existing regulatory mechanisms; (v) and other natural or man-made factors affecting the species' continued existence.

Climate variability and change can affect any of the five factors, but here we focus on effects of climate change on habitat. Climate change can reduce habitat volume and reduce habitat suitability leading to decreased population growth rates and increased extinction risk (Thomas et al., 2004; Deutsch et al., 2008). Habitat volume and habitat suitability may also increase, leading to increases in population abundance (Hare et al., 2010). Climate change will cause shifts in habitat, and extinction probability will be influenced by species dispersal characteristics, as well as the speed and direction of habitat shifts (Loarie et al., 2009). Climate change also can lead to the fragmentation of habitats and affect the spatial dynamics and resilience of populations (Wilson et al., 2009).

Most studies of endangered species and climate change pertain to terrestrial systems and, overall, the effect in marine systems is less well known (Richardson and Poloczanska, 2008). Changes in habitat and associated species distributions have been projected in marine systems for marine mammals (Kaschner et al., 2011), fishes (Cheung et al., 2009), zooplankton (Reygondeau and Beaugrand, 2010) and phytoplankton (Bopp et al., 2005) (see also Fulton, 2011). Key structural components of marine ecosystems will also be affected by both warming and ocean acidification (Hoegh-Guldberg et al., 2007). In addition, broad scale shifts and changes in the volume of marine biomes are likely (Polovina et al., 2011). Thus, climate change has the potential to dramatically affect the status of marine species, both negatively and positively, and these effects need to be incorporated into the ESA determination process.

Marine fishery species present unique challenges to species conservation efforts. Overexploitation is prevalent through both directed fisheries (Myers et al., 1997) and bycatch (Davies et al., 2009). Fisheries are typically regulated under different legislation than endangered species creating regulatory complexity. In the USA, marine fisheries management is regulated by the federal Magnuson–Stevens Fishery Conservation and Management Reauthorization Act of 2006 (MSFCMA, 2006) and a variety of state laws and cooperative acts (e.g. the Atlantic Coastal Fisheries Cooperative Management Act, 1993). In the past two decades, concern has increased over the extinction risk of marine, estuarine, and diadromous fish species, many of which are impacted by fishing activities (Musick et al., 2000). Species status as determined under fisheries management and species conservation frameworks can be dramatically different. For example, in European waters many fishery species have declined in abundance by more than 70%, which is the criteria for an endangered classification under the International Union for Conservation of Nature framework (World Conservation Union, 2001), yet remain above the biomass level for an overfishing determination (Rice and Legacè, 2007). In some instances, the implementation of fisheries regulations have failed to prevent collapses and failed to provide for species recovery (Myers et al., 1997). In addition, habitats and incidental takes (e.g. bycatch), which are a major focus of endangered species legislation, are negatively impacted by many fishing practices (Auster and Langton, 1999; Barnes and Thomas, 2005). In this context, concern has increased over the status of marine fishery species, and there have been recommendations to classify more marine fishery species as either threatened or endangered (Musick et al., 2000; Dulvy et al., 2003).

Our purpose here is to consider the effect of climate change on cusk (Brosme brosme), a marine fishery species found in moderately deep water throughout the North Atlantic. The US National Oceanic and Atmospheric Administration, National Marine Fisheries Service (NMFS) initiated a status review for cusk due to concerns over population declines and ongoing threats in US and Canadian waters (Federal Register, 2007). Climate change has been identified as a possible factor for threatened and endangered species (Ruhl, 2008); the US portion of the Northwest Atlantic Ocean represents the southern-most distribution of cusk and includes the Gulf of Maine and Georges Bank region. In analysing historical data (1960s to the present), Nye et al. (2009) found that of 36 species studied, more than half shifted distributions poleward or to deeper water, including cusk. To evaluate the effect of climate change on cusk, we develop a species niche model and then couple this model with projections of bottom temperature derived from downscaling global climate models. We use the coupled niche-climate model to project distributions of cusk habitat in the future and to quantify change in habitat area and fragmentation. We then consider the projected changes in habitat in the context of the ESA status review for cusk. An ancillary objective is to develop and document an approach for projecting future changes in species distribution in the northeast US shelf ecosystem.

Material and methods

Cusk biology and status

Cusk are distributed across the North Atlantic from the Northeast US Continental Shelf to the European Shelf (Knutsen et al., 2009). Recent population genetic studies indicate the Northwest Atlantic cusk are distinct from cusk in the Northeast Atlantic. Cusk from Rockall Bank and the Mid-Atlantic Ridge are also distinct, suggesting that deepwater (>1000 m) is a barrier to gene flow (Knutsen et al., 2009). Based on these results, the Northwest Atlantic population should be considered a distinct population segment (DPS) (Lea Harris, Fisheries and Oceans Canada, pers. comm.).

Our study addresses cusk in the Northwest Atlantic only, and specifically the Gulf of Maine, Georges Bank, and the Scotian Shelf, which represents the southern extent of the species range (Figure 1a). In this region, where sampling has been consistent over the last 40 years (Azarovitz, 1981), cusk are found in moderately deep waters and have a preference for waters 120–240 m deep. Bigelow and Schroeder (1953) reported that cusk were found in waters as shallow as 20 m; the cause of this apparent loss of cusk from shallower habitats is not clear, but may be due to declining population sizes (Davies and Jonsen, 2011) or changing environmental conditions (Nye et al., 2009). Throughout their range, cusk use complex habitats including rough areas of rocks and boulders and areas of gravel or pebble bottom (Oldham, 1972; Collette and Klein-MacPhee, 2002). They also have been reported associated with deepwater corals (Husebø et al., 2002; Auster and Lindholm, 2005).

Figure 1.

(a) Map of the occurrence of cusk (Brosme brosme) in the Northeast US Continental Shelf Large Marine Ecosystem. Data were obtained from the Northeast Fisheries Science Center spring and fall trawl survey, which sample the area using a stratified random design (Azarovitz, 1981). Place names used in the text are also provided. (b) Abundance time series from three fishery-independent trawl surveys. (c) Mean length-at-capture from two fishery-independent trawl surveys. (d) US commercial and recreational landings data.

The abundance of cusk in the Gulf of Maine-Georges Bank–Scotian Shelf region has decreased over the past 40 years in both fishery-independent and fishery-dependent surveys. Fisheries-independent trawl surveys in the USA indicate a decline of 75–80% over the past 50 years (Figure 1b). Mean length has also declined in these surveys from the 1960s to the 2000s (Figure 1c). This decline in length has not been continuous and shows a sharp decrease in the early 1990s coincident with relatively steep declines in cusk abundance. Commercial landings in the USA also have decreased (Figure 1d) from more than 2000 metric tons in 1985 to less than 100 metric tons since 2004. Cusk is not managed in the USA so there is no limit on landings, but their value is relatively low so directed effort is minimal. Recreational landings are low compared to commercial landings, but in recent years recreational landings have exceeded commercial landings. Similar trends have been observed in Canada with decreases in both fishery-independent and fishery-dependent indices (Harris and Hanke, 2010). Cusk are managed in Canadian waters with limits on landings (Harris and Hanke, 2010). Some caution must be used in interpreting fishery-dependent data because of changes in fishing effort and regulations over time, but fisheries-dependent and fisheries-independent data from the USA and from Canada indicate declines in cusk over the past 30 years.

The declines in cusk abundance have prompted both the US and Canadian governments to consider whether cusk is threatened or endangered. The USA determination is underway. In Canada, the Committee on the Status of Endangered Wildlife in Canada (COSEWIC) determined that cusk is threatened (Harris and Hanke, 2010; COSEWIC, 2003) and a decision to list cusk under the Species at Risk Act (SARA) is pending. Decreases in survey indices met the criteria for endangered status under the Canadian legislation (>90% decrease in abundance), but the recommendation was for threatened status owing to cusk's widespread distribution and reductions in fishing, which is the main source of mortality. Recently, Davies and Jonsen (2011) indicate some of the decline in cusk in Canadian trawl surveys can be attributed to a change in catchability. They argue that at lower population size the proportional abundance of cusk in non-trawlable areas may be higher, thereby lowering catchability in the Canadian fishery-independent surveys; using a surplus production state-space mode they estimate declines of ∼60%. They point out that their estimated declines would likely still result in a listing as threatened, but the perception of population status is very different compared to a 93% decline. The argument of Davies and Jonsen (2011) likely applies to the US trawl survey data (Figure 1b) and thus, these data may also overestimate the decline in cusk.

General approach

Our goal was to evaluate the potential change in cusk habitat in the Gulf of Maine, Georges Bank, and the Scotian Shelf region as a result of climate change. This goal was accomplished through a four-part effort: climate model downscaling, species niche modelling, projections of future cusk habitat, and analysis of projected habitat distributions (Figure 2).

Figure 2.

Schematic of modelling approach used in this study. First, downscaling estimated the response of regional water temperatures to large-scale climate change forcing. Second, a species-niche model was created for cusk using generalized additive models with bottom temperature and bottom ruggedness as niche dimensions. Third, downscaled bottom water projections and bottom ruggedness maps were used as inputs to the species niche model to project future distributions of cusk habitat. Fourth, projected distributions of cusk habitat were analysed to evaluate possible changes in habitat area and fragmentation with future climate change.

Climate downscaling

The climate downscaling consisted of three parts: (i) development of a bottom temperature climatology from observed data (⁠⁠), (ii) calculation of change between present and future temperatures from global climate models (Δ T), and (iii) estimate of projected bottom temperature in the future by adding change in temperature to temperature climatology (⁠⁠).

Bottom temperature climatology

A bottom temperature climatology from Cape Hatteras, North Carolina to the Laurentian Channel, Canada was developed using two regional hydrographic data sources (Gregory, 2004; Fratantoni et al., 2011) with more than 33000 observations (Figure 3). Near-bottom temperature observations were obtained from hydrographic temperature data recorded within 10 m of the bottom, as verified by 2-minute gridded global relief seafloor terrain data (ETOPO2v2, 2006).

Figure 3.

Bottom temperature climatologies derived from more than 33 000 observations collected from 1977–2009. The six panels show the bottom temperature in two-month periods through the year.

Bottom temperatures were averaged by two-month periods and 0.25° bins. Smaller bin sizes resulted in too many empty cells and larger bin sizes produced much larger variances in the averaged observations. A simple bin average was chosen over more complex gridding and objective analyses for ease, and because the bin size still captures the primary scales of hydrographic variability in the study area; future iterations of regional downscaling should consider more complex methods of calculating a regional climatology. There were fewer observations earlier in the series, so a weighted mean was used to calculate the climatological value of bottom temperature; equal weights were assigned per decade (1977–1986; 1987–1996; 1997–2009) to ensure that the climatology represented observations over the entire time series and not disproportionately for the latter period owing to the greater number of observations.

Temperature change calculation

The “delta” method was used to develop future bottom temperature projections. This approach has been widely used for climate projections in both terrestrial and marine systems (e.g. Akhtar et al., 2008; Fogarty et al., 2008; Anandhi et al., 2011). The method uses the difference between a climate variable in the future (e.g. the mean ocean temperature from 2060–2100, T_2060−2100) and an historical period (e.g. the mean ocean temperature from 1977–2009, T_1977−2009) as estimated from an Atmosphere-Ocean Global Circulation Model (AOGCM). The model-derived difference (ΔT = T_2060−2100 – T_1977−2009) is then added to an observation-derived climatology for the same historical time period to produce a projection of the climate variable.

The delta method removes the mean climate model bias from the projection and, if certain assumptions are met (Stock et al., 2011), provides an assessment of climate-change. The assumptions include: (i) differences between the observed and modelled ocean temperatures arise primarily from model biases in the mean climate state and not differences in the phase of climate variability, (ii) the mean climate state and the projected change are not strongly correlated, and (iii) changes in ocean temperature from the broad-scale changes in radiative forcing and ocean dynamics resolved by AOGCMs are not strongly counteracted by unresolved changes in local shelf-scale dynamics.

Projections based on the delta method were derived from an ensemble of eight AOGCM's to account for inter-model spread in projections. The models were used in the Fourth Assessment Report of the International Panel on Climate Change (IPCC, 2007); the Fifth Assessment Report is expected in 2014. AOGCMs depict the climate using a three dimensional grid over the globe; the resolution of these models is coarse and subgrid scale processes (e.g. turbulence in the boundary layer, thunderstorms and ocean eddies) are parameterized based on large-scale conditions, i.e. variables that are simulated on the model's coarse grid. AOGCMs with ocean resolutions of less than 2° were chosen to maximize the resolution of shelf bathymetry and dynamics, and thus reduce the potential for violating assumption (iii) (Table 1). The selected models had between 29 and 50 vertical levels and depths of between 150–350 m in the Gulf of Maine region. One model was subsequently removed from the ensemble due to a large cold bias and excessive sea-ice estimated for the present period. While the cold bias could be removed by the delta method, the initial presence of ice cover and its subsequent melting under global warming creates an unrealistically strong response in temperature, thereby violating assumption (ii).

Table 1.

Atmosphere-Ocean Global Circulation Models used in this study.

Model	Ocean Resolution [latitude, longitude, depth levels (L)]	B1	A1B	A2
BCCR CM1	1.9° × 1.9° L30	X	X	X
CGCM3.1 T 63	1.4° × 0.94° L29	X	X
MPI ECHAM5	1.5° × 1.5° L40	X	X	X
IAP FGOALS	1° × 1° L33	X	X
GFDL 2.0	1° × 1° L50	X	X	X
GFDL 2.1	1° × 1° L50	X	X	X
UKMet HadCM3	1.25° × 1.25° L30	X	X	X
MIROC Med	1.4° × 0.5° L43	X	X	X

The ocean resolution of each model, and the emissions scenarios available for each model are listed.

Changes in temperature were calculated for three greenhouse gas emission scenarios: low (B1), moderate (A1B) and high (A2); these scenarios are standards used by the IPCC to project future climate change (Nakicenovic and Swart, 2000). These scenarios were developed for the Third Assessment Report in 2001 and also used for the Fourth Assessment Report. All the models that were used here were used in the Fourth Assessment Report. Output was available for only seven models for the A2 scenario. The 30-year climatology and 40-year average projections were used to minimize the effect of assumption (i). Temperature changes were calculated between the period over which the observed climatologies were constructed (1977–2009) and two future periods (2020–2060 and 2060–2100), at 20-m depth intervals. Deltas were calculated over three regions to account for spatial differences in temperature changes and for the six bimonthly periods to minimize the effect of assumption (iii) (see Supplemental Materials). In general, regional and seasonal differences in delta T's were minimal.

Projected bottom temperatures

Regional depth-specific changes in temperature (ΔT) were added to the bottom temperature climatology (⁠⁠) to develop bottom temperature projections (⁠⁠). Depth-specific ΔT's (20-m resolution) were averaged across all climate models to form an ensemble mean depth-specific estimate. These ensemble depth-specific climate projections were mapped to the bottom temperature climatology using linear interpolation. For example, for a grid cell with a 50-m bottom depth, the ΔT was calculated from linear interpolation of the 40- and 60-m projection. When depths were deeper than the ΔT estimates, the deepest ΔT was used.

Analysing the projected bottom temperatures shows that most of the projected variability in warming is associated with the time period and the specific model (Table 2). The ensemble method averages across the inter-model variability and the examination of distinct time periods accounts for the strong influence of time on changes in temperature. There is relatively little variability associated with region and season. The relatively large scale of the climate models (1–2° latitude) limits the resolution of the delta method. But the minimal effect of region and season on the projected temperatures implies that the signal of warming is coherent over the scale of the northeast US Shelf ecosystem.

Table 2.

Sources of variability in the change in bottom temperature from the present to a time in the future.

Factor	Proportion variance explained	Number of levels
Region	0.002	3
Scenario	0.052	3
Time period	0.396	2
Depth	0.012	12
Season	0.000	6
Model	0.183	7
Error	0.354
Total	1.000

Calculations based on a fully orthogonal ANOVA of the estimated changes in bottom temperature (ΔT) including the effects of region (Scotian Shelf, Gulf of Maine, Southern New England), emission scenario (B1, A1B, A2), time period (2020–2060, 2060–2100), depth (0–360 in 20-m intervals), season (six bimonthly intervals), and model. This analysis only included ΔT's from the six Atmosphere-Ocean Global Circulation Models with results for all three emission scenarios (see Table 1).

Species niche model

The development of the species-niche model for cusk consisted of three components: (i) data assembly, (ii) model construction, and (iii) model evaluation.

Data assembly

Response variable

Cusk is considered a data-poor species in terms of stock assessment; there is not enough information to use stock assessment methods to estimate population abundance or fishing mortality (e.g. Quinn and Desiro, 1999). However, there are numerous data sets pertaining to cusk that could potentially be used to develop a species niche model (Table 3). We chose to use cusk presence/absence from two fishery-independent trawl surveys as the response variable. The two surveys have been conducted over 40 years from Cape Hatteras, North Carolina to New Brunswick, Canada (Azarovitz, 1981; Shackell and Frank, 2003). In all, there were 39858 trawl samples, and cusk were present in 2256. In addition, these trawl samples have coincident environmental observations, including temperature and salinity, which can be used in the species niche model. Data from a number of other surveys were considered but were not used here because of limited spatial coverage, lack of corresponding environmental data, or relatively short-time period.

Table 3.

List of potential data sources reviewed for inclusion in the species niche model. All weblinks last accessed date 26 September 2012.

Program	Lead	Years	Area of coverage	More information
Fishery independent longline survey	Maine Department of Marine Resources	2007 and 2008	Coastal Maine waters	2007 report available on request; 2008 report online: http://www.maine.gov/dmr/rm/halibut/08halibutcusk.pdf (2008 report)
Maine/New Hampshire inshore trawl survey	Maine Department of Marine Resources, New Hampshire Fish and Game Department	Fall 2000 through present	Coastal waters of Maine and New Hampshire	http://www.maine.gov/dmr/rm/trawl/index.htm
Lobster sea sampling program	Maine Department of Marine Resources	1985 through present	Coastal Maine waters	http://www.maine.gov/dmr/rm/lobster/lobstersamplingprograms.htm
Cusk mortality study	Maine Department of Marine Resources	2011	Coastal Maine waters	Results pending analysis. Information available on request.
Massachusetts trawl survey	Massachusetts Division of Marine Fisheries	1978 through present (spring and fall)	Massachusetts' territorial waters	http://www.mass.gov/dfwele/dmf/programsandprojects/resource.htm#resource
Massachusetts industry-based survey for Gulf of Maine cod	Massachusetts Division of Marine Fisheries	November 2003 through February 2007	Inshore Gulf of Maine	http://www.mass.gov/dfwele/dmf/programsandprojects/ibs_final_report.htm
SMAST Study Fleet	The School for Marine Science and Technology (SMAST), UMASS Dartmouth	November 2000 through present	Georges Bank	http://www.smast.umassd.edu/Fisheries/Trawler/index.php
Northeast Fishery Observer Program	Fisheries Sampling Branch, NMFS Northeast Fisheries Science Center	1989 through present	Maine through North Carolina	http://www.nefsc.noaa.gov/fsb/

After comparison with our inclusion criteria, these datasets were not used, but there is certainly valuable information regarding cusk ecology which could be useful to future studies.

Predictor variables

Our purpose was to develop a species-niche model defined by environmental variables that could be downscaled from AOGCMs, or that could be assumed constant over the next 90+ years. Initially, variables were considered for inclusion in the species niche model if the physiological or ecological mechanisms behind the species responses were well understood. This initial selection was culled so that only those variables for which field measurements were available at the spatial scale (Cape Hatteras to the Scotian Shelf) and spatial resolution (e.g. grain) of the trawl surveys, served as the primary source of cusk occurrence data. These selection criteria were developed a priori and the resulting list of variables included bottom ruggedness, bottom temperature, bottom salinity, and solar elevation. Solar elevation was included since capture efficiency in trawls can vary as a result of changes in the behavior of animals with photoperiod and time of day (Casey and Myers, 1998).

Importantly, bottom depth was excluded from the analysis to allow depth distribution to change in the future as a result of changes in other habitat variables. Cusk are found from 10–700 m (Collette and Klein-MacPhee, 2002) but have a preference for 120–240 m. However, a mechanism for this depth-dependent distribution is uncertain. There is evidence for depth-dependent physiological differences within marine fishes (Sullivan and Somero, 1980), but these differences may be more related to energy availability and not strictly depth (Vetter et al., 1994; Drazen and Seibel, 2007). Further, several studies have found fishes shifting to deeper waters over time, possibly as a response to increasing temperatures (Dulvy et al., 2003; Nye et al., 2009). Thus, depth was excluded from the species-niche model to allow depth distribution to change in the future as a result in changes in other habitat variables (e.g. temperature). We recognize at the outset that excluding depth from the model will likely result in an overprediction of cusk distribution in shallow waters where adult cusk rarely occur, and this bias was observed in the model (Figure 1a).

The model selection process is described in the Supplemental Material and only methods related to the final model are provided here. Two variables were included in the final species niche model: bottom temperature and bottom ruggedness. Bottom temperature was measured with all the trawl samples included in the analysis. A terrain ruggedness index (TRI, sensuRiley et al., 1999) was selected for each trawl sample location from a grid of TRI estimates, which was calculated from a 15-arc second (∼350 m East–West, 430 m North–South) bathymetric grid that merged Canadian and US soundings. TRI was defined as the square root of the sum of the difference of squared elevations between a grid cell and the neighboring eight cells. The spatial resolution (grain) of the TRI grid was ∼1.3 km² (Figure 4).

Figure 4.

Terrain ruggedness index (TRI) for the Northeast US Continental Shelf and the Scotian Shelf. The upper panel shows the entire region and the lower panel shows the Gulf of Maine region in more detail.

Model construction

Generalized additive modelling (GAM) using the mgcv package in R software was used for the species niche modelling (Wood, 2006). A binomial link function was used in the model and cusk presence/absence was used as the response variable. An iterative approach was used to develop the final species-niche model. The response variable—presence/absence of cusk derived from fishery-independent trawl surveys—was examined combining all trawl surveys and for each survey separately to evaluate any potential differences in cusk occurrence among surveys. All initially selected variables were included iteratively (Stage 1). Response curves of models built using all the data were compared to those from a more conservative dataset of fall trawl data from the USA and summer trawl data from Canada to determine whether seasonal and geographical variation was present in the response of cusk to the environment (Stage 2). Finally, a density-dependent effect in the niche model was evaluated by evaluating the model for two subsets of data (Stage 3): a high abundance period (1970 to 1990) and a low abundance period (1991 to 2008) (Figure 1b).

Initially, the GAM models were constrained to have relatively simple responses of cusk occurrence to the predictor variables (e.g. ≤3 knots). The argument was that species response to environmental factors should be relatively simple (e.g. a bell-shaped curve). However, constraining the partial effects function created oscillations in the error envelop around the function. Thus, an iterative process determined the optimal number of knots to use. Gamma was set to 1.4, increasing the penalty for models of greater complexity (greater degrees of freedom), and a backward selection procedure and analysis of nested models was used to select habitat covariates for the final model (Wood, 2006). Models were compared using analysis of deviance (∼log-likelihood ratio test) and Akaike information criterion (AIC).

Model evaluation

The performance of the species-niche model was quantified using Receiver Operator Characteristic (ROC) analysis and a summary confusion matrix developed from 10-fold cross validation of the model (Fielding and Bell, 1997). For the 10-fold cross validation the 39858 observations were randomly partitioned into 10 nearly equally sized segments that served as independent test sets (Refaeilzadeh et al., 2009). Each test set was then used to predict probabilities of cusk presence using a model trained with remaining observations and habitat data. ROC was then applied to the test set predictions and observations to find an optimal probability threshold for classification of predictions of the niche model to construct a summary confusion matrix and binary habitat maps. To classify predictions, the minimum difference threshold was selected (see Figure S2 in the Supplemental Material); the probability at which the sensitivity (the true positive rate) and specificity of the model (the true negative rate) are equivalent (Jimenez-Valverde and Lobo, 2007; Lobo et al., 2008).

The difference between the sensitivity (the true positive predictions; i.e. presences) and specificity (the true negative predictions; i.e. absences) of the final statistical niche model was minimized at a probability of 0.075 (Figure S2). This minimum difference threshold was therefore used to construct the confusion matrix summarizing model performance (Table 4) and to classify predictions and construct binary cusk habitat maps for analysis of seascape structure. False positives are not necessarily errors but can have an ecological basis (e.g. unoccupied habitat in the case of a depleted species). False negatives are a greater concern – not predicting habitat that is in fact used. This second form of error was quite low (mean = 0.015, 95% quantile = 0.011, 0.017). As indicated above, the minimum difference threshold minimizes the overall error rate. For species conversation, an argument can be made for a threshold that minimizes the occurrence of false negatives only. The method of determining the optimal probability threshold is an area of active debate (Meynard and Kaplan, 2012). Our goal here is to examine the change in habitat relative to current conditions. If future studies aim to quantify habitat area, the sensitivity to thresholding methods should be evaluated.

Table 4.

Confusion matrix showing proportion of test samples successfully and unsuccessfully predicted with the niche model defined by bottom ruggedness and density-dependent response to temperature.

		Observed
		Present	Absent
Predicted	Present	0.04 (0.039, 0.050)	0.23 (0.229, 0.244)
	Absent	0.02 (0.011, 0.017)	0.71 (0.70, 0.72)

Estimates (median, 95% quantile) were derived from 10-fold cross validation tests and a minimum difference threshold probability of 0.075.

Classified predictions were compared with observations to populate a confusion matrix quantifying proportions of successes and failures in the cross-validated test sets. The species niche model was also evaluated qualitatively by comparing predicted habitat distributions to cusk distribution measured in the surveys. To make habitat distribution maps, the range of predictor variables that had “positive effects” on cusk occurrence were used to classify raster maps and to project the individual environmental niche dimensions in space and time. These “positive effects” ranges were defined as the range over which the lower confidence band was greater than zero (i.e. a significant positive effect). Presence/absence data were then overlaid on these maps to evaluate concurrence.

Analysis of cusk occurrence projections

Two approaches were used to evaluate the effect of climate change on cusk distribution. The first approach involved increasing temperatures incrementally starting from the climatology, thereby evaluating the response of cusk presence/absence to a range of temperature increases. The second used ΔT estimated from the AOGCM ensemble projections for the periods 2020–2060 and 2060–2100 and three emission scenarios (B1, A1B, and A2) coupled to the November–December temperature climatology (⁠⁠). This approach allowed for the evaluation of the response of cusk to AOGCM projected climate change.

For both approaches, binary habitat maps were constructed with probabilities of cusk occurrence projected from the niche model using the minimum difference threshold described above. Habitat area and fragmentation metrics were calculated from binary maps using the SDMTools package in R (VanDerWal et al., 2011) to implement methods of FRAGSTATS (http://www.umass.edu/landeco/research/fragstats/fragstats.html) (Mcgarigal et al., 2002). The following metrics were calculated: (i) total surface area of potential cusk habitat, (ii) the total number of distinct habitat patches, (iii) the mean surface area habitat patches and, (iv) patch cohesion, which measured the physical connectedness of habitat patch types. Patch cohesion is an important indicator of the effects of habitat fragmentation on population connectivity, genetics and spatial population dynamics in terrestrial systems (Cushman et al., 2012). In marine systems, increases in fragmentation are linked to decreases in abundance, biomass, and species richness (Godet et al., 2011).

Results

Climate downscaling

The ensemble ΔT calculations project warming in the bottom waters of the Northeast US Continental Shelf that is dependent on time and emission scenario (Figure 5). In the 2020–2060 time period, bottom temperatures are projected to increase by about 1°C across all three emission scenarios. In the 2060–2100 time period, under the B1 scenario (lower emissions), bottom temperatures are projected to increase by ∼1.8°C. Under the A1B and A2 scenario (higher emissions), bottom temperatures are projected to increase by ∼ 2.4°C (Figure 5). The inter-model spread is smaller than the differences between the B1 and A1B/A2 climate scenarios. Analysis of climate model projections hereafter thus focuses on the ensemble mean patterns for each scenario, along with analysis of a range of fixed temperature changes spanning the scenarios.

Figure 5.

Projected changes in the Northwest Atlantic US Continental and Scotian Shelves bottom temperature by time period and emission scenario. Projections were based on an ensemble of eight Atmospheric-Ocean General Circulation Models for two time periods (2020–2060 and 2060–2100) and three emission scenarios (B1, A1B, and A2). Ensemble means and standard deviations are provided.

Species niche model

The statistical niche model for cusk included high and low abundance temperature responses and one bottom complexity response (Figure 6). The temperature response differed among time periods, which was interpreted as a population density effect. Model parameters from the latter time period were used in projections to represent the distribution in the future in low abundance conditions. The model explained approximately 18% of the deviance in the cusk presence/absence data (Table 5) with most of the error associated with projection cusk at locations where they were not observed (e.g. overprediction). This result is expected for a species that is depleted in abundance and for a model that only partially defines the niche space. The modes of the temperature response curves were nearly identical during high and low-density periods (low 7.9°C; high 7.4°C), but breadth along the niche dimension was wider during the early period when cusk were more common (1970–1990: 2.5–12°C; 1991–2008: 5–10.4°C). Bottom complexity had a positive effect on cusk occurrence at values >3 and the response was not different between high and low-density periods. Probabilities of cusk occurrence declined and became variable at TRI values greater 9. Based on the literature, cusk prefer complex habitats (Collette and Klein-MacPhee, 2002), so the variable nature of the response surface and the decline after approximately TRI = 80 are questionable and likely result from the very low capture efficiency of bottom trawls in highly complex habitats (Davies and Jonsen, 2011) and the relative low number of samples in these habitats. For the projections, we used the low abundance model under the assumption that bycatch is continuing.

Figure 6.

Deviance plots derived from the statistical niche model for cusk showing (a) response to temperature during high abundance, (b) response to temperature during low abundance, and (c) response to bottom ruggedness. High and low abundance periods are not shown for bottom ruggedness because there was no significant effect found in the GAM model. The confidence interval bands (2 s.e.) are also shown and where these do not overlap with zero, the effect of the variable has either a positive (above zero) or negative (below zero) effect. Vertical lines in the plots mark the boundaries of positive effects and define “preferred” ranges of the specific variable.

Table 5.

Results showing deviance, partial deviance, and estimated degrees of freedom (DoF) of smoother statistic from the generalized additive model analysis used to construct the statistical niche model for cusk.

Variable	Deviance	Partial deviance	Estimated DoF
Bottom temperature: cusk abundance high	2847.1	2662.5	7.5
Bottom temperature: cusk abundance low			4.4
Bottom complexity	549.1	363.7
Bottom temperature: cusk abundance + bottom complexity	3211.3
Deviance explained %	18.5

Analysis of projections of cusk habitat

Cusk occurrences from trawl collections were confined to deep waters of the Gulf of Maine and the Scotian Shelf, where bottom complexity was relatively high and bottom temperatures remained within the “preferred” range (see Figure 6). The overlap between predicted distribution and observed distribution was greatest for the November–December period (Figure 7). During the other months of the year, potential thermal habitat extended from this region into shallow water and the Mid-Atlantic Bight where cusk are rarely collected. Further, distributions were projected based on the climatologies for all six two-month periods, and the most limited distribution was for the November–December period; temperatures are greatest in the season cycle during this period at the depths cusk prefer. Thus, we chose the November–December period as the basis for our subsequent climate projections on the assumption that habitat distribution was most limiting during this period.

Figure 7.

Distribution of cusk habitat estimated from the statistical niche model compared to actual distribution of cusk from fishery independent trawl surveys. The green area shows distribution of areas with a “positive” effect for cusk occurrence from the statistical niche model. The black dots show the location of trawls where cusk were captured. The surface area of modelled habitat was largest for the (a) May–June climatology and smallest for the (b) November–December climatology. Future projections used the November–December period based on the argument that habitat is most limiting during this season.

Projections of the niche model indicated that habitat distributions can change dramatically with increasing temperatures. Based on fixed temperature increases, half of the surface area classified as cusk habitat disappeared with increases >1.5°C, and patch cohesion exhibited a marked decreased at approximately this value (Figure 8). Temperature increases between 1 and 1.5°C resulted in broad scale habitat fragmentation. With a 2°C increase, patch cohesion and average patch area decrease by 50% (Figure 8).

Figure 8.

Changes in indices for cusk habitat based upon classified maps of statistical niche model projection and fixed temperature increases relative to the November–December bottom temperature climatology. Bars represent the range of values from the analyses of cusk distribution based on the species niche model using the ± 2 standard error GAM predictions. Four indices of habitat are provided: (a) habitat area, (b) patch number, (c) mean patch area, and (d) patch cohesion.

Projections of the niche model based on the downscaling of climate models were comparable to the appropriate fixed temperature increases. The projections for the 2020–2060 period (ΔT = 0.7–1.5° for all scenarios) produced declines in total habitat area and a mean patch size of 40 to 50 percent (Figure 9). Projected changes in bottom temperature were greater for the 2060–2100 period but with greater difference between emission scenarios; thus, the projected changes in habitat were greater than the 2020–2060 period but more variable among scenarios (Figure 10). Under the B1 scenario, the total area of project cusk habitat was 46% of the amount available for present day climatology, and ∼35% of habitat patches were lost (Figure 9). As with the simple fixed temperature increase approach, cusk habitat became less contiguous and confined to the western Gulf of Maine and the eastern Scotian Shelf. Under the A1B and A2 scenarios, habitat area and mean habitat patch area declined to 15% and 32% from present day climatology and patch cohesion also was substantially lower. Under these scenarios, cusk habitat is projected to disappear from the central part of the Gulf of Maine by the end of the century but remain in the western Gulf of Maine.

Figure 9.

Changes in indices for cusk habitat based upon classified maps of statistical niche model projection using the projected future bottom temperatures derived from an ensemble of Atmospheric Ocean General Circulation Models. Bars represent the range of values from the analyses of cusk distribution based on the species niche model using the ± 2 standard error GAM predictions. Four indices of habitat are provided: (a) habitat area, (b) patch number, (c) mean patch area, and (d) patch cohesion.

Figure 10.

Binary maps of potential habitat for adult cusk based upon classified maps of statistical niche model projection using the projected future bottom temperatures derived from an ensemble of Atmospheric Ocean General Circulation Models.

Discussion

The future of cusk in the Gulf of Maine region is in part dependent on future greenhouse gas emissions. Under the higher emission scenarios (A1B and A2), the habitat distribution of cusk will be greatly restricted by the end of the century (∼80% reduction, Figure 9 and 10). In addition, habitat will be more fragmented, potentially decreasing local population viability. Under the lower emission scenario (B1), the decrease in cusk habitat is less dramatic (∼50% reduction, Figure 9). The effects of changing climate increase with time; in the 2020–2060 time period, there are reductions in cusk habitat (30–40% reduction), and in the 2060–2100 time periods these reductions continue to increase (50–80%). There is no quantitative threshold for habitat loss under the ESA, but qualitatively, based on the analyses here, climate change poses a significant threat to cusk in the 50–100 year time period and at temperature increases >1.5°C. That said, the generally coarse resolution of present climate models, regional model biases, and the potential aliasing of low-frequency climate variability into the climate change signal, suggests that these analyses be viewed as initial estimates subject to future refinement.

The consistency of results using a large ensemble of climate models, multiple climate change scenarios, and multi-decadal averages provides confidence that our projections reflect current best estimates of the magnitude and uncertainty of regional changes in bottom temperature over the next century. Our use of an ensemble of climate models is designed to ameliorate the uncertainties associated with individual models. Several modes of climate variability may impact the projections presented here (North Atlantic Oscillation [NAO], Atlantic Multidecadal Oscillation [AMO]). NAO is a primary driver of variations in ocean temperature in the North Atlantic and exhibits variability across a wide spectrum of time-scales, including the multi-decadal scales analysed herein (Hurrell, 1995). AMO is an approximately 70-year pattern of sea surface temperatures across the North Atlantic (Kerr, 2000). The climatology, even though a 30-year average, could be biased by these natural climate modes, which could potentially result in an underestimate or overestimate of change in our projections.

The potential impact of unresolved physics on the bottom temperature projections is more difficult to assess. The seasonal heat budget of the study region is largely determined by atmosphere-ocean heat exchanges (i.e. short- and long-wave radiation, sensible and latent heating) (Umoh and Thompson, 1994; Mountain et al., 1996). Climate models capture broad-scale patterns in these fluxes and the impact of greenhouse gas accumulation in the atmosphere is primarily reflected in a change in the long-wave flux. However, dynamic fluctuations in bottom temperatures on inter-annual to decadal time-scales can be driven by changes in heat transport on the shelf and between the shelf and surrounding waters (Petrie and Drinkwater, 1993; Mountain, 2003). Our understanding of these fluctuations is limited. The near equal regional-scale bottom temperature increase between the A1B and A2 scenario is notable; global warming is greater under the A2 scenario. One explanation is changes in regional physics (e.g. Labrador Current transport, Gulf Stream transport) act to partially counteract the effect of higher greenhouse gas emissions under the A2 scenario. Additional research is needed to assess whether unresolved changes in heat transport can significantly alter the strong broad-scale warming signal. This will require continued development of high-resolution climate models that can better resolve shelf-scale processes and high-resolution regional simulations that can be forced by climate projections.

Our projections are based on a simple model estimating the response of cusk to two-niche dimensions and based on trawl survey data, which does not effectively sample cusk in their preferred complex bottom habitat (see Davies and Jonsen, 2011). Based on the niche concept, a species distribution is the spatial realization of an n-dimensional hypervolume, with the dimensions defined by environmental factors or ecological processes (Morin and Lechowicz, 2008; Holt, 2009). We considered only bottom temperature and bottom complexity explicitly. Our results support the concept of density-dependence in species distribution (the temperature response differed between periods of high and low abundance), which implies ecological processes such as competition for structural refugia within regions of optimal thermal habitat may also influence distribution (MacCall, 1990; Auster and Lindholm, 2005). Predation, competition, consumption, growth, recruitment and other processes may also interact with climate change to affect the distribution of cusk (see Kordas et al., 2011). For example, the effect of increased habitat fragmentation may be small as larval cusk are planktonic and may be exchanged among isolated groups of adults (Cowen and Sponaugle, 2009). Thus, our results should be viewed as a first order approximation of changes in the distribution of cusk habitat, and subject to future refinement. For a more complete list of research needs see the Supplemental Material.

The projected effect of climate change on bottom temperature and its effects on cusk habitat extent and fragmentation raises serious challenges to the NOAA mandate under the ESA, as well as under the MSFCMA. As has been argued for the US Fish and Wildlife Service's responsibilities under the ESA (Ruhl, 2008), it is highly likely that climate change will be a factor causing extirpation, and possibly extinction, of some species under NOAA's jurisdiction, even with the protections of the ESA. By analogy, some species will remain overfished even with regulations under the MSFCMA. Countries such as Canada (SARA), Australia (Commonwealth's Endangered Species Protection Act) and the European Union (Species Directive) with similar legislation to conserve biodiversity will also need to consider how to confront these challenges scientifically, legally, economically, and socially. Cusk is a good example of many of these challenges. Cusk are mainly taken as bycatch in the New England groundfish and American lobster trap fisheries. Decreases in cusk habitat availability and increases in habitat fragmentation as a consequence of climate change may exacerbate population decreases resulting from fishing activities (see Mieszkowska et al., 2009 for an overview of the effects of climate and fishing on Atlantic cod). We are not proposing remedies to these challenges; rather we hope that identifying the challenges will spur an evaluation of current management and regulations in the context of climate variability and change.

Cusk also is a transboundary species, occurring in both US and Canadian waters of the Northwest Atlantic Ocean. In Canada, the Committee on the Status of Endangered Wildlife in Canada (COSEWIC) assessed cusk as threatened in 2003 (COSEWIC, 2003) and the decision to list cusk under the Species at Risk Act (SARA) is currently pending. Cusk is also taken in several fisheries including longline and lobster trap (Harris and Hanke, 2010), and Powles (2011) noted a disconnect in the treatment of cusk under fishery and conservation legislation. Without knowing the degree of connectivity between cusk in US and Canadian waters, differential regulation under conservation or fishing legislation could hamper the management goals of one or both countries. That said, our modelling effort projects that cusk habitat will increase in the Scotian Shelf waters; if cusk occupy this new habitat, then abundances in these areas could increase.

The issues related to cusk (e.g. fishing, climate change, habitat, international jurisdictions) have implications for NOAA's responsibilities to many marine fish species under both ESA and MSFCMA. Although our projections are only a first order approximation of changes in the distribution of cusk habitat, and are not directly coupled with extinction risk models or fishery assessment models, the direction and magnitude of projected changes indicate that climate change is going to have a negative effect on the status of cusk in the future. Incorporating this finding into regional and international fisheries management and species conservation efforts is necessary and will be challenging.

Funding

PJA was supported, in part, by funds from the National Marine Sanctuary Foundation and NOAA Stellwagen Bank National Marine Sanctuary. EH was supported by a grant (to PJA) from The David and Lucile Packard Foundation Conservation and Science Program. Support for AC was provided by National Science Foundation grant 1544901 as part of the synthesis phase of the GLOBEC Program.

Supplementary Data

Supplementary data are available at ICES Journal of Marine Science online.

Acknowledgements

We thank the participants of the National Marine Fisheries Service Endangered Species Act/Climate Change Working Group for their comments and insights during the development of this manuscript, including Michele McClure, Erin Sweeney, Jeff Jorgensen, and others. We also thank James Scott (ESRL) for his assistance accessing AOGCM model output. We thank the following people for commenting on earlier drafts of the manuscript: Roger Griffis, Matthew Poach, Harvey Walsh, and two anonymous reviewers. Acknowledgement of the above individuals does not imply their endorsement of this work; the authors have sole responsibility for the content of this contribution. Meetings among co-authors were supported, in part, by NOAA's North Atlantic Regional Team, a collaborative network of NOAA employees and close partners representing the agency with their diverse capabilities from Virginia through Maine.

Author notes

‡

Present address: 65 Winterberry Rd, Saunderstown, RI 02874, USA