Decadal-scale phenology and seasonal climate drivers of migratory baleen whales in a rapidly warming marine ecosystem.

Species' response to rapid climate change can be measured through shifts in timing of recurring biological events, known as phenology. The Gulf of Maine is one of the most rapidly warming regions of the ocean, and thus an ideal system to study phenological and biological responses to climate change. A better understanding of climate-induced changes in phenology is needed to effectively and adaptively manage human-wildlife conflicts. Using data from a 20+ year marine mammal observation program, we tested the hypothesis that the phenology of large whale habitat use in Cape Cod Bay has changed and is related to regional-scale shifts in the thermal onset of spring. We used a multi-season occupancy model to measure phenological shifts and evaluate trends in the date of peak habitat use for North Atlantic right (Eubalaena glacialis), humpback (Megaptera novaeangliae), and fin (Balaenoptera physalus) whales. The date of peak habitat use shifted by +18.1 days (0.90 days/year) for right whales and +19.1 days (0.96 days/year) for humpback whales. We then evaluated interannual variability in peak habitat use relative to thermal spring transition dates (STD), and hypothesized that right whales, as planktivorous specialist feeders, would exhibit a stronger response to thermal phenology than fin and humpback whales, which are more generalist piscivorous feeders. There was a significant negative effect of western region STD on right whale habitat use, and a significant positive effect of eastern region STD on fin whale habitat use indicating differential responses to spatial seasonal conditions. Protections for threatened and endangered whales have been designed to align with expected phenology of habitat use. Our results show that whales are becoming mismatched with static seasonal management measures through shifts in their timing of habitat use, and they suggest that effective management strategies may need to alter protections as species adapt to climate change.

INTRODUCTION

Species' adaptation to changing environmental and ecological conditions is imperative for their survival in a rapidly changing climate. One way to measure species' adaptive capacity to novel conditions is through shifts in the timing of recurring biological events, known as phenology; including changes in migration, peaks in productivity, and breeding or spawning events (Parmesan & Yohe, 2003). Shifting phenology has been well studied in terrestrial animals, with the majority of studies showing earlier timing of spring events including leaf‐out and flower blooming, insect emergence, and songbird migration (Fu et al., 2015; Mayor et al., 2017; Parmesan & Yohe, 2003). Relatively, fewer studies have shown definitive shifts in phenology in marine systems, but where responses have been documented, earlier timing is most prevalent, particularly for small‐bodied, lower‐trophic level species (e.g., phytoplankton, zooplankton) during the growing season (Cohen et al., 2018; Poloczanska et al., 2013; Staudinger et al., 2019). Body size and the ability to thermoregulate are key traits that affect the rate and magnitude of phenological responses to changes in seasonal temperature (Cohen et al., 2018). Disparate responses among taxa and trophic levels are widely predicted to decouple historical ecological relationships and cause resource mismatches that affect growth, survival, and population dynamics (Anderson et al., 2013; Bewick et al., 2016; Both et al., 2006; Cushing, 1990).

A major gap in our understanding of shifts in marine phenology exists for higher trophic levels, particularly species that make long‐distance seasonal migrations (Edwards & Richardson, 2004; Kellermann & van Riper, 2015; Staudinger et al., 2019). Long‐distance seasonal migrants are thought to be disproportionately at risk due to differential climate change impacts in departure and destination habitats (Cherry et al., 2013; Middleton et al., 2013; Otero et al., 2014; Wood & Kellermann, 2015; Youngflesh et al., 2021).

Understanding climate‐induced changes in marine phenology is becoming increasingly urgent in the Gulf of Maine region, which is one of the most rapidly warming marine ecosystems in the world (Mills et al., 2013; Pershing et al., 2015). Sea surface and bottom temperatures of the continental shelf have increased at rates three times faster than the global average over recent decades (Saba et al., 2016). Warming has been strongest during the summer and autumn seasons, with onset of the growing season for primary and secondary producers (e.g., phytoplankton and zooplankton) occurring earlier in the year, leading to longer summer and shorter winter seasons (Thomas et al., 2017). Future projections for the region predict that warming trends will continue (Pershing et al., 2021; Saba et al., 2016), with additional climate impacts on coastal ecosystems from increased stratification and changes in ocean circulation (Brickman et al., 2021). The Gulf of Maine is known for its high seasonal productivity (Townsend, 1991), which attracts a diversity of migratory fish and wildlife; these animals time their arrival into specific habitats to take the advantage of seasonally occurring nutrient‐rich prey and other resources to support growth and reproduction (Staudinger et al., 2019). Predictable timing of habitat use is especially important for endangered, threatened, and otherwise vulnerable species including large migratory whales, specifically, North Atlantic right whales (Eubalaena glacialis, hereafter “right whales”), humpback whales (Megaptera novaeangliae), and fin whales (Balaenoptera physalus).

Right whales are critically endangered (NOAA, 2008a), and the most recent median estimate of their population size is 336 (95% confidence range ± 14) (Pettis et al., 2022). They were recently uplisted from Endangered to Critically Endangered on the International Union for Conservation of Nature (IUCN) Red List of Threatened Species (Cooke, 2020). The uplisting was due to a decrease in population size (Pace et al., 2017), an increase in vessel strikes and entanglements in fixed fishing gear (Daoust et al., 2017; Davies & Brillant, 2019; Knowlton et al., 2012; Sharp et al., 2019), and a decrease in birth rate (Pettis et al., 2021). Right whale migratory phenology within the Gulf of Maine—documented over several decades (Brown et al., 2007)—describes large aggregations of right whales in Cape Cod Bay during winter and early spring (Ganley et al., 2019; Mayo et al., 2018; Pendleton et al., 2009), the Great South Channel from late spring to summer (Kenney et al., 1995; Pendleton et al., 2009), and the Bay of Fundy and Roseway Basin from late summer through autumn (Kraus et al., 1982; Patrician & Kenney, 2010; Winn et al., 1986). In response to changes in temperature, circulation, and the supply of right whales' primary copepod prey, late‐stage Calanus finmarchicus. Except for Cape Cod Bay (Ganley et al., 2019; Mayo et al., 2018), right whales have reduced their use of most of these important feeding habitats (Meyer‐Gutbrod et al., 2021; Record et al., 2019). While right whales' reduced use of most traditional habitats has been observed, there has been an increased use of geographically marginal habitats such as the Gulf of St. Lawrence (Crowe et al., 2021; Davies & Brillant, 2019), southern New England waters (O'Brien et al., In Review; Quintana‐Rizzo et al., 2021), and the mid‐Atlantic waters of the United States (Davis et al., 2017). Based on high‐resolution climate models, the redistribution of right whales across the western North Atlantic Ocean is projected to continue at least through the year 2050 (Ross et al., 2021).

Humpback and fin whales also utilize Cape Cod Bay as a seasonal habitat (Clapham et al., 1993; Clapham & Seipt, 1991; Stamieszkin et al., 2010). All marine mammals in United States waters are protected under the Marine Mammal Protection Act. However, under the United States Endangered Species Act, humpback whales in this region are not listed as threatened or endangered, while fin whales are listed as endangered throughout their range. Winter use of Cape Cod Bay by humpback whales is not completely understood, but it likely reflects the habitat‐use patterns of the wider Gulf of Maine; sightings of humpback whales are sparse from January through March, as the majority of the population undertakes a southward migration (Clapham et al., 1993; Katona & Beard, 1990; Kenney et al., 1981). The migratory pattern and winter distribution of fin whales is largely unknown, and they are rarely seen during the winter months in the Gulf of Maine (Hain et al., 1992; Seipt et al., 1990). Both whale species feed on capelin (Mallotus villosus), herring (Clupea harengus), sand lance (Ammodytes spp.), and krill (Meganyctiphanes norvegica) (Payne et al., 1990; Simard et al., 2002; Weinrich et al., 1997), but most observations point to sand lance as being their primary prey in the southwestern Gulf of Maine, especially in the Stellwagen Bank area, adjacent to Cape Cod Bay (Overholtz & Nicolas, 1979; Silva et al., 2021; Staudinger et al., 2020). Sand lance are thought to determine the local distribution and abundance of both species in this region (Kenney et al., 1981; Payne et al., 1990).

Better information on the strength and direction of phenological changes in long‐distance seasonal migrants, particularly those of high conservation concern, is needed to mitigate threats from climate change and conflicts with human activities (Mills et al., 2013; Peer & Miller, 2014). The timing of large whale occurrence in and around Cape Cod Bay is of high concern because this region is heavily used by commercial shipping, fishers, whale watch vessels, and recreational users. In particular, Cape Cod Bay has remained an important and persistent habitat for right whales, even as their abundance has declined and their distribution has shifted in recent years (Ganley et al., 2019; Pace et al., 2017). Vessel speed limits and restrictions on certain fishing practices within Cape Cod Bay (MADMF, 2021; NOAA, 2008b, 2013) were designed to coincide with the historical peak in right whale occurrence while benefiting other species such as humpback and fin whales. However, shifts in the timing of right and other whales could reduce the effectiveness of protective measures if not accounted for in planning and management (Davies & Brillant, 2019; Meyer‐Gutbrod et al., 2018).

Our study sought to improve the understanding of the magnitude, direction, and drivers of large whale phenology in Cape Cod Bay. We did this by testing the overall hypothesis that the phenology of large whale habitat use in Cape Cod Bay has changed, and that it is related to regional‐scale shifts in seasonal temperature indicators due to climate change. Using a novel multi‐season occupancy model parameterized to measure phenological shifts, we tested for trends in peak habitat use in Cape Cod Bay for right, fin, and humpback whales over 21 years, from 1998–2018. We then evaluated interannual variability in peak habitat use relative to thermal spring transition dates in the Gulf of Maine. We hypothesized that right whales, which are planktivorous specialist feeders, would exhibit a stronger response to bottom‐up indicators of phenology (e.g., temperature and changes in primary and secondary production) than fin and humpback whales, which are more generalist and piscivorous feeders in the western Gulf of Maine. Based on previous studies showing earlier arrival and departure times for humpback and fin whales in the Gulf of St. Lawrence (Ramp et al., 2015), and suggestive trends of earlier presence of right whales in Massachusetts Bay (Charif et al., 2020), we tested for evidence of temporal shifts in habitat use for these three species in this important seasonal habitat.

MATERIALS AND METHODS

Marine mammal aerial surveys

Between 1998 and 2018, the right whale researchers at the Center for Coastal Studies (Provincetown, Massachusetts) flew line‐transect aerial surveys over Cape Cod Bay, Massachusetts (Figure 1). They used a twin engine, high wing, Cessna Skymaster 336 or 337 airplane. Surveys occurred as weather allowed between December or January and mid‐May, and were aborted if observers determined the weather conditions deteriorated to a Beaufort sea state greater than 5 or visibility less than 3.7 km. In all, 15 east‐west tracklines, spaced 2.8 km apart, were flown at a minimum altitude of 229 m, and a ground speed of approximately 185 km/h.

FIGURE 1

Cape Cod Bay study region located in Massachusetts, USA. East‐west aerial survey tracklines are shown as thin black lines. Tracklines are separated by 2.8 km. Shaded contour regions, from light to dark, are shown for 0–10, 10–25, 25–50, 50–75, 75–100, and >100 m depths. Inset shows the Gulf of Maine with our study area outlined in black. Green and orange areas show the western and eastern Gulf of Maine regions, respectively

Between 1998 and 2002, the survey crew consisted of one pilot, one data recorder, and two observers. In 2003, the crew configuration was altered for increased flight safety to two pilots, but the two observers maintained the same protocols during surveys. Observers, one on either side of the airplane, scanned the water surface outward in an arc‐shaped swath beginning from the forward limit of visibility to directly abeam of the aircraft and from the downward limit of visibility out to at least 3.6 km (Brown et al., 2007; Brown & Marx, 1998). Observers recorded the time, location, and behaviors of all marine mammals in their viewing swath (Mayo et al., 2018). Data from these surveys were curated by and obtained through the North Atlantic Right Whale Consortium (NARWC, 2019).

Environmental data

We obtained bathymetric data from the United States Geological Survey at 15‐arc seconds (~0.46 km) spatial resolution (Roworth & Signell, 1998), and we resampled those data to a 4.6 km spatial resolution. We obtained sea surface temperature (SST) and chlorophyll‐a data from multiple sources to span the 21 years of our study. We used AVHRR Pathfinder Version 5.2 (PFV5.2) SST data for years 1998–2002, which we obtained from the United States National Oceanographic Data Center and GHRSST (https://pathfinder.nodc.noaa.gov). We used SST data from the Moderate Resolution Imaging Spectrometer (MODIS) instrument on the Aqua satellite for years 2002–2018. We used chlorophyll‐a data from the Sea‐viewing Wide Field‐of‐view Sensor (SeaWiFS) for years 1998–2002. To avoid periodic gaps in coverage caused by SeaWiFS outages, we used chlorophyll‐a data from the MODIS instrument on the Aqua satellite, which came online in mid‐2002, for years 2003–2018. SeaWiFS and MODIS‐Aqua chlorophyll‐a data sources are comparable (Zhang et al., 2006). All SeaWiFS and MODIS data were downloaded as monthly 4.6 km resolution Level 3 coverages from the NASA Ocean Color Web (https://oceancolor.gsfc.nasa.gov). We calculated thermal spring transition dates for the western and eastern Gulf of Maine according to the methods of Friedland et al. (2015). Spring thermal transition date is the day‐of‐year that regional ocean temperatures, which are smoothed with a five‐point moving average, surpasses the mean annual temperature for that location. Thermal transition dates were calculated using a 0.5° grid and NOAA's 0.25° daily Optimum Interpolation Sea Surface Temperature (or daily OISST; https://www.ncdc.noaa.gov/oisst).

Data processing

Aerial survey data were converted to detection/non‐detection data for our analysis. We imposed a 4.6 km spatial resolution to define a set of grid cells, or sites, within our study area (Figure 1; Figure S1). The spatial gridding we used is consistent with that used by many earth‐orbiting satellite data collection and dissemination programs (e.g., NOAA Pathfinder, NASA Ocean Color Web). For each aerial survey, we recorded non‐detections (0) at every grid cell visited by the aircraft and for which no animals were observed. We recorded a detection (1) at every grid cell visited by the aircraft in which one or more whales (of a given species) were observed. We followed this procedure for each of our target species: right, humpback, and fin whales. We also compiled data on Beaufort sea state and the ordinal day‐of‐year associated with each detection and non‐detection.

We binned aerial survey detection/non‐detection data into sampling “periods” which were one half month in duration (hereafter, “half‐month sampling period”) during which whale use of Cape Cod Bay was assumed to be stationary (i.e., “closed” to immigration/emigration). To do this, we divided each month into two sampling periods representing the first and second half of each month. Our full annual span of study (i.e., covering all half‐month sampling periods) lasted from 1 January to 15 May for right whales (9 half‐month sampling periods per year), and from 16 February to 15 May for humpback and fin whales (6 half‐month sampling periods per year). The species‐specific beginning and end dates were chosen to capture the typical annual temporal dynamics of primary habitat use by each species. Annual surveys and their start and end dates were designed to monitor right whales, and this dictated the seasonal boundaries we selected for the right whale model. Although surveys were conducted during December in some years, there were not enough surveys to justify inclusion of December in our study. Surveys usually ended during the first half of May, after the survey team confirmed that right whales were no longer present in Cape Cod Bay (Brown et al., 2007; Stamieszkin et al., 2010). We selected a shorter modeling time span for fin and humpback whales after exploratory data analyses revealed sporadic presence/absence of these species within Cape Cod Bay during January, and more consistent buildup of their presence during February.

Modelling approach

We tested for temporal trends in whale phenology in the Cape Cod Bay region by building a hierarchical model of intra‐ and inter‐annual whale usage of our study area for each of our three target species. We followed the general approach of a “multi‐season occupancy model” (MacKenzie et al., 2006), which allowed us to study changes in occurrence at sites over time while accounting for imperfect detection (i.e., the probability that a whale may have occurred at a location during a specific sampling period but went undetected). In particular, our model estimated the local “colonization” and “persistence” of whales within 4.6 km grid cells and half‐month sampling periods, decomposing changes in occurrence status of whales in grid cells into whether whales entered an empty grid cell (colonization) or remained in an occupied grid cell (persistence). The decomposition of occurrence into colonization and persistence dynamics provides a more mechanistic approach to modeling spatial dynamics, which can lead to superior predictive ability (Rushing et al., 2019). Whale behaviors manifest differently as seen through colonization and persistence probabilities: high colonization with high persistence indicates whales entering the study system, low colonization with low persistence indicates whales leaving the study system, high colonization with low persistence indicates higher rates of movement within the study area, and low colonization with high persistence indicates lower rates of movement within the study area and whales not entering the study system. Additionally, the combination of colonization and persistence together manifest as overall occupancy (i.e., the proportion of occupied sites at a given time), and overall occupancy can be estimated through time, allowing the precise estimation of phenometrics, such as the peak in whale use of the study area in a given year. Our general modeling approach thus built annual models of colonization and persistence dynamics for each whale species within our study area, but hierarchically linked parameters within species across all years. We used this approach to improve model fit and because we did not expect species' mechanistic relationships to covariates to change over time. The specific model formulation is as follows.

For a given whale species, we modeled the observations of that species in grid cell j (1, …, 109), on replicate k (1, …, 7), of half‐month sampling period l (1, …, 9), of year t (1, …, 21) as a Bernoulli random variable, observed with a probability of , where represents the probability of detecting a whale in a grid cell during a survey given its occurrence and is a binary indicator representing the true occurrence status of a whale species in a cell during that half‐month sampling period. The probability of detection, , is modeled as a logit‐linear function following the general form, where is a vector of fitted parameters, including intercepts and slopes, and X is a design matrix including observed covariates. We modeled detection as a function of two a priori covariates: the ordinal day of the year and Beaufort sea state (0–6) as measured during aerial surveys.

Occurrence status for a whale species was modeled differently depending on whether it was the first half‐month sampling period (l = 1) of the year or any subsequent half‐month sampling period (l > 1). The true occurrence status of a grid cell, , was modeled as a Bernoulli random variable with occupancy probability . In the first half‐month sampling period, was logit‐linearly derived from a single year‐specific intercept . In subsequent half‐month sampling periods, occupancy probability was decomposed into dynamic processes, such that, where and are the probabilities of persistence and colonization, respectively, for grid cell j during half‐month sampling period l of year t. Persistence and colonization probabilities were modeled using identical parameterizations, as we did not a priori expect them to differ. Consequently, we provide a full description of the parameterization of persistence, as follows, but do not present the same parameterization for colonization for the sake of redundancy. Persistence probability was modeled as a logit‐linear function of the form,

where is a vector of fitted parameters, including an intercept and slope for three covariates (SST, chlorophyll‐a, and bathymetry), X is a design matrix including the covariate values, and is a term to quantify a phenological shift across years taking the form, where and are slopes, day is the ordinal day of the year, is an offset term, and is the year of the survey (1–21). This parameterization, derived from Socolar et al. (2017), specifies a quadratic relationship between persistence probability and day of year, but it allows this phenological relationship to shift horizontally over time via the offset parameter . Estimating phenological shifts as a temporal translation of a dynamic occurrence process allows inference from the full annual whale‐visiting season across Cape Cod Bay, thus avoiding common pitfalls in phenological studies such as the bias that comes from focusing on imperfectly estimated dates such as “firsts” (Miller‐Rushing et al., 2008), or confounding of phenology and spatial variation (de Keyzer et al., 2017).

All slope and intercept covariates for detection (), initial occupancy (b), persistence (), and colonization () were modeled as year specific but were drawn from Gaussian hyper‐distributions, akin to a “random slopes” model formulation (Bolker et al., 2009). This hierarchical structure allowed the model to share information across years on how covariates impacted different parameters (i.e., the mean effect), but unlike the formulation for phenological shifts, year‐to‐year variation in covariate effects and intercepts (i.e., inter‐annual variation) was assumed to be random and uncorrelated.

A key feature of the phenological structure of the model is that it assumes a linear phenological shift over the 21‐year time frame. While rates of phenological shifts may not be constant over time, we were most interested in estimating a trend over the two‐decade time frame, and a linear form was the simplest assumption.

The model, like most occupancy‐type models, assumes that the system is “closed” within sampling periods (MacKenzie et al., 2006); that is, that the true occurrence status of a cell does not change over the three or so aerial surveys within each sampling period. Violations of this assumption can bias parameters (Rota et al., 2009), specifically through the deflation of detection probabilities. We sought to avoid closure violations by specifically defining sampling periods to as short a sampling period as possible while still maintaining approximately three surveys within each sampling period, as is needed to robustly estimate detection probability (Mackenzie & Royle, 2005).

We fit the model to the data with JAGS (Plummer, 2003) using the R statistical programming language version 3.6.3 (R Core Team, 2020) and the package “R2jags” (Su & Yajima, 2021). We used vague priors (i.e., normal with  = 0,  = 0.1). We ran nine chains of 240,000 iterations thinned by 200 with a burn‐in of 120,000, yielding a posterior sample of 4,500 across all chains. We checked convergence visually with traceplots and we ensured that all fixed effects and hyper‐parameters had a Gelman–Rubin statistic < 1.1 (Gelman et al., 2004). Inference on parameters was made using 95% Bayesian credible intervals (95 CrI). Distributions of parameter estimates that did not overlap 0 at the 95 CrI were considered to be extremely likely to have an effect, while distributions of estimates that did overlap 0 at the 95 CrI were not considered further. Posterior predictive checks of model fit were conducted by calculating Bayesian p‐values (Gelman et al., 1996) for two test statistics representing the average number of occupied half‐month sampling periods for each site and year, and the standard deviation of occupied half‐month sampling periods across all sites and years. Our models converged and Bayesian p‐values from posterior predictive checks were greater than 0.05 for each species, indicating no lack of model fit (Figures S2–S4). Full JAGS code and data for our model are provided in our online data archive.

Indicators of phenological change and their responses

We derived quantities to directly estimate trends in the magnitude and direction of phenological shifts of whale occurrence in our study area. We determined the posterior estimate of the true occurrence status for grid cells in each half‐month sampling period, , and from this we determined the day of each year on which the maximum number of occupied cells occurred. This day of maximum occupancy, ψ max, is an estimate of a direct phenological quantity of interest, specifically the day of year on which whale use of Cape Cod Bay peaked. We determined the exact day of year as follows: when there were two or more surveys in a sampling period, we used the median day of year between the first and last survey, when there was only one survey in a sampling period, we used that date, and when there were no surveys within a sampling period we used the median between the first and last day of the sampling period. A trend in days of year, within half‐month sampling periods, could bias our results. Thus, we tested for linear trends in the number of days between the start of each half‐month sampling period and the day of year (computed as described above); we found no significant trends (α = 0.01). Second, we calculated the mean proportion of occupied cells within each half‐month sampling period and year, . This provided a relative heuristic on the mean occupancy within each half‐month sampling period by which to understand how half‐month sampling period‐specific occupancy has changed over the 21 years of our study. To examine general within‐year habitat‐use dynamics, we plotted the average proportion of occupied sites as a function of day‐of‐year, for each year, for each species. Finally, to place phenological changes into the context of regional climate change that has occurred across our study area, we compared maximum occupancy () to the thermal spring transition date in the western and eastern Gulf of Maine. Propagating posterior uncertainty using inverse‐variance weighting of and for each year, we then estimated linear trends in and over the 21 years of our study, and in and the thermal spring transition date in the western and eastern Gulf of Maine. As we used a frequentist framework to estimate the relationships between the thermal spring transition date and maximum and seasonal mean occupancy, we describe the results from these models using frequentist terminology (i.e., statistical significance).

RESULTS

Aerial surveys

After spatially binning the study area, there were 109 sites. For the years 1998–2018, between 1 January and 15 May (the annual time frame of study for right whales), the aerial survey team conducted 576 surveys, with a mean of 3.0 (sd of 1.5) surveys per half‐month sampling period. The mean number of detections and non‐detections per year (maximum of one per site per survey) for right whales was 135 (sd of 63) and 1869 (sd of 574), respectively. Between 16 February and 15 May (the annual time frame of study for humpback and fin whales), the aerial survey team conducted 432 surveys with a mean of 3.4 (sd of 1.4) surveys per half‐month sampling period. The mean number of detections and non‐detections per year (maximum of one per site per survey) for humpback whales was 29 (sd of 19) and 2,225 (sd of 688), and for fin whales it was 77 (sd of 32) and 2,180 (sd of 677), respectively.

General temporal dynamics of habitat use

Modeled habitat use by right whales (Figure 2a) was consistent with historically known patterns in Cape Cod Bay (Ganley et al., 2019; Mayo et al., 2018). Occupancy was low in January, increased in February, peaked in March and April, and declined in May from its peak every year except 2018 when it increased—the final year included in this study. Humpback whale occupancy (Figure 2b), considerably lower than that of right whales, was usually near zero during the first half‐month sampling period we modeled (1–15 February) and peaked during April of every year—except in 2004, 2008, and 2018, when peak occupancies occurred in May. Fin whale occupancy (Figure 2c) was generally higher than that of humpback whales, but lower than that of right whales. Within‐year occupancy dynamics show a unimodal pattern similar to that for right and humpback whales. Across all years, fin whale occupancy was low in the second half of February and increased to a peak during April of every year.

FIGURE 2

Modeled changes in occupancy within each year, 1998–2018 for right whales (a), humpback whales (b), and fin whales (c). Occupancy was estimated for each half‐month sampling period; only posterior means are plotted for clarity. Earlier years in the time series are represented by cooler colors and later years are represented by warmer colors and. Whale illustration credit: NOAA Fisheries

Interannual period‐specific trends in habitat use

Across the 21 years of our study, period‐mean occupancy, , of right whales increased significantly during 16–29 February (p = .033) and in the four subsequent half‐month sampling periods: 1–15 March (p = .001), 16–31 March (p = .001), and 1–15 April and 16–30 April (p < .001; Figure 3d–h). Period‐mean occupancy during 1–15 May was highly variable, having near‐zero values in several years, and being above 0.5 in other years (Figure 3i). To better understand the trend in occupancy in years, when there was an appreciable presence of right whales in the 1–15 May half‐month sampling period, we tested for a linear trend after removing years in which was less than 0.005 (i.e., 2000, 2009, 2010, 2014, and 2016). After doing this, there was a significant increasing linear trend in period‐mean occupancy for the 1–15 May half‐month sampling period (p = .008). There were no significant trends in right whale period‐mean occupancy during the first three half‐month sampling periods, 1 January–15 February. There were significant negative trends in period‐mean occupancy of humpback whales during the 16–29 February (p = .004) and 16–31 March (p = .004) half‐month sampling periods, but there were no trends in the 1–15 March, or in the last three, half‐month sampling periods, 1 April–15 May, when overall humpback whale occupancy was generally higher than it was during the first three half‐month sampling periods (Figure S5). There was a significant positive trend in period‐mean occupancy for fin whales in the 1–15 March half‐month sampling period (p = .005), but no significant trends were seen in any other half‐month sampling periods (Figure S6).

FIGURE 3

Mean occupancy of right whales in Cape Cod Bay from 1998 to 2018 within sampling periods: 1–15 January (a), 16–31 January (b), 1–15 February (c), 16–29 February (d), 1–15 March (e), 16–31 March (f), 1–15 April (g), 16–30 April (h), and 1–15 May (i). Vertical lines show 50% CrI. Solid and dashed fitted lines are inversely weighted by the variance in each estimate of mean occupancy. Dashed fitted line in (h) indicates a significant trend in mean occupancy as a function of year, for all years in which mean occupancy for this sampling period was ≥0.005. Mean occupancy < 0.005 is indicated by hollow circles. Shaded regions represent 95% CIs. Bold p‐values indicate statistical significance (p < .05)

Phenological shifts in habitat use

Significant shifts in the day of maximum occupancy, , were detected over the 21 years of our study for two out of the three whale species (Figure 4). The day of maximum occupancy for right whales shifted later in time (p < .001) with an estimated phenological change of +18.1 days (+0.90 days/year), or approximately 1.2 half‐month sampling periods. The day of maximum occupancy for humpback whales shifted later in time (p = .009) with an estimated phenological change of +19.1 days (+0.96 days/year), or approximately 1.2 half‐month sampling periods. The day of maximum occupancy for fin whales suggested by the linear trend was −5.8 d earlier in time, but this trend was not statistically significant (p = .219).

FIGURE 4

Day of maximum occupancy for right whales (a), humpback whales (b), and fin whales (c) in Cape Cod Bay from 1988 to 2018. Vertical lines show 50% CrIs. Fitted lines are inversely weighted by the variance in each estimate of mean occupancy. Shaded regions represent 95% CIs. Bold p‐values indicate statistical significance (p < .05)

Relationship to thermal spring transition date

There was a significant negative effect of the western Gulf of Maine spring transition date on the right whale date of maximum occupancy (p < .001; Figure 5a), and a significant positive effect of the eastern Gulf of Maine spring transition date on the fin whale date of maximum occupancy (p = .044; Figure 5f). Humpback whales showed decreasing trends with spring transition date in both the western and eastern Gulf of Maine (Figure 5c,d), but these relationships were not statistically significant (p = .217 and p = .447, respectively).

FIGURE 5

Relationship between the day of maximum occupancy (day‐of‐year) for right (a, b), humpback (c, d), and fin (e, f) whales with the spring thermal transition date (day‐of‐year) for the western (a, c, e) and eastern (b, d, f) Gulf of Maine. Error bars show the 50% CrIs for dates of maximum occupancy. Fitted lines are linear model fits inversely weighted by the variance in each estimate of mean occupancy. Shaded regions represent 95% CIs. Bold p‐values indicate statistical significance (p < .05). Note that y‐axes differ between species

Detection probabilities

Large whales spend a substantial amount of time below the surface, and during this time they are normally unavailable to be visually detected by survey crew observers. Therefore, our estimates of detection probability (Table 1) combine observer bias and availability bias. The average probability of detection per aerial survey for right whales was 0.142 (95 CrI: 0.105–0.181), for humpback whales was 0.056 (95 CrI: 0.029–0.092) and for fin whales was 0.093 (95 CrI: 0.058–0.145). There was an extremely likely negative effect of Beaufort sea state on detection probability for right whales such that right whales were less likely to be detected when sea state was high. There was an extremely likely positive effect of day‐of‐year on the detection probability for all three species such that all species were more likely to be seen as the year progressed.

TABLE 1

Posterior estimates and their 95% credible intervals for models of right, humpback, and fin whales in Cape Cod Bay from 1998 to 2018. Estimates whose 95% credible interval did not cross zero are shown in bold as they are considered to have strong support and thus extremely likely to have had an effect

Parameter name	Parameter interpretation	Right whales	Humpback whales	Fin whales
		Mean (2.5% 97.5%)	Mean (2.5% 97.5%)	Mean (2.5% 97.5%)
Phenological shift
∆_g	Phenological offset term for colonization	0.005 (−0.026 0.033)	−0.046 (−0.097 0.005)	0.019 (−0.118 0.112)
∆_e	Phenological offset term for persistence	−0.025 (−0.051 0.000)	−0.005 (−0.061 0.054)	−0.023 (−0.060 0.012)
Detection
mu.a.0	Intercept	−1.798 (−2.146 −1.442)	−2.818 (−3.527 −2.195)	−2.282 (−2.785 −1.776)
mu.a.bft	Beaufort sea state	−0.170 (−0.268 −0.072)	−0.151 (−0.366 0.058)	−0.092 (−0.220 0.032)
mu.a.jday	Day‐of‐year	0.333 (0.044 0.630)	0.726 (0.108 1.333)	0.521 (0.092 0.903)
Initial occupancy
mu.b.0	Intercept	−1.326 (−2.262 −0.693)	−4.445 (−6.409 −2.989)	−2.520 (−3.711 −1.473)
Colonization
mu.g.0	Intercept	−2.911 (−4.487 −1.623)	−4.207 (−6.160 −2.914)	−4.237 (−7.423 −1.788)
mu.g.bat	Bathymetry	−1.661 (−2.503 −1.021)	−2.004 (−2.834 −1.400)	−1.918 (−3.257 −1.059)
mu.g.chl	Chlorophyll‐a	1.882 (0.693 3.315)	0.092 (−0.807 0.809)	0.423 (−0.126 1.122)
mu.g.sst	Sea surface temperature	−2.222 (−3.876 −0.778)	−1.561 (−4.048 0.619)	−3.149 (−6.208 −0.317)
mu.g.jday	Day‐of‐year	−0.639 (−3.107 1.736)	0.207 (−2.805 2.994)	1.636 (−0.816 3.857)
mu.g.jday2	Day‐of‐year2	−2.426 (−5.103 −0.053)	−2.261 (−3.839 −1.127)	−1.019 (−2.052 0.141)
Persistence
mu.e.0	Intercept	8.326 (5.749 11.221)	3.410 (−0.958 7.665)	4.566 (2.401 7.483)
mu.e.bat	Bathymetry	1.234 (0.032 2.460)	−2.062 (−5.202 −0.087)	0.850 (−0.701 2.310)
mu.e.chl	Chlorophyll‐a	−0.113 (−0.732 1.078)	0.010 (−3.728 4.266)	−0.540 (−2.068 0.916)
mu.e.sst	Sea surface temperature	0.305 (−2.143 2.218)	0.230 (−3.340 3.763)	0.244 (−1.687 1.883)
mu.e.jday	Day‐of‐year	−2.455 (−6.248 1.085)	−1.138 (−6.382 4.218)	−0.071 (−3.874 3.597)
mu.e.jday2	Day‐of‐year2	−5.781 (−7.735 −3.948)	−1.255 (−5.252 3.750)	−3.130 (−5.355 −1.105)

Colonization and persistence

For all species, bathymetry had an extremely likely negative effect on colonization (Table 1). Chlorophyll‐a, a proxy for primary productivity, had an extremely likely positive effect on colonization for right whales. SST had an extremely likely negative effect on colonization for right and fin whales. Bathymetry had an extremely likely positive effect on the persistence of right whales, but had an extremely likely negative effect on the persistence of humpback whales.

DISCUSSION

Our study suggests that highly migratory marine mammals can and do adapt the timing of their habitat use in response to climate‐driven changes in their environment, and this emphasizes the need for improved understanding of habitat‐specific linkages between physical drivers and lower trophic level food‐web dynamics. Using a novel multi‐season occupancy model, we found strong evidence for later timing of peak seasonal habitat use for right and humpback whales in Cape Cod Bay but no detectable change in timing for fin whales. Understanding the magnitude and direction of phenological changes for highly migratory species provides insight into how climate change affects seasonal habitat use, and provides crucial information to assist managers and decision‐makers tasked with balancing protective measures for species and their habitats with diverse stakeholder needs (e.g., shipping and fishing industries). Of particular concern are increases in occupancy during May (the end of the seasonal foraging time in Cape Cod Bay) in comparison to stable phenology (no detectable changes) at the beginning of the seasonal foraging time (January–February). Regional thermal indicators of the onset of spring show that right whales are sensitive to long‐term interannual variability in temperature conditions in the western Gulf of Maine, while fin whales appear to be more sensitive to conditions in the eastern Gulf of Maine; the relationship for humpback whales was inconclusive. Overall, we found increased habitat use by right whales in Cape Cod Bay during winter and spring and relative to thermal indicators. Given recent trends of decreased use of other seasonal habitats in the Gulf of Maine (Davies et al., 2019; Record et al., 2019), our results emphasize the need for regional protections to integrate and keep pace with new information on climate‐induced shifts in whale habitat use in both space and time.

Phenological shifts

Only a few studies have measured phenological changes for large migratory whales in the western Atlantic Ocean. Ramp et al. (2015) documented earlier shifts in the timing of arrival and departure of fin and humpback whales in the Gulf of St. Lawrence over a 27‐year time span (1984–2010). Specifically, they estimated an arrival date shift of >1 day per year earlier for fin and humpback whales, and a departure date shift of 1 day per year earlier for humpback whales and 0.4 day per year earlier for fin whales. Our study used the day of maximum occupancy to represent the time of year when peak habitat use occurred in Cape Cod Bay by each species. We found later timing of maximum habitat use for right and humpback whales of 0.9 and 0.96 days per year, respectively. Although modeling approaches differed, our estimated rates of change were similar in magnitude to the shifts in arrival and departure dates estimated by Ramp et al. (2015), but opposite in direction. In contrast, Charif et al. (2020) analyzed the acoustic presence of right whales in Massachusetts Bay, including part of Stellwagen Bank which adjoins Cape Cod Bay on its northern boundary. Their results were limited to 6 years of data and were not statistically significant, but they suggested that peak presence for right whales had shifted earlier by approximately 40 days (five 8‐day periods) from 2008 to 2011. However, this pattern reversed in the final 2 years of their study, from 2012 to 2013, when the date of peak presence shifted later in time. Detections of presence/absence from acoustic and visual observations are problematic (Clark et al., 2010), making estimation of an annual rate and comparisons with Charif et al. (2020) difficult to fully ascertain.

Physical drivers have been shown to affect timing of habitat use by marine mammals in other regions, particularly high‐latitude environments. For example, the autumn migration timing of the Eastern Chukchi Sea population of beluga whales (Delphinapterus leucas) has shifted later by 2 to > 4 weeks over 15 years (1998–2012; Hauser et al., 2017). These shifts were strongly related to the onset of regional sea‐ice freeze‐up, which limits whale access to northern habitats in the Chukchi and Beaufort seas. In the same region, the autumn migration of Western Arctic bowhead whales (Balaena mysticetus) was found to occur 7 days per year later from 2008 to 2017, though the physical or biological reasons for this were unclear (Stafford et al., 2021).

Determination of prey‐driven habitat use is not always straightforward because biological factors such as zooplankton or finfish (whale prey) abundance can be difficult to measure. As a result, many studies resort to using chlorophyll‐a as a proxy for primary production. Visser et al. (2011) found that fin and humpback whales track secondary production generated by the North Atlantic spring bloom, but that their abundance was more strongly associated with the timing of the onset of the phytoplankton spring bloom than with later stages of bloom development. Szesciorka et al. (2020) found that blue whales (Balaenoptera musculus) arrived (as determined by acoustic detections) in the Southern California region 42 days earlier over a 10‐year time span, for a phenological shift of 4.2 days per year. Blue whale arrival was earlier when the SST anomaly from the previous year was colder, and later when the SST anomaly from the previous year was warmer. Furthermore, colder years had lower krill biomass and warmer years had higher krill biomass. Thus, day‐of‐year trends were related to a combination of physical and biological cues, but the phenological shifts in habitat use manifested themselves in a context (day‐of‐year) that was meaningful for resource managers who usually develop management plans based upon a calendar rather than biological cues.

Right whales are specialist planktivorous feeders that meet most of their energetic needs by consuming late‐stage C. finmarchicus; however in Cape Cod Bay, right whales target smaller copepods that occur in greater numbers, primarily Pseudocalanus spp. and Centropages typicus (Mayo & Marx, 1990; Pendleton et al., 2009). Humpback and fin whales are generalist piscivorous feeders that target small fishes (primarily sand lance) in the southwestern Gulf of Maine (Overholtz & Nicolas, 1979; Payne et al., 1986, 1990; Staudinger et al., 2020). Our results, similar to Ganley et al. (In Review), showed a significant relationship between right whales and the western Gulf of Maine spring transition date, but not for fin or humpback whales and the western Gulf of Maine spring transition date. This result supports our hypothesis that baleen whales which feed at a relatively lower trophic level are more tightly coupled to the physical environment than higher trophic level feeders. Surprisingly, the right whale maximum occupancy date moved later while the spring transition date moved earlier. This is the opposite of our expectation based on global trends in marine phenology (Cohen et al., 2018; Friedland et al., 2018; Poloczanska et al., 2013) and regional trends in bloom timing of spring primary productivity (Friedland et al., 2015).

A possible explanation for the inverse relationship between the right whale maximum occupancy date and the spring transition date is that climate‐induced reductions in prey (late‐stage C. finmarchicus) in other right whale habitats within the Gulf of Maine have caused Cape Cod Bay to become a “waiting room” for right whales. In this scenario, right whales wait in a habitat with adequate prey (Pseudocalanus spp. and C. typicus) while richer C. finmarchicus resources develop in northerly habitats such as the Gulf of St. Lawrence. Indeed, decades of observations of right whales (Brown et al., 2007; CETAP, 1982; Watkins & Schevill, 1982) have revealed that after departing Cape Cod Bay, right whales typically moved to the Great South Channel to feed on C. finmarchicus (Beardsley et al., 1996; Kenney et al., 1995; Wishner et al., 1988) before moving to other deep‐water Gulf of Maine habitats such as the Bay of Fundy and Roseway Basin in search of late‐stage C. finmarchicus. In recent years, right whales have not followed this migratory pattern (Davies et al., 2019; Record et al., 2019), lending support to the hypothesis that right whales use Cape Cod Bay as a waiting room while richer prey resources develop elsewhere.

Little is known about fin whale migration in the Gulf of Maine. The Western North Atlantic population studied in this paper ranges from the United States mid‐Atlantic waters to Newfoundland (NOAA, 2020). The positive relationship we found between fin whale day of maximum occupancy and the eastern Gulf of Maine spring transition date suggests that fin whales are tracking changing production dynamics in the eastern Gulf of Maine. More research is needed on fin whale distribution and diet in the western North Atlantic Ocean to better understand their phenological patterns and responses to climate change.

Colonization and persistence relationships

Bathymetry had an extremely likely negative effect on colonization for all three species (Table 1), indicating whales are more likely to occupy deeper areas, but results differed for persistence. For right whales, bathymetry had an extremely likely positive effect on persistence, indicating that, once occupied, shallower sites were likely to remain occupied. This finding is consistent with repeated observations of right whales along a shallow isobath, a southwesterly arc from the southwest to northeast quadrants of Cape Cod Bay (Mayo et al., 2018; Nichols et al., 2008). The general circulation pattern in Cape Cod Bay is counter clockwise, with the Western Maine Coastal Current circulation flowing in a southerly direction and penetrating into Cape Cod Bay before the flow pattern moves water around the eastern arm of Cape Cod (Bigelow, 1927; Geyer et al., 1992; Jiang et al., 2007). These circulation patterns may concentrate zooplankton along Cape Cod Bay's shallow eastern isobaths. For humpbacks, bathymetry had an extremely likely negative effect on persistence, further reinforcing our finding that humpback whales are more likely to occupy relatively deeper waters, which are located centrally and toward the northern boundary of our study area, off of Stellwagen Bank which is known habitat for their primary prey, sand lance (Silva et al., 2021).

Chlorophyll‐a is a proxy for primary productivity, and the extremely likely positive relationship between right whale colonization and chlorophyll‐a suggests a close linkage between primary and secondary production in Cape Cod Bay. The relationship between primary production and right whale colonization suggests that right whales use Cape Cod Bay to feed, and aligns with field observations and previous studies that demonstrate the importance of, and strong linkage between, right whales and their prey in Cape Cod Bay (Mayo & Marx, 1990; Pendleton et al., 2009). Humpback and fin whales feed at higher trophic levels, so their relationship with chlorophyll‐a may be temporally lagged or decoupled with colonization or persistence depending on how well synchronized their prey are with primary and secondary producers.

Sea surface temperature had an extremely likely negative effect on right and fin whale colonization, indicating that these species are more likely to use areas with lower temperatures and/or at times of year with relatively low temperatures. Given the rapid rates of seasonal and annual warming observed and projected in the region (Pershing et al., 2021; Saba et al., 2016; Thomas et al., 2017), these species may reduce their use of the Cape Cod Bay habitat in the future if preferred thermal conditions are not maintained.

Cohen et al. (2018) suggest that phenological delays in timing may arise due to insufficient time series and sampling error. We acknowledge that our dataset is limited to only two decades, and detecting whale occurrence is inherently problematic due to the ability of whales to dive and intermittent aerial survey coverage (e.g., Ganley et al., 2019). Our modeling framework accounted for imperfect detection to address the issue of false non‐detections resulting from subsurface whales that are present but not visible. Since our ultimate inference is based on temporal trends in phenology, and there is no indication that closure violations would have increased or decreased over time, our inference is robust to any violations that may still exist in our period‐aggregated data. We also note that estimates of occupancy do not represent density. Rather, occupancy is a spatial measure of habitat use, with greater occupancy indicating that a larger area of the habitat is occupied. Due to a lack of availability of fine‐scale data on key zooplankton species in Cape Cod Bay, we did not directly analyze zooplankton productivity or any potential prey shifts, even though changing prey fields are likely an important determinant of the timing of whale presence in our study area (Mayo & Marx, 1990; Pendleton et al., 2009). A detailed understanding of sand lance and C. finmarchicus dynamics, and the importance of each to right whales versus fin and humpback whales in both Cape Cod Bay and Stellwagen Bank (e.g., Payne et al., 1990) would greatly improve our understanding of when to expect large whales in and around Cape Cod Bay.

Management implications

Globally, fishing gear entanglements and vessel strikes are the major source of human‐caused mortality of marine mammals (Read, 2008; Thomas et al., 2016). In United States waters and in Cape Cod Bay, protections for large whales have evolved over the years. In 2008, to protect right and other large whales from vessel strikes, federal speed limits for large vessels (>19.8 m) in and around Cape Cod Bay were implemented between 1 January and 15 May (NOAA, 2008b), and in 2013 they were extended in perpetuity (NOAA, 2013; Laist et al., 2014). In 2018 (the last year of our study), the Commonwealth of Massachusetts implemented speed limits for small vessels (<19.8 m) within Cape Cod Bay. A changing patchwork of federal and state fixed fishing gear restrictions intended to protect right whales from entanglement have been in place since the late 1990s. In 2015, a seasonal trap/pot gear closure was implemented from 1 February–30 April within state and federal waters in and around Cape Cod Bay. In 2021, state protections on fishing and shipping were expanded to 15 May (MADMF, 2021). Recent state protections have been written to allow dynamic extension of protections if whales remain in the area. Our study shows that right whale occupancy has increased significantly from 16 February to 30 April, and that 1–15 May often has very high right whale occupancy but also high inter‐annual variability. Thus, our results support the extension of fishing and vessel speed restrictions to 15 May, and later extension of protections may be warranted.

Our understanding of the Cape Cod Bay habitat has benefited greatly from over two decades of systematic, high resolution, monitoring. The resultant dataset has allowed resource managers to adapt regulations in response to observations of whales. The patterns revealed in our study provide definitive evidence of long‐term phenological changes by large whales in Cape Cod Bay. Results can be used to support nimble resource management plans that adapt to climate‐driven ecological changes in the greater Gulf of Maine and other rapidly warming marine ecosystems. For example, projections of potential future phenological shifts could be used in management scenario plans which attempt to anticipate a range of outcomes and design corresponding management actions (Borggaard et al., 2020). Additionally, Climate Vulnerability Assessments (Hare et al., 2016; Lettrich et al., 2019) could also incorporate results from phenological studies like ours to help managers anticipate future risks and identify conservation priorities.

A major conservation challenge for human society and regulatory agencies charged with protecting wildlife populations is to adapt rules in a flexible and timely fashion in response to environmental and biological changes. In addition, it is important to set predictable expectations regarding when and where potentially harmful human activities can take place. Our oceans are facing increasing pressure from a combination of commercial uses, recreation, and renewable energy development. Continued long‐term monitoring data and quantitative tools to measure and predict phenological changes, such as those presented here, are critical to understanding how wildlife populations are responding to climate change. Our results provide needed information to resource managers who are challenged with balancing demands from multiple stakeholders, with the ultimate goal of advancing sustainable use of natural resources.

CONFLICT OF INTEREST

The authors declare no conflicts of interest.

Appendix S1 Supporting information

Click here for additional data file.

Click here for additional data file.