Can a key boreal Calanus copepod species now complete its life-cycle in the Arctic? Evidence and implications for Arctic food-webs

The changing Arctic environment is affecting zooplankton that support its abundant wildlife. We examined how these changes are influencing a key zooplankton species, Calanus finmarchicus, principally found in the North Atlantic but expatriated to the Arctic. Close to the ice-edge in the Fram Strait, we identified areas that, since the 1980s, are increasingly favourable to C. finmarchicus. Field-sampling revealed part of the population there to be capable of amassing enough reserves to overwinter. Early developmental stages were also present in early summer, suggesting successful local recruitment. This extension to suitable C. finmarchicus habitat is most likely facilitated by the long-term retreat of the ice-edge, allowing phytoplankton to bloom earlier and for longer and through higher temperatures increasing copepod developmental rates. The increased capacity for this species to complete its life-cycle and prosper in the Fram Strait can change community structure, with large consequences to regional food-webs. Supplementary Information The online version contains supplementary material available at 10.1007/s13280-021-01667-y.

: Temperature-Salinity curves for pairs of stations sampled in 2018 (red) and 2019 (black) within the Nordic seas at a) S1, b) S2 and c) S3. Water mass boundaries are as described in Rudels et al (2005) with Atlantic Water (AW), Arctic Deep Waters (ArDW), Polar Surface Water (PSW and warm Polar Surface Water (PSWw) noted.

Occurrence data
Six online repositories (OBIS, PANGAEA, NSF Arctic Data Center, BODC, COPEPOD global plankton database, NOAA NODC), comprising of 46 individual datasets, were used to compile 65,037 georeferenced occurrence records of Calanus finmarchicus lifestages CIadult. Duplicate records were removed, and the remaining records were spatially thinned to remove the fewest records necessary to substantially reduce the effects of sampling bias, while simultaneously retaining the greatest amount of useful information (Aiello-Lammens et al. 2015). Thinning distance was equal to environmental data resolution (0.25° x 0.25°) and the final number of records used in analyses was 27,726 ( Figure S2a).

Environmental data
Ten sea-surface environmental predictors were identified as candidate variables for the niche model (Table S1). For each predictor, seasonal climatologies from 30° -90°N were obtained for two eras of approximately 30 years; 1955-1984 and 1985-2017. Based on data availability, the seasonal partitions represent the months; Jan-Feb-Mar (JFM), Apr-May-Jun (AMJ), Jul-Aug-Sep (JAS), and Oct-Nov-Dec (OND). Occurrence data were then matched to environmental predictor data from the corresponding season and time period in which they were collected from by adapting the method and code from (Duffy & Chown 2017).

MaxEnt niche model
Data were fitted to the presence-only ecological niche modelling algorithm MaxEnt v. 3.4.1 (Phillips & Dudik 2008, Phillips et al. 2017) using the SDMtune R package (Vignali et al. 2020). MaxEnt estimates the conditional probability of 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 environments at locations the species is known to occupy.
A kernel density surface of zooplankton data was 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 S2b). Zooplankton data were taken as all OBIS-stored occurrence records of the phylum "arthropoda" from the upper 200m collected between 1960 and 2017 in the region of interest. These points were randomly assigned to a season and time period which was proportionate to the temporal distribution of zooplankton data. As with occurrence records, background points were subsequently matched to corresponding environmental data.
The spatialBlock function of the blockCV package in R (Valavi et al. 2019) was used to create fold partitions (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 assigned to a fold. The MaxEnt model is then run several times, withholding a different fold for evaluation each time.
To prevent model overfitting, several tuning procedures were followed and 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. First, the gridsearch function of SDMtune package was used to find optimal combination of MaxEnt hyperparameters. Models were run with varying combinations of regularisation parameter (0.2 -3) and iteration number parameter (300 -900). Feature class settings did not vary with only linear and quadratic transformations allowed. Average AUCTEST and TSSTEST metrics were stored for each parameter combination and the one with highest values was considered optimal. Secondly, the varSel function was used to remove any correlated environmental predictors, choosing to remove whichever reduced model permutation importance the least. Thirdly, the reduceVar function was used to find and remove environmental predictors with low model contribution (<3% permutation importance).

Predicting habitat change
The optimised MaxEnt 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 . To assess changes in the spatial pattern of habitat suitability between the two eras, we subtracted model predictions of habitat suitability for the most recent era from the former and repeated this calculation for each season. Here, we focus on the change in habitat suitability for the early productive season (AMJ) and late productive season (JAS) clipped to the Fram Strait region.

2.5.
Results summary The final model, optimised in terms of MaxEnt hyperparameters and environmental predictors, retained five predictors (in order of permutation importance: sea-ice concentration, temperature, salinity, chlorophyll a, bathymetry) and had good discriminatory ability (mean AUCTEST = 0.72, mean TSSTEST = 0.41). The model predicts that C. finmarchicus habitat is characterized by optimal surface temperatures between 4 and 12°C, with a peak at 9°C. This corresponds to a geographic distribution spanning 40-80°N, with high suitability across the North Atlantic and Norwegian Sea, and southern and northern range edges that vary seasonally ( Fig. S3a-b). This biogeography is consistent with previous model findings (Helaouet & Beaugrand 2007, Helaouet & Beaugrand 2009, Beaugrand et al. 2013, Albouy-Boyer et al. 2016) and observations (Bonnet et al. 2005, Strand et al. 2020. Moreover, our results demonstrate that suitable habitat for C. finmarchicus has increased at Arctic latitudes in the last 30 years while temperate regions have declined in suitability ( Fig  S3c-d). This is in line with observations of boreal species entering, and becoming more dominant within, Arctic ecosystems (Kortsch et al. 2012, Fossheim et al. 2015, Aarflot et al. 2018, Dalpadado et al. 2020, Moller & Nielsen 2020, Polyakov et al. 2020. Full assessment of model results are within Freer et al. (In press).