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.