Range edge analysis
We quantified species range edges as the 0.01 and 0.99 quantiles of density along spatial axes running the length of each study region. We chose these quantiles to capture the extremes of each species’ distribution; because edges were calculated from VAST’s spatiotemporal biomass estimates and not from the raw data, they were unlikely to be influenced by the rare, high biomass observations that are common in shoaling species such as fishes (Thorson et al. 2011). Species range edges are conventionally measured in units of degrees latitude along a north-south axis (e.g., Hickling et al. 2005). However, in marine regions with complex coastline topographies and/or coastlines that are not oriented parallel to lines of longitude, coastal distance is a more accurate metric of range edge position than latitude (Hare et al. 2010, Bell et al. 2015, Fredston‐Hermann et al. 2020). We therefore developed a coastal distance metric for the West Coast and the Northeast (Figure 1A, 1B; methods in Appendix 2). We then associated points along the coastline with the grid of VAST knots (see Appendix 2) by finding the points with the minimum Euclidean distances. Matching points along the coastline to the VAST knots enabled us to estimate density, and thus range edge position, along the coastal distance axis. In the Eastern Bering Sea, the coastal distance axis was less applicable because the shelf is so wide that many species fall quite far from the coast and also because the presence of islands makes the delineation of a smoothed coastline more complex. In the Eastern Bering Sea, we therefore estimated density along a line drawn from the Aleutian Islands (56ºN, 161ºW) to the edge of the US Exclusive Economic Zone (62ºN, 176.5ºW) through the Middle Domain, a hydrographic region with similar bathymetry defined as lying between two oceanographic mixing zones that partition the middle from inner and outer domains (Coachman 1986, Ortiz et al. 2016).
To ensure that the species analyzed had at least one range edge in the study region, we eliminated range edges with mean positions over time that fell close to the boundary of the study region, defined as less than 10% from the end of the axis of measurement in either direction (Northeast axis length = 1368 km, West Coast axis length = 2037 km, Eastern Bering Sea axis length = 1102 km). This removed 18 Northeast species, 24 West Coast species, and 13 Eastern Bering Sea species. After this filter, we proceeded with 153 range edges—55 in the Northeast, 25 on the West Coast, and 73 in the Eastern Bering Sea—across 145 fishes and invertebrates from 17 taxonomic classes (eight had both equatorward and poleward range edges; see Appendix 3). We note that for almost all species, only one range edge fell within the study region (see Appendix 3). Thus, our analysis evaluated the evidence for our different hypotheses by evaluating many isolated range edges, not by evaluating both range edges of a single species.
We tested whether range edge positions had shifted significantly over time with single-species Bayesian linear regressions of range edge position on time (n = 151 models; two models were removed that did not estimate standard errors around range edge position due to parameter reduction to aid convergence in VAST). Single-species models were fitted using the “rstanarm” package (Goodrich et al. 2018) with four chains, 40,000 iterations including 10,000 burn-in draws, and a target average proposal acceptance probability of 0.99. We selected a normally distributed vague prior with a mean of 0 and standard deviation of 50 km/year; this standard deviation was chosen to exceed the upper range of climate velocities in the oceans (Burrows et al. 2011). Range positions were weighted by VAST-estimated standard errors (\(\frac{1}{\text{SE}^{2}}\)) for each year so that estimated edge positions with higher associated uncertainties were less influential (Thorson et al. 2016b). All of these models converged, as evaluated by Gelman-Rubin convergence statistic below 1.1.