Impacts of marine heatwaves on top predator distributions are variable but predictable

Domain: earth_science
Published: 2023-09-05
Editorial outcome: minor_revision

Abstract
Marine heatwaves cause widespread environmental, biological, and socio-economic impacts, placing them at the forefront of 21st-century management challenges. However, heatwaves vary in intensity and evolution, and a paucity of information on how this variability impacts marine species limits our ability to proactively manage for these extreme events. Here, we model the effects of four recent heatwaves (2014, 2015, 2019, 2020) in the Northeastern Pacific on the distributions of 14 top predator species of ecological, cultural, and commercial importance. Predicted responses were highly variable across species and heatwaves, ranging from near total loss of habitat to a two-fold increase. Heatwaves rapidly altered political bio-geographies, with up to 10% of predicted habitat across all species shifting jurisdictions during individual heatwaves. The variability in predicted responses across species and heatwaves portends the need for novel management solutions that can rapidly respond to extreme climate events. As proof-of-concept, we developed an operational dynamic ocean management tool that predicts predator distributions and responses to extreme conditions in near real-time.

Introduction
Long-term climate trends (e.g., global warming) and short-term extreme events (e.g., heatwaves) have global impacts on ecosystem structure and functioning, and human well-being1–3. The impacts of long-term climate trends have received considerable attention through the examination of the warming signal in both historical observations and future climate projections4–9. However, mounting evidence indicates that episodic events like fires, floods, and heatwaves can have catastrophic ecosystem and socio-economic impacts1,10–12.

In particular, heatwaves on land13 and at sea1 adversely impact an additional 157 million people today compared to the turn of the century14. The largest and warmest marine heatwaves (MHWs) on record have occurred in the last decade and include the 2012 Northwest Atlantic15, the 2015–16 Tasman Sea16, and 2013–16 Northeast Pacific17 events, exerting considerable physical, ecological, socio-economic, and human health impacts1.

Faced with a changing climate, mobile species’ first responses are often to shift their geographic ranges to remain within suitable environmental conditions2. This phenomenon is well-described within the context of long-term warming4–6,8 and El Niño events18,19, yet there is a paucity of information on species redistribution in response to unprecedented warming during recent MHWs20. Temperatures observed during MHWs can be similar to projected mean future conditions at the end of the 21st century21, providing valuable insight into species redistribution in the coming decades. Previous investigations on MHW-driven species range shifts have largely relied on observational data from opportunistic sightings, surveys, or tagging programs, and have demonstrated poleward or vertical distributional shifts toward cooler temperatures22–26. However, these data are often patchy, offering only snapshots of impacts from a single MHW on a single species22,23,27 or on several species24,25,28. Comparisons of multi-species responses across multiple MHWs are rare despite ample evidence that there is high variability in the evolution, drivers, and physical characteristics of MHW events21,29,30. Statistical models provide a means of interpolating across space, time, and taxa, providing information on MHW-driven redistribution by offering inferences on unobserved locations, MHWs, and individuals. Furthermore, statistical models can relate species distributions to multiple environmental drivers, thereby accounting for the complex physical and biogeochemical changes beyond the increased temperature that occurs during MHWs29.

Here, we quantify the impacts of four major North Pacific MHWs (2014, 2015, 2019, 2020; Supplementary Fig. 1) on the spatial distributions of 14 marine top predators, spanning several major guilds: seabirds, mammals, turtles, tunas, and sharks (Fig. 1). The impacts to most of these species during the MHWs were previously unexamined, and thus, we provide an unprecedented look at how responses vary across a group of iconic species and across MHWs.

The Northeast Pacific Ocean provides an ideal testbed for examining the effects of MHWs. This region has experienced some of the longest, largest, and most intense MHWs on historical record31,32 (Fig. 1A), with thermal displacements of >1000 km during the strongest events21. The Northeast Pacific is also a biodiversity hotspot and major foraging ground33, attracting a range of top predator species from across the broader Pacific Ocean, many of which are ecologically and commercially valuable and/or threatened with extinction33 (Fig. 1B). These predators are highly mobile and conduct large-scale movements throughout the Pacific basin34, and thus have the physical ability to actively redistribute in response to MHWs. Many of the most extreme temperature anomalies have occurred where predator density is highest (Fig. 1), providing an opportunity to examine MHW impacts across several higher trophic levels that are particularly sensitive to climate perturbations35. Furthermore, we examine how MHWs shift predators across political boundaries, revealing which exclusive economic zones (EEZs) lost and gained predator habitat and where adaptable governance strategies are likely to be important. By implementing a multi-species, multi-MHW analytical framework, our work provides the basis for informing future management of the impacts of these extreme climate events.

Results
We used boosted regression tree models fit extensive telemetry datasets33,36,37 to predict the redistribution of species’ preferred habitats during each MHW event. Our models accurately predicted distribution shifts during MHW years captured by an extensive (>one million records) independent top predator dataset collated from public, private, and government sources (Supplementary Table 2 & Supplementary Fig. 7). This independent dataset included records from tagging programs, shipboard surveys, opportunistic sightings, and fisheries observer programs, and allowed us to explore predictive performance in the telemetry-based models. For some species, independent datasets were more closely aligned with those used to build models due to a similar spatial domain and/or data type (e.g., telemetry data for blue whales), which resulted in higher predictive performance. However, for others, opportunistic sightings and shipboard surveys covered different spatial domains than the telemetry datasets and typically had lower predictive performance (e.g., Laysan albatross, Supplementary Fig. 8).

Predicted MHW impacts were highly variable within and among MHW events (Fig. 2). For example, responses among MHWs (rows, Fig. 2) were highly variable for Coastal species. While all Coastal species were displaced to the northwest during the 2014 and 2015 events, the 2019 and/or 2020 events drove southeastward displacement for bluefin and albacore tuna, and blue and mako sharks (Figs. 2B and 3C; Supplementary Fig. 10A, C, I). Southeastward displacement may be related to the emergent cool water refugia along the southern US and Mexican coasts during these two events (Fig. 1A). Sooty shearwaters and sea lions underwent a reversal in range shifts between the 2014–15 and 2019–20 events, experiencing range compression during the 2014 and 2015 events (22–64% decrease), and range expansion during the 2019 and 2020 events (335–377% increase) (Figs. 2C and 3A; Supplementary Fig. 10E). Similarly, elephant seals and sea lions lost habitat during the 2014 MHW (29% and 65% decrease, respectively), but gained habitat during the 2019 and 2020 MHWs (71–46% and 117–158% increase) (Fig. 3B; Supplementary Fig. 10E).

In addition to among-MHW differences, we found differences in predicted impacts within-MHWs (columns, Fig. 2). While most species were predicted to experience relatively low displacement distances during the 2019 and 2020 events, sooty shearwaters were displaced large distances (536 and 721 km, respectively) during these two MHWs. Displacement direction during the 2020 event ranged from north (blue shark), east (salmon shark), south (Laysan albatross), and west (blue whale) (Fig. 2B; Supplementary Fig. 10). The 2019 and 2020 events caused large range expansions for sooty shearwaters and California sea lions whereas most other species experienced range compression (Figs. 2C and 3A; Supplementary Fig. 10E). While most species lost habitat in 2019 and 2020, elephant seals, salmon sharks, California sea lions, and yellowfin tuna gained habitat (Fig. 2D).

Among- and within-MHW differences were also apparent in the predicted redistribution of predators across political boundaries (Fig. 4). The US EEZ was predicted to gain predator habitat during each MHW, with the largest increase during the 2015 event (10%), and the smallest increase during the 2020 event (2%; Fig. 4A). The Mexican EEZ was predicted to lose predator habitat during 2014, 2015, and 2020 events (with the largest decrease in 2015 of 8%). The Canadian EEZ and high seas were predicted to gain and lose predator habitat during each MHW, respectively. Coastal species experienced the largest jurisdictional redistributions during the 2014 and 2015 events, particularly yellowfin and albacore tunas, with 31% and 22% of their respective habitats predicted to shift into US waters (Fig. 4B). Southern species experienced the largest redistributions during the 2019 event, with 39% of white shark habitat and 8% of leatherback turtle habitat predicted to shift into the US EEZ. Northern species experienced the largest redistributions during the 2020 event, with 14% of salmon shark habitat predicted to shift into the Canadian EEZ.

Methods
Telemetry data used to build the species distribution models (SDMs) were acquired for 14 top predators tagged from 2000–2010, including data from the Tagging of Pacific Predators Project33 and private datasets. For all species except albatrosses, detailed methods on the number of individuals tagged, handling of tagging bias, and state space modeling are included in Block et al.33, Winship et al.36, and Jordan et al.37. Further details on albatross telemetry data are available in Supplementary methods section 1.3.

Species were assigned to one of three geographical groups—Northern, Coastal, and Southern—based on the locations of telemetry data collected from August–October across all years. Coastal species had over 60% of observations located within the California Current and Gulf of California Large Marine Ecosystems (blue [Prionace glauca] and mako [Isurus oxyrinchus] sharks; albacore [Thunnus alalunga], yellowfin [T. albacares], and bluefin [T. thynnus] tunas; blue whales [Balaenoptera musculus], sooty shearwaters [Ardenna grisea], and California sea lions [Zalophus californianus]), despite some of these species undertaking long migrations across the North Pacific. Northern species had observations to the northwest of Coastal species (black-footed [Phoebastria nigripes] and Laysan [P. immutabilis] albatross; elephant seals [Mirounga angustirostris] and salmon sharks [Lamna ditropis]); and Southern species had observations to the southwest (white sharks [Carcharodon carcharias] and leatherback turtles [Dermochelys coriacea]). These regional groupings were useful to describe broad patterns, as MHW impacts were most similar across species within the same geographical group. The 75th percentile kernels for each geographical group were calculated by randomly subsampling each species’ locational data to the same number of records, and then calculating kernel densities and the 75th percentile kernel across all species within each grouping (Fig. 1B). Kernel densities and 75th percentile kernels were also calculated for each species individually and used to quantify environmental conditions within each species’ 75th percentile kernel during the MHW events (Supplementary Fig. 3).

For each year, MHWs were assessed from August–October: the months in which the highest temperature anomalies were observed across the Pacific (Supplementary Fig. 1, Supplementary methods section 1.2). In addition to SST anomalies (Fig. 1A), anomalies of oxygen at 200m, mean primary productivity within the upper 200m, and surface chlorophyll-a conditions were also assessed during MHWs (Supplementary Table 1). Anomalies during each MHW were calculated relative to mean conditions across August–October 2000–2020. We were precluded from using a longer baseline by the availability of satellite observations for chlorophyll-a (a covariate in the SDMs), which came online as a science-quality product in 1998.

The boosted regression tree (BRT) models used the multi-year animal telemetry data to describe species-environment relationships from a suite of environmental variables (Supplementary Table 1, Supplementary Methods section 1.1). Dynamic variables used to fit the BRTs included primary productivity, oxygen, sea surface temperature (SST) and its spatial standard deviation, sea level anomaly, eddy kinetic energy, mixed layer depth, chlorophyll-a, and day of the year. All environmental variables were resampled from their native resolutions to 0.25 degrees to match the coarsest resolution of the environmental datasets (sea surface height products). All analyses were performed in R version 4.0.4.

Background pseudo-absences were generated at a 1:1 ratio of presences for the telemetry datasets58. For each presence point, a pseudo-absence was generated for the same date (Supplementary methods section 1.5.1). Presence and pseudo-absence data were matched to the environmental datasets in space and time. We used a BRT framework to model the probability of species presence as a function of the environment. BRTs are a common machine learning model, popularized by their ability to fit complex nonlinear relationships and their robustness to wide varieties of data types and distributions. For all 14 species, BRTs with a binomial distribution were used to model the probability of species presence as a function of the environment. BRTs were built with a bag fraction of 0.6, and a tree complexity of three, and a learning rate that varied between 0.0001 and 0.00001 to ensure at least 2,000 trees were fit for each model. In addition to predicting suitable habitats, the relative importance of variables in each SDM was identified (Supplementary Fig. 5).

Independent datasets were used to validate the temporal extrapolation of telemetry-derived models beyond the time-series of the training data (i.e., post-2010), particularly during MHW years (Supplementary Table 2, Supplementary methods sections 1.4. and 1.5.4). Datasets were acquired from a diverse range of public and government sources including fisheries observer programs, animal sightings from citizen science databases, and tagging data from dedicated programs. For each species, all data sources were pooled for validation, e.g., the blue shark model was validated on an independent dataset composed of both California drift gillnet observer data and survey data from the North Pacific Pelagic Seabird Database (Supplementary Table 2). In addition to validation on independent data, several methods of cross-validation were performed (Supplementary Table 3, Supplementary methods section 1.5.3).

The BRT models were predicted over the daily environmental data from 2000–2020, and spatially constrained within a minimum convex hull of the training data. The resultant habitat suitability predictions ranged from 0–1, with low and high values indicating unsuitable and suitable habitats, respectively. These continuous predictions were averaged to create an August–October 2000–2020 mean, and August–October anomalies during each MHW year were calculated relative to the long-term mean (Supplementary Fig. 9).

The continuous habitat suitability predictions were reclassified into daily binary core habitat using species-specific thresholds (see Supplementary methods section 1.6.2 for a comparison of binary versus continuous core habitat). Core habitat was defined as pixels with habitat suitability greater than or equal to the top 50% of predicted values at true presences (i.e., models were predicted back on the telemetry data used in model fitting, and 50% prediction quantile was used as the threshold). The 50% threshold was selected over the more conservative thresholds used in the climate projection literature (e.g., 25%4) as these conservative thresholds sometimes resulted in a complete loss of core habitat during extreme events (e.g., 2015).

We quantified MHW impacts on species’ core habitats using four metrics (Fig. 2). Displacement direction and distance captured how the center of gravity of core habitat changed during MHWs, and was measured as a vector of both cardinal direction (degrees) and distance (kilometers). Range extent captured the distance between the leading and trailing edges of the core habitat in the north-south and east-west directions, providing a metric of range compression or expansion due to MHWs (negative and positive percent change range, respectively). Core habitat area captured the total available core habitat and its loss or gain due to MHWs (negative and positive percent change area, respectively). These metrics are independent of one another; for example, during a MHW, a species can gain new habitat beyond its typical range (i.e., range expansion), while overall having less available habitat throughout its range (i.e., habitat loss).

Three metrics were calculated from the daily binary core habitat rasters: the center of gravity of core habitat (mean latitude/longitude coordinate pairs), the interquartile range of core habitat (latitude/longitude coordinate pairs for the north-south and east-west interquartile ranges), and the total amount of core habitat (km2). To quantify species’ range extents, the daily north-south interquartile range was multiplied by the daily east-west interquartile range (both measured in decimal degrees). The interquartile range was selected over the absolute range to remove the effect of outliers (e.g., individual pixels of core habitat with extreme distributions), and to ensure that expansion of range was not artificially truncated due to proximity to the convex hull boundary.

The three metrics were averaged across each MHW event (August–October of 2014, 2015, 2019, 2020), and deviations from the 2000–2020 August–October mean was calculated for each habitat metric and each species, producing estimates of habitat displacement, changes to the range extent, and changes to the habitat area. Habitat displacement was assessed as a vector (i.e., cardinal direction and distance) joining the long-term mean position of the core habitat center of gravity to its MHW-specific location. Range extent change reflected the MHW-specific compression or expansion of range extent compared to the mean range, measured by area, while habitat area change assessed the MHW-specific loss or gain of total core habitat, measured by area. In addition to species-specific habitat metrics (Fig. 2), the mean and standard deviation of each metric was calculated across the regional species groups (Supplementary Table 4). To capture model sensitivity to input data, 20 replicate BRTs were fit to 20 different 75% subsamples of the telemetry and pseudo-absence data. Replicate models were predicted across the time-series and habitat metrics were calculated as above. Mean, standard error and coefficient of variation were then calculated for each species, MHW, and habitat metric across the 20 models (see Supplementary methods section 1.6.2, Supplementary Table 5).

The absolute magnitude of change to species’ range extent and habitat area varied greatly across species, due in part to differences in occurrence extent (e.g., seabirds have large ranges relative to California sea lions). In order to make MHW impacts comparable across species, range compression/expansion and habitat area loss/gain were expressed as a percent change from mean conditions.

We quantified how species core habitats redistributed across national jurisdictions during MHW events. Version 11 of Exclusive Economic Zone (EEZ) boundaries were used to define national jurisdictions (marineregions.org)59 (Supplementary Fig. 1). The amount of each species core habitat within the EEZs of Canada, the US, and Mexico, and within the high seas was calculated for each day in August–October of 2000–2020. For each species, the average core habitat area in each jurisdiction was calculated for August–October of each MHW year, and the August–October 2000–2020 baseline was subtracted to obtain the anomaly. Anomalies were expressed as a percent change relative to baseline conditions. Habitat redistribution was summarized as both a percent change for each jurisdiction/MHW, and as the largest percent loss and gain for each species. Mean anomaly, standard error, and coefficient of variation were calculated for each species, MHW, and jurisdiction across the 20 replicate BRT models (see Supplementary methods section 1.7, Supplementary Table 7).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "General comments:\nThe manuscript \"Idiosyncratic effects of marine heatwaves on top predator predicted habitat within\nthe Northeast Pacific\" presents novel findings on the redistribution of top predators during extreme\nenvironmental events – marine heat waves (MHW). The analysis is based on a large and valuable\ndataset, producing important results for our understanding of the responses of marine species to\nabrupt changes in the environment. The main noteworthy result of the manuscript is the variable\neffect of MHW on core habitat (in terms of direction, distance, percent change range, or area)\nacross species and MHW events, which has further consequences for variable cross-jurisdictional\nredistributions of these important organisms. The original results presented in this manuscript are\nbased on the unprecedented data set, which considers multiple species and multiple MHW in the\nstandardized analysis. In my opinion, this is the main strength of this work. The analytical\napproach applied in this analysis is sound, meets the expected standards in the field, and I have\njust moderate comment to the methodological part. The conclusions and claims are supported by\nthe data and the results. The main results of this manuscript point out significant challenges to the\nfuture management of marine fauna in the periods of MHW. Importantly, the Authors extend the\nwork with the presentation of the potential way forward – an operational dynamic ocean\nmanagement tool that can potentially be used for the prediction of the distribution of the top\npredators in the near real-time. I have to congratulate the authors for their work which can be a\nvaluable contribution to the existing literature. My comments are mainly composed of suggestions\nthat can improve the presentation of the results or the manuscript’s clarity. I would like also to ask\nfor a few clarifications. My only critical points are i) lack of appreciation of uncertainty in the\npredictions of the core habitat and the calculated metrics of habitat change, ii) the location of\ncentroids of some core habitats on land, which is an artefact of the method, and iii) the operational\ndynamic ocean management tool is only briefly introduced and lacks serious evaluation in my\nopinion. Specifically,\ni) Prediction errors are present in all habitat models (Barry and Elith, 2006). Booster regression\ntrees (BRT) are powerful and handy tools to predict species’ habitats, but they are known to\nproduce different results with different subsets of data. While it could be a demanding task, a\ndeeper evaluation of the source, magnitude, and pattern of errors and uncertainties is essential if\nthe models are to be used transparently in decision-making (Barry and Elith, 2006). It would add\nsignificant strength to this manuscript.\nOne of the potential solutions would be to use replicate BRT models for each species, using\nindependent random samples data, predict the core habitat (and metrics) for each replicate and\nthen calculate uncertainty based on the results of these replicates. The models developed during\nthe cross-validations can be used for this purpose. A similar approach was applied by (Barker et\nal., 2014) to assess model performance and provide a map of the mean prediction uncertainty\nacross species. Another example can be found in (Woodman et al., 2019). The among model\nstandard error (SE) of the ensemble predictions and the associated coefficient of variation (CV)\n‐\nwere used to calculate map cell-specific uncertainty. But here among model variance can be used\nto assess uncertainty in the final metrics.\n‐\nii) I understand that some predicted habitat centroids (and ranges) are located on land due to the\nspecific distribution of the species along the coast and in the Gulf of California. But is it valid? I\nfound them unrealistic in some cases. This is my main reservation which might have consequences\nfor the calculated metrics and general results.\niii) I argue that the leave-one-year framework is inadequate for the assessment of the prediction\nmodel that supposes to work in a near real-time framework. These are two different time\nresolutions, and I think that some patterns in the daily distributions can be poorly recreated by the\nmodel, even though relatively good performance is observed for year averages.\nSpecific comments:\n-Line 1 (title): While I find “Idiosyncratic” an accurate term, it might be difficult to grasp for some\nreaders. Please consider a change to a more common and straightforward word.\n-lines 33-35: Why results for 2015 and the US zone specifically are reported in the abstract? I\nwould suggest something more general\n-Line 179 (Figure 1): I assume that areas under the silhouettes of the studied species are\ncharacterized by low densities and low variability, but my preference would be to move the\nsilhouettes out of the maps (or exclude them completely, since the species are listed in the\ncaption) to show undisturbed information. Instead, consider showing EEZ borders and indicating\neach country/high seas here.\n-Lines 108-110: In my understanding, this analysis showed irregular responses across species and\nMHW, indicating challenges and difficulties of future predictions, rather than serving directly as a\nbasis for informing future management. Please reconsider this sentence.\n-Line 127 (Figure 2). It would be easier to match groups if they are indicated on the left side of the\nfigure (next to the species names). In B you might consider adding arrows (actual directions, with\nN up, S – down, etc.) on top of the colored cells. If uncertainty can be measured for these metrics\n(see my general comment), it can be presented using transparency (e.g. alpha in ggplot2 package\nof R), where more transparent cells indicate more uncertain results.\n-Line 135: I’m wondering if a second-order analysis can be conducted to assess the variance\nassociated with each component (species/event). It could be a model with two random effects:\nspecies (14 levels factor) and MHW event (4 levels factor) fitted to each metric. Modeling\ndirections can be more problematic.\n-Line 143: It is not clear if to 22-64% of the initial area or by 22-64%. The same with 335-377%\nand other percentages in this paragraph.\n-Lines 149-153: it is not clear to me if the centers of crosses indicate the center of gravity. I argue\nthat crosses do not indicate process (expansion/compression), but simply a range of core habitat.\nI suggest a change (here and in extended data) to: “crosses indicate a range of core habitat” and\nmentioning that the black cross is the climatological mean, while colored crosses are MHW-specific\ncore habitat ranges. And “Inset: percent change in core habitat area relative to climatological\nconditions.”.\n-Line 167 (Figure 4): I suggest dropping silhouettes in Figure 4 A as they are not especially\ninformative. If the color legend can be moved to the bottom part (below 0), the upper part of this\nsubplot can be used to create an inset with a map of the study area and indicated EEZ/high seas.\nThen, arrows can be drawn to link the groups of bars and specific areas on the map. If uncertainty\nin the prediction can be assessed (see my main comment), I would add confidence intervals to the\nbars (both in Figure 4A and Figure 4B).\n-Lines 187-222: Figures and extended data figures are intensively referred to in this part of the\nmanuscript. These results, if essential, should be placed in the Results section and the only main\ninterpretation/conclusion presented in the Discussion section.\n-Lines 181-185: Coastal and Southern species are described in this paragraph, but there is no\ninformation about Northern species.\n-Line 191: Indicate that you refer to the results of this study. “Here,” at the beginning?\n-Lines 209-212: Is there a reference for this claim?\n-Line 278: Interesting observation based on the presented results is that there are irregular\nresponses of different species across MHW, but they are not unpredictable. It is quite opposite -\nusing a range of environmental variables, including variables other than temperature alone, it is\npossible to predict species redistributions. I am not sure if this perspective is strong enough in the\ncurrent version.\n-Line 325: why do you refer to extended data figure 2?\n-Lines 353-361: day of the year is missing in this paragraph while it was included in extended data\ntable 1. The second sentence of this paragraph repeats what is in the first sentence – delete\n“including temperature, chlorophyll-a, oxygen, and primary productivity”.\n-Lines 363-365: it would be good to provide at least one reference for the methods of pseudo-\nabsence (as in the supplementary materials).\n-Line 449: I have the impression that the Authors refer to different figures in different parts of the\nmanuscript, but in some places, it is not essential, while distracting from the main message.\n-Line 597: day of the year is not mentioned in the main text.\n-Line 598: I suggest keeping the same order of species consequently throughout the manuscript\nfor easier reading. There are two species (white shark, and leatherback turtle) that are missing in\nthis table. If there are no observations in the independent data sets – they should be included in\nthe table with zeros.\n-Line 620 (Extended data figure 2): It is difficult to assess the size of points because it is\n“dominated” by the numerous observations on the elephant seal. I suggest using a log scale to\nsize points accordingly (log-transformed number of observations).\nAlso, since colors are not easy to distinguish (e.g. tunas), I strongly suggest using faceting by\nspecies, without grouping.\nI encourage authors to test the differences in the accuracy between MHW years and non-MHW\nyears. It would provide important information on how well the model predicts species distribution\nunder MHW in comparison to “normal” years. For example, in albacore tuna, the accuracy obtained\nfor the MHW years seems to be above the long-term average.\n-Line 633 (Extended Data Figure 3): Some predicted habitat centroids and ranges (e.g. California\nsea lion) are located on land. I can imagine that it is some kind of artefact of the calculation\nprocedure, e.g. California sea lions are distributed along the coastline and are present in the Gulf\nof California, thus calculated centroid and range falls on land. The example species selected for the\nmain body of the manuscript are free of these issues. Shouldn’t these centroids be truncated to\nthe marine areas before calculating metrics? For example, ArcGIS tools allow the creation of\ncentroid points located inside or contained by the bounds of the input feature (“contained by\ninput features” option), while GDAL/Shapely for Python allows calculating representative point\n(representative_point() function) which is located inside the polygon. Most importantly, how this\naffects the calculation of metrics? How this affects the general results?\n-Line 643: It would be helpful to indicate in the caption that these anomalies and baseline\nconditions are calculated for August-October. Please add units for all variables.\n-Line 649: bathymetry gained high importance in the models. This variable can be interpreted as a\nproxy of other conditions (temperature, oxygen, etc.), thus taking over the significance. A short\ncomment might be valuable in the appropriate section.\n-Line 652: Overall, it looks OK, but I am not able to properly assess this figure – it has too low\nresolution.\n-Supplementary materials:\n-line 56: day of the year – static?\n-line 172-174: according to the extended data table 2, there were numerous observations on\nelephant seal – only yellowfin tuna had inadequate data, but some species were missing in the\ntable, as indicated in the comment above\n-line 199: please indicate which measure was used for the assessment of the variable importance\n-line 244-248: multivariate/multi-covariate - multivariate methods are not the same as\nmultivariable methods. Multivariate methods have more than one dependent variable, while\nmultivariable methods have one dependent variable and more than one independent variable or\ncovariates. Please correct and use one term consequently.\n-line 300: Just a thought: MHW years are included in the climatological mean (2000-2020). If\nthese years with abnormal environmental conditions (2014, 2015, 2019, 2020) are excluded, the\nredistribution would be even more pronounced.\n-line 310: Since the Southern group is represented by two species, the standard deviation is\ncalculated with two values. Is it meaningful? At least, indicate the number of species included in\neach group.\n-line 349: The Authors provided information here and in the Reporting Summary that R code is\navailable through the GitHub repository, but it was not available during the review (404-page not\navailable). It would help understand some of the steps of the analysis, however, the study doesn’t\nintroduce novel algorithms or methods, thus R code review was not crucial.\nReferences:\nBarker, N. K. S., Cumming, S. G., and Darveau, M. 2014. Models to predict the distribution and\nabundance of breeding ducks in Canada. Avian Conservation and Ecology, 9.\nBarry, S., and Elith, J. 2006. Error and uncertainty in habitat models. Journal of Applied Ecology,\n43: 413–423.\nWoodman, S. M., Forney, K. A., Becker, E. A., DeAngelis, M. L., Hazen, E. L., Palacios, D. M., and\nRedfern, J. V. 2019. esdm: A tool for creating and exploring ensembles of predictions from species\ndistribution and abundance models. Methods in Ecology and Evolution, 10: 1923–1933.", "reviewer_2": "Given the increasing occurrence and prevalence of marine heatwaves - MHWs - it has been\nsurprising that their impact on the distribution of marine animals has received so little attention.\nThis is a very timely submission, that I believe will be of significance to the field, especially given\nthat it is the first look at how the responses of marine megafauna might vary in response to short-\nterm extreme events.\nThe results are, hence, quite novel, with the authors predicting distribution shifts during MHW\nyears. The predicted impacts were, unsurprisingly, variable (range expansions/contractions) for\nthe different species considered; but importantly, different MHWs were predicted to have impacted\ncertain species (e.g., shearwaters; sea lions) differently.\nIt would have been ideal if tracking data spatio-temporally overlapped the MHWs, but the authors\nbuilt a robust predictive framework to work around the gap in data. The authors did a great job\nexploring and presenting their results, especially when linking redistribution patterns in relation to\npolitical boundaries.\nThe manuscript is very well written and structured. I have no major comments and am supportive\nof its publication.\nMINOR COMMENTS:\n#1 I understand it is quite hard to make figures fit the limited space that journals provide, but\nFigure 1 should be made bigger since it is especially hard to see the kernel densities. The large\nnumber of species outlines for the coastal plot also does not help. This figure should also include a\nlegend for the different species (like Figure 2).\n#2 The insets in Figure 3 should indicate what the axis refer to (percent change in area).\n#3 L192-194; In the first paragraph of the Discussion the summary of the results is sometimes\nvague. For example \"[...] some predators were predicted to experience near total loss of habitat\nand range compression, while others were predicted to experience a two-fold habitat increase and\nsignificant range expansion.\" It would be better if authors mention which predators were predicted\nto experience range expansions/contractions.\n#4 L209. \"Sooty shearwaters were predicted to expand their range during the 2019 and 2020\nevents, likely associated with elevated primary productivity in coastal Alaskan waters\". It is very\nminor comment, but Alaska is not shown in Figure 4 under the US EEZ.\nI can only congratulate the authors on an excellent manuscript. Also, the supplementary materials\nare quite thorough, and the authors clearly aware of the caveats of the approach they used.", "reviewer_3": "Review of “Idosyncratic effects of marine heatwaves on top predator predicted habitat within the\nNortheast Pacific” by Welch et al.\nThis study investigates the effects on top predator species of multiple large-scale marine heatwaves,\nassociated with ‘The Blob’, in the Northeast Pacific from 2014-2019. Telemetry data is used together\nwith selected environmental variables to fit a species distribution model to predict distributions during\nmarine heatwave events. The authors present different responses amongst species but also amongst\nheatwave events and potential consequences for ecosystem management due to shifting distributions\ninto or out of US waters. The results highlight the urgent need to develop operational ocean\nmanagement tools.\nThe manuscript is well written and results presented in a clear way. For full disclosure, this review\nfocuses on the oceanographic aspects and marine heatwave definition of the manuscript. I recommend\nthe manuscript for publication after some (potentially minor) revisions. Detailed comments are outlined\nbelow.\nComments\nI have a few more general comments about the definitions of marine heatwaves and associated\nterminology.\nThe authors do not specifically detect MHWs here but build on existing literature that has identified the\n4 big MHW events in the specific years that are referred to here. For completion it might be helpful to\nadd some MHW metrics about those 4 events in the supplementary, e.g. start, end and peak date since\nmarine heatwaves are defined as discrete events. I understand that these might not be as relevant for\nthis study as the focus is on the months with the warmest total temperatures but if the term MHW is\nused it could be helpful to state them specifically or add a couple of sentences along the lines of what I\nmentioned above.\nFurthermore, the authors speak a few times about averages across each MHW event (e.g. in L299) which\ncan be misleading. Instead, I would say something like Aug-Oct average during each MHW year.\nThe years 2000-2020 are used as baseline to derive anomalies and the authors refer to it as climatology. I\nwould not call it climatology as that term typically refers to a 30-year average, in particular for MHW\nstudies. Furthermore, the authors justify the use of 2000-2020 with the limitation of the dataset, in\nparticular because of the chlorophyll data, however, SST is available earlier and one could use a more\ncommon baseline (e.g. 1982-2010) to derive the anomalies to be more consistent with previous MHW\nwork. Of course, it is under current debate how to best define a MHW\n(https://pubmed.ncbi.nlm.nih.gov/37012469/) and species’ responses will be very different depending\non the thermal tolerance and potentially depend stronger on the absolute temperatures instead of\nanomalies? So, one overall question is, how sensitive are the results to the magnitude of anomalies and\nit was not clear to me if the habitat model was fit with SST anomalies or absolute SST values?\nThe authors make the point several times in the text that species impacts will vary across MHW events. It\nis worth highlighting in the discussion (e.g. L199-230) that MHWs are just a statistical construct, but they\ncan be caused by very different physical mechanisms, which likely is reflected in other variables too. In\nthe supplementary the authors mention that for example the events in 2019 and 2020 were much\nshallower than 2014 and 2015, which likely drives different responses. In my opinion it is crucial to\ndiscuss the different regional ocean processes associated with different types of marine heatwaves and\nhow these can for example impact oxygen at depth though for example variability in the coastal\nupwelling. I am aware that a detailed discussion of the physical processes is likely out of scope for this\nstudy but still think that it would be important to have a short paragraph in the discussion dedicated to\nthis.\nLine-based comments\nL135ff: Before talking about the effects of MHWs, I think it would be useful to briefly mention that this\nstudy focuses on Aug-Oct only in each MHW year. I am aware that the specifics of each MHW are not the\nfocus of the study but I believe it will help the reader with the interpretation of the results.\nFigures\nFig1: This might be personal taste, but I believe maps should always include some meridians and\nparallels for reference. The authors are likely very familiar with the geography of the region but not all\nreaders will be.\nFig3: Add axis labels or change tick labels to include °N etc.\nFigure S1: The panels are very small and have low resolutions which makes it impossible to read the\nnumbers on the colorbar."}, "author_response": "Response Reviewer 1, Comment 1. Thank you very much for your thorough read of our manuscript and\npositive and constructive comments. We particularly appreciate your suggestions to improve Figure 2,\nand for helping us improve our messaging on the predictability of MHWs.\n\nResponse Reviewer 1, Comment 2. Thank you for bringing up this point and providing guidance for\nhow to evaluate our models’ prediction errors. We followed the methods in Barker et al 2014 and created\n20 replicate BRT models using independent 75% samples of the telemetry and pseudo-absence data.\nThese 280 models (20 replicates x 14 species) were predicted from 2000-2020, and then habitat metrics\nwere calculated (displacement distance, displacement direction, % change range, % change area). For\neach metric, species and marine heatwave year, the mean, standard error, and coefficient of variation were\ncalculated across the 20 model outputs to understand the sensitivity of metrics to the subset of data used\nin model fitting.\nWe have created a new supplementary section on habitat metric sensitivity which contains these results in\naddition to results from our sensitivity analysis on binary vs continuous centroid calculation. Our added\ntext is as follows:\n“In addition, we tested the sensitivity of habitat metrics to the data used in model fitting. For each\nspecies, telemetry and pseudo-absence data were randomly subsampled into 20 different subsets\ncontaining 75% of the data and maintaining the 1:1 ratio of presences to pseudo-absences. A BRT model\nwas fit to each subset (n=280 models; 14 species x 20 data subsets), predicted, and summarized to\ncalculate habitat metrics following methods in Sections 1.5.2., 1.5.6., and 1.6.1. Then, habitat metrics\nwere summarized across the 20 models by mean value, standard error, and coefficient of variation (Table\nS2). Coefficient of variation was included in addition to standard error because values are standardized\nby the mean, controlling for larger values having larger standard errors. We found low sensitivity of\nhabitat metrics to data used in model fitting. The average standard errors across species and MHW years\nwere: displacement distance: 9.81 (±19.2km 95% confidence interval); displacement direction: 2.75\n(±5.39° 95% confidence interval); percent change range: 3.52 (±6.89% 95% confidence interval); and\npercent change area: 1.68 (±3.29% 95% confidence interval). Displacement distance had the largest\nstandard error, but its confidence interval indicates that error (19.2km) was less than the size of one\n25km pixel. The mean coefficient of variations across species and MHW years were also low:\ndisplacement distance: 0.24; displacement direction: 0.18; % change range: 1.14; and % change area:\n0.31.”\nFurther, we added a table to visualize the results and enable comparisons:\nTable S2. Sensitivity of habitat metrics to data used in model fitting. For each habitat metric, species, and\nmarine heatwave (MHW), the mean, standard error (SE), and coefficient of variation (CV) were\ncalculated across 20 models fit to random 75% subsets of the telemetry and pseudo-absence data. SE and\nCV columns have been shaded to ease interpretation such that darkest reds indicate the largest value in\neach column.\n\nResponse Reviewer 1, Comment 3. While a mean centroid is coarse and does not necessarily represent\nwhere species actually are, it still shows whether the relative proportions of habitat in different areas\nchanges over time or not. Our method of centroid calculation is common in the literature for marine\nspecies1–3. In particular, Pinsky et al. (2018 Science)1 Fig. 1A shows individual species centroids on land\natop the US west coast and the Aleutians due solely to the coastal orography.\nTo address this point in the text, we have added an explanation in the supplementary materials:\n“Due to the convex curvature of the US west coast, centroids for some species were located on land\n(Extended Data Figure 3). While a centroid is coarse and does not necessarily represent where species\nactually are, it still captures the relative relocations of preferred habitat among species.”\nFor clarification, predictions of habitat suitability were truncated to marine areas prior to centroid\ncalculation. In addition, our calculations of species jurisdictional shifts (which require absolute locations\nof habitat as opposed to relative shifts of habitat) were directly calculated from core habitat as opposed to\ncentroids. There is low variability in centroid location when models are fit on different subsets of data\n(Reviewer 1, Comment 3), or when centroids are calculated from binary versus continuous habitat\nsuitability (supplementary methods section 1.6.2.):\n“We tested the effect of calculating core habitat centroids on binary core habitat (all core habitat pixels\nassigned a value of 1) versus continuous core habitat (all core habitat pixels maintain values of\ncontinuous habitat suitability predictions). This comparison only affects values of habitat displacement\n(direction and distance) as range extent change and habitat area change do not utilize centroids. We\nfound minimal differences between the two methods of centroid definition, with displacement distance\nvarying on average by 3km, and displacement direction varying on average by 1 degree (Fig. S2).”\nFigure S2. Sensitivity of displacement distance (A,B,C) and direction (D,E,F) to core habitat centroids\ncalculated on binary core habitat (all core habitat pixels assigned a value of 1 (A,D)) versus continuous\ncore habitat (all core habitat pixels maintain values of continuous habitat suitability predictions (B,E)).\nC,F show binary core habitat minus continuous core habitat for habitat displacement distance and\ndirection, respectively.\n\nResponse Reviewer 1, Comment 4. While we are suggesting these models could be used in a daily\ndynamic tool, the focus of this paper is instead the general annual patterns in responses to marine\nheatwaves. For this reason the leave one year out approach is useful to understand how these predictions\nwould remain skillful in years that do not have model-training data, and has become the standard for\ntesting operational models4,5. We have included information about how collecting and assimilating daily\nobservations could be an important next step in the discussion :\n“Real-time predictions of species distributions during MHWs could be used to guide observational\nprograms (e.g. field surveys and tagging programs) during anomalous conditions. Importantly, these new\nobservations could be assimilated with distribution model outputs to improve real-time model predictions,\na common process in oceanography which is not yet a standard practice in ecology”\nIn light of this comment, we have added a supplementary analysis of the daily ratio of observed to\npredicted presences from the telemetry datasets used in model fitting to test how well the model performs\non withheld data at this higher temporal resolution. We found that the models did perform well and thus\ncould be operationalized in a daily time step:\n“Additionally, we calculated the daily ratio of observed to predicted values at presences from the\ntelemetry datasets used in model fitting (Supplementary methods Figure S3) to evaluate the predictive\ncapabilities of models in a near-real time framework. Each presence was assigned a 1 and then summed\nto calculated observed values at presences. To calculate predicted values at presences, model-predicted\nvalues at each presence were summed. Ratios of 1 indicate perfect agreement between observed and\npredicted presences (i.e. all model predictions at presences were 1s), while ratios above and below one\nindicate under and over prediction, respectively6. Ratios showed high alignment between daily observed\nand predicted presences (mean 1.5 +/- 1.3 SD), indicating strong capacity to generate predictions in a\nnear-real time framework.”\nFigure S3. The daily ratio of observed to predicted presences for each of the boosted regression tree\nmodels. Ratios of 1 indicate perfect agreement between observed and predicted presences, while ratios\nabove and below one indicate under and over prediction, respectively.\n\nResponse Reviewer 1, Comment 5. Thank you for the suggestion to make it more accessible. We have\nchanged the title to “Impacts of marine heatwaves on top predator distributions are variable but\npredictable”.\n\nResponse Reviewer 1, Comment 6. Those results were reported because they represent the greatest\nredistribution, but we agree that this level of detail is not needed in the abstract. We’ve changed this\nsentence to read:\n“MHWs rapidly altered political bio-geographies, with up to 10% of predicted habitat across all species\nshifting jurisdictions during individual MHWs.”\n\nResponse Reviewer 1, Comment 7. Thank you for this feedback. We received several comments\nsuggesting edits to this figure, including adding EEZ borders and country names, adding meridians and\nparallels, and that it was too small to see the kernel densities and contours (Reviewer 1, Comment 7;\nReviewer 2, Comment 2; Reviewer 3, Comment 7). However, even making the maps larger, we feel it is\nstill too small to show EEZ borders and still convey the kernels and contours (which is the primary\ninformation in this figure):\nIn lieu of adding EEZs to this figure, we have added an inset map showing EEZs to Figure 4, as suggested\nby Reviewer 1, Comment 13. As requested, we have made the maps in Figure 1 bigger, and have removed\nthe species outlines:\n\nResponse Reviewer 1, Comment 8. Thank you for flagging this lack of clarity in our messaging. We\nhave rewritten the discussion section on early warning systems to better communicate how our analysis\nserves as a basis for informing future management (in response to Reviewer 1, Comment 18.):\n“Our results show a wide range of species impacts across MHW events and jurisdictions. This diversity\nof responses poses a daunting management challenge: how to plan ahead and respond swiftly to MHW-\ndriven species redistribution. Our results indicate that species responses to MHWs are highly variable yet\npredictable: our models performed well through extensive validation across space, time, and on novel\ndata (Extended Data Figure 2 and Table 3). The high variability of species responses to MHWs suggests\nthat we cannot assume future MHWs will impact species in the same way as past events. However, high\npredictability indicates that species responses to future MHWs could be predicted in real-time to provide\naccurate information on impacts. This combination of high variability and predictability is also seen in\nhurricanes – future hurricane paths cannot be extrapolated from past paths, yet they can be predicted in\nreal-time to provide accurate information on risk.”\n\nResponse Reviewer 1, Comment 9. Thank you, it's a good idea to add direction arrows and move\nspecies groups to the left side of the figure, we have made both of these changes:\nFigure 2. Predicted impacts on top predator habitat within (columns, e.g. 2014) and among (rows, e.g.\nWhite shark) marine heatwave events measured using four metrics: A. displacement distance\n(kilometers), B. displacement direction (degrees, where 0/360 is north (N), 90 is east (E), 180 is south (S),\nand 270 is west (W)), C. range compression or expansion (percent change relative to baseline\nconditions), D. habitat area gain or loss (percent change relative to baseline conditions). All metrics\nwere calculated from August-October in each MHW year relative to baseline conditions (August-October\n2000-2020), see Table S2 for an analysis of metric uncertainty. Northern, Coastal, and Southern regional\ngroupings indicate the geographies where the majority of the species telemetry data occurs.\nHowever, we have elected to not use alpha to indicate uncertainty simply because it was extremely\ndifficult to discern the original intent of the figures. Visualizing uncertainty using alpha masks the\nrelationship between the tile and legend colors, making it challenging to communicate magnitudes. For\nexample, in the plot below, the left panel doesn’t use alpha, while the right panel maps alpha to\nuncertainty. In addition, because the habitat metrics vary across different magnitudes (e.g. distance ranges\n(0-700 km) versus percent change range (-100% to 370%), the standard errors aren’t comparable across\nhabitat metrics (larger numbers will have larger standard errors and lower alphas, penalizing habitat\nmetrics that vary over larger magnitudes). Lastly, to include uncertainty in this figure, we would need to\nchange the habitat metric values currently displayed (outputs from the full model containing 100% of the\ndata) to the mean habitat metrics averaged across the 20 model replicates of 75% data (the mean\ndescribed by standard error). However, we feel that outputs from the full model containing 100% of the\ndata are the most accurate, given that they were derived from the most data.\nInstead, we have added Table S2 that shows, for each species, MHW and habitat metric, the mean,\nstandard error, and coefficient of variation across 20 replicate models, and point to it in the caption of\nFigure 2 (see Response Reviewer 1, Comment 2.).\n\nResponse Reviewer 1, Comment 10. We have added a new section and Table S3 to the supplementary\ninformation following this suggestion (section 1.6.3.). As you mention, it is challenging to build models\nfor direction, which is a circular variable and must be handled differently than the other metrics. There are\nseveral functions we explored for this purpose (lm.circular() and bpnr()), but we had concerns about how\nto appropriately structure and interpret circular linear models, and how to compare these models with\nmodels for the other habitat metrics. As a result, we have elected to not model direction, but have added\nthe following text and table for the other habitat metrics following your suggestions:\n“We explored the extent to which variability in continuous habitat metrics was associated with species vs\nMHWs using linear mixed effect models. For each habitat metric (distance, percent change range,\npercent change area), we built two models with one fixed effect and one random effect, e.g. for percent\nchange area:\nMod1: lmer(percent change area ~ MHW + (1|species), data = dat)\nMod2: lmer(percent change area~ species + (1| MHW), data = dat)\nIn Mod1, MHW is fixed and species is random, and in Mod2, MHW is random and species is fixed. Next,\nwe compared the marginal R2 values for each model, i.e. the R2 considering only the fixed effect and\nignoring the random effect. This comparison allowed us to determine if MHW explains more variance\nwhile controlling for species (Mod1), or if species explains more variance while controlling for MHW\n(Mod2) (Table S3).\nTable S3. Linear mixed models for percent change range, percent change area, and distance habitat\nmetrics. Each model contained one random effect (either species or marine heatwave (MHW)) and one\nfixed effect. Marginal R2 values indicate the variance explained by the fixed effect while controlling for\nthe random effect.\nResults indicate that species explains more of the variance than MHW for percent change range, percent\nchange area, and distance habitat metrics.”\n\nResponse Reviewer 1, Comment 11. The percentages in this paragraph are percent increases and\ndecreases, i.e. the difference between the final and initial values, expressed as a percentage of the initial\nvalue. We have clarified this in the caption of Fig. 2 which is referenced in this paragraph, and now report\npercentages as “% increase” and “% decrease”.\n\nResponse Reviewer 1, Comment 12. We agree that our wording was not clear, and have edited the\ncaption of Figure 3 and Extended Data Figure 3 to better distinguish pattern from process:\n“Large maps: arrows indicate predicted habitat displacement from the center of gravity in each MHW\n(center of colored crosses) relative to the center of gravity during baseline conditions (center of black\ncross). Crosses indicate predicted longitudinal and latitudinal range extents during each MHW and\nduring baseline conditions (colored and black crosses, respectively); an increase in cross size during a\nMHW compared to baseline conditions indicates range expansion with the converse representing\ncompression. Inset: percent change habitat area relative to baseline conditions.”\n\nResponse Reviewer 1, Comment 13. We appreciate your suggestion to drop the silhouettes from Figure\n4A and include a map of the study area indicating jurisdictional boundaries, and have made both changes.\nWe elected to not draw arrows from each bar group to each map jurisdiction because it felt too busy.\nFigure 4. Cross-jurisdictional shifts in predators’ predicted habitats during marine heatwaves. A. Total\nloss and gain of species habitat area across the US, Mexican (MEX), and Canadian (CAN) Exclusive\nEconomic Zones (EEZs), and the high seas during each MHW. B. Largest loss and gain of habitat for\neach predator in any MHW event (x-axis text color indicates species regional groupings: Southern\n(green), Coastal (black), Northern (blue)). Percent change in habitat area is calculated relative to\nbaseline conditions (2000-2020), see Table S4 for an analysis of uncertainty in cross-jurisdictional shifts.\nIn lieu of adding uncertainty bars to this figure, we have created a new supplementary Table S4 and now\npoint to it in the Figure 4 caption. To include uncertainty in this figure, we would need to change the\nanomaly values currently displayed (outputs from the full model containing 100% of the data) to the mean\nanomaly values across the 20 model replicates of 75% data (the mean described by standard error).\nHowever, we feel that outputs from the full model containing 100% of the data are the most accurate,\ngiven that they were derived from the most data.\nHere is our new supplementary text and table describing the analysis of jurisdictional uncertainty:\n“Similar to the habitat metrics, we tested the sensitivity of species core habitat redistribution across\njurisdictions to the data used in model fitting by fitting and predicting 20 replicate BRT models (Section\n1.6.2.). Habitat redistributions were calculated as above, and summarized across the 20 models by mean\nanomaly, standard error, and coefficient of variation (Table S4). The average standard errors and\ncoefficients of variations across species, MHW years, and jurisdictions were low: 0.0031 (±.0061% 95%\nconfidence interval) and 1.47 for standard error and coefficient of variation, respectively.\nTable S4. Sensitivity of species core habitat redistribution across jurisdictions to data used in model\nfitting. For each jurisdiction, species, and marine heatwave (MHW), the mean anomaly, standard error\n(SE), and coefficient of variation (CV) were calculated across 20 models fit to random 75% subsets of the\ntelemetry and pseudo-absence data. SE and CV columns have been shaded to ease interpretation such\nthat darkest reds indicate the largest value in each column.”\n\nResponse Reviewer 1, Comment 14. We have removed redundant extended data figure references, and\nlimited our references to main text figures in the Discussion.\n\nResponse Reviewer 1, Comment 15. We now include the following text on northern species:\n“Northern species experienced the largest redistributions during the 2020 event, with 14% of salmon\nshark habitat predicted to shift into the Canadian EEZ.”\n\nResponse Reviewer 1, Comment 16. Good point, we have now added this for clarification..\n\nResponse Reviewer 1, Comment 17. We have added a new supplementary figure showing the conditions\nof variables beyond temperature to support this claim:\nFigure S1. Anomalies of sea surface temperature (°C; A), oxygen concentration at 200m (mmol m-3; B),\nmean primary productivity in the upper 200m (mg m-3; C), and log chlorophyll-a (mg m-3; D) during each\nmarine heatwave event. Anomalies were calculated from August-October in each year relative to a 2000-\nWe have also included a couple sentences on regional oceanographic processes in the section to help\nexplain sooty shearwater and elephant seal responses:\n“However, temperature anomalies were fairly consistent across the 75th kernels of sooty shearwaters and\nelephant seals during the four MHW events, yet each species experienced markedly different predicted\nimpacts among MHWs. Sooty shearwaters were predicted to expand their range during the 2019 and\n2020 events, likely associated with elevated primary productivity in coastal Alaskan waters (Fig. S1).\nElephant seals were predicted to experience large habitat loss during the 2014 and 2015 events, likely\ndue to a reduction of primary productivity and subsurface oxygen in the north-central Pacific (Fig. S1).\nThe deep warm water anomalies of the 2014 and 2015 events decreased the solubility of oxygen and led\nto a near-surface decline in oxygen concentration7. The 2019 and 2020 events, in contrast, had much\nshallower temperature anomalies, and likely had a weaker impact on local biogeochemical signals and\nproductivity. In general, variations in winds, air-sea gas exchange, circulation patterns and water column\nstratification can lead to very different physical and biogeochemical conditions among MHW events. The\nidiosyncratic responses of species to MHWs reflect these variable conditions and differing species-\nenvironment relationships (Extended Data Figure 5).”\n\nResponse Reviewer 1, Comment 18. We appreciate this comment, we’ve had some trouble articulating\nthat while it is not possible to extrapolate future MHW impacts based on historical impacts, it is possible\nto predict impacts in near real-time. We’ve rewritten this section using a more familiar example -\nhurricanes, which are highly variable but predictable:\n“Our results show a wide range of species impacts across MHW events and jurisdictions. This diversity\nof responses poses a daunting management challenge: how to plan ahead and respond swiftly to MHW-\ndriven species redistribution. Our results indicate that species responses to MHWs are highly variable yet\npredictable: our models performed well through extensive validation across space, time, and on novel\ndata (Extended Data Figure 2 and Table 3). The high variability of species responses to MHWs suggests\nthat we cannot assume future MHWs will impact species in the same way as past events. However, high\npredictability indicates that species responses to future MHWs could be predicted in real-time to provide\naccurate information on impacts. This combination of high variability and predictability is also seen in\nhurricanes – future hurricane paths cannot be extrapolated from past paths, yet they can be predicted in\nreal-time to provide accurate information on risk.”\n\nResponse Reviewer 1, Comment 19. Thank you, this was an error caused by a reordering of the figures.\nWe have removed this figure reference.\n\nResponse Reviewer 1, Comment 20. Thank you, we have now included day of year, and removed the\nrepetitive text.\n\nResponse Reviewer 1, Comment 21. This is a good point. We have included the following reference:\nBarbet‐Massin, Morgane, et al. \"Selecting pseudo‐absences for species distribution models: How, where\nand how many?.\" Methods in ecology and evolution 3.2 (2012): 327-338.\n\nResponse Reviewer 1, Comment 22. We have removed these figure references.\n\nResponse Reviewer 1, Comment 23. We now include day of year in the Online Methods.\n\nResponse Reviewer 1, Comment 24. We have added white shark and leatherback to the table and\nindicated that they have no records. We have also reordered Extended Data items Table 2, Table 3,\nFigure 2, Figure 3, Figure 5, and Figure 6 to list species alphabetically. Main text Figures 2, 4 and\nExtended Data Figure 4 have different orders to better communicate the regional species groups and the\nmagnitudes of cross-jurisdictional shifts that are presented in the figures.\n\nResponse Reviewer 1, Comment 25. We have made both of these changes:\n\nResponse Reviewer 1, Comment 26. We have performed a t-test comparing AUC during MHW vs\nnormal years and added the results to the figure caption:\n“Across all species, there was not a significant difference in AUC between marine heatwave years\n(mean= 0.78) and normal years (mean=0.78); t(52.20) = -0.27, p-value = 0.791.”\n\nResponse Reviewer 1, Comment 27. Please see our full response in Reviewer 1, Comment 3 - our\nmethod is commonly used in the scientific literature. In addition, we investigated the two tools suggested\nhere as well as a similar R tool and found that none were suitable for our analysis. Our analysis includes\n(for a given species and day): a raster of core habitat patches (A), and optionally the minimum bounding\npolygon which we used to constrain the species predictions (B).\nTool 1. “Contained by input features option” in ArcGIS.\nWe could convert our raster (A) to points to work with this tool, which requires point, line or polygon\ndata. If we checked the “inside” parameter to ensure the output point fell within the spatial domain of the\ninput (note this wouldn’t be (B), it would be the points in (A)), the input point that is nearest to the\ncentroid would be selected as the output: “The output point will be coincident with one of the points in\nthe multipoint”. In other words, the pixel containing species core habitat that is closest to the centroid\nwould be selected, which is not the metric we are aiming to calculate. https://pro.arcgis.com/en/pro-\napp/latest/tool-reference/data-management/feature-to-point.htm.\nTool 2. “representative__point()” in GDAL/Shapely for Python\nThis tool will find the representative point located within a polygon by shifting the true centroid to a\nlocation inside a polygon that balances polygon complexity and shape. This means that for polygons that\nare irregularly shaped such as the minimum bounding polygons (B), portions of the polygon that are\nlarger will pull the representative point towards them. This functionality would change our results from\ndescribing the center of gravity of species core habitat (A) to the center of gravity of the polygon used to\nconstrain predictions (B). Furthermore, this tool finds the representative center of a polygon, as opposed\nto the representative center of a point cluster that is contained within a polygon.\nTool 3. st_point_on_surface() in R\nThis tool works similarly to Tool 2 - for a line or polygon, it will find the center of gravity that falls\nwithin or on the input surface. However, we are aiming to find the centroid of (A) that falls within (B). If\nwe use points as an input (by converting the raster (A) to point), it returns all input points (which are the\ncenter of gravity of themselves) as opposed to one point summarizing all points. Like tool 2, this tool\nfinds the representative center of a polygon or line, as opposed to the representative center of a point\ncluster that is contained within a polygon.\nIn summary, none of these tools achieve the desired functionality, which is to summarize the center of\ngravity of (A) while constraining it within polygon (B). Even if the desired functionality was possible, we\nfeel it could bias results for species for which centroids appear on land and need to be shifted vs species\nfor which centroids appear on water and don’t need to be shifted. We feel the simplest and most\ndefensible (given its prevalence in the literature) way forward is to find geometric centroids for all\nspecies, with the understanding that values do not describe the absolute locations of habitat but rather how\nproportions of habitat in different locations change over time.\n\nResponse Reviewer 1, Comment 28. Thank you, we have updated the caption to read:\n“Anomalies of Sea surface temperature (°C; A), log chlorophyll-a (mg m-3; B), oxygen at 200m (mmol m-\n3; C), and mean primary productivity in the upper 200m (mg m-3; D) measured over the 75% kernel of\neach species. Anomalies were calculated from August-October in each year relative to a 2000-2020\nAugust-October baseline.”\n\nResponse Reviewer 1, Comment 29. It’s true that many of our variables are correlated (correlation plot\nbelow). We’ve included a comment in the supplements that BRTs are robust to collinearity effects8, but\nthat collinearity can obscure relative importance:\n“BRTs are a common machine learning model, popularized by their ability to fit complex nonlinear\nrelationships, automatic handling of collinearity effects, and their robustness to wide varieties of data\ntypes and distributions. Relative importance captures each variable’s influence on the response8, and\nreflects the number of times each variable is selected for splitting across all trees. Importance is scaled\nacross all variables such that the sum adds to 100, producing relative importance, where higher values\nindicate stronger contributions. However, relative importance may not accurately reflect variable\ninfluence when variables are highly correlated.”\n\nResponse Reviewer 1, Comment 30. Apologies, we have remade it at a higher resolution.\n\nResponse Reviewer 1, Comment 31. We have removed the word “static”.\n\nResponse Reviewer 1, Comment 32. Thank you for catching this, we received the elephant sea data last\nand failed to update this sentence to reflect the new data. We did not have adequate validation data for\nyellowfin tuna, white shark, or leatherback turtle and have updated this sentence and extended data table 2\naccordingly.\n\nResponse Reviewer 1, Comment 33. We have added the following text:\n“Relative importance captures each variable’s influence on the response8, and reflects the number of\ntimes each variable is selected for splitting across all trees. Importance is scaled across all variables\nsuch that the sum adds to 100, producing relative importance, where higher values indicate stronger\ncontributions.”\n\nResponse Reviewer 1, Comment 34. Thank you for noticing this inconsistency. We’ve replaced\n“multivariate/multi-covariate” with “multivariable”.\n\nResponse Reviewer 1, Comment 35. While it is true that excluding the warmest years from the\nclimatology would separate them even more from mean conditions, we include all years in the\nclimatology for consistency with the MHW literature (and with calculation of anomalies more\ngenerally)9–11.\n\nResponse Reviewer 1, Comment 36 We have indicated the number of species in each group:\n\nResponse Reviewer 1, Comment 37. Our apologies for this oversight. We have changed the repository\nfrom private to public https://github.com/HeatherWelch/MHW_impacts_top_predators\n\nResponse Reviewer 2, Comment 1. Thank you very much for your encouraging remarks.\n\nResponse Reviewer 2, Comment 2. We’ve done our best to increase the size of this figure by removing\nwhite space, and moving the map labels onto land to free up more space for the maps. We’ve also\nremoved the species outlines, as this comment and and Reviewer 2, Comment 7 noted they obscured the\nkernel densities:\n\nResponse Reviewer 2, Comment 3. Done.\n\nResponse Reviewer 2, Comment 4. We have edited this sentence to read:\n“Predicted responses were highly variable across species and MHW (Fig. 2): some predators were\npredicted to experience near-total loss of habitat and range compression, e.g. bluefin tuna during the\n2015 event, while others were predicted to experience a two-fold habitat increase and significant range\nexpansion, e.g. California sea lion during the 2019 event.”\n\nResponse Reviewer 2, Comment 5: True, we consolidated US west coast, Alaskan, and Hawaiian waters\ninto one group to simplify Figure 4. We have added a new supplementary figure to support this statement:\n“Sooty shearwaters were predicted to expand their range during the 2019 and 2020 events, likely\nassociated with elevated primary productivity in coastal Alaskan waters (Fig. S1).”\nFigure S1. Anomalies of sea surface temperature (°C; A), oxygen concentration at 200m (mmol m-3; B),\nmean primary productivity in the upper 200m (mg m-3; C), and log chlorophyll-a (mg m-3; D) during each\nmarine heatwave event. Anomalies were calculated from August-October in each year relative to a 2000-\n\nResponse Reviewer 3, Comment 1. This is a good point. We have now edited this discussion paragraph\nto include start, end, and peak aerial extent data for each MHW:\n“Prominent marine heatwaves (MHWs) in the North Pacific occurred in 2014, 2015, 2019, and 2020\n(Extended Data Figure 1). Warm water conditions persisted across 2014 and 2015, however we elected\nto examine these years as distinct events due to their distinct warming patterns. Strongly positive\ntemperature anomalies in the Gulf of Alaska appeared in October 2013, expanding to reach the US west\ncoast in September 2014, penetrating to over 100m in depth across large portions of the region12 (Fig.\n1A). In fall of 2015, this pre-existing heatwave was amplified by a strong El Niño event, extending the\nspatial extent of the MHW westward towards Hawaii, increasing the depth penetration of warm\nanomalies, and producing the warmest tropical SST anomalies on record13,14. This MHW dissipated in\nlate 2015. The 2019 and 2020 MHWs originated in the Gulf of Alaska during late spring (April-May) and\nreached the US west coast in September, though there were differences in their seasonality and drivers15.\nThe 2020 event was the second largest MHW on record (following the 2014 event), and at its maximum\nextent, the 2019 event was only 6% smaller than the 2014 event16. Despite their comparable sizes, the\ndepth penetration of warm anomalies was much shallower in 2019 and 2020, reaching maximum depths\nof 40-50m12. Both of these events rapidly retracted from the US west coast during October and November\nof their respective years, and dissipated below thresholds for marine heatwave classification11 during\nDecember and January. Maximum aerial extents for each event were reached in September 2014, July\n2015, September 2019, and September 2020.”\n\nResponse Reviewer 3, Comment 2. Thank you, we have clarified that averages are calculated from Aug-\nOct in each MHW year throughout the text, and replaced all instances of “climatology” with “baseline”.\n\nResponse Reviewer 3, Comment 3. The habitat models were fit with absolute SST values as opposed to\nSST anomalies, thus none of our biological response metrics are sensitive to baseline length. While the\navailable SST time-series is longer than that of the other variables, we have elected to maintain all\nanomaly calculations (e.g. Extended Data Figure 4) across the same time-series to allow anomalies to be\ncomparable across variables and with biological responses (e.g. Fig 2).\n\nResponse Reviewer 3, Comment 4. We have added some regional oceanographic context to explain the\npatterns seen in sooty shearwaters and elephant seals:\n“However, temperature anomalies were fairly consistent across the 75th kernels of sooty shearwaters and\nelephant seals during the four MHW events, yet each species experienced markedly different predicted\nimpacts among MHWs. Sooty shearwaters were predicted to expand their range during the 2019 and\n2020 events, likely associated with elevated primary productivity in coastal Alaskan waters. Elephant\nseals were predicted to experience large habitat loss during the 2014 and 2015 events, likely due to a\nreduction of primary productivity and subsurface oxygen in the north-central Pacific. The deep warm\nwater anomalies of the 2014 and 2015 events decreased the solubility of oxygen and led to a near-surface\ndecline in oxygen concentration7. The 2019 and 2020 events, in contrast, had much shallower\ntemperature anomalies, and likely had a weaker impact on local biogeochemical signals and productivity.\nIn general, variations in winds, air-sea gas exchange, circulation patterns and water column\nstratification can lead to very different physical and biogeochemical conditions among MHW events. The\nidiosyncratic responses of species to MHWs reflect these variable conditions and differing species-\nenvironment relationships (Extended Data Figure 5).”\nAnd an additional supplementary figure showing conditions of oxygen, primary productivity, and\nchlorophyll-a during the MHW years:\nFigure S1. Anomalies of sea surface temperature (°C; A), oxygen concentration at 200m (mmol m-3; B),\nmean primary productivity in the upper 200m (mg m-3; C), and log chlorophyll-a (mg m-3; D) during each\nmarine heatwave event. Anomalies were calculated from August-October in each year relative to a 2000-\nWe are in full agreement that variables beyond temperature are driving species response to MHWs, and\nhave a paragraph devoted to this in the discussion:\n“MHW-driven species redistributions are most commonly ascribed to the relocation of preferred\ntemperature conditions17–19. However, our findings suggest that a multi-variable approach allows for\nadditional ecological inferences with regard to MHW biodiversity impacts. Many of the predators\nexamined here have broad thermal tolerances, and are distributed in warmer waters elsewhere in the\nPacific than those encountered during the MHWs20. Thus it is likely that anomalous conditions of\nvariables beyond temperature (Extended Data Figure 4) are driving their responses. Indeed, a\ncomparable suite of temperature-only models had significantly worse predictive performance on novel\nvalidation data than the multivariate models (median Area Under the Receiver Operator Characteristic\nCurve of 0.6 in the temperature-only models vs 0.8 in the multivariate models, t-test p-value <0.05). In\nfollowing with evidence that temperature alone cannot account for species responses to climate\nchange2,21, we suggest that a multi-variable approach is critical to capture species responses to short-\nterm warming. Several programs exist to monitor and forecast extreme ocean warming based on\nobserved and predicted ocean temperatures10,22, and these results indicate the utility of concurrently\ntracking ocean conditions such as oxygen and productivity for a more nuanced understanding of possible\nspecies impacts.”\n\nResponse Reviewer 3, Comment 5. Thank you, we have indicated that metrics are calculated from Aug-\nOct in the caption of Fig. 2, which displays the metrics discussed in this section:\n“All metrics were calculated from August-October in each MHW year relative to baseline conditions\n(2000-2020).”\n\nResponse Reviewer 3 Comment 6. This is a good idea, we have made this change.\n\nResponse Reviewer 3, Comment 7. Done:\n\nResponse Reviewer 2, Comment 1. Thank you - we have increased the size of this figure and improved\nthe resolution"}, {"round": 2, "reviewer_comments": {"reviewer_1": "I re-evaluated the manuscript after the first round of revisions. The authors did an excellent job\naddressing the reviewers’ comments and suggestions. I appreciate the effort put to conduct\nadditional tests, including sensitivity runs. In a few cases, when authors did not follow my\nsuggestions, they provided very reasonable rebuttals, which I agree with. Added or modified parts\nof the text improved the clarity and provided interesting perspectives to the discussion. The\nparallel to the hurricanes seems very useful and helps the reader understand the issue. The quality\nof the figures has been improved. Overall, after the last changes, the manuscript gained a lot of\nstrength. I am satisfied with all the changes made.\nMy only suggestion would be to consider mentioning the general results of the new analyses\npresented in the supplementary materials in the main body of the text. For example, it is reported\nin the supplementary materials that “Results indicate that species explains more of the variance\nthan MHW for percent change range, percent change area, and distance habitat metrics”, or\n“Ratios showed high alignment between daily observed and predicted presences (mean 1.5 +/- 1.3\nSD), indicating strong capacity to generate predictions in a near-real time framework”, but these\nresults are not mentioned in the main body of text. What is the general interpretation of the\nresults of the sensitivity test? You may consider including very brief information on each result in\nappropriate sections with reference to the supplementary materials.\nI don't have any more comments. If authors find my suggestion useful, very minor modifications\ncan be done to integrate summaries of the new results. Otherwise, I believe this manuscript can\nbe considered for publication in the journal.", "reviewer_3": "The authors have addressed all my comments – I recommend the manuscript for publication as is.\nI also want to commend and thank the authors for the thorough response to all reviewers; it was\nalso very helpful to have the changed/added text and figures within the document.\nThis study is an important contribution as to the impacts of marine heatwaves on top predators\nand highlights the need for effective multi-species ocean management as well as multi-variable\nanalyses and the consideration of the depth-extent for marine heatwave analyses.\nI again commend the authors for a very well written manuscript and clear figures to support the\npresented results."}, "author_response": "Response Reviewer #1. Thank you for your encouragement, and for improving the quality and clarity of\nour manuscript during two rounds of review. In following with the above suggestion, we now reference\nour analysis of daily predictive capacity in the discussion:\n“Our results indicate that species responses to MHWs are highly variable yet predictable: our models\nperformed well through extensive validation across space, time, and on novel data ((Supplementary Figs.\n6,7; Supplementary Table 3)”.\nWe have elected to not report results from the sensitivity analysis to data used in model fitting, but these\nsupplementary tables are referenced in the relevant figures (Fig. 2 and Fig. 4). We have also elected to not\ninclude the analysis of variance between vs among marine heatwaves in the main text, as we were not\nable to complete this analysis for direction. We feel there is value in maintaining a simple and\nstraightforward narrative in the main text to increase accessibility for a broad readership; while including\nmore detailed and nuanced analyses/results in the supplements for subject experts.\n\nResponse Reviewer #3. We appreciate your help improving the clarity and messaging of our manuscript."}]

Decision letter