Modelling the biogeographic boundary shift of Calanus finmarchicus reveals drivers of Arctic Atlantification by subarctic zooplankton

Biological communities in the Arctic are changing through the climate‐driven encroachment of subarctic species. This “Atlantification” extends to keystone Calanoid copepods, as the small‐bodied Calanus finmarchicus increases in abundance in areas where it overlaps with larger Arctic congeners. The environmental factors that are facilitating this shift, whether related to optimal conditions in temperature or seasonality, remain unclear. Assessing these drivers at an Arctic‐wide scale is necessary to predict future ecosystem change and impacts. Here we have compiled range‐wide occurrences of C. finmarchicus and a suite of seasonal biophysical climatologies to build a boreo‐Arctic ecological niche model. The data set was divided into two eras, 1955–1984 and 1985–2017, and an optimized MaxEnt model was used to predict the seasonal distribution of the abiotic niche of C. finmarchicus in both eras. Comparing outputs between eras reveals an increase in habitat suitability at the Arctic range edge. Large and significant increases in suitability are predicted in the regions of the Greenland, Labrador, and Southern Barents Seas that have experienced reduced sea‐ice cover. With the exception of the Barents Sea, these areas also show a seasonal shift in the timing of peak habitat suitability toward an earlier season. Our findings suggest that the Atlantification of Arctic zooplankton communities is accompanied by climate‐driven phenology changes. Although seasonality is a critical constraint to the establishment of C. finmarchicus at Arctic latitudes, earlier sea‐ice retreat and associated productivity is making these environments increasingly favorable for this subarctic species.

Vital to the Arctic food web are three coexisting species of Calanoid copepods; Calanus hyperboreus, C. glacialis, and C. finmarchicus. This Calanus complex dominates herbivorous mesozooplankton throughout the Arctic and northern seas (Mauchline, 1998).
Although they display differences in growth, development, and reproduction as imposed by different optimal environments, all three species are adapted to graze on the seasonal phytoplankton bloom, converting low-energy carbohydrates and proteins into high-energy wax esters (lipids), which they store to aid survival over-winter and fuel reproduction (Falk-Petersen et al., 2009). In doing so, they are a fundamental source of energy for higher trophic levels, sustaining vast fish stocks, seabird colonies, and marine mammal populations (Wassmann et al., 2006).
The Calanus life-cycle is complex. For C. finmarchicus at its northern range, it lasts 1-year. During spring and summer it develops from eggs via six naupliar and four copepodite stages to its major overwintering stage (CV), accumulates lipid reserves, and then descends to an overwintering depth below 500 m. Molting to the final adult stage and mating occur at the end of the winter. C. finmarchicus is generally an income breeder, thus while it accumulates enough lipid reserves to survive the winter, it relies on the nutritional input of the spring bloom upon resurfacing to fuel gonad maturation, egg production, and nauplii development (Hirche et al., 1997;Niehoff et al., 2002).
In the last 30 years, C. finmarchicus-a species typically associated with Atlantic water masses-has undergone a poleward distribution shift in the North Atlantic . This small, subarctic species has also increased in contribution to overall Calanus biomass in several Arctic regions (Aarflot et al., 2018;Hop et al., 2019;Moller & Nielsen, 2020;Weydmann et al., 2014). The fate of its larger congeners endemic to Arctic water masses is less clear due to a lack of long-term data from the high Arctic. There are indications of a decline in the southern margins of their distribution (Aarflot et al., 2018;Chust et al., 2014) while coastal populations of C. glacialis show stable population size (Hop et al., 2019;Moller & Nielsen, 2020). Overall, there is growing evidence that Atlantification is favoring a shift toward an Arctic Calanus community with smaller body size and less lipid content (Renaud et al., 2018). Shifts in the composition of Arctic Calanus and their functionally important traits (body size and lipid content) are of significant concern as a reduction in lipid production may impact the energy available to higher trophic levels (Kwasniewski et al., 2010;Renaud et al., 2018).
Any effect of C. finmarchicus' increasing presence is likely to be exacerbated under future climate warming. To be better equipped to predict this, we must move from local observations of encroachment to a regional understanding of the underlying mechanisms facilitating it. The thermal tolerance range of C. finmarchicus was found to be important in determining recent and future distribution shifts within the North Atlantic and subarctic (Beaugrand et al., 2008Chust et al., 2014;Helaouet & Beaugrand, 2007;Reygondeau & Beaugrand, 2011;Villarino et al., 2015), yet it is the length of the growing season, that is, the period of phytoplankton availability, which is thought to be a critical limitation in successful recruitment of C. finmarchicus at their Arctic range edge (Hirche & Kosobokova, 2007;Ji et al., 2012). Here, the timing of the phytoplankton bloom is strongly dictated by local light conditions related to seasonal sea-ice cover (Falk-Petersen et al., 2009). Recent studies have highlighted a link between the biogeography of C. glacialis and C. hyperboreus and changing sea-ice characteristics (Ershova et al., 2021). If the biogeographic boundary of C. finmarchicus is also limited more by seasonality and food availability than by temperature per se, areas which have increased in suitability for C. finmarchicus in recent decades would have also experienced the greatest change in seasonality due to reduced sea-ice cover.
Ecological niche models (ENMs) and allied species distribution models are valuable tools for understanding environmental correlates of Calanus biogeography (Albouy-Boyer et al., 2016;Beaugrand et al., 2013;Helaouet & Beaugrand, 2007;Helaouet et al., 2011;Record et al., 2018;Villarino et al., 2015;Wilson et al., 2016). These have largely taken advantage of long-term monitoring data sets within the North Atlantic yet often lack the seasonal resolution necessary to uncover seasonal as well as decadal shifts in suitable habitat. Here, we extend these efforts by compiling rangewide collections of C. finmarchicus and a suite of seasonal biophysical climatologies to build the first boreo-Arctic ENM for this species.
We assess how the distribution of their ecological niche, particularly at the poleward boundary, has shifted seasonally between two eras (1955-1984 and 1985-2017), which correspond to cool and warm thermal regimes in the region, respectively (Beaugrand, 2009).
We aim to (a) characterize the spatial and temporal patterns of C. finmarchicus distribution on a pan-Arctic scale and (b) determine the abiotic drivers, including the relative importance of thermal tolerance limits and seasonality, in facilitating the Atlantification of zooplankton communities in the Arctic. After identifying overlapping data sets between repositories, remaining records were thinned to retain only one occurrence per season per grid cell (resolution: 0.25° × 0.25°) using the "spThin" R package (Aiello-Lammens et al., 2015). This removes the fewest records necessary to substantially reduce the effects of sampling bias, while simultaneously retaining the greatest amount of useful information.

| Environmental predictors
Ten environmental predictors were identified as candidate variables for the niche model. These include: temperature, bathymetry, slope, chlorophyll a, sea-ice concentration, salinity, silicate, pH, photosynthetically active radiation (PAR), and current velocity. Representing the sea surface only, these data were obtained for the full study region (longitude: −180°E, 180°W; latitude: 30°N-90°N) using a combination of empirical observations and model re-analyses (see Table S1 in File S1 for details on all data sources). As C. finmarchicus is absent from the Pacific Arctic, these regions were removed from environmental raster data to limit model outputs to the known spatial distribution range of the species. Temperature and PAR were the only variables to be highly correlated (Pearson's r = .71, Figure S1).
Methods used to deal with correlation are described in the section "model evaluation and tuning." Occurrence data were matched to the most appropriate environmental data in relation to the season and year in which they were collected by adapting the method from Duffy and Chown (2017).
For example, a record collected in May 2002 was assigned the environmental conditions from the AMJ 1985-2017 climatology. This accounts for variation in environmental conditions (particularly sea-ice extent and primary productivity) between seasons and multidecadal time periods, and, as only the most appropriate environmental data are paired with each occurrence record, the accuracy of the ENM can, in principle, be improved (Duffy & Chown, 2017).
MaxEnt estimates the conditional probability of the presence of a species relative to locations where the species has been observed by sampling the environment at a range of "background" locations across the study region and discriminating these from locations where species is known to be present. MaxEnt assumes background locations adequately cover areas accessible to the species and that the presence of localities are unbiased and cover important environmental gradients (Jarnevich et al., 2015). Although a lack of absence data prevents probability estimates of a species presence and predictions of a species' realized distribution, presence-only outputs more closely represent the existing, fundamental niche of a species (Soberon & Nakamura, 2009).

| Sampling bias and background data selection
Sampling in this region is skewed in favor of sea-ice free areas, leading to spatial and environmental sampling bias. To avoid generating distribution maps which overly reflect sampling effort (Botella et al., 2020), the geographic and temporal distribution of background data were selected based on an estimate of zooplankton sampling effort across the region (Fourcade et al., 2014;Phillips et al., 2009). All occurrence records within the phylum "arthropoda" occurring in the upper 200 m and collected between 1960 and 2017 were downloaded from OBIS.

F I G U R E 1
Annual mean sea surface temperature and sea-ice extent within the North Atlantic and European Arctic (−75°W, 75°E, 30°N, 90°N) between 1955 and 2017. Dashed vertical line denotes transition between the two eras compared within the study (1955-1984 and 1985-2017). Data plotted from Hadley Centre Sea Ice and Sea Surface Temperature data set (Rayner et al., 2003) and A kernel density surface of these records was generated ( Figure S2) and used to weight the selection of 10,000 random points across the region (i.e., more points taken from areas with higher-density values; Figure S2). These points were also randomly assigned to a season and era proportionate to the temporal distribution of zooplankton data.
Background points were matched to the most appropriate environmental data as described for the occurrence records.

| Model evaluation and tuning
Cross-validation approaches partition the data in to K number of folds. The model is then run K times, withholding a different fold for model evaluation each time (Araujo et al., 2019). The spatialBlock function of the "blockCV" package in R (Valavi et al., 2019) was used to create folds (K = 5) that account for spatial autocorrelation in the environmental data. Spatial block size was determined by fitting isotropic variogram models using 5000 random points from each environmental predictor raster. This finds the effective range of spatial autocorrelation and the spatial block size was based on median of these ranges. The occurrence and background data within each block were then allocated to a fold. To ensure block-to-fold allocations achieved an even spread of data, 100 iterations were ran with the most even allocation of data being used. See Figure S3 for map of spatial blocks and fold assignment.
The area under the receiving operator curve (AUC) and true skill statistic (TSS) metrics were used to evaluate model discriminatory performance on the evaluation (test) fold. The AUC score is a widely used rank-based measure of predictive accuracy that can be interpreted in the context of MaxEnt as the probability that a randomly chosen presence location is ranked higher than a randomly chosen background point (Merow et al., 2013). A model with no discriminatory power will have an AUC value equal to 0.5 (no better than random), whereas a model with perfect fit would have an AUC value of 1.0. TSS values range from −1 to 1, with values of 0 or less reflecting a model that is no better than random and values closer to 1 being better at discerning presence and background points (Allouche et al., 2006).
To prevent model overfitting, the gridsearch function of "SDMtune" package was used to find the optimal combination of MaxEnt hyperparameters. Varying combinations of regularization parameter (0.2-3) and iteration number parameter (300-900) were tested, and the combination with the highest AUC TEST and TSS TEST scores was considered optimal. Only linear and quadratic feature class settings were used in all combinations. Additional, more complex transformations such as the "hinge" feature class did not improve model performance and generated less realistic response curves to environmental parameters.
The model with optimal hyperparameters was further tuned using (1) the varSel function to remove any correlated variables, removing the one that results in the best performing model when removed and (2) the reduceVar function to find and remove environmental predictors with low model contribution (<3% permutation importance) when their removal did not decrease model performance based upon the model's mean AUC TEST value.

| Model prediction
The final, optimized model was used to predict the habitat suitability of C. finmarchicus across the region of interest. Separate predictions were made for each season (JFM, AMJ, JAS, OND) and era . For each prediction, outputs from each crossvalidation fold as well as the mean across all folds were retained to explore between-fold variation.

| Model output analyses
To assess changes in the spatial pattern of habitat suitability between the two eras, we subtracted model outputs for the most recent era from the former. To assess seasonal changes between the eras, we used ArcGIS Cell Statistics tool to determine, for both eras, which season gave the highest habitat suitability value in each grid cell. Subtracting these two outputs from each other allowed an assessment of where timings in peak suitability had shifted between eras. To assess the significance of habitat suitability change, the 95% confidence interval (CI) was calculated for each grid cell and era based on variation from the cross-validation folds. Grid cells that had no overlap in CI range between eras were classed as being significantly different.

| Sensitivity analyses
Sensitivity analyses were carried out to account for potential misidentification of occurrences between C. finmarchicus and its congener, C. glacialis (Choquet et al., 2017). In regions where their ranges overlap, a random subset (10% and 20%) of C. finmarchicus occurrences were replaced with a corresponding number of C. glacialis records. These models are named 10%_c.glacialis and 20%_c.glacialis, respectively, and Figure S4 shows localities of dropped/replaced occurrences in each model. To check that the seasonal assignment process of occurrences to environmental data did not influence model results, we carried out a third sensitivity analysis (model name: 30%_season_shift) whereby 50% of records collected in months 1, 4, 7, and 10 (first month of each climatology) were re-assigned to one season earlier. Similarly, 50% of records collected in months 3, 6, 9, and 12 (last month of each climatology) were re-assigned to one season later. This led to 34% of all occurrence records being assigned to a different season and the remaining 66% of occurrences remained unchanged. A final analysis (model name: original_no_ice) was carried out to test the sensitivity of model outputs to the inclusion of seaice concentration. These ENMs were optimized following the same procedure as for the original model. Further details of these analyses are given in Table S2.

| Spatial change in habitat suitability
The predicted geographic distribution of C. finmarchicus for each season and era are shown in Figure 2 with the variation between cross-validation folds given in Figure S5. Habitat suitability is predicted to be highest in the East Atlantic during Jan-Feb-Mar, whilst it is predicted to be highest during JAS at high latitudes including the Barents Sea (Figure 2).
Across its Arctic range edge, the model predicts there to have been large increases in C. finmarchicus habitat suitability between eras. This has predominantly occurred where sea-ice cover has declined (Figure 3). In these regions of increased suitability, the 95% CIs in predictions do not overlap between eras, suggesting that the predicted increase in suitability is significant ( Figure   S6). Regions and levels of increased suitability are consistent in the outputs from sensitivity analyses 10%_c.glacialis, 20%_c.glacialis, and 30%_season_shift ( Figure 4). When sea-ice concentration is removed (original_no_ice), predicted increases are present throughout the Arctic range edge of C. finmarchicus, but values are lower ( Figure S8). In all but the last sensitivity analysis, the vast majority of cells predicted to have increased in suitability remain significant ( Figure S6).

| Temporal change in habitat suitability
We find a seasonal shift in the timing of peak habitat suitability, most notably in the Labrador Sea, northern North Atlantic, and Greenland Sea regions ( Figure 5). Within these areas, the timing of peak habitat suitability has advanced forward by one season, that is, changing from being optimal during JAS in the former era to being optimal in AMJ more recently. All cells predicted to have advanced in timing of peak suitability were found to be significant, based on a lack of 95% CI overlap between eras. These regions coincide with areas of retreating sea-ice cover, and this result is consistent across sensitivity analyses with the exception of model original_no_ice ( Figure S7).

| Model performance and predictor importance
The optimized MaxEnt model retained five environmental predictors; sea-ice concentration, temperature, chlorophyll a, bathymetry, and salinity. Model performance metrics indicate strong discrimination ability with mean AUC TEST = 0.73, and TSS TEST = 0.41. On average, sea-ice concentration had the highest permutation importance at 45.5%, followed by temperature, salinity, chlorophyll a, and bathymetry (Table 1) ( Figures S9-S13) show a strong negative relationship between C. finmarchicus habitat suitability and sea-ice concentration ( Figure   S9) and indicate an optimal temperature of 8.9°C ( Figure S10).
Sea-ice concentration and temperature ranges within the upper quartile of suitable habitat values are 0%-24.2% and 4.9-12.7°C, respectively. F I G U R E 3 Regions predicted to have increased in habitat suitability (>0.1) for Calanus finmarchicus between eras (1985-2017 and 1955-1984) and for each season (a-d). Green and black lines denote the mean position of the sea-ice edge for the older era and the recent era, respectively [Colour figure can be viewed at wileyonlinelibrary.com] F I G U R E 4 From left to right, predictions of increased habitat suitability (>0.1) between eras (1985-2017 and 1955-1984) using (a) the original model and (b-e) four different sensitivity analyses (see Section 2 for details). Green and black lines denote the mean position of the sea-ice edge for the older era and the recent era, respectively [Colour figure can be viewed at wileyonlinelibrary.com] Sensitivity analyses to account for potential misidentification of C. finmarchicus and C. glacialis (Choquet et al., 2017) resulted in minimal change to AUC and TSS performance metrics, a small decline in the contribution of sea-ice concentration (Table 1), and a less severe negative relationship between habitat suitability and sea-ice concentration ( Figure S14). Sensitivity analyses to account for seasonal assignment of occurrence records resulted in no change to model performance or variable contribution. A further test that withheld the inclusion of sea-ice concentration (model: original_no_ice) was found to have decreased model performance and increased importance of salinity and temperature predictors (Table 1).

| DISCUSS ION
Large-scale changes in the abundance and distribution of marine species are omnipresent and consistent with ocean warming over the last century (Hastings et al., 2020). Robust to the sensitivity analyses tested, our results reveal that suitable habitat for C. finmarchicus has increased at Arctic latitudes in the last 30 years, extending the previously known range shift within the North Atlantic . This is consistent with regional observations of boreal plankton and benthic species becoming more dominant within Arctic ecosystems (Aarflot et al., 2018;Dalpadado et al., 2020;Fossheim et al., 2015;Kortsch et al., 2012;Moller & Nielsen, 2020;Polyakov et al., 2020).
Predictions of "suitable habitat" from presence-only ENMs can be interpreted as showing the potential distribution of a species. This is typically broader than their realized distribution because the model does not incorporate absence records and does not account for biological interactions such as competition which further constrain where a species can persist (Soberon & Nakamura, 2009). In this study, results show a similar pattern to the core distribution of C. finmarchicus described by Choquet et al. (2017). Known zones of expatriation, where the species is found through advection but cannot successfully complete its life-cycle, such as in the Arctic Ocean basin, were not predicted to be highly suitable. Thus, our outputs represent regions where surface conditions are suitable for copepodite survival and where population recruitment may occur locally. Moreover, our model predicts that C. finmarchicus habitat is characterized by optimal surface temperatures between 4 and 12°C, with a peak at 9°C.
These are consistent with previous regional model estimates (Albouy-Boyer et al., 2016;Beaugrand et al., 2013;Helaouet & Beaugrand, 2007, 2009) and observations (Bonnet et al., 2005;Strand et al., 2020) for this species. Although temperature has a strong influence on the biogeography of C. finmarchicus, our findings suggest that other factors-in addition to temperature-may have influenced the opening up of suitable habitat at their Arctic range edge, as we detail below.
We find a strong overlap between regions of sea-ice retreat and regions predicted to have undergone: (a) an increase in suitable habitat and, (b) a seasonal advancement in suitable habitat.
Although this may, in part, be influenced by greater sampling effort in low sea-ice conditions, the importance of sea-ice parameters in determining the biogeography of other Arctic Calanus (C. glacialis and C. hyperboreus) has recently been demonstrated (Ershova et al., 2021;Feng et al., 2016Feng et al., , 2018. Our findings suggest that F I G U R E 5 Solid filled areas represent the season containing highest suitability value for Calanus finmarchicus during era 1 . Hatching denotes areas where the highest suitability value advanced forward by one season during era 2 (1985-2017). White and black lines denote the mean position of the sea-ice edge for the older era and the recent era, respectively [Colour figure can be viewed at wileyonlinelibrary.com] TA B L E 1 Percent contribution (%) of each environmental variable to model performance (fivefold mean ± 1 SD). Model names are Original = final Calanus finmarchicus model; 10% and 20%_c.glacialis = replaced 10% and 20% of occurrences with C. glacialis; 30%_season_ shift = reassigned 30% of occurrences to earlier/later seasonal climatology; original_no_ice = original model without sea-ice concentration as a variable phenological changes caused by the retreating ice-edge may also be an important driver of Arctic Atlantification of zooplankton: with sea-ice loss, the seasonal conditions necessary for C. finmarchicus to succeed at Arctic latitudes have started to emerge. This gives cause and context to recent empirical observations as areas of earlier and/or increased suitability from this model correspond with Atlantic-Arctic gateway areas known to have experienced biomass increases of C. finmarchicus including the Barents Sea (Aarflot et al., 2018), Disko Bay (Moller & Nielsen, 2020), and the Fram Strait (Weydmann et al., 2014).
Our study is limited to inferring indirectly the link between sea-ice concentration and suitable seasonality for C. finmarchicus. However, this premise is supported by studies showing a strong correlation between decreased summer sea-ice concentration and an earlier Arctic phytoplankton bloom (Kahru et al., 2011;Song et al., 2021), the peak of which has advanced by up to 50 days in the Baffin Sea (Kahru et al., 2011) and over a month in areas of the Barents Sea (Dalpadado et al., 2020). Studies from the Bering Sea have also found zooplankton community shifts linked to combined effects of temperature and primary production between ice-covered and ice-free years (Kimmel et al., 2018).
Although C. finmarchicus does also use some capital resources, it is predominantly an income breeder, that is, requiring food provided by the spring bloom to contribute to facilitate maturation and reproduction, and is, thus, less equipped to deal with high inter-annual variability in the bloom phenology (Falk-Petersen et al., 2009). Its inability to cope with a short growing season in addition to low temperatures is regarded as the main limiting factors to allow it to survive and reproduce in the polar basin (Hirche & Kosobokova, 2007;Ji et al., 2012). Thus, as long as environmental conditions such as timing of ice break up and onset of the spring bloom remain highly variable in the areas defined as newly suitable, C. finmarchicus may struggle to reproduce and survive on a year-to-year basis, and successful establishment will be dependent on the constant replenishment of the population from the south. habitat earlier in the year, there is a limit as to how far these factors will aid the northward extension of C. finmarchicus (Ji et al., 2012;Slagstad et al., 2011). Even if the Arctic Ocean becomes ice free, the light climate at extreme high latitudes will limit primary production for long periods each year, and high inter-annual variability in sea-ice extent leads to unpredictable bloom phenology, providing conditions in which the Arctic congener species are still better adapted than C. finmarchicus (Daase et al., 2013;Falk-Petersen et al., 2009). In areas of co-existence such as the Barents Sea and in the fjords of Svalbard and Greenland, there is no conclusive evidence that an increase in one Calanus species is directly detrimental to another (Hop et al., 2019;Moller & Nielsen, 2020), although the degree to which their respective niches overlap remains poorly understood. Furthermore, difficulties in distinguishing between C. finmarchicus and C. glacialis morphologically (Choquet et al., 2018), and a lack of long-term data with reliable species identification, limits our ability to assess any changes to their co-occurrence.
Understanding the adaptability and resilience of each Calanus species will require life-stage specific estimates of suitable habitat and fitness under a range of Arctic conditions. Although the correlative model applied here is the first step in examining broad-scale patterns of change, potential sources of bias remain (e.g., geographic sampling bias) because of its data-driven approach. Mechanistic approaches that estimate an organism's energetic budget at a fine temporal resolution and that incorporate the Arctic's extreme light environment, will have an important, complementary role in predicting the success of subarctic species at Arctic latitudes (Ljungstrom et al., 2021). Novel observation methods (Vilgrain et al., 2021) and increased sampling within sea-ice environments will also be important in overcoming data limitations in the Arctic and for parameterizing and validating model outputs.
Results from this study suggest that the seasonal conditions necessary for C. finmarchicus to survive have emerged at their Arctic range edge in recent decades. In these Arctic gateway regions, encroachment of C. finmarchicus is likely to alter the overturning and availability of energy in the pelagic ecosystem due to their smaller size, lower lipid content, and shorter life-cycle durations compared to Arctic congeners. An unprecedented warm and ice-free year in the Bering Sea saw an increase in small, low-lipid zooplankton, concurrent poor catches of pelagic fish, and low reproductive success and mass mortality at seabird colonies (Duffy-Anderson et al., 2019).
Yet trait-based models have demonstrated that a C. finmarchicus-like life-history also brings a shorter generation time and faster population turnover, which may compensate or even enhance the transfer of energy to predators (Renaud et al., 2018). The ecological implications of changes in the Calanus complex remain uncertain. These will depend on the dynamics of sea-ice decline and associated phenology shifts, as well as the adaptability of, and interactions between, Arctic pelagic species.

ACK N OWLED G M ENTS
We thank David Pond as DIAPOD principal investigator. JJF Deep Impact project (300333).

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
A list of the environmental data used within this study is provided in Table S1 within File S1. A list of the occurrence record data sets used within this study is provided in File S2. The data and R code that support the findings of this study are available from the UK Polar Data Centre at https://doi.org/10.5285/FC660 BC3-09AB-4C1A-9D2A-40269