ARTICLE

Methods, Tools, and Technologies

# An ecological niche model that considers local relationships among variables: The Environmental String Model

Grégory Beaugrand

CNRS, University of Lille, Université du Littoral Côte d'Opale, UMR 8187, LOG, Laboratoire d'Océanologie et de Géosciences, Wimereux, France

**Correspondence**

Grégory Beaugrand

Email: gregory.beaugrand@univ-lille.fr

**Funding information**

French National Research Agency,
Grant/Award Number:
ANR-21-EXES-0011

Handling Editor: Samniqueka J. Halsey

**Abstract**

Many methods have been proposed to model the spatial distribution of a spe-
cies. While some methods have been specifically designed for this purpose,
others are well-known statistical tools that can be used in many scientific
fields. In this paper, I propose a new ecological niche model, called the
Environmental String Model (ESM), that is based on the concept of environ-
mental string, which is defined as being a combination of environmental vari-
ables, with as many nodes as environmental variables. There are two types
of environmental strings: (1) the abundance-known string and (2) the
abundance-unknown string (or target string) for which an estimation of abun-
dance is searched. The novelty of the model is that it assesses the abundance
associated with a target string from nearby abundance-known strings, which
preserve the local multidimensional relationships with the target string. The
model does not provide an abundance estimate in the absence of data from a
similar environment and it can therefore deal with truncated spatial distribu-
tions or niches. It is tested in the North Atlantic Ocean on two key copepod
species, *Calanus finmarchicus* and *Calanus helgolandicus*, which have been
monitored by the Continuous Plankton Recorder (CPR) survey for decades.
I investigate the influence of variables on model performance. I show that the
model reconstructs the mean spatial distribution and seasonal fluctuations in
both *Calanus* well. When compared with generalized linear models (GLMs),
generalized additive models (GAMs) and generalized regression neural net-
work (GRNN), the ESM gives the best performance. I propose a number of
indicators to evaluate the robustness of estimated abundance in space and
time and show how the model may be extended to presence/absence or
presence-only data. I think that the ESM could be used to fill gaps in any
sampling program such as the CPR survey and many satellite databases
(e.g., ocean color and photosynthetically active radiation).

**KEYWORDS**

Continuous Plankton Recorder, ecological theory, model performance, modeling, niche

This is an open access article under the terms of the [Creative Commons Attribution](https://creativecommons.org/licenses/by) License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

© 2024 The Author(s). *Ecosphere* published by Wiley Periodicals LLC on behalf of The Ecological Society of America.
# INTRODUCTION

Many models can be used to predict or project the spatial distribution of a species (Helaouët & Beaugrand, 2009; Hirzel et al., 2002; Melo-Merinoa et al., 2020; Philips et al., 2006; Schickele et al., 2020; Thuiller et al., 2009; Valavi et al., 2021). According to the types of data on which such models can be applied, three families can be distinguished: (1) models that work on presence-only data such as the ecological niche factor analysis (ENFA; Hirzel et al., 2002), MaxEnt (Philips et al., 2006), and the nonparametric probabilistic ecological niche (NPPEN) model (Beaugrand et al., 2011); (2) models that work on presence/absence data such as generalized linear models (GLMs) (McCullagh & Nelder, 1983), generalized additive models (GAMs) (Hastie & Tibshirani, 1990), and random forests (Zhang et al., 2019); and (3) models working on abundance data such as generalized regression neural network (GRNN) (Specht, 1991), GLMs, and GAMs. Because many datasets have only occurrence data, the pre-requisite to the application of many models working on presence/absence data is to generate pseudo-absence (Barbet-Massin et al., 2012; Valavi et al., 2021). Some models have been specifically designed for species distribution modeling such as MaxEnt, ENFA, and NPPEN, whereas others such as GLMs, GAMs (Guisan et al., 2002; Hastie & Tibshirani, 1990), quantile regression (Cade et al., 2005; Cade & Noon, 2003; Koenker & Bassett, 1978; Vaz et al., 2008), flexible discriminant analysis (FDA) (Hastie et al., 1994), or GRNN have been applied in many scientific fields. Some are primarily species distribution models (SDMs) such as MaxEnt, whereas others are ecological niche models (ENMs) such as NPPEN. Although SDMs and ENMs are frequently considered to be synonymous terms, their meaning is not quite the same because SDMs do not always rely on the concept of the ecological niche in contrast to ENMs (Peterson & Soberon, 2012).

ENMs are based on the concept of ecological niche sensu Hutchinson (1957, 1978). The Hutchinsonian niche represents the combination of environmental tolerances and resources required by an organism to grow and reproduce. Hutchinson (1957) conceptualized this notion with the so-called *n*-dimensional hypervolume, in which *n* ideally corresponds to all environmental (biotic and abiotic) factors. This way to define the niche was considered to be revolutionary because it was operational, enabling a straightforward quantification of a species' niche (Kléparski & Beaugrand, 2022). The concept of the ecological niche remains discussed in the scientific literature. In the context of ecological niche modeling, discussions have taken place on the niche concept itself (Chase & Leibold, 2003; Hirzel & Le Lay, 2008; Holt, 2009), the relationships between fundamental and

realized niches (Bruno et al., 2003; Pulliam, 1988, 2000), the niche-distribution duality (Colwell & Rangel, 2009; Pulliam, 2000), organism mobility and accessibility (Soberon & Peterson, 2005), and niche conservatism (Crisp et al., 2009; Peterson, 2011; Peterson et al., 1999). Another issue, perhaps less discussed in the scientific literature (Frontier & Pichot-Viale, 1993), is the difficulty of considering the dependency among environmental variables in the models. This dependency is frequently considered globally by means of linear or nonlinear correlation or regression coefficients, the use of the generalized Mahalanobis distance (e.g., NPPEN), or pair interaction in GLMs and GAMs (Beaugrand et al., 2011; Durant et al., 2005; Hastie & Tibshirani, 1990; McCullagh & Nelder, 1983).

Here I propose a simple comprehensive ENM, based on the concept of environmental strings. In this ENM, the niche is envisioned as a combination of strings, defined as being a combination of interrelated environmental variables associated with a certain level of species abundance, nil or positive. These environmental strings are used to estimate species abundance in space and time, providing that they are sufficiently close and numerous to the environment in which a value of abundance is searched (see the precise definition of a string in the next section). The assessment is performed by applying a weighted average using a negative exponential function. The concept of environmental string is interesting because it considered locally the relationships among variables, which is only accounted for globally in many approaches (Beaugrand et al., 2011). I first explain the Environmental String Model (ESM) using a simple graphical example. I then apply the model to two key zooplankton species sampled by the Continuous Plankton Recorder (CPR) survey in the North Atlantic Ocean and its adjacent seas. I propose a series of indices that helps understanding when and where modeled abundance may be less reliable. I compare the results with other approaches (i.e., GLM, GAM, and GRNN). I show how to identify the most important variables in the model and finally demonstrate that the method can be used on presence/absence data and probably on presence-only data after generating pseudo-absence.

# MATERIALS AND METHODS

## The model

### Rationale

My model is an ENM based on the concept of environmental string (Figure 1). I defined an environmental string as a combination of environmental variables where
<figure><img src="image_1.png" /><figcaption>**FIGURE 1** Schematics showing how the model works. An environmental string is a combination of standardized environmental variables associated (or not) with a value of abundance (or alternatively any measure of concentration, e.g., chlorophyll concentration). TES, One abundance-unknown environmental string or target environmental string (TES) and five abundance-known environmental strings (ESs) are considered in this example. The node of a string corresponds to the value of an environmental variable (standardized between 0 and 1) corresponding to one environmental dimension (or variable). The gray zone represents the target zone (TZ) where each node of an ES must be in order to take part of the calculation of the abundance associated with the TES. The average distance between the node of an ES and the TES is calculated for all ecological variables. The average distance of each ES within TZ is then used to estimate the abundance of the TES.</figcaption></figure>

all the values are interdependent; each node of a string corresponds to the value of an environmental variable (all variables being standardized between 0 and 1). The term string is preferred to vector because it emphasizes the fact that all environmental values of a string are dependent upon one another, an important property that enables the ENM to consider locally relationships among variables (and not globally as most ENMs). As an example, Figure 1 shows a number of environmental strings characterized by six nodes (i.e., six hypothetical environmental variables). Some environmental strings have an associated level of species abundance, called abundance-known environmental string (ES) hereafter. In Figure 1, the blue string corresponds to an abundance level of 0, the red, green, black, and yellow, 50, 60, 5, and 0, respectively. Other strings have no associated abundance estimation; they are called abundance-unknown environmental string, or target environmental string (TES). In Figure 1, the white string (i.e., TES) has no associated abundance.

The gray zone in Figure 1 represents the target zone (TZ) where each node of an ES must occur to be included in the calculation of the abundance associated with the TES. This explains why the term string is preferred to the term vector because if only one node of a string is outside TZ, the string is not used to estimate the abundance associated with TES. The distance between the node of an ES and the TES is calculated for all environmental variables. The total distance between an ES and TES within TZ is subsequently used to weight their influence in the estimate of the abundance of the TES (here the white TES). In the example, the blue, yellow, and black ESs cannot be used to assess the white ES because at least one node is outside TZ. In most models, the values of the black string (i.e., at least those within TZ) may not only be used but may also have a strong weight in the assessment of the value of TES, despite the fact that the value of the first environmental variable is far from TZ (Figure 1). Only the green and red ESs can be used in this example. When the number of local ESs (i.e., ESs within TZ) is high enough, the abundance of a TES can be estimated; a threshold of ESs can be implemented in the ENM to assess the level of abundance of a species.

## Mathematical description

### Step 1: Representing the environmental string in a two-dimensional Euclidean space

Representing the different environmental strings (i.e., ES and TES) in a two-dimensional Euclidean space is accomplished by performing a standardization of the environmental dataset between 0 and 1. Each vector (i.e., environmental variable) is standardized as follows:

$$ \mathbf{x}^* = \frac{\mathbf{x} - \min(\mathbf{x})}{\max(\mathbf{x}) - \min(\mathbf{x})} \quad (1) $$

where x is the column vector, with all values of a given environmental variable (e.g., temperature, salinity); note that the minimum and maximum values can be determined from large-scale environmental datasets (e.g., the whole spatial and temporal range of a variable at a biogeographical scale) or based on physiology (e.g., minimum or maximum values associated to a given physiological state). After standardization, the value of each environmental variable is between 0 (minimum value of the variable) and 1 (maximum value of the variable). This stage is important because it enables the representation of the multidimensionality of the environmental strings in a two-dimensional space (Kléparski et al., 2021) and allows an equal contribution of each environmental variable to the estimates, as in neural network (Kim, 2017).
Each sample with an associated value of abundance (positive or nil) can therefore be represented by an environmental string. All strings associated with a value of abundance (i.e., ES, e.g., blue, red, green, black and yellow strings in Figure 1) are in matrix $C_{np}$ and their associated abundance (or biomass) values are in $A_{ns}$, where n is the number of samples, p is the total number of environmental variables or dimensions, and s is the number of species.

Environmental strings (i.e., TES, e.g., white string in Figure 1) for which there are no associated values of abundance are stored in matrix $U_{mp}$, with m being the number of observations where an estimation of abundance is sought. Observations can correspond to the environment observed in a geographical cell at a given day, month, and year.

### Step 2: Selection of ESs

To estimate the abundance of a species, the algorithm first searched ESs (matrix $C_{np}$) in the vicinity of a given TES (i.e., within the TZ); the string may be represented by the environment within a geographical cell for a specific time. To do so, the distance between each environmental value between the ESs and a TES is first calculated. Only ESs that are within a given threshold of distance $d_j$, that is, the TZ, which can vary from one variable to another, are kept. The vector of thresholds of distance is $D_p = [d_j]$. Because each environmental value is standardized between 0 and 1 for each environmental variable, this number should be between ~0.01 and 0.5 (i.e., distance above or below a point of a string). Note that, in practice, the threshold ranges between ~0.01 and ~0.2 depending on the number of available ESs and the degree of smoothing considered for the estimation of abundance (see below). Because a threshold of distance exists for each environmental dimension, this enables the consideration of the local relationships among all environmental variables in contrast to most ENMs that only considered global relationships.

### Step 3: Estimation of abundance by weighted average

The estimated abundance $\hat{A}$ of a TES for species s is estimated as follows:

$$ \hat{A} = \frac{\sum_{i=1}^{t} w_i A_i}{\sum_{i=1}^{t} w_i} \qquad (2) $$

where t is the total number of ESs in the TZ of TES, i is the ith ES within the distance imposed by the

vector of thresholds of distance $D_p$. $w_i$ is the ith weight associated to $A_i$. The resulting m estimations of abundance, corresponding to $U_{mp}$, are stored in $\hat{A}_{ms}$. Therefore, an estimated abundance $\hat{A}$ of a TES for species s is based on ESs that are close in the Euclidean space of a TES (not a geographical space). Equation (2) is similar to many equations applied in geostatistics (Wackernagel, 1995), including the inverse distance weighting (IDW) method (Chen & Liu, 2012). The weights $w_i$ are, however, as follows:

$$ w_i = v e^{-\lambda k_i} \qquad (3) $$

where v is the maximum weight for a nil distance between TES and ES, that is, $k_i = 0$, with $v \ge 1$. When $v = 1$, all $w_i = 1$ and Equation (2) is a unweighted arithmetic mean. Note that v cancels in Equation (2) but its influence still remains through $\lambda$ (Equation 6). When $v > 1$, Equation (2) becomes a weighted arithmetic mean. An increase in v gives more weight in the calculation of the average to a string that is close to a TES. The multivariate distance k (Manhattan metric) between a TES and an ES is calculated as follows (Legendre & Legendre, 1998):

$$ k_i = \sum_{j=1}^{p} |h_j - q_j| \qquad (4) $$

with $h$ ($0 \le h \le 1$) and $q$ ($0 \le q \le 1$) being the values of each environmental variable (standardized between 0 and 1) between a TES ($h$) and an ES ($q$).

The parameter $\lambda$ is fixed, so that $w_i = 1$ at $k_i = d_{tot}$, with $d_{tot}$ assessed as follows:

$$ d_{tot} = \sum_{j=1}^{p} d_j \qquad (5) $$

with p being the total number of environmental dimensions (or variables) and $d_j$ being the threshold of distance fixed for each dimension ($j = 1, ... p$). When $|h_j - q_j| > d_j$, the ES is discarded and does not contribute to the estimation of a TES (see below).

$\lambda$ is determined from Equation (3) to scale $w_i$ to 1 when $k_i = d_{tot}$:

$$ \frac{1}{v} = e^{-\lambda d_{tot}} \quad \text{so} \quad \ln\left(\frac{1}{v}\right) = -\lambda d_{tot}, \quad \text{which gives} \\ \lambda = \frac{-\ln(1/v)}{d_{tot}}, \quad \text{so:} $$

$$ \lambda = \frac{\ln(v)}{d_{tot}} \qquad (6) $$

From Equation (3), when $k_i = 0$ (i.e., the TES and ES are identical), $w_i = v$, and when $k_i = d_{tot}$, $\lambda$ is calculated
from Equation (6), so that $w_i = 1$. For example, for a nine-dimensional ES, $v = 100$, and $d_{tot} = 0.9$ (e.g., $d_j = 0.1$ for each of the nine environmental variables), it is easy to check that $\lambda = 5.1169$.

Therefore, this model depends upon the selection of two thresholds: (1) $v$ is the maximum weight for a distance $k_i = 0$. When $v = 100$, there is a maximum factor of 100 in the weight between ES and TES (i.e., $k(i) = 0$, $w_i = 100$ and $k(i) = d_{tot} w_i = 1$); and (2) thresholds of distance $D_p = [d_j]$ with $0 < d_j \le 0.5$ for each environmental variable $j$.

### Step 4: Calculation of some indices on the quality of the estimated abundance

I used four indicators to assess the robustness of the abundance estimation in space and time.

First, the total number of ES used to estimate a TES is an important indicator to have a first idea of the robustness of the estimated abundance. When there is no ES within TZ around TES, there is no estimation. The absence of estimation due to lack of abundance-known ES within the interval fixed by $D_p$ is coded -9999. This is important as most ENMs/SDMs would provide an estimate, whereas ecological theory suggests that the estimate may be far from actual abundance.

Second, for each estimation of abundance $\hat{A}(u_0)$, which is a weighted arithmetic mean, the weighted SD, $s_w$, can be calculated as follows:

$$s_w = \sqrt{\frac{\sum_{i=1}^{t} w_i (A_i - \hat{A}_i)^2}{\frac{t-1}{t} \sum_{i=1}^{t} w_i}} \quad (7)$$

with $t$ being the number of ESs within the TZ.

Third, from Equations (2) and (7), the $CV_w$ at TES can be assessed:

$$CV_w = \frac{s_w}{\bar{A}_i} 100 \quad (8)$$

For such indices to be reliable, a sufficient number of ESs in the TZ is needed. This can be controlled from the knowledge of the total number of ESs within TZ (see first indicator above).

Four, the average multivariate distance $\phi$ between ES and TES is assessed as follows:

$$\phi = \frac{\sum_{i=1}^{t} k_i}{t} \quad (9)$$

with $k_i$ calculated from Equation (4). A small average distance suggests that estimations are likely to be more accurate.

## Data

### Biological data

The CPR survey is a long-term, sub-surface marine plankton monitoring program consisting of a network of CPR transects towed monthly across the major geographical regions of the North Atlantic and its adjacent seas (Batten et al., 2003; Reid et al., 2003). It has been operating in the North Sea since 1931, with some standard routes existing with a virtually unbroken monthly coverage back to 1946. The CPR survey is recognized as the longest sustained and geographically most extensive marine biological survey in the world (Reid et al., 2003). The dataset comprises a uniquely large record of marine biodiversity covering ~1000 taxa over multidecadal periods. The survey determines the abundance and distribution of phytoplankton and zooplankton (including fish larvae) in the oceans and shelf seas. Using ships of opportunity from ~30 different shipping companies, it obtains samples at monthly intervals on ~50 transocean routes. The CPR uses a high-speed plankton recorder that is towed behind “ships of opportunity” through the surface layer of the ocean (~10 m depth) (Warner & Hays, 1994). Water passes through the recorder and plankton are filtered by a slow-moving silk (mesh size 270 µm). A second layer of silk covers the first and both are reeled into a tank containing 4% formaldehyde. Upon returning to the laboratory, the silk is unwound and cut into sections corresponding to 10 nautical miles and approximately 3 m³ of filtered sea water (Jonas et al., 2004). In this paper, I used data on the abundance of *Calanus finmarchicus* and *Calanus helgolandicus* because these species are well sampled by the CPR survey (Helaouët & Beaugrand, 2007). These species (>2 mm) are enumerated off-silk (Warner & Hays, 1994). Abundance data corresponding to the period 1958–2016 were used here.

### Environmental data

Sea surface temperature (SST) originated from Centennial in situ Observation Based Estimates (COBE) SST2 dataset provided by the National Oceanic and Atmospheric Administration (NOAA), Ocean and Atmospheric Research (OAR), Earth System research Laboratory's (ESRL), Physical Sciences Division (PSD), Boulder Colorado, USA (https://www.esrl.noaa.gov/psd/). It is a 180 × 360 global grid of 1° latitude by 1° longitude from 89.5° N to 89.5° S and 0.5° E to 359.5° E covering the time period from 1850 to 2016. In this study, I considered the period 1958–2016, which corresponded to the CPR timescale. Units were in degrees Celsius.
Bathymetry data originated from the BIO-ORACLE V2.0 dataset. It encompassed 23 ecologically relevant surface and benthic marine layers tailored for mechanistic and correlative modeling approaches, at a monthly scale with a spatial resolution of 0.08° or 9.2 km at the equator (http://www.bio-oracle.org/). Bathymetry was expressed in meters.

Monthly mean downward solar radiation flux at surface (DSRF) originated from the National Centers for Environmental Prediction (NCEP), National Center for Atmospheric Research (NCAR) reanalysis monthly means, and other derived variables, provided by the NOAA/OAR/ESRL PSD, Boulder Colorado, USA (https://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.html). The temporal coverage was from 1948 to present at a monthly scale. The dataset had a global spatial coverage with a grid of 94 latitudes (from 88.542° N to 88.542° S) by 192 longitudes (from 0° E to 358.125° E). I used this variable as a proxy of PAR (photosynthetically active radiation) because it covered the northern part of the studied zone (above 56° N, there are no PAR data in winter). Unit was in watts per square meter.

Nutrients data originated from the World Ocean Database (WOD) (https://www.nodc.noaa.gov/OC5/WOD/pr_wod.html). I used three variables: (1) nitrate, (2) phosphate, and (3) silicate concentration. Data were expressed in micromoles per kilogram.

Distance to nearest coastline was provided by the University Corporation for Atmospheric Research (UCAR), community programs, Boulder Colorado, USA (http://oos.soest.hawaii.edu/thredds/ncss/dist2coast_1deg_ocean/dataset.html). The data were on a grid from 90° N to 89.9899° S and from 179.99° E to 180° W, with a spatial resolution of 0.01°. Units were in kilometers.

Density mix layer depth (DMLD) originated from Monthly Isopycnal and Mixed-layer Ocean Climatology (MIMOC) provided by NOAA (https://www.pmel.noaa.gov/mimoc/). MIMOC is a set of NetCDF maps, based on Argo CTD data, shipboard data, and ice-tethered profiler CTD data, with a covering from 90° N to 80° S and at a 0.5° resolution. DMLD is the pressure at the bottom of the mixed layer, thus very similar to the mixed-layer depths in meter. Units were in decibar.

Sea surface salinity (SSS) originated from NCPE Global Ocean Data Assimilation System (GODAS) dataset, provided by NOAA ESRL, PSD, Boulder Colorado, USA. The dataset contains the mean monthly salinity from 1980 to present. Units were in kilograms per kilogram (https://www.esrl.noaa.gov/psd/cgiin/DataAccess.pl?DB_dataset=NCEP+GODASandDB_variable=salinityandDB_statistic=Monthly+MeanandDB_tid=72048andDB_did=98andDB_vid=1914).

All environmental data were arranged on a spatial grid 20.5° N-79.5° N and 99.5° W-19.5° E, every 1° of latitude and 1° of longitude, and on a monthly basis for 1958-2016. When environmental data were not available for the period, it was replaced by the long-term monthly average. Bathymetry data were duplicated for all months of the period 1958-2016.

## Test of the model

### Sensitivity analysis

I ran the model for two calanoid copepod species well sampled by the CPR survey: *C. finmarchicus* and *C. helgolandicus*. I performed two types of sensitivity analyses. The first type of sensitivity analysis was carried out to examine the influence of two parameters *d* and *v* on the model: (1) *v* the maximum weight at distance *k* = 0 and (2) *D*<sub>*p*</sub> = [*d*<sub>*j*</sub>] the vector of distance to the target string. Here, *v* was fixed to 10, 10<sup>2</sup>, and 10<sup>3</sup> for *d*<sub>*j*</sub> = 0.05. Different thresholds, *d*<sub>*j*</sub> = 0.025, *d*<sub>*j*</sub> = 0.05, and *d*<sub>*j*</sub> = 0.01, were used for all environmental dimensions (or variables), with *v* = 10<sup>2</sup> (Table 1). An additional sensitivity analysis was carried out for *d*<sub>*j*</sub> = 0.1 and *v* = 10<sup>9</sup>. For all sensitivity analyses, I used CPR data of the period 1958-1999 (i.e., model calibration) to model the abundance of both *Calanus* for each month for the period 2000-2016. An order-2 bidimensional simple moving average was applied to smooth the abundance in space for each month because of the strong noise evident in the monthly CPR data (Beaugrand, 2004; Beaugrand & Edwards, 2001). All tests to assess the validity of the model against observations were therefore performed on the period 2000-2016. This period was chosen because of a major warming that substantially altered the spatial distribution of copepods in the North Atlantic prior to 2000 (Beaugrand et al., 2009).

### Observed abundance data of *Calanus*

I examined the spatial patterns in *C. finmarchicus* and *C. helgolandicus* on a monthly basis from January 2000 to December 2016. To do so, I interpolated the abundance data on the same spatial grid than environmental data (20.5° N-79.5° N and 99.5° W-19.5° E, every 1° of latitude and 1° of longitude) and therefore the modeled abundance. The method of interpolation was the inverse squared distance method, with a search radius of 250 km and a number of neighbors comprised between 10 and 20 CPR samples (Beaugrand et al., 2000; Wackernagel, 1995). An order-2 bidimensional simple moving average was applied to slightly smooth the abundance estimates in space for each month of 2000-2016 (Kléparski et al., 2021).
## Visual comparison between observed and modeled data (average of all months of 2000–2016)

I mapped the average observed and modeled abundance based on all months of the period 2000–2016 to visually compare results from the ESM using different values of distance to target string *d* and maximum weight *v*. I used six values of *d*<sub>j</sub>/<i>v</i> for both *Calanus* (Figure 2): *d*<sub>j</sub> = 0.05 for *v* = 10, *v* = 10<sup>2</sup>, and *v* = 10<sup>3</sup>, *d*<sub>j</sub> = 0.025 and *d*<sub>j</sub> = 0.01 for *v* = 10<sup>2</sup>, and *d*<sub>j</sub> = 0.1 for *v* = 10<sup>9</sup>. This last pair of values was chosen to see how results differed for such an extreme maximum weight in comparison with other more moderate values of *d*<sub>j</sub>/<i>v</i>.

## Numerical comparison between observed and modeled data for all months of 2000–2016

I compared observed and modeled abundance for all months of 2000–2016. Here also, I used six values of *d*<sub>j</sub>/<i>v</i> for both *Calanus* (Table 1): *d*<sub>j</sub> = 0.05 for *v* = 10, *v* = 10<sup>2</sup>, and *v* = 10<sup>3</sup>, *d*<sub>j</sub> = 0.025 and *d*<sub>j</sub> = 0.01 for *v* = 10<sup>2</sup>, and *d*<sub>j</sub> = 0.1 for *v* = 10<sup>9</sup>.

Comparisons were carried out by means of five indices:

First, the Pearson coefficient *r* of linear correlation was calculated (Sokal & Rohlf, 1995). This analysis was performed to examine whether patterns of abundance varies in space or time in the same way than observations.

Second, the bias was estimated as follows:

$$ \text{Bias} = \frac{\sum_{i=1}^{\beta} (\hat{A}_i - A_i)}{\beta} \qquad (10) $$

with *β* being the number of pairs of points involved in the comparison. When *bias* = 0, it means that there is no systematic underestimation (when *bias* < 0) or overestimation (when *bias* > 0) of abundance.

Third, the mean absolute error (MAE) that measures the average magnitude of the errors in a set of predictions was used, without considering their direction.

$$ \text{MAE} = \frac{\sum_{i=1}^{\beta} |\hat{A}_i - A_i|}{\beta} \qquad (11) $$

Equation (11) represents the absolute differences between predictions $\hat{A}_i$ and observations $A_i$, divided by the number of differences to be tested (with all individual differences having equal weight). It is a negatively

oriented score, which means that the lower values are related to the best correlations. This analysis gives an idea of the mean absolute difference between predictions and observations. However, values also depend upon the way observations are averaged in space and time, so that this is only a rough indicator of model efficiency.

The root mean square error (RMSE) is calculated as follows:

$$ \text{RMSE} = \sqrt{\frac{\sum_{i=1}^{\beta} (\hat{A}_i - A_i)^2}{\beta}} \qquad (12) $$

The RMSE is more sensitive to extreme differences than the MAE. Last, the RMSE was also expressed as a percentage of the average observed abundance. Percentage RMSE (PRMSE) is calculated as follows:

$$ \text{PRMSE} = \frac{\text{RMSE}}{\frac{\sum_{i=1}^{\beta} A_i}{\beta}} \qquad (13) $$

## Mapping of the indices of modeled abundance quality

For all months, I counted the number of strings and averaged the distance to target string *φ* used to estimate the abundance in each geographical cell of the North Atlantic Ocean. I also assessed the weighted SD (*s*<sub>w</sub>) and the weighted CV(CV<sub>w</sub>) of both *Calanus* (see previous section, Figure 3).

## Seasonal changes in the spatial distribution of observed and modeled abundance (2000–2016)

I compared the seasonal changes in observed and modeled abundance for both *Calanus*. For this analysis, I averaged the data by 3-month periods for 2000–2016: (1) January to March, (2) April to June, (3) July to September, and (4) October to December (Figure 4).

## Monthly comparisons in the North Sea

Monthly comparisons were carried out in the North Sea for 2000–2016. All observed and modeled abundances were averaged monthly for the period 2000–2016 in the area 51° N–60° N and 3° W–10° E (Figure 5). The coefficient *r* of linear correlation and the bias were calculated (see previous section).
## Comparison with three existing models

I also compared the ESM with three regression models: GLMs (McCullagh & Nelder, 1983), GAMs (Hastie & Tibshirani, 1990), and GRNN (Specht, 1991). As with GAMs, GRNN is a method that can be used when the assumption of linearity is not justified. The method provides smooth transitions from one observation to another, even when data are sparse in the multidimensional space (Specht, 1991). In contrast to the GRNN, GLMs and GAMs are statistical techniques that have been frequently used in niche modeling (Austin & Meyers, 1996; Hirzel & Le Lay, 2008; Qiao et al., 2019; Thuiller et al., 2009).

I applied two different GLMs (MATLAB function “fitglm”), one with linear responses (Poisson distribution), with no interaction among environmental variables, and a second with quadratic responses (Poisson distribution), with interaction among variables (i.e., linear and squared terms for each predictor, and all products of pairs of distinct predictors).

I used two GAMs with and without interactions among independent variables. The MATLAB function (“fitrgam”) fits a GAM by means of a gradient boosting algorithm (least-squares boosting) (Yin et al., 2012, 2013). It uses a boosted tree as a shape function for each predictor to capture the nonlinear relation between each explanatory and the response variable. The second GAM included all possible pairs of interaction terms (by adding bivariate shape functions), which was $p(p-1)/2$ number of terms in total; because $p=9$, the number of possible interactions was 36.

I applied the GRNN with eight values of smoothing: spread = 0.001, spread = 0.005, spread = 0.01, spread = 0.05, spread = 0.1, spread = 0.2, spread = 0.3, and spread = 0.7. When the spread is larger, the smoothing of the approximation function is greater, and vice versa (Specht, 1991). I only kept the value of spread that gave the best GRNN in terms of correlations with the observations.

Comparisons for all months of the period 2000–2016 were done by means of the same metrics than the ones used to compare the ESM with observed data (Table 1). I also mapped the spatial distribution of the average modeled abundance of both *Calanus* for the different models (Figure 6).

## Investigation of the most relevant variables in the model

My ESM was applied on nine environmental variables. To apply the model on each combination of variables

from $p=1$ to $p=9$, that would have involved $C=427$ sensitivity analyses:

$$C = \sum_{i=1}^{p} \binom{p}{k_i} = \sum_{k=1}^{p} \frac{p!}{k!(p-k)!} = 427 \quad (14)$$

Such a high number of analyses would have required a long time. I therefore reduced the number of dimensions (or variables) by first making analyses using all variables but one (i.e., a total of nine analyses, e.g., all but distance to coast, all but bathymetry and so on). On the basis on these first results, I chose to always keep three variables in the combinations, bathymetry, SST, and DSRF, which are key variables in marine ecology (Raybaud et al., 2013; Rombouts et al., 2012; Schickele et al., 2021). This decreased the number of analyses from $C=427$ to 63 (see Equation 14, with $n=6$ variables). I finally retained the variables or combinations of variables that gave the best results from a combination of four (3 being already known: bathymetry, SST, and DSRF) to nine variables with observed data, considering all months of the period 2000–2016 in the comparisons (Table 3). Analyses were done for $d_j = 0.05$ and $v = 10^2$.

Because the increase in the number of descriptors in a model may raise artificially the percentage of explained variance, I also calculated the adjusted $R^2$ and Akaike information criterion (AIC) (Akaike, 1974; Ostertagová, 2012). Adjusted $R^2$ is more a measure of statistical quality (percentage of explained variance), whereas AIC is more indicative of the predictive capability of a model (Harrell, 2001; Ostertagová, 2012; Tredennick et al., 2021).

## Use of the ESM on presence/absence data

To investigate whether or not the ESM could be applied on presence/absence data, and therefore on presence-only after the generation of pseudo-absence, I converted the abundance data of both *Calanus* in presence/absence and performed an analysis with $d_j = 0.05$ and $v = 10^2$. I investigated the visual relationships between the average annual and seasonal maps of observed *Calanus* abundance (Figure 7) and linear correlations were calculated for all months of the period 2000–2016. In addition, I used the true skill statistic (TSS) and the kappa statistic to evaluate the relationships between (Allouche et al., 2006). These calculations were done using a threshold of 0.5 to separate unsuitable from suitable areas (Hirzel et al., 2006).
# RESULTS

## Influence of the distance to target string and maximum weight

I first examined the influence of the choice of the distance to target string $d_j$ and maximum weight $v$ by comparing modeled and observed abundances (Figure 2 and Table 1); recall that only 1958–1999 observations were considered to assess the abundance of the species for 2000–2016.

For both copepods, there was a good correspondence between observed and modeled abundance averaged for the period 2000–2016 (Figure 2). A visual inspection

reveals that modeled abundances were close to observations for all values of $d_j$ and $v$. No changes were observed for the different values of $v$ but an alteration was detected for small values of $d_j$ ($d_j = 0.01$) because the number of strings within the TZ diminishes when $d_j$ decreases (i.e., more cases with no estimates because of an absence of string; Figure 2F).

This visual inspection was confirmed by the performance metrics. Note that for this analysis, performance metrics were applied for all geographical cells (1° latitude × 1° longitude), all months and years of the period 2000–2016 (i.e., no average as in Figure 2). The choice of different values of $d_j$ and $v$ had only a minor influence on the results (Table 1). For *C. finmarchicus*,

<table><caption>TABLE 1 Comparisons between observed and modeled abundance for all months of the period 2000–2016, based on the proposed Environmental String Model with different values of distance to target string $d_j$ and maximum weight $v$, generalized linear model (GLMs; Poisson distribution, with and without interaction), generalized additive models (GAMs; with and without interaction), and generalized regression neural networks (GRNNs; different degree of smoothing).</caption><thead><tr><th>Distance to target string, d<sub>j</sub></th><th>Maximum weight, v</th><th>m</th><th>r</th><th>Bias</th><th>MAE</th><th>RMSE</th><th>PRMSE</th></tr></thead><tbody><tr><td colspan="8"><em>Calanus finmarchicus</em></td></tr><tr><td>0.05</td><td>10</td><td>121,999</td><td>0.78<sup>a</sup></td><td>0.18</td><td>0.31</td><td>0.21</td><td>46.09</td></tr><tr><td>0.05</td><td>10<sup>2</sup></td><td>121,999</td><td>0.78<sup>a</sup></td><td>0.18</td><td>0.31</td><td>0.21</td><td>45.24</td></tr><tr><td>0.05</td><td>10<sup>3</sup></td><td>121,999</td><td>0.79<sup>a</sup></td><td>0.17</td><td>0.30</td><td>0.20</td><td>44.54</td></tr><tr><td>0.025</td><td>10<sup>2</sup></td><td>110,354</td><td>0.77<sup>a</sup></td><td>0.11</td><td>0.28</td><td>0.20</td><td>44.79</td></tr><tr><td>0.01</td><td>10<sup>2</sup></td><td>49,684</td><td>0.71<sup>a</sup></td><td>0.04</td><td>0.27</td><td>0.21</td><td>52.60</td></tr><tr><td>0.1</td><td>10<sup>9</sup></td><td>122,812</td><td>0.79<sup>a</sup></td><td>0.21</td><td>0.32</td><td>0.21</td><td>46.42</td></tr><tr><td>GLM (linear, Poisson distribution, no interaction)</td><td></td><td>122,854</td><td>0.64<sup>a</sup></td><td>0.46</td><td>0.58</td><td>0.46</td><td>100.24</td></tr><tr><td>GLM (quadratic, Poisson distribution, interaction)</td><td></td><td>122,854</td><td>0.76<sup>a</sup></td><td>0.21</td><td>0.34</td><td>0.23</td><td>49.33</td></tr><tr><td>GAM (no interaction)</td><td></td><td>99,150</td><td>0.66<sup>a</sup></td><td>0.63</td><td>0.72</td><td>0.67</td><td>119.17</td></tr><tr><td>GAM (all possible pairs of interaction)</td><td></td><td>111,367</td><td>0.64<sup>a</sup></td><td>0.76</td><td>0.82</td><td>0.84</td><td>172.12</td></tr><tr><td>GRNN (spread = 0.05)</td><td></td><td>122,854</td><td>0.78<sup>a</sup></td><td>0.23</td><td>0.33</td><td>0.23</td><td>49.28</td></tr><tr><td colspan="8"><em>Calanus helgolandicus</em></td></tr><tr><td>0.05</td><td>10</td><td>121,999</td><td>0.73<sup>a</sup></td><td>-0.02</td><td>0.17</td><td>0.09</td><td>32.66</td></tr><tr><td>0.05</td><td>10<sup>2</sup></td><td>121,999</td><td>0.73<sup>a</sup></td><td>-0.03</td><td>0.17</td><td>0.09</td><td>32.76</td></tr><tr><td>0.05</td><td>10<sup>3</sup></td><td>121,999</td><td>0.73<sup>a</sup></td><td>-0.03</td><td>0.17</td><td>0.09</td><td>32.92</td></tr><tr><td>0.025</td><td>10<sup>2</sup></td><td>110,354</td><td>0.70<sup>a</sup></td><td>-0.07</td><td>0.17</td><td>0.10</td><td>36.11</td></tr><tr><td>0.01</td><td>10<sup>2</sup></td><td>49,684</td><td>0.61<sup>a</sup></td><td>-0.12</td><td>0.22</td><td>0.16</td><td>46.24</td></tr><tr><td>0.1</td><td>10<sup>9</sup></td><td>122,812</td><td>0.72<sup>a</sup></td><td>-0.01</td><td>0.18</td><td>0.09</td><td>33.40</td></tr><tr><td>GLM (linear, Poisson distribution, no interaction)</td><td></td><td>122,854</td><td>0.55<sup>a</sup></td><td>0.07</td><td>0.27</td><td>0.13</td><td>49.33</td></tr><tr><td>GLM (quadratic, Poisson distribution, interaction)</td><td></td><td>122,854</td><td>0.72<sup>a</sup></td><td>-0.02</td><td>0.18</td><td>0.09</td><td>32.24</td></tr><tr><td>GAM (no interaction)</td><td></td><td>99,150</td><td>0.54<sup>a</sup></td><td>0.18</td><td>0.36</td><td>0.18</td><td>55.65</td></tr><tr><td>GAM (all possible pairs of interaction)</td><td></td><td>111,367</td><td>0.47<sup>a</sup></td><td>0.29</td><td>0.41</td><td>0.25</td><td>87.36</td></tr><tr><td>GRNN (spread = 0.05)</td><td></td><td>122,854</td><td>0.72<sup>a</sup></td><td>-0.00</td><td>0.18</td><td>0.09</td><td>33.39</td></tr></tbody></table>

*Note*: When the bias is positive, modeled abundances are on average higher than observed abundance, and vice versa. *m* is the total number of geographical cells considered in the comparisons, and *r* is the coefficient of linear correlation. For the meaning of other performance metrics, see *Materials and methods*.

Abbreviations: MAE, mean absolute error; PRMSE, percentage RMSE; RMSE, root mean square error.

<footnote><sup>a</sup>Significant coefficient of correlation.</footnote>
<figure><img src="image_2.png" /><figcaption>FIGURE 2 Legend on next page.</figcaption></figure>
actually observed in Figure 2, the best results were achieved for $d_j = 0.05$ and $v = 10^3$ for *C. finmarchicus* and for $d_j = 0.05$ and $v = 10$ for *C. helgolandicus* (Table 1). There were low differences among performance metrics, however. That is why I chose to consider $d_j = 0.05$ and $v = 10^2$ for both *Calanus* in all analyses below.

## Indicators on the quality of abundance estimates

The ESM provides some key indicators on the quality of the abundance estimates (Figure 3). The first indicator is the number of strings that contributed to the abundance estimates for each geographical cell (Figure 3A); the more strings available, the better the estimation is likely to be. This number was high around the British Isles and much lower in the oceanic part south of Greenland and Iceland. In the present study, the number of strings had a low influence on the large-scale spatial distribution of both *Calanus*, and all indicators were virtually not affected by the removal of geographical cells whose abundance estimates were based on less than 10 strings for

the whole period 2000–2016 (e.g., correlations rose from 0.78 to 0.80 for *C. finmarchicus* and did no change for *C. helgolandicus*, and MAE decreased from 0.31 to 0.30 for *C. finmarchicus* and from 0.169 to 0.166 for *C. helgolandicus*). The second indicator is the mean distance to target string φ (Figure 3B). The examination of this indicator showed that the mean distance to target string was generally in the range of 20%–30% in the area sampled by the CPR survey. This indicator was slightly higher over the mid-Atlantic ridge, especially to the south (~30%–40%), and therefore one might expect less reliable abundance estimates in this area. The third indicator is the weighted SD of the abundance estimates (Figure 3C,E). Because the higher the abundance, the greater the weighted SD ($s_w$; Figure 3C,E vs. Figure 2C,J), I used a fourth indicator: the weighted CV ($CV_w$; Figure 3D,F). This latter indicator was high south of the spatial range of *C. finmarchicus* and to a lesser extent south of the North Sea and to the west of Ireland. In this latter area, the estimates might be less reliable. For *C. helgolandicus*, weighted CV was high in the subarctic gyre and to a lesser extent south of NewFoundand. In this latter area, the estimates might also be less reliable.

<figure><img src="image_3.png" /><figcaption>**FIGURE 4** Seasonal changes in the spatial distribution of observed and modeled abundance of *Calanus finmarchicus* and *Calanus helgolandicus* for 2000–2016. The Environmental String Models based on $d_i = 0.05$ and $v = 10^2$ were used here for both *Calanus*.</figcaption></figure>
## Seasonal changes in the spatial distribution of species (2000–2016)

Although regional differences were sometimes detected, the ESM reproduced the seasonal changes in the spatial distribution of both *Calanus* globally well (Figure 4). Some differences between observations and modeled data were found south of Greenland for January–March and in the northern part of the North Atlantic Drift Province (Longhurst, 1998) in April–June for *C. finmarchicus*, as well as in the North Sea for most seasons for *C. helgolandicus*.

## Seasonal changes in the abundance of *Calanus* in the North Sea

Seasonal changes in the abundance of both *Calanus* were well reconstructed by the model for the period 2000–2016 (Figure 5); the linear correlation coefficient was *r* = 0.86 (*p* < 0.01) and *r* = 0.84 (*p* < 0.01) for *C. finmarchicus*.

<figure><img src="image_4.png" /><figcaption>**FIGURE 5** Seasonal changes in the observed (blue) and modeled (red) abundance of *Calanus finmarchicus* and *Calanus helgolandicus* for 2000–2016 in the North Sea. The chosen Environmental String Models were based on *d*ᵢ = 0.05 and *v* = 10² for both *Calanus*.</figcaption></figure>

and *C. helgolandicus*, respectively. The model overestimated the observed abundance for *C. finmarchicus* (bias = 0.27, i.e., ~0.86 individual per month on average) and underestimated (bias = -0.31, i.e., ~1.04 individual per month on average) for *C. helgolandicus*.

## Comparisons with other ENMs

I compared results obtained from the ESM with GLMs, GAMs, and GRNN (Table 1 and Figure 6). A visual comparison between observed and modeled average annual abundance of both *Calanus* revealed that best results were achieved with the ESM and GRNN, followed closely by GLM (with interaction) (Figure 6 vs. Figure 2A,H). GLMs (without interaction) and GAMs (with and without interactions) gave poor results compared with observations (Figure 6 vs. Figure 2A,H). This first visual inspection was confirmed by the performance metrics performed on observed and modeled data from all months of the period 2000–2016 (Table 1). The ESM and the GRNN had similar performance metrics, although the former was slightly better in terms of correlations and MAE/RMSE for both *Calanus*. In terms of bias, the GRNN performed better than the ESM for *C. helgolandicus*, but not for *C. finmarchicus* (Table 1). The GLM with interaction did not perform better than the ESM and the GRNN for both *Calanus*, with the exception of the PRMSE metric for *C. helgolandicus* (Table 1).

## Identification of the most important environmental variables

To identify the most important variables in the model, I first carried out sensitivity analyses based on all but one variable (Table 2). This first series of analyses revealed that SST and DSRF (for both *Calanus*), as well as bathymetry (for *C. helgolandicus*), had a high influence. Indeed, excluding these variables from the ESM diminished its predictive power for the period 2000–2016 (Table 2).

I then performed 63 sensitivity analyses, which represented the combinations of all possible variables, but keeping all the time SST, DSRF, and bathymetry, which were important in terms of model performance for at least one species of *Calanus* (Table 2). Therefore, all combinations were tested based on four to nine simultaneous variables (Table 3). I only reported the best fit with the observed data for all months of the period 2000–2016.

Model performance was evaluated by many metrics (Table 3). Using adjusted *R*² and AIC that both penalize the addition of a variable, model performance was best
<figure><img src="image_5.png" /><figcaption>**FIGURE 6** Legend on next page.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**FIGURE 7** Modeled spatial distribution in the annual and seasonal spatial distribution of *Calanus finmarchicus* and *Calanus helgolandicus* for the period 2000–2016 ($d_i$ = 0.05 and $v$ = 10²). (A, F) Annual abundance. (B, G) Average abundance for January to March. (C, H) Average abundance for April to June. (D, I) Average abundance for July to September. (E, J) Average abundance for October to December. Comparison of the mean occurrence of both *Calanus* can be made with average abundance shown in Figures 2 and 4.</figcaption></figure>

for *C. finmarchicus* when either 8 or 7 variables (adjusted $R^2$, MAE, PRMSE, and RMSE as indicative of model quality) or 6 variables (AIC as indicative of predictive capability) were considered simultaneously (Table 3).

Distance to coast, silicate, and SSS were not particularly relevant as the ESM also performed well when these variables were not considered. Bias was the smallest when the model was based on all variables.

<caption>**FIGURE 6** Mapping of the mean annual spatial distribution in the abundance of *Calanus finmarchicus* and *Calanus helgolandicus* modeled by generalized linear models (GLMs), generalized additive models (GAMs), and generalized regression neural network (GRNN). (A, G) Environmental String Model ($d_i$ = 0.05 and $v$ = 10²). (B, H) GLM (linear, no interaction). (C, I) GLM (quadratic, interaction). (D, J) GAM (no interaction). (E, K) GAM (with interaction). (F, L) GRNN (spread = 0.05). Comparisons can be made with observed annual abundance of both *Calanus* in Figure 2A,H.</caption>
For *C. helgolandicus*, model performance was achieved when only seven variables were taken simultaneously (Table 3), revealing that silicate and DMLD were not relevant. Nitrate was also not so important as the removal of this variable did not affect substantially model performance. According to the evaluation metrics, different conclusions can be reached for both *Calanus*. PRMSE was smallest when the model was based on seven variables. Bias was nil when the model was based on five or four variables. Using the AIC that penalizes the addition of a variable, best model was composed of only four variables. Such differences have been regularly reported, depending on the emphasis of the comparison (e.g., exploration, inference, prediction) (Harrell, 2001; Tredennick et al., 2021).

For both *Calanus*, nutrients did not have a substantial influence on model performance, which was expected as these variables have only an indirect effect on their ecology (Helaouët & Beaugrand, 2007). DMLD was significant for *C. finmarchicus* (oceanic species) but not for *C. helgolandicus* (pseudo-oceanic species; species occurring in both oceanic and neritic domains but being more abundant over the shelf-edge). As expected, distance to coast and SSS were more important for *C. helgolandicus* than its congeneric more oceanic species (Table 3).

## Application of the ESM on presence/absence data

The application of the ESM on presence/absence data gave results that corresponded well with observed annual and seasonal patterns of changes in the abundance for both *Calanus* (Figure 7 vs. Figure 2A,H and Figure 4). Correlations between mean observed abundance and modeled mean occurrence for all months of the period 2000–2016 were highly significant (*p* < 0.01): 0.76 (vs. 0.79 for modeled abundance) for *C. finmarchicus* and 0.69 (vs. 0.73 for modeled abundance) for *C. helgolandicus* (see also Table 1). In addition, the kappa statistic and TSS were relatively high: 0.88 and 0.87 for *C. finmarchicus* and 0.85 and 0.82 for *C. helgolandicus*, respectively. These results therefore suggest that the model is applicable to presence/absence data, and therefore to presence-only data after generating pseudo-absences.

## DISCUSSION

### Advantages of the model

Many ENMs/SDMs can be applied to model the spatial distribution of a species (Helaouët & Beaugrand, 2009;

Hirzel et al., 2002; Melo-Merinoa et al., 2020; Philips et al., 2006; Schickele et al., 2020; Thuiller et al., 2009; Valavi et al., 2021). Model predictions may vary substantially (Raybaud et al., 2015; Schickele et al., 2020, 2021), however, and it is therefore important to investigate model performance (Araújo et al., 2019; Boria et al., 2014; Qiao et al., 2019; Saupe et al., 2012; Tredennick et al., 2021; Zhang et al., 2015) and continue to propose new ENMs that may improve modeled species' spatial distribution. The ESM has been specifically designed for this purpose and is deeply rooted in niche theory; the niche is here envisioned by means of environmental strings, which are associated with a level of abundance, nil or positive (i.e., ESs).

The ESM is easy to understand and to interpret. It simply averages locally (in the Euclidean space) the abundance values associated with combinations of variables within a TZ, an average that is weighted by the total distance *d*<sub>tot</sub> between each environmental string and the TES by means of a negative exponential function standardized between maximum weight *v* and 1 (Equations 3 and 6); the two parameters, distance to target string *d*<sub>j</sub> and maximum weight *v* needed to run the model, are therefore easy to comprehend. Although the difference was small with GRNN and to a lesser extent with GLM (Poisson distribution, with interaction), the ESM performed better than models considered here (Table 1 and Figure 6). This was especially true for GLMs (without interaction) and GAMs (Guisan et al., 2002), which are frequently applied in spatial distribution modeling, perhaps because they are available in many packages or software (e.g., R, Splus or MATLAB) (Lehmann et al., 2002; Thuiller et al., 2009). It should be noted, however, that the GRNN, less frequently used in spatial distribution modeling, also reproduced spatial changes in the monthly abundance of both *Calanus* well (Table 1 and Figure 6).

Here, I examined, in detail, the influence of *d*<sub>j</sub> and *v* on the ESM performance (Table 1). In future versions of the model, a procedure based on a cost function or performance metrics could be implemented to select the optimal combination of parameters as in machine learning (Hastie et al., 2009; Wlodarczak, 2020). In contrast to most SDMs/ENMs, the ESM can be applied on species for which there is only partial information on their spatial distribution or niche. In such a case, many models (e.g., GLM and GAM) would provide an abundance assessment even if the niche is truncated, which leads to biased estimations (Chevalier et al., 2022; Mateo et al., 2019; Scherrer et al., 2021). This is not the case with the ESM, which does not assess the abundance (or the mean occurrence) of a species where or when a biotope is composed of environmental conditions for which there is no ES in the TZ (called hereafter new environmental
<table><caption>TABLE 2 Influence of the removal of one variable on Environmental String Model performance.</caption><thead><tr><th>Variables</th><th>Ref</th><th>D2Coast</th><th>Bathy</th><th>SST</th><th>DSRF</th><th>Nitrate</th><th>Phosphate</th><th>Silicate</th><th>DMLD</th><th>SSS</th></tr></thead><tbody><tr><td colspan="11"><em>Calanus finmarchicus</em></td></tr><tr><td>j</td><td>121,999</td><td>122,210</td><td>122,286</td><td><strong>122,514</strong></td><td><strong>122,475</strong></td><td>122,431</td><td>122,377</td><td>122,055</td><td>122,119</td><td>122,160</td></tr><tr><td>r</td><td>0.78<sup>a</sup></td><td>0.79<sup>a</sup></td><td>0.79<sup>a</sup></td><td><strong>0.76<sup>a</sup></strong></td><td><strong>0.76<sup>a</sup></strong></td><td>0.78<sup>a</sup></td><td>0.78<sup>a</sup></td><td>0.78<sup>a</sup></td><td>0.78<sup>a</sup></td><td>0.78<sup>a</sup></td></tr><tr><td>Bias</td><td>0.18</td><td>0.19</td><td>0.19</td><td><strong>0.26</strong></td><td><strong>0.23</strong></td><td>0.19</td><td>0.19</td><td>0.18</td><td>0.19</td><td>0.18</td></tr><tr><td>MAE</td><td>0.31</td><td>0.31</td><td>0.31</td><td><strong>0.37</strong></td><td><strong>0.35</strong></td><td>0.32</td><td>0.31</td><td>0.31</td><td>0.32</td><td>0.31</td></tr><tr><td>RMSE</td><td>0.21</td><td>0.20</td><td>0.21</td><td><strong>0.27</strong></td><td><strong>0.24</strong></td><td>0.22</td><td>0.21</td><td>0.21</td><td>0.21</td><td>0.21</td></tr><tr><td>PRMSE</td><td>45.24</td><td>43.47</td><td>45.64</td><td><strong>57.99</strong></td><td><strong>52.74</strong></td><td>47.62</td><td>46.65</td><td>45.45</td><td>46.58</td><td>45.63</td></tr><tr><td colspan="11"><em>Calanus helgolandicus</em></td></tr><tr><td>j</td><td>121,999</td><td>122,210</td><td><strong>122,286</strong></td><td><strong>122,514</strong></td><td><strong>122,475</strong></td><td>122,431</td><td>122,377</td><td>122,055</td><td>122,119</td><td>122,160</td></tr><tr><td>r</td><td>0.73<sup>a</sup></td><td>0.73<sup>a</sup></td><td><strong>0.71<sup>a</sup></strong></td><td><strong>0.69<sup>a</sup></strong></td><td><strong>0.70<sup>a</sup></strong></td><td>0.73<sup>a</sup></td><td>0.72<sup>a</sup></td><td>0.73<sup>a</sup></td><td>0.73<sup>a</sup></td><td>0.72<sup>a</sup></td></tr><tr><td>Bias</td><td>-0.03</td><td>-0.01</td><td><strong>-0.01</strong></td><td><strong>-0.01</strong></td><td><strong>-0.00</strong></td><td>-0.02</td><td>-0.01</td><td>-0.02</td><td>-0.02</td><td>-0.02</td></tr><tr><td>MAE</td><td>0.17</td><td>0.17</td><td><strong>0.18</strong></td><td><strong>0.18</strong></td><td><strong>0.19</strong></td><td>0.17</td><td>0.18</td><td>0.17</td><td>0.17</td><td>0.17</td></tr><tr><td>RMSE</td><td>0.09</td><td>0.09</td><td><strong>0.09</strong></td><td><strong>0.09</strong></td><td><strong>0.09</strong></td><td>0.09</td><td>0.09</td><td>0.09</td><td>0.09</td><td>0.09</td></tr><tr><td>PRMSE</td><td>32.76</td><td>32.69</td><td><strong>33.81</strong></td><td><strong>35.34</strong></td><td><strong>34.86</strong></td><td>32.54</td><td>33.67</td><td>32.79</td><td>32.39</td><td>33.42</td></tr></tbody></table>

Note: All analyses were performed for $d_j = 0.05$ and $v = 10^2$ by removing only the variable indicated at the top of the column. Ref: Reference for $d_j = 0.05$ and $v = 10^2$ for all variables (see Table 1). Boldface highlights variables that are important in the model.

Abbreviations: Bathy, bathymetry; D2Coast, distance to nearest coasts; DMLD, density mixed-layer depth; DSRF, downward solar radiation flux at surface; MAE, mean absolute error; PRMSE, percentage RMSE; RMSE, root mean square error; SSS, sea surface salinity; SST, sea surface temperature.

<sup>a</sup>Significant correlation.

conditions). Therefore, the ESM also solves the issue of SDM/ENM transferability (Fernández & Hamilton, 2015; Qiao et al., 2019; Randin et al., 2006). Indeed, when new environmental conditions occur in a biotope, transferability may imply extrapolation in the niche (i.e., Euclidean) space (Fernández & Hamilton, 2015). The ESM only provides abundance (or mean occurrence) estimates if there are ESs available in the environmental regime corresponding to the new area, either invaded by an exotic species or after a climate-induced biogeographical movement has taken place. The number of strings available for each geographical cell is known (Figure 3A). It is therefore possible to evaluate its influence on the estimates. The user may also choose to remove estimates associated with a low number of samples. In the present study, more than one ES within the TZ were used to assess the abundance of a species.

The model preserves the local relationships among environmental variables. Most models use a global assessment of the potential interactions among variables; these models are frequently based on correlation or regression coefficients assessed globally (Beaugrand et al., 2011; Hastie & Tibshirani, 1990; McCullagh & Nelder, 1983). However, the strength and sign of interactions among variables may vary in the Euclidean space of the niche. The ESM only considers the local relationships within the TZ using only ESs that are close to a TES (Figure 1). Note that there exist other methods that also locally weight estimations such as kernel

regression (Cleophas & Zwinderman, 2022; Wolberg, 2000) or Geographically Weighted Regression (GWR; Brunsdon et al., 1996). Kernel regression uses a bandwidth (i.e., a smoothing parameter) (Horová et al., 2012), which should not be strictly compared with parameter $d_j$ in my model because the former is a smoothing factor, whereas in my model, the smoothing is controlled by both $d_j$ and $v$. In the ESM, an increase in $d_j$ raises the smoothing, but this can be counterbalanced by a large weight $v$ in contrast to kernel regression (Horová et al., 2012). Other important parameters in kernel regression are (1) the kernel $K$, which is a weight function (e.g., Gaussian or triangle), and (2) the order of the kernel (Horová et al., 2012). The ESM could be viewed as similar to GWR, although GWR assesses changes in the local relationships among variables in the geographical and not the Euclidean space (Brunsdon et al., 1996), which differs if the relationships among variables change rapidly in space. In the ESM, the amount of changes in the relationships among variables is controlled by $d_j$ in the Euclidean space of the niche (Figure 1).

The ESM can deal with a high number of environmental dimensions (i.e., variables). To know the contribution of each variable, it is possible to run the ESM for all combinations of variables and retain the combination that best explains the observations; this can lead to high calculation times, but this way allows one to gain more confidence about the influence of variables into the model (Tables 2 and 3). In this paper, I have considered
<table><caption>TABLE 3 Best combinations of variables for the Environmental String Model applied from four to nine simultaneous variables.</caption><thead><tr><th>Performance metric</th><th>All variables</th><th>Eight variables</th><th>Seven variables</th><th>Six variables</th><th>Five variables</th><th>Four variables</th></tr></thead><tbody><tr><td><em>Calanus finmarchicus</em></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Variable</td><td>1-9</td><td>All but 1</td><td>2-8</td><td>2-6, 8</td><td>2-5, 8</td><td>2, 3, 4, 8</td></tr><tr><td>j</td><td>121,999</td><td>122,210</td><td>122,262</td><td>122,395</td><td>122,692</td><td>122,825</td></tr><tr><td>r</td><td>0.78<sup>a</sup></td><td>0.80<sup>a</sup></td><td>0.80<sup>a</sup></td><td>0.79<sup>a</sup></td><td>0.78<sup>a</sup></td><td>0.77<sup>a</sup></td></tr><tr><td>Bias</td><td>0.18</td><td>0.19</td><td>0.19</td><td>0.19</td><td>0.21</td><td>0.23</td></tr><tr><td>MAE</td><td>0.31</td><td>0.31</td><td>0.31</td><td>0.32</td><td>0.33</td><td>0.35</td></tr><tr><td>RMSE</td><td>0.21</td><td>0.20</td><td>0.20</td><td>0.20</td><td>0.22</td><td>0.23</td></tr><tr><td>PRMSE</td><td>45.24</td><td>43.47</td><td>43.68</td><td>44.36</td><td>46.68</td><td>49.95</td></tr><tr><td>Adj. R<sup>2</sup></td><td>0.61</td><td>0.64</td><td>0.64</td><td>0.62</td><td>0.61</td><td>0.59</td></tr><tr><td>AIC (×10<sup>5</sup>)</td><td>-3.81</td><td>-3.93</td><td>-3.93</td><td>-3.94</td><td>-3.71</td><td>-3.61</td></tr><tr><td><em>Calanus helgolandicus</em></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Variable</td><td>1-9</td><td>All but 8</td><td>1-6, 9</td><td>1-4, 6, 9</td><td>1-4, 9</td><td>2-4, 9</td></tr><tr><td>j</td><td>121,999</td><td>122,119</td><td>122,189</td><td>122,612</td><td>122,804</td><td>122,834</td></tr><tr><td>r</td><td>0.73<sup>a</sup></td><td>0.73<sup>a</sup></td><td>0.73<sup>a</sup></td><td>0.73<sup>a</sup></td><td>0.72<sup>a</sup></td><td>0.71<sup>a</sup></td></tr><tr><td>Bias</td><td>-0.03</td><td>-0.02</td><td>-0.02</td><td>-0.01</td><td>0.00</td><td>0.00</td></tr><tr><td>MAE</td><td>0.17</td><td>0.17</td><td>0.17</td><td>0.17</td><td>0.18</td><td>0.18</td></tr><tr><td>RMSE</td><td>0.09</td><td>0.09</td><td>0.09</td><td>0.09</td><td>0.09</td><td>0.09</td></tr><tr><td>PRMSE</td><td>32.76</td><td>32.39</td><td>32.35</td><td>32.44</td><td>33.00</td><td>33.00</td></tr><tr><td>Adj. R<sup>2</sup></td><td>0.53</td><td>0.53</td><td>0.53</td><td>0.53</td><td>0.52</td><td>0.50</td></tr><tr><td>AIC (×10<sup>5</sup>)</td><td>-5.87</td><td>-5.88</td><td>-5.88</td><td>-5.90</td><td>-5.91</td><td>-5.92</td></tr></tbody></table>

Note: All analyses were performed for *d*<sub>j</sub> = 0.05 and *v* = 10<sup>2</sup> with bathymetry, sea surface temperature, and downward solar radiation flux at surface.

Ref: Reference for *d*<sub>j</sub> = 0.05 and *v* = 10<sup>2</sup> and for all variables (see Table 1). 1, distance to nearest coasts; 2, bathymetry; 3, Sea Surface Temperature; 4, Downward solar radiation flux at surface; 5, nitrate; 6, phosphate; 7, silicate; 8, density mixed-layer depth; 9, sea surface salinity. In contrast to adjusted *R*<sup>2</sup>, a model is best when the AIC is the smallest.

Abbreviations: AIC, Akaike information criterion; MAE, mean absolute error; PRMSE, percentage RMSE; RMSE, root mean square error.

<sup>a</sup>Significant correlation.

that bathymetry, SST, and DRSF were key variables that influence the spatial distribution of the two *Calanus* after removing all variables individually from the model (Table 2) and by expert knowledge (Beaugrand et al., 2013, 2014; Helaouët & Beaugrand, 2007). In the second analysis, I investigated the influence of the variables by testing all potential combinations (i.e., 63 combinations; see also *Materials and methods*, Equation 14). A penalty algorithm could be implemented to select the best number of environmental descriptors using shrinkage methods such as the ridge or lasso (least absolute shrinkage and selection operator) regression (Harrell, 2001; Hastie et al., 2009; Tibshirani, 1996).

The ESM can also be applied on presence/absence data (Figure 7). A mean is performed on 0/1 values, each value being weighted by the distance to target string. Results should therefore be interpreted in terms of mean occurrence and not probability of occurrence, as in other models such as NPPEN (Beaugrand et al., 2011). Because the number of strings needed to have an assessment may be higher, the spatial coverage might be lower. The ESM

could therefore be applied on presence-only data after the generation of pseudo-absence (Barbet-Massin et al., 2012; Hefley & Hooten, 2016), although techniques such as ENFA (Hirzel et al., 2002) or NPPEN (Beaugrand et al., 2011) that specifically deal with presence-only data may be more appropriate. The data used here, which are available for testing, could be used to investigate model performance for presence-only data in depth because the CPR survey (Batten et al., 2003; Reid et al., 2003) is among rare large-scale monitoring programs that provide clear information on abundance, nil or positive, and could be used as a reference data against which all SDMs/ENMs could be correctly evaluated in the marine realm.

## Limitations, caveats, and caution when applying the ESM

As most SDMs/ENMs, my model has some limitations or potential caveats that may originate from conceptualization, data quality and preprocessing, collinearity, and
interpretation, to name a few (Araújo et al., 2019; Jarnevich et al., 2015; Leroy, 2023). The niche concept itself has been discussed (Chase & Leibold, 2003; Holt, 2009), and it is clear that it remains difficult to identify and implement biotic interactions that affect species abundance locally and the influence of dispersal/mobility that is rarely taken into account in ecological niche modeling (Soberon & Nakamura, 2009; Soberon & Peterson, 2005). These issues have been well summarized by the so-called BAM diagrams, where B stands for biotic, A stands for abiotic (scenopoetic variables), and M stands for mobility (or movement) (Soberon & Peterson, 2005). In the pelagic realm, it is possible that plankton species be expatriated well beyond their suitable area (Angel, 1993, 1997; van der Spoel, 1997), which may affect abundance and therefore model predictions regionally. This potential limitation, as well as biotic interaction, influences all ENMs/SDMs, however.

My approach remains dependent upon the quality of the environmental variables that form an ES. In the study, many environmental variables were not available on a monthly basis for the whole period covered by the CPR survey. I had to assume that some variables did not change on a year-to-year or decadal scale, which is obviously a wrong assumption. That is why I did not examine year-to-year changes in modeled species abundance. Fortunately, temperature, which is a key variable in biology (Beaugrand, 2015; Brown et al., 2004), was available at a good spatial resolution for the whole time period covered by the CPR sampling. Many environmental data measured from satellites have excellent spatiotemporal resolutions and could be used to improve the resolution of the CPR dataset. However, most satellite data are available from the end of 1997 (Alvain et al., 2008; Druon et al., 2019).

An important limitation in using the ESM is the standardization procedure of the environmental data between 0 and 1, which is applied prior to the calculation of the ESM (Equation 1). A user needs to use a large range of environmental gradients in the (multidimensional) niche space to be sure to run the model inside the range of environmental conditions experienced by the species in the past, present, and future. Two parameters also need a careful examination: distance to target string $d_j$ and maximum weight $v$ (Equations 3, 5, and 6). Although the values chosen here did not alter substantially the estimates (Table 1), it remains important to select the most optimal values for these parameters; the same was done for the spread coefficient in GRNN (Specht, 1991). A high distance to target string $d_j$ allows more ESs to be considered but increases the smoothing. In such a case, an increase in the weight gives more influence to the ES closer to the

TES and might also prevent underfitting. It is also recommended to increase $v$ or/and diminish distance to target string $d_j$ when the number of ESs is high within the TZ (Figure 3A). A small distance to target string, in addition to a strong weight, is not advised when the number of ESs is small because it may lead to overfitting. It is possible to choose a distance to target string $d_j$ specific for each dimension (or variable). In such a case however, the dimension with the higher distance will have more weight in the calculation of the abundance (Equation 5). This can be suitable in some circumstances but not in all cases; a standardization of individual distances $d_j$ prior to the calculation of $d_{tot}$ would ensure an equal contribution of all individual distances $d_j$. Last, weights are calculated by means of a negative exponential function ranging from $v$ (with $v \ge 1$, i.e., maximum weight when the ES and TES are identical) to 1 (i.e., when the multivariate distance between the ES and TES is the highest possible) within the TZ. The variation in the weights is therefore meaningful. However, other functions are possible, for example, linear or inverse squared distance method (Wackernagel, 1995).

Even if a large range of environmental variables was tested in this particular study, the study of the influence of the removal of one variable on model performance was influenced by collinearity (Table 2), for example, the removal of SST had a limited influence in the model because it was compensated by DRSF, both of which being highly correlated positively at a large scale. This issue has been well identified in ecological niche modeling (Dormann et al., 2013; Sillero & Barbosa, 2021), and some authors have advocated to keep only uncorrelated variables (Sillero & Barbosa, 2021). However, this might be a problem in situations where globally correlated variables diverge locally in some circumstances. For example, oxygen concentration and temperature are globally highly correlated negatively (Beaugrand et al., 2019). One could choose to remove oxygen from the model. However, if a new situation emerges and that oxygen diminishes due to eutrophication or an upwelling bringing anoxic waters in shallow coastal regions, for example (Diaz et al., 2010; Grantham et al., 2004), the absence of consideration into the model of the variable may strongly affect model prediction locally. DRSF and SST are correlated positively, but in the context of global climate change, SST is likely to increase in most regions while DRSF may respond regionally in different ways (Intergovernmental Panel on Climate Change, 2019). DRSF is a key variable (i.e., proxy of the photoperiod), perhaps involved in the diapause of *C. finmarchicus* (Häfker et al., 2018), therefore controlling the timing of its annual cycle, whereas SST influences the physiology and the ecology of active individuals at the surface.
These two examples illustrate the risk of removing variables on the basis of global correlation indices. This might also be the case with some other variables (e.g., distance to coast and bathymetry).

In the present study, temporal data splitting (i.e., twofold cross-validation using a block in 1958–1999 for training and 2000–2016 for testing) was applied to test and compare the performance of my model with others. The central assumption of data splitting is that training and evaluation data are independent (Roberts et al., 2017). If this assumption is violated, error estimates are too minimized. Although this might indeed be the case in this study because of spatial autocorrelation, this has unlikely altered the comparison of the performance of the different models. However, if the ESM is applied in ecological niche modeling, the method of *k*-fold cross-validation (with *k* > 2) should be preferred (Roberts et al., 2017).

## Application of the model to CPR data

Since the beginning of the CPR survey, an important amount of biological samples has been collected on a monthly basis. The number of samples collected for the period 1958–2016 was 234,486 in the spatial domain 99.5° W–19.5° E and 20.5° N–79.5° N. However, at the scale of the North Atlantic Ocean, this represents only an average of 234,486/(12 months × 59 years) = 331.19 samples per month for a surface that covers million square kilometers. Because sampling is also irregular in space and time, most analyses performed on CPR data have regularized species abundance by calculating arithmetic or geometric average in a given geographical cells or by spatial interpolation (Beaugrand et al., 2003; Fromentin & Planque, 1996; Helaouët & Beaugrand, 2007; Kléparski et al., 2021; Planque et al., 1997). At a monthly scale, regularization or spatial interpolation is only valid for some areas with a sufficient spatial coverage in CPR samples such as the North Sea. In many areas of the North Atlantic Ocean, many gaps remain at a monthly scale (Beaugrand et al., 2019). My approach could be used to better fill gaps in the dataset at a monthly scale with a better control on the quality of the estimates. Finally, I think that the ESM could also be used to fill gaps in large remote sensing datasets such as ocean color or PAR (Ruddick et al., 2014).

## ACKNOWLEDGMENTS

I thank the owners and crews who have towed the CPRs on a voluntary basis for over 80 years, contributing to one of the world's largest and longest ongoing ecological experiments. Without these early pioneers of citizen science and broadscale volunteer monitoring projects,

this unique ecological dataset would never have been financially or logistically viable. The author thanks the National Centre for Scientific Research (CNRS) and the Continuous Plankton Recorder (CPR) survey. The author is indebted to Roy Elhourany for advice on generalized regression neural network, Jérôme Castant for performance metrics, and Eric Goberville for discussion about ecological niche models. The author also thanks the two anonymous reviewers who contributed to improving the quality of the research paper. This work is part of the Graduate school IFSEA that benefits from grant ANR-21-EXES-0011 from the French National Research Agency, under the Investments for the Future Programme France 2030, and the French Region Hauts-de-France and Ifremer in the framework of the project CPER IDEAL 2021–2027.

## CONFLICT OF INTEREST STATEMENT

The author declares no conflicts of interest.

## DATA AVAILABILITY STATEMENT

Data (Beaugrand, 2023a) and code (Beaugrand, 2023b) are available from Figshare: https://doi.org/10.6084/m9.figshare.22638808 and https://doi.org/10.6084/m9.figshare.22638766.v1, respectively.

## ORCID

Grégory Beaugrand https://orcid.org/0000-0002-0712-5223

## REFERENCES

Akaike, H. 1974. "A New Look at the Statistical Model Identification." *IEEE Transactions on Automatic Control* 19: 716-723.

Allouche, O., A. Tsoar, and R. Kadmon. 2006. "Assessing the Accuracy of Species Distribution Models: Prevalence, Kappa and the True Skill Statistic (TSS)." *Journal of Applied Ecology* 43: 1223-32.

Alvain, S., C. Moulin, Y. Dandonneau, and H. Loisel. 2008. "Seasonal Distribution and Succession of Dominant Phytoplankton Groups in the Global Ocean: A Satellite View." *Global Biogeochemical Cycles* 22:GB3001.

Angel, M. V. 1993. "Biodiversity of the Pelagic Ocean." *Conservation Biology* 7: 760-772.

Angel, M. V. 1997. "Pelagic Biodiversity." In *Marine Biodiversity: Patterns and Processes*, edited by R. F. G. Ormond, J. D. Gage, and M. V. Angel, 35-68. Cambridge: Cambridge University Press.

Araújo, M. B., R. P. Anderson, A. M. Barbosa, C. M. Beale, C. F. Dormann, R. Early, R. A. Garcia, et al. 2019. "Standards for Distribution Models in Biodiversity Assessments." *Science Advances* 5: eaat4858.

Austin, M. P., and J. A. Meyers. 1996. "Current Approaches to Modelling the Environmental Niche of Eucalypts: Implication for Management of Forest Biodiversity." *Forest Ecology and Management* 85: 95-106.
Barbet-Massin, M., F. Jiguet, C. H. Albert, and W. Thuiller. 2012. "Selecting Pseudo-Absences for Species Distribution Models: How, Where and How Many?" *Methods in Ecology and Evolution* 3: 327-338.

Batten, S. D., R. Clark, J. Flinkman, G. Hays, E. John, A. W. G. John, T. Jonas, J. A. Lindley, D. P. Stevens, and A. Walne. 2003. "CPR Sampling: The Technical Background, Materials, and Methods, Consistency and Comparability." *Progress in Oceanography* 58: 193-215.

Beaugrand, G. 2004. "Monitoring Marine Plankton Ecosystems (1): Description of an Ecosystem Approach Based on Plankton Indicators." *Marine Ecology Progress Series* 269: 69-81.

Beaugrand, G. 2015. *Marine Biodiversity, Climatic Variability and Global Change*. London: Routledge.

Beaugrand, G. 2023a. "Data Matrices for Ecological String Model." Figshare. https://doi.org/10.6084/m9.figshare.22638808.v1.

Beaugrand, G. 2023b. "MATLAB Functions Ecological String Model." Figshare. https://doi.org/10.6084/m9.figshare.22638766.v1.

Beaugrand, G., and M. Edwards. 2001. "Comparison in Performance among Four Indices Used to Evaluate Diversity in Pelagic Ecosystems." *Oceanologica Acta* 24: 467-477.

Beaugrand, G., M. Edwards, and P. Helaouët. 2019. "An Ecological Partition of the Atlantic Ocean and Its Adjacent Seas." *Progress in Oceanography* 173: 86-102.

Beaugrand, G., F. Ibañez, and J. A. Lindley. 2003. "An Overview of Statistical Methods Applied to the CPR Data." *Progress in Oceanography* 58: 235-262.

Beaugrand, G., S. Lenoir, F. Ibañez, and C. Manté. 2011. "A New Model to Assess the Probability of Occurrence of a Species Based on Presence-Only Data." *Marine Ecology Progress Series* 424: 175-190.

Beaugrand, G., C. Luczak, and M. Edwards. 2009. "Rapid Biogeographical Plankton Shifts in the North Atlantic Ocean." *Global Change Biology* 15: 1790-1803.

Beaugrand, G., D. Mackas, and E. Goberville. 2013. "Applying the Concept of the Ecological Niche and a Macroecological Approach to Understand How Climate Influences Zooplankton: Advantages, Assumptions, Limitations and Requirements." *Progress in Oceanography* 111: 75-90.

Beaugrand, G., P. C. Reid, F. Ibañez, and P. Planque. 2000. "Biodiversity of North Atlantic and North Sea Calanoid Copepods." *Marine Ecology Progress Series* 204: 299-303.

Boria, R. A., L. E. Olson, S. M. Goodman, and R. P. Anderson. 2014. "Spatial Filtering to Reduce Sampling Bias Can Improve the Performance of Ecological Niche Models." *Ecological Modelling* 275: 73-77.

Brown, J. H., J. F. Gillooly, A. P. Allen, V. M. Savage, and G. B. West. 2004. "Toward a Metabolic Theory of Ecology." *Ecology* 85: 1771-89.

Bruno, J. F., J. J. Stachowicz, and M. D. Bertness. 2003. "Inclusion of Facilitation into Ecological Theory." *Trends in Ecology & Evolution* 18: 119-125.

Brunsdon, C., A. S. Fotheringham, and M. E. Charlton. 1996. "Geographically Weighted Regression: A Method for Exploring Spatial Nonstationarity." *Geographical Analysis* 28: 281-298.

Cade, B. S., and B. R. Noon. 2003. "A Gentle Introduction to Quantile Regression for Ecologists." *Frontiers in Ecology and the Environment* 1: 412-420.

Cade, B. S., B. R. Noon, and C. H. Flather. 2005. "Quantile Regression Reveals Hidden Bias and Uncertainty in Habitat Models." *Ecology* 86: 786-800.

Chase, J. M., and M. A. Leibold. 2003. *Ecological Niches: Linking Classical and Contemporary Approaches*. Chicago: The University of Chicago.

Chen, F.-W., and C.-W. Liu. 2012. "Estimation of the Spatial Rainfall Distribution Using Inverse Distance Weighting (IDW) in the Middle of Taiwan." *Paddy and Water Environment* 10: 209-222.

Chevalier, M., A. Zarzo-Arias, J. Guélat, R. G. Mateo, and A. Guisan. 2022. "Accounting for Niche Truncation to Improve Spatial and Temporal Predictions of Species Distributions." *Frontiers in Ecology and Evolution* 10: 944116.

Cleophas, T. J., and A. H. Zwinderman. 2022. *Kernel Ridge Regression in Clinical Research*. Cham: Springer.

Colwell, R. K., and T. F. L. V. B. Rangel. 2009. "Hutchinson's Duality: The Once and Future Niche." *Proceedings of the National Academy of Sciences of the United States of America* 106: 19651-58.

Crisp, M. D., M. T. K. Arroyo, L. G. Cook, M. A. Gandolfo, G. J. Jordan, M. S. McGlone, P. H. Weston, M. Westoby, P. Wilf, and H. P. Linder. 2009. "Phylogenetic Biome Conservatism on a Global Scale." *Nature* 458: 754-56.

Diaz, R. J., M. Selman, and C. Chique. 2010. *Global Eutrophic and Hypoxic Coastal Systems. Eutrophication and Hypoxia: Nutrient Pollution in Coastal Waters*. Washington D.C: World Resources Institute.

Dormann, C. F., J. Elith, S. Bacher, C. Buchmann, G. Carl, G. Carré, J. R. G. Marquéz, et al. 2013. "Collinearity: A Review of Methods to Deal with it and a Simulation Study Evaluating Their Performance." *Ecography* 36: 27-46.

Druon, F.-N., P. Helaouët, G. Beaugrand, J.-M. Fromentin, A. Palialexis, and N. Hoepffner. 2019. "Satellite-Based Indicator of Zooplankton Distribution for Global Monitoring." *Scientific Reports* 9: 4732.

Durant, J. M., D. Hjermann, T. Anker-Nielssen, G. Beaugrand, A. Mysterud, N. Pettorelli, and N. C. Stenseth. 2005. "Timing and Abundance as Key Factors Affecting the Climate-Trophic Interaction Relationship." *Ecology Letters* 8: 952-58.

Fernández, M., and H. Hamilton. 2015. "Ecological Niche Transferability Using Invasive Species as a Case Study." *PLoS One* 10: e0119891.

Fromentin, J. M., and B. Planque. 1996. "Calanus and Environment in the Eastern North Atlantic. II. Influence of the North Atlantic Oscillation on *C. finmarchicus* and *C. helgolandicus*." *Marine Ecology Progress Series* 134: 111-18.

Frontier, S., and D. Pichot-Viale. 1993. *Ecosystèmes: Structure, Fonctionnement, Évolution*. Paris: Masson.

Grantham, B. A., F. Chan, K. J. Nielsen, D. S. Fox, J. A. Barth, A. Huyer, J. Lubchenco, and B. A. Menge. 2004. "Upwelling-Driven Nearshore Hypoxia Signals Ecosystem and Oceanographic Changes in the Northeast Pacific." *Nature* 429: 749-754.

Guisan, A., T. C. Edwards, Jr., and T. Hastie. 2002. "Generalized Linear and Generalized Additive Models in Studies of Species Distributions: Setting the Scene." *Ecological Modelling* 157: 89-100.

Häfker, N. S., M. Teschke, K. S. Last, D. W. Pond, L. Hüppe, and B. Meyer. 2018. "Calanus finmarchicus Seasonal Cycle and
Diapause in Relation to Gene Expression, Physiology, and Endogenous Clocks." *Limnology and Oceanography* 63: 2815-38.

Harrell, F. E., Jr. 2001. *Regression Modeling Strategies with Applications to Linear Models, Logistic Regression, and Survival Analysis*. New York: Springer.

Hastie, T. J., R. Tibshirani, and A. Buja. 1994. "Flexible Discriminant Analysis by Optimal Scoring." *Journal of the American Statistical Association* 89: 1255-70.

Hastie, T. J., R. Tibshirani, and J. Friedman. 2009. *The Elements of Statistical Learning: Data Mining, Inference, and Prediction*. New York: Springer.

Hastie, T. J., and R. J. Tibshirani. 1990. *Generalized Additive Models*. London: Chapman and Hall.

Hefley, T. J., and M. B. Hooten. 2016. "Hierarchical Species Distribution Models." *Current Landscape Ecology Reports* 1: 87-97.

Helaouët, P., and G. Beaugrand. 2007. "Macroecology of *Calanus finmarchicus* and *C. helgolandicus* in the North Atlantic Ocean and Adjacent Seas." *Marine Ecology Progress Series* 345: 147-165.

Helaouët, P., and G. Beaugrand. 2009. "Physiology, Ecological Niches and Species Distribution." *Ecosystems* 12: 1235-45.

Hirzel, A. H., J. Hausser, D. Chessel, and N. Perrin. 2002. "Ecological-Niche Factor Analysis: How to Compute Habitat-Suitability Maps without Absence Data?" *Ecology* 83: 2027-36.

Hirzel, A. H., and G. Le Lay. 2008. "Habitat Suitability Modelling and Niche Theory." *Journal of Applied Ecology* 45: 1372-81.

Hirzel, A. H., G. Le Lay, V. Helfer, C. Randin, and A. Guisan. 2006. "Evaluating the Ability of Habitat Suitability Models to Predict Species Presences." *Ecological Modelling* 199: 142-152.

Holt, R. D. 2009. "Bringing the Hutchinsonian Niche into the 21st Century: Ecological and Evolutionary Perspectives." *Proceedings of the National Academy of Sciences of the United States of America* 106: 19659-65.

Horová, I., J. Koláček, and J. Zelinka. 2012. *Kernel Smootinh in Matlab: Theory and Practice of Kernel Smoothing*. Hackensack, NJ: World Scientific.

Hutchinson, G. E. 1957. "Concluding Remarks." *Cold Spring Harbor Symposia on Quantitative Biology* 22: 415-427.

Hutchinson, G. E. 1978. *An Introduction to Population Ecology*. New Haven, CT: Yale University Press.

Intergovernmental Panel on Climate Change. 2019. "Summary for Policymakers." In *IPCC Special Report on the Ocean and Cryosphere in a Changing Climate*, edited by H. O. Pörtner, D. C. Roberts, V. Masson-Delmotte, P. Zhai, M. Tignor, E. Poloczanska, K. Mintenbeck, M. Nicolai, A. Okem, J. Petzold, B. Rama, and N. Weyer, 35. Cambridge: Cambridge University Press.

Jarnevich, C. S., T. J. Stohlgren, S. Kumar, J. T. Morisette, and T. R. Holcombe. 2015. "Caveats for Correlative Species Distribution Modeling." *Ecological Informatics* 29: 6-15.

Jonas, T. D., A. Walne, G. Beaugrand, L. Gregory, and G. C. Hays. 2004. "The Volume of Water Filtered by a CPR: The Effect of Ship Speed." *Journal of Plankton Research* 26: 1499-1506.

Kim, P. 2017. *MATLAB Deep Learning: With Machine Learning, Neural Networks and Artificial Intelligence*. Seoul: APRESS.

Kléparski, L., and G. Beaugrand. 2022. "The Species Chromatogram, A New Graphical Method to Represent,

Characterise and Compare the Ecological Niches of Different Species." *Ecology and Evolution* 12: e8830.

Kléparski, L., G. Beaugrand, and M. Edwards. 2021. "Plankton Biogeography in the North Atlantic Ocean and Its Adjacent Seas: Species Assemblages and Environmental Signatures." *Ecology and Evolution* 11: 5135-49.

Koenker, R., and G. Bassett. 1978. "Regression Quantiles." *Econometrica* 46: 33-50.

Legendre, P., and L. Legendre. 1998. *Numerical Ecology*, 2nd ed. Amsterdam: Elsevier Science B.V.

Lehmann, A., J. M. Overton, and J. R. Leathwick. 2002. "GRASP: Generalized Regression Analysis and Spatial Prediction." *Ecological Modelling* 157: 189-207.

Leroy, B. 2023. "Choosing Presence-Only Species Distribution Models." *Journal of Biogeography* 50: 247-250.

Longhurst, A. 1998. *Ecological Geography of the Sea*. London: Academic Press.

Mateo, R. G., M. J. Aroca-Fernandez, A. Gaston, V. Gomez-Rubio, S. Saura, and J. I. Garcia-Vinas. 2019. "Looking for an Optimal Hierarchical Approach for Ecologically Meaningful Niche Modelling." *Ecological Modelling* 409: 108735.

McCullagh, P., and J. A. Nelder. 1983. *Generalized Linear Models*. London: Chapman and Hall.

Melo-Merinoa, S. M., H. Reyes-Bonillab, and A. Lira-Noriegac. 2020. "Ecological Niche Models and Species Distribution Models in Marine Environments: A Literature Review and Spatial Analysis of Evidence." *Ecological Modelling* 415: 108837.

Ostertagová, A. 2012. "Modelling Using Polynomial Regression." *Procedia Engineering* 48: 500-506.

Peterson, A. T. 2011. "Ecological Niche Conservatism: A Time-Structured Review of Evidence." *Journal of Biogeography* 38: 817-827.

Peterson, A. T., and J. Soberon. 2012. "Species Distribution Modeling and Ecological Niche Modeling: Getting the Concepts Right." *Natureza & Conservação* 10: 1-6.

Peterson, A. T., J. Soberon, and V. Sanchez-Cordero. 1999. "Conservatism of Ecological Niches in Evolutionary Time." *Science* 285: 1265-67.

Philips, S. J., R. P. Anderson, and R. E. Shapire. 2006. "Maximum Entropy Modeling of Species Geographic Distributions." *Ecological Modelling* 190: 231-259.

Planque, B., G. C. Hays, F. Ibañez, and J. C. Gamble. 1997. "Large Scale Spatial Variations in the Seasonal Abundance of *Calanus finmarchicus*." *Deep Sea Research, Part I* 44: 315-326.

Pulliam, H. R. 1988. "Sources, Sinks, and Population Regulation." *The American Naturalist* 132: 652-661.

Pulliam, R. 2000. "On the Relationship between Niche and Distribution." *Ecology Letters* 3: 349-361.

Qiao, H., X. Feng, L. E. Escobar, A. T. Peterson, J. Soberón, G. Zhu, and M. Papeş. 2019. "An Evaluation of Transferability of Ecological Niche Models." *Ecography* 42: 521-534.

Randin, C. F., T. Dirnböck, S. Dullinger, N. E. Zimmermann, M. Zappa, and A. Guisan. 2006. "Are Niche-Based Species Distribution Models Transferable in Space?" *Journal of Biogeography* 33: 1689-1703.

Raybaud, V., G. Beaugrand, J.-M. Dewarumez, and C. Luczak. 2015. "Climate-Induced Range Shifts of the American Jackknife Clam *Ensis Directus* in Europe." *Biological Invasions* 17: 725-741.
Raybaud, V., G. Beaugrand, E. Goberville, G. Delebecq, C. Destombe, M. Valero, D. Davoult, P. Morin, and F. Gevaert. 2013. "Decline in Kelp in West Europe and Climate." *PLoS One* 8: e66044.

Reid, P. C., J. M. Colebrook, J. B. L. Matthews, J. Aiken, R. Barnard, S. D. Batten, G. Beaugrand, et al. 2003. "The Continuous Plankton Recorder: Concepts and History, from Plankton Indicator to Undulating Recorders." *Progress in Oceanography* 58: 117-173.

Roberts, D. R., V. Bahn, S. Ciuti, M. S. Boyce, J. Elith, G. Guillera-Arroita, S. Hauenstein, et al. 2017. "Cross-Validation Strategies for Data with Temporal, Spatial, Hierarchical, or Phylogenetic Structure." *Ecography* 40: 913-929.

Rombouts, I., G. Beaugrand, and J.-C. Dauvin. 2012. "Potential Changes in Benthic Macrofaunal Distributions from the English Channel Simulated under Climate Change Scenarios." *Estuarine, Coastal and Shelf Science* 99: 153-161.

Ruddick, K., G. Neukermans, Q. Vanhellemont, and D. Jolivet. 2014. "Challenges and Opportunities for Geostationary Ocean Colour Remote Sensing of Regional Seas: A Review of Recent Results." *Remote Sensing of Environment* 146: 63-76.

Saupe, E. E., V. Barve, C. E. Myers, J. Soberon, N. Barve, C. M. Hensz, A. T. Peterson, H. L. Owens, and A. Lira-Noriega. 2012. "Variation in Niche and Distribution Model Performance: The Need for a Priori Assessment of Key Causal Factors." *Ecological Modelling* 237-238: 11-22.

Scherrer, D., M. Esperon-Rodriguez, L. J. Beaumont, V. L. Barradas, and A. Guisan. 2021. "National Assessments of Species Vulnerability to Climate Change Strongly Depend on Selected Data Sources." *Diversity and Distributions* 27: 1367-82.

Schickele, A., E. Goberville, B. Leroy, G. Beaugrand, T. Hattab, P. Francour, and V. Raybaud. 2021. "Redistribution of Small Pelagic Fish in Europe and Climate Change." *Fish and Fisheries* 22: 212-225.

Schickele, A., B. Leroy, G. Beaugrand, E. Goberville, T. Hattab, P. Francour, and V. Raybaud. 2020. "Modelling European Small Pelagic Fish Distribution: Methodological Insights." *Ecological Modelling* 416: 108902.

Sillero, N., and A. M. Barbosa. 2021. "Caveats for Correlative Species Distribution Modeling." *International Journal of Geographical Information Science* 35: 213-226.

Soberon, J., and M. Nakamura. 2009. "Niches and Distributional Areas: Concepts, Methods, and Assumptions." *Proceedings of the National Academy of Sciences of the United States of America* 106: 19644-50.

Soberon, J., and A. T. Peterson. 2005. "Interpretation of Models of Fundamental Ecological Niches and Species' Distributional Areas." *Biodiversity Informatics* 2: 1-10.

Sokal, R. R., and F. J. Rohlf. 1995. *Biometry*. New York: W.H. Freeman and Company.

Specht, D. F. 1991. "A General Regression Neural Network." *IEEE Transactions on Neural Networks* 2: 568-576.

Thuiller, W., B. Lafourcade, R. Engler, and M. B. Araújo. 2009. "BIOMOD – A Platform for Ensemble Forecasting of Species Distributions." *Ecography* 32: 369-373.

Tibshirani, R. 1996. "Regression Shrinkage and Selection Via the Lasso." *Journal of the Royal Statistical Society: Series B: Methodological* 58: 267-288.

Tredennick, A. T., G. Hooker, S. P. Ellner, and P. B. Adler. 2021. "A Practical Guide to Selecting Models for Exploration, Inference, Andprediction in Ecology." *Ecology* 102: e03336.

Valavi, R., G. Guillera-Arroita, J. J. Lahoz-Monfort, and J. Elith. 2021. "Predictive Performance of Presence-Only Species Distribution Models: A Benchmark Study with Reproducible Code." *Ecological Monographs* 92: e01486.

van der Spoel, S. 1997. "Pelagic Biogeography and Temperature." *Annales de l'Institut Océanographique* 73: 185-194.

Vaz, S., C. S. Martin, P. D. Eastwood, B. Ernande, A. Carpentier, G. J. Meaden, and F. Coppin. 2008. "Blackwell Publishing Ltd Modelling Species Distributions Using Regression Quantiles." *Journal of Applied Ecology* 45: 204-217.

Wackernagel, H. 1995. *Multivariate Geostatistics. An Introduction with Applications*. Berlin: Springer-Verlag Berlin Heidelberg.

Warner, A. J., and G. C. Hays. 1994. "Sampling by the Continuous Plankton Recorder Survey." *Progress in Oceanography* 34: 237-256.

Wlodarczak, P. 2020. *Machine Learning and Its Applications*. Boca Raton: CRC Press.

Wolberg, J. R. 2000. *Expert Trading Systems: Modeling Financial Markets with Kernel Regression*. New York: John Wiley & Sons, Inc.

Yin, L., R. Caruana, and J. Gehrke. 2012. "Intelligible Models for Classification and Regression." In *Proceedings of the 18th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (KDD'12)*, 150-58. Beijing: ACM Press.

Yin, L., R. Caruana, J. Gehrke, and G. Hooker. 2013. "Accurate Intelligible Models with Pairwise Interactions." In *Proceedings of the 19th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (KDD'13)*, 623-31. Chicago: ACM Press.

Zhang, L., F. Huettmann, S. Liu, P. Sun, Z. Yu, X. Zhang, and C. Mi. 2019. "Classification and Regression with Random Forests as a Standard Method for Presence-Only Data SDMs: A Future Conservation Example Using China Tree Species." *Ecological Informatics* 52: 46-56.

Zhang, L., L. Shirong, P. Sun, T. Wang, G. Wang, X. Zhang, and L. Wang. 2015. "Consensus Forecasting of Species Distributions: The Effects of Niche Model Performance and Niche Properties." *PLoS One* 10: e0120056.

**How to cite this article:** Beaugrand, Grégory. 2024. "An Ecological Niche Model that Considers Local Relationships among Variables: The Environmental String Model." *Ecosphere* 15(10): e70015. https://doi.org/10.1002/ecs2.70015