Classic or hybrid? The performance of next generation ecological models to study the response of Southern Ocean species to changing environmental conditions

In the context of intensifying threats of climate change on marine communities, ecological models are widely applied for conservation strategies, though polar studies remain scarce given the limited number of datasets available. Correlative (e.g. species distribution models, SDM) and mechanistic (e.g. dynamic energy budget models, DEB) modelling approaches are usually used independently in studies. Using both approaches in integrative, hybrid models could help to better estimate the species potential ecological niche, as mechanistic and correlative models complement each other very well, giving more insights into species potential response to fast‐changing environmental conditions.


| INTRODUC TI ON
For the last two decades, an ever-growing number of ecological studies have used modelling approaches to highlight the main ecological drivers of species distribution and evaluate the response of species to changing environmental conditions and anthropogenic stressors (Elith et al., 2006;Franklin, 2009).
The overall tendency is to use these models across groups of organisms and regions (Gutt et al., 2012) to inform stakeholders and conservation policies (Mouquet et al., 2015;Singer et al., 2016;Thuiller et al., 2013).
Current developments are focussed on the integration of distinct modelling methods (i.e. hybrid modelling) that has long been considered as a way to improve the understanding of ecosystem functioning (Benito Garzón et al., 2019;Dormann et al., 2018;Guillaumot, Fabri-Ruiz, et al., 2018;Gutt et al., 2012). For instance, combining correlative methods, which rely on the spatial relationship between species occurrence records and the environment (e.g. Species Distribution Models, SDMs), with ecophysiological approaches (e.g. mechanistic models) was shown to improve the modelling performance compared with single correlative methods (Elith et al., 2010;Pertierra et al., 2019;Schouten et al., 2020;Singer et al., 2016). Correlative models statistically assess the main drivers of species distribution (Elith et al., 2006;Peterson et al., 2011) and are used to estimate the potential ecological niche Soberón, 2010). As a consequence, SDMs perform well when species distribution and the environment are in equilibrium, in static systems, a prerequisite that is not verified in highly dynamic ecosystems subject to environmental changes or in studies addressing environmental rapid changes (Fabri-Ruiz et al., 2021;Loehle & Leblanc, 1996;Schouten et al., 2020).
Mechanistic models can evaluate the effect of environmental conditions on the physiological performance of individuals or populations (Kearney & Porter, 2009). Such models typically require a greater level of biological knowledge, but, in contrast to static, correlative approaches, they explicitly include dynamic processes, offering the opportunity to describe process-based causes of species distribution change (Dormann et al., 2012;Kearney & Porter, 2009), even in nonequilibrium systems (Kearney et al., 2008;Keith et al., 2008).
They include a set of mathematical functions relating to species' functional traits (morphology, behaviour, physiology) or associated life history (development, growth, reproduction) and then evaluate the effect of environmental drivers on species physiological traits (Dormann et al., 2012;Kearney & Porter, 2009), which leads to estimating the species' fundamental niche (Kearney & Porter, 2009).
Several methods have been developed to integrate correlative and mechanistic models. For instance, mechanistic models can be spatially-projected and used as a input predictor in SDMs (Buckley et al., 2011;Elith et al., 2010;Mathewson et al., 2016;Rodríguez et al., 2019). Other close approaches consist in defining absence records from the mechanistic model and use the set of presenceabsence records to implement SDMs (Elith et al., 2010;Feng & Papeş, 2017) or to fine-tune thresholds for lethal conditions from the mechanistic approach and associate uncertainty estimates to SDM predictions accordingly (Woodin et al., 2013). Bayesian inference methods have also been widely used (Brewer et al., 2016;Ellison, 2004;Feng et al., 2020;Gamliel et al., 2020;Talluto et al., 2016), following the development and better accessibility of high-performance computers and programs (Van Dongen, 2006) and the development of more complex models (de Rivera et al., 2019).
They were proved interesting to optimize the estimation of species habitat suitability (Zurell et al., 2016), to better assess the effect of seasonality in predictions and highlight critical tipping points in changing ecosystems (Oberle et al., 2016;Zhao et al., 2019) providing accurate uncertainty estimates (Zhao et al., 2019). Bayesian methods combine the information of a prior belief (i.e. the prior distribution, for instance, our knowledge of species physiology) with new information (i.e. the conditional probability given the data) to produce a posterior estimation (Van Dongen, 2006). These two steps therefore update the probability of the hypothetical distribution as more evidence or information on species physiology is available (Van Dongen, 2006).
Many regions of the Southern Ocean, either in Antarctic or sub-Antarctic zones (Convey et al., 2009;Féral et al., 2017), are currently exposed to fast environmental changes (Convey & Peck, 2019;Cook et al., 2016;Turner et al., 2016), including increasing seawater temperatures and shifting seasonality (Bers et al., 2013;Henley et al., 2019;Schofield et al., 2017); glacier melting, changing wind speed (Cook et al., 2016;Meredith & King, 2005), which in turn have an impact on food chains, organic matter production and processes of the benthopelagic coupling (see Convey & Peck, 2019;Henley et al., 2019 as reviews). Climate changes together with the everincreasing maritime traffic (i.e. fisheries, tourism and science) boost the introduction of non-native species in Southern Ocean coastal areas, a major threat to polar ecosystems usually characterized by high levels of endemic species (Hughes et al., 2020;McCarthy et al., 2019). These combined issues strongly urge the need to fill the gaps in our knowledge of ecological processes and ecosystem dynamics (Kennicutt et al., 2015). a broad panel of natural examples, noteworthy for decision-making and conservation management purposes.

K E Y W O R D S
Bayesian inference, data-poor systems, integrated approaches, Kerguelen Islands, sea urchin, species distribution modelling Due to remoteness and harsh weather conditions, above all in winter, access to the field and data collection in the Southern Ocean are strongly limited (De Broyer et al., 2014), resulting in missing data, spatial and temporal aggregations of observations and difficulties to conduct biological experiments (see Guillaumot et al., 2021 as a review). However, research on marine life of the Southern Ocean has recently benefited from a significant coordinated and international effort with the emergence of oceanographic campaigns and international scientific programs such as the International Polar Year (IPY 2007(IPY -2008, the Census of Antarctic Marine Life (CAML 2005(CAML -2010 or the Scientific Committee on Antarctic Research, Evolution and Biodiversity in Antarctica (SCAR-EBA 2006 (nonexhaustive list) (De Broyer et al., 2014;Schiaparelli et al., 2013).
Several studies have used correlative approaches to characterize the relationship between environmental conditions and the distribution of Southern Ocean species (Bombosch et al., 2014;Fabri-Ruiz et al., 2019;Freer, 2018;Pinkerton et al., 2010) or used physiological models to evaluate the influence of environmental conditions on organisms' physiological performances (Agüera et al., 2015;Agüera et al., 2017;Jager & Ravagnan, 2015) and population dynamics (Arnould-Pétré et al., 2020;Goedegebuure et al., 2018;Groeneveld et al., 2015). However, surprisingly, no study has used integrated modelling approaches despite their considerable potential for analysing dynamic, complex and ill-known systems.
In this study, we used data from ongoing research on a sea urchin species, Abatus cordatus (Verrill, 1876), in the Baie du Morbihan, the most visited area of the otherwise highly remote archipelago of the Kerguelen Islands (French sub-Antarctic islands). We tested the performance of integrated modelling approaches to deal with a study on a Southern Ocean marine species and compared model outputs with other 'classic' correlative (SDM) and mechanistic (Dynamic Energy Budgets) approaches. In addition, we assessed the effect of environmental changes (i.e. related to two periods with contrasting conditions), a fundamental feature of ecosystem functioning in high latitudes and a key to understand the functioning of marine life in the Southern Ocean. Dealing with two different dates was here chosen to test the performance of different modelling procedures in an environmentally dynamic context.

| Overview/conceptualization
In this study, we aim at improving correlative models ('Classic' Species Distribution Modelling, SDM) with physiological information for a sub-Antarctic marine species example. This physiological information is integrated as an a priori knowledge in the SDM and represented by a physiological submodel that describes the relationship between growth performance and food availability. This physiological submodel was calculated from the DEB model of the species.
The a priori information was integrated inside the SDM using two approaches: (1) by integrating the spatial projection of the DEB model as an environmental predictor in the SDM or (2) using a Bayesian inference procedure. For the latter, the parameters that fit the physiological submodel were added to the initial matrix of parameters that compose the SDM.
Methods are detailed step by step in the following paragraphs and are completed with the summarizing flowchart in Figure 1. R codes developed for this study are available at https://github.com/ charl enegu illau mot/THESIS.

| Study species
In this study, the heart urchin A. cordatus (Verrill, 1876) was selected as the study example as it constitutes a relatively well-documented species compared with other Southern Ocean species. A. cordatus is a shallow deposit-feeder and sediment swallower restrained to soft sediment habitats (De Ridder & Lawrence, 1982;Poulin, 1996) (Appendix S1). Endemic to the Kerguelen Plateau, the species is distributed from shallow subtidal (<2 m depth) to deep shelf areas exceeding 500 m depth (Poulin, 1996). In coastal zones, populations of A. cordatus can locally reach densities of up to 280 individuals per square meter (Magniez, 1980;Poulin, 1996). High population densities along with the species endemicity were interpreted as a consequence of the species reproduction strategy and direct development that includes no dispersal larval stage (Mespoulhé, 1992;Poulin & Féral, 1995). Females brood their young on the aboral side of the test, inside four brood chambers formed by the sunken paired ambulacra, until juveniles exit the pouch and reach the sea bottom in the proximity of their mothers (Appendix S1, Magniez, 1983). In most places of the Baie du Morbihan, individuals invest energy into the growth of gonads in March, when food is the most abundant (Magniez, 1983). Once fertilized, the eggs are brooded in the female incubating chambers for almost 9 months (a period of low food availability and low temperature) before the young are released and settle on the seabed (Schatt & Féral, 1996) or live sheltered between holdfasts of the giant kelp Macrocystis pyrifera (Poulin, 1996).
The reproduction cycle of A. cordatus is constant across years for a given place (Magniez, 1983). However, it was observed that the reproduction period can shift from a few months between sites (Mespoulhé, 1992;Poulin, 1996;Schatt & Féral, 1991), which was explained by spatial and temporal variations in food availability and sediment enrichment in nutrients (Schatt & Féral, 1991).
Depth, temperature and primary production were identified as major environmental drivers of the distribution of A. cordatus (Poulin, 1996). Sediment granulometry and hydrodynamics were also shown to be important drivers of population densities in A. cordatus (Poulin & Féral, 1995). These two key factors cannot, however, be included in our models as they are not available for the Baie du Morbihan. In shallow-water areas, the species was shown to be tolerant to environmental stressors induced by high variations in salinity, as a result of fresh-water runoffs (Guille & Lasserre, 1979), and sudden temperature shifts including heat waves in the austral summer (Motreuil et al., 2018).

| Biodiversity data and data partitioning
A set of 26 presence-only records of A. cordatus sampled from 1898 to 2015 in the Baie du Morbihan was compiled by Guillaumot et al. (2016). Most of the data were collected after 1975 (Guillaumot et al., 2016; and the temporal heterogeneity of data sampling was proved to barely influence the results of SDM predictions . No absence data could have been gathered due to the difficulties in frequently accessing the area. Presence-only records from Guillaumot et al. (2016) were checked for georeferencing errors and complemented with data from Poulin and Féral (1995) (Figure 2a).
Data are homogeneously distributed in the area with a Moran's I score of −0.01 (p-value = .15). Consequently, background records were randomly sampled in the area without any targeted sampling approach as the effect of spatial autocorrelation was not significant Phillips et al., 2009). In order to sample environmental conditions prevailing in the study area as pre-

| Study area
The study area focusses on the Baie du Morbihan, a 700 km 2 silled basin 50 m deep on average, located in the east of the Kerguelen Islands (sub-Antarctic) (Appendix S1 and Figure 2a). Since the 1960s, the area F I G U R E 1 Flowchart of the methodological framework. The DEB model developed for A. cordatus is used both for spatially projecting species physiological performances on the Baie du Morbihan and for assessing the relationship between species growth performance and food availability. These products are then used to enrich the classic species distribution model (SDM) by (1) addition of the spatial DEB layers to the initial environmental layers (integrated SDM-DEB model) or (2) addition of the a priori physiological parameters into the SDM equation (integrated Bayesian procedure). has been recurrently studied by marine ecologists who conducted research programs in biological oceanography including studies of micro-and macrobenthic communities (Delille et al., 1996;Poulin, 1996).
Depth, sea surface temperature and primary production were used as environmental predictors of the distribution of A. cordatus.
'Seasonality' was assumed by focussing on environmental contrasts between the austral summer and the austral winter. Monthly values over the 2002-2021 period were studied (Appendix S2). Two dates were selected to represent seasonal conditions: 2017/02/09 for the summer period (warm temperatures and medium chlorophyll-a concentration) and 2017/08/20 for the winter period (colder temperatures and low chlorophyll-a concentration). These dates have also suitable satellite images that could be processed for the study and belong to the same year, which is appreciable to limit uncertainties associated with temporal heterogeneities in species distribution modelling .
By selecting these two dates, we aim at comparing the contrasts between model predictions as a mean of evaluation of the influence of environmental conditions on species likelihood of distribution. In any case, we could assume that generated model projections represent the true distribution of the species at the two considered dates, given the lack of precision in our datasets.

| Environmental data: Bathymetry, chlorophyll-a concentration and seawater temperature
The bathymetric chart was obtained from Beaman and O'Brien (2011), available at https://resea rchda ta.edu.au/kergu elen- Water is colder in August (temperatures range between 2.7 and 3.3°C) and food availability much lower than in February, with the richest environments located nearshore. 100 m). It was updated by Sexton (2005) using new single beam echosounder data from commercial fishing and research voyages, and some new multibeam swath bathymetry data. Satellite-derived datasets were used to provide island topography and to fill in no data areas (see Beaman & O'Brien, 2011).
As a deposit-feeder, A. cordatus feeds upon organic grain coatings and particles present in sediments (Pascal et al., 2021). Seawater chlorophyll-a concentration was used as a proxy of food availability because data on the exact organic content of sediments are not available at the scale of the entire bay (Arnould-Pétré et al., 2020). Values were retrieved using imagery from Operational Land Imager (OLI) and Thermal InfraRed Sensor (TIRS) of Landsat 8 obtained from USGS (United States Geological Survey, 2019, https://earth explo rer.usgs. gov/, accessed on May 2020). Chlorophyll-a concentration was derived from OLI data using the Case-2 Regional Coast Colour processor (C2RCC) (Brockmann et al., 2016) for the SentiNel Application Platform

| Dynamic energy budget (DEB) model
The DEB theory defines individuals as dynamic systems and provides a mathematical framework for the life cycle of an organism, from the start of the embryo development to the death. It describes the physiological processes with four primary state variables: reserve, structure, maturity and reproduction buffer (the latter for adults only), directly linked to mass and energy flows and influenced by two forcing environmental variables: temperature and food resources availability (Appendix S4, Kooijman, 2010). DEB theory relies on key concepts such as first laws of thermodynamics for conservation of mass, energy and time (Jusup et al., 2017) and assumes that the various energetic processes, such as assimilation and maintenance rates are dependent either on surface area or on body volume (van der Meer, 2006, more details given in Appendix S4).
The model was specifically built for A. cordatus using zerovariate (single data) and uni-variate (x ~ y relationship data) datasets According to DEB theory, the somatic maintenance ṗ M has priority over growth and reproduction to ensure survival. Maturity maintenance ṗ J has priority over reproduction (Kooijman, 2010). These conditions imply that if the energy available in the reserve compartment ṗ C is not sufficient to pay for the required maintenance costs (ṗ C < ṗ M + ṗ J), the organism cannot reproduce, and will progressively starve and die.

| Physiological submodel based on food availability
Using DEB equations and parameters (Equation 1), average growth rates were calculated for individuals measuring from 2.5 to 4.5 cm, according to food availability (for all values available in the projection area, Figure 2a-e) and a random selection of temperatures within the range of values of the considered date (Appendix S6). This constitutes the 'physiological submodel' that therefore takes into account both food availability and temperature. Twenty-five replicates of individual sizes and temperature selection were performed. The growth rate was calculated with the following DEB equation (Kooijman, 2010): with kap being the fraction of energy directed towards complex- M the somatic maintenance rate coefficient (time −1 ) and T C the temperature correction factor (−) (see Appendix S4 for more details).
Concretely, the physiological submodel was built by generating a Bayesian beta regression, with food availability as a predictor and growth performance probability as a response. A total of 4000 MCMC (Markov Chain Monte Carlo) samples were used for burn-in and the posterior distribution was estimated using 4000 additional samples. The physiological submodel coefficients were initiated with Gaussian priors, with the mean taken from the maximum likelihood estimation to improve convergence and a vague prior set on the variance (set at 1000).

| 'Classic' species distribution modelling (SDM)
A Generalized linear model (GLM) was used to relate species occurrences with the three environmental predictors previously described (depth, food availability, sea surface temperature and their square forms, Figures 1 and 2), in order to fit the equation with the equation y = b0 + b1*depth + b2*food + b3*temperature + b4*tempe rature 2 + b5*f 2 . In this approach, presence and background data are

| Integrated 'SDM-DEB' model
Integrating correlative and mechanistic models were first tested by using the spatial projection of the DEB model as an environmental predictor in the SDM (Buckley et al., 2011;Elith et al., 2010;Mathewson et al., 2016;Rodríguez et al., 2019). The procedure is similar to the 'classic' SDM model approach, except that the DEB layer (i.e. 'ṗ C > (ṗM + ṗ J)?') was added to the initial set of environmental predictors (depth, temperature, food availability). Similarly, the procedure was replicated for 50 samplings of background records, and the average relative likelihood of occurrence was predicted on a map.

| Integrated Bayesian model
The method developed by Talluto et al. (2016), and applied by Gamliel et al. (2020)  trarily fixed at 100, as we considered them as vague priors (Gamliel et al., 2020). The detail of prior values is given in Appendix S10.

| MODEL A SS E SS MENT
The DEB model was validated by estimating the goodness of fit using Model relative likelihood of occurrence for all approaches was evaluated by measuring the Area Under the Curve (AUC) (Allouche et al., 2006;Elith et al., 2006;Fielding & Bell, 1997) using the R package ROCR (Sing et al., 2005). In complement, the percentage of correctly classified presence data was measured by extracting likelihood values over the position of each presence data and compared with the MaxSSS threshold (Maximum Sensitivity plus Specificity threshold), highlighted to be the best threshold to characterize predicted suitable (>MaxSSS value) and unsuitable areas (<MaxSSS value) for presence-only models (Liu et al., 2013). Standard deviations of model replicates was used as uncertainty maps (Buisson et al., 2010;Swanson et al., 2013).

| MODEL PREDIC TIONS
In all cases, 50 model replicates were generated to represent model variability in link with background data resampling. Each group of 50 model replicates were averaged and plotted for comparison.
Partial dependence plots were used to represent the relationship between model predictions and environmental values and compared between models. They are built by plotting model likelihood values of each grid-cell pixel (y axis) against the value of the environment at the same pixel (x axis; each partial dependence plot is specific to a single environmental layer). Partial dependence curves were also used as a mean of model evaluation based on expert knowledge.

| Spatial projection of the DEB model
Spatial projections of DEB model outputs show important contrasts between the two dates (Figure 3a,b). In February, when temperatures are higher than 6°C and food availability homogeneously higher than 0.5 over the entire bay area (Figure 2), high species survival and reproduction are predicted almost everywhere (Figure 3a), except in some areas where food availability is very low (Figure 2).
Nearly four times more energy is predicted to be contained in the reserve compartment of A. cordatus in February compared with August (Appendix S7), energy available for individuals' maintenance and development.
In contrast, in August, the DEB model predicts maintenance costs of up to three times higher than in February while the energetic load available is lower (Appendix S7), leading to reduced reproduction and survival abilities in the majority of the study area.
Individual survival is modelled to be higher closer to the shoreline due to higher food availability (Figure 3b).

| Classic species distribution model (SDM)
The overall likelihood of occurrence predicted by 'classic SDMs' are low (<0.5, Figure 3c,d) for the entire area and both dates, and standard deviations are comparatively high (homogeneously close to 0.45 for February and more contrasted in space but coastal areas reaching 0.45 too for August, Appendix S8), stressing an important variability between model replicates.
Average relative likelihood scores are more contrasting in August than in February (Figure 3c,d). For August, the model predicts the highest relative likelihood of occurrence (around 0.5) near the shoreline, in shallow-water areas, and the lowest relative likelihood of occurrence (around 0.2) in the center of the bay and in Northwestern Fjord characterized by deep waters (Figure 3d). In February, the relative likelihood of occurrence is homogeneous in all the areas and close to 0.4 (Figure 3c).
Areas where model extrapolation occurs correspond from 36 (in February) to 37.8% (in August) of the total surface of the projection area and are mainly to be related to depth and to temperature in large patches for February (black and light grey patches, respectively, Appendix S8). Extrapolation in link to food availability (dark grey patches) is almost fully absent for February, whereas patchy and as frequent as extrapolation linked to temperature (light grey patches) for August (Appendix S8).

| 'Integrated SDM-DEB' model
Model relative likelihood scores are highly contrasting between

| 'Integrated Bayesian' model
'Integrated Bayesian' models were implemented using the following set of parameters as priors (Appendix S10). The coefficient values of f and f 2 are high compared to that of the other parameters (average and tau scores), increasing the influence of food availability in final model outputs (Appendix S10). In August, the coefficient value of the f parameter is eight times higher than in February (8.43 compared with −0.89), but f 2 is twice lower (11.38 compared with 27.78) (Appendix S10).
In 'integrated Bayesian' models, the relative likelihood of occur-

| Contribution of predictors and model performance
Model performance ( Partial dependence plots ( Figure 5) were generated to evaluate the influence of each environmental predictor (depth, food availability and temperature) on model predictions. Overall, a comparison between models shows that integrated modelling approaches ('integrated SDM-DEB and 'integrated Bayesian) provide more contrasting response curves for all three predictors compared with the 'classic SDM' approach, both for February and August ( Figure 5).
The 'integrated Bayesian' model results ( Figure 5) suggest a more substantial influence of environmental values on predicted relative likelihood of occurrence, with higher temperatures, higher food availability and lower depths associated with higher predicted habitat suitability. This sensitivity of environmental influence on model predictions is confirmed by the higher performance metrics observed for the 'integrated Bayesian' approach, noteworthy in February ( Table 1).

| Potential and main limitations to the different modelling approaches
Correlative approaches ('classic SDMs') are aimed at describing the correlation between species occurrence records and environmental conditions. SDM outputs can provide knowledge on the main environmental factors that drive species distribution (Elith et al., 2006;Peterson et al., 2011). Because presence records are used as input data, SDMs also indirectly integrate the influence of other factors such as the effect of biotic interactions (either competition, exclusion or facilitation between species) and the biogeographic context (barriers or dispersal vectors) on species distribution, thereby simply and explicitly assessing the species potential ecological niche (Soberón, 2010). However, the relevance of niche estimation often constitutes the main limitation to 'classic SDMs', because their predictive performance strongly relies on sampling completeness (Araújo et al., 2005;Broennimann et al., 2007;Holt, 2009;Loehle & Leblanc, 1996;Randin et al., 2006;Vaughan & Ormerod, 2003). The heterogeneity of presence sampling induces statistical artefacts that can bias model predictions (Bahn & McGill, 2007;Currie, 2007), a substantial limitation that has already been stressed in former works on the Southern Ocean .
Compared with SDMs, mechanistic models require more data (and require a good knowledge of species ecology or physiology) for parameter estimation and model implementation (Kearney & Porter, 2009). However, if the model can be built, the approach is powerful to evaluate the survival capacity of individuals in given environmental conditions (Arnould-Pétré et al., 2020;Fabri-Ruiz et al., 2021) and can estimate the species fundamental niche (Kearney & Porter, 2009).
Combining the merits of both correlative and mechanistic approaches to fine-tune the estimation of the species potential ecological niche can provide important benefits (Dormann et al., 2012), as prior information on the influence of the environment on species metabolism, given by physiological models, can be used to improve correlative models (Feng et al., 2020). This combined approach is also valuable to assess the effect of fast-changing environmental conditions (e.g. seasonality or future predictions), which generate non-equilibrium states (Kearney et al., 2008;Keith et al., 2008) that cannot be accurately modelled by static, correlative approaches (Fabri-Ruiz et al., 2021;Loehle & Leblanc, 1996;Schouten et al., 2020).
In the present study, the comparison of 'classic SDM', the most commonly used approach in ecological studies of Southern Ocean species, with 'integrated' approaches, was performed. All approaches have good performance statistics (Table 1), except for the 'spatial DEB' model. Spatial projections of the 'spatial DEB' approach are strongly F I G U R E 4 Spatial projections of the 'integrated SDM-DEB' models for February (a,c) and August (b,d), averaged 50 model replicates. Average distribution (a,b) and associated standard deviations (c,d). The available energy after paying off the somatic and maturity maintenances is integrated with the model as a predictor that assesses for each pixel the value of ṗ C − (ṗM − ṗ J), with ṗ C the amount of energy contained in the reserve compartment, ṗ M the amount of energy required for somatic maintenance and ṗ J the amount of energy required for maturity maintenance. Spatial projections of the 'integrated Bayesian' models for February (e,g) and August (f,h), averaged 50 model replicates. Average distributions (e,f) and associated standard deviations (g,h) TA B L E 1 Comparison of model performances (percentage of presence data correctly classified and area under the curve, AUC, metric) for the two seasons   (Table 1), as 'low food' areas are simply and systematically predicted as unsuitable to the species survival, with no consideration for the influence of the other environmental drivers. However, the model is interesting because it stresses the link between energetic costs and one major environmental driver (Appendix S7), a good complement to physiological submodels, and is interesting to assess the environmental conditions that drive species distribution.
'The classic SDM' is characterized by good validation scores (AUC > 0.71 and percentage of correctly classified presence data >77.8%) ( Table 1), but the relative likelihood of occurrence is contrasting for August compared with February (Figure 3c, d) when food concentration is high and evenly distributed in the all bay area ( Figure 2). As a consequence, the contribution of this variable to model predictions is low (Figure 5), an unrealistic prediction that contrasts with results obtained with the integrated approaches ('integrated SDM-DEB' and 'integrated Bayesian') ( Figure 5). Using a physiological submodel to inform an SDM has been applied in recent works by directly adding a physiological layer to the SDM (Buckley et al., 2011;Elith et al., 2010;Mathewson et al., 2016;Rodríguez et al., 2019) or by generating absence data from the modelled physiological information (Elith et al., 2010;Feng & Papeş, 2017). Model outputs are easy to interpret, but the approach requires the combination of several models, as in any hybrid approach, and implies a risk inherent in the addition of biased estimations of each individual model (Feng & Papeş, 2017). be used also influences model extrapolation (Appendix S9) (Rodríguez et al., 2019), which must be taken into consideration when interpreting model results (Buckley et al., 2011;Elith et al., 2010), and increases the complexity of the model calibration. Therefore, the real benefits of adding modelled physiological information to SDMs are case-dependent, and the improvement of modelling performances is not certain (Buckley et al., 2011;Rodríguez et al., 2019). However, the method can prove helpful for future predictions and analyses of non-equilibrium states, which constitutes the main limitation of the SDM approach (Buckley et al., 2011;Elith et al., 2010;Martínez et al., 2015;Mathewson et al., 2016). When there are few data available and the causal relationship between organism physiology and environment drivers difficult to model in a robust way, using the 'integrated SDM-DEB' approach can be problematic, and model outputs must be interpreted with caution.
Bayesian methods are increasingly used in marine sciences (Colloca et al., 2009;Gamliel et al., 2020;Muñoz et al., 2013;Pennino et al., 2014;Roos et al., 2015). They were proved to have several advantages compared with other methods, including (1) a more accurate and realistic estimation of uncertainty as observations and model parameters are both used as random variables in model predictions (Robert, 2007) and (2) the possibility to integrate information from different sources, scales or nature Hobbs & Ogle, 2011;Peters et al., 2004) with the inclusion of a priori knowledge to improve model goodness of fit and more accurate uncertainty estimates (Van Dongen, 2006).
In the present work, the highest AUC scores and correctly classified presence data were obtained with the 'integrated Bayesian' approach. Models performed well in representing uncertain areas, compared with other approaches (Figures 3 and 4), as the areas pre-

| Environmental changes (seasonality) and predicted distribution of A. cordatus
In all model predictions, the relative likelihood of occurrence is the highest in coastal areas, where populations of A. cordatus were known to be the most abundant (Poulin, 1996;Poulin & Féral, 1995).
Predicted suitable areas for the species perfectly match these conditions of high food availability and high temperature that prevail in coastal areas.
Important contrasts, however, were obtained in model predictions between February and August, suggesting that 'seasonal' variations significantly affect the metabolism of A. cordatus as organisms face different conditions in terms of food availability and temperature. By considering two dates with contrasting environmental conditions, we aimed at evaluating whether such differences between model outputs could be highlighted, as a mean of evaluation of the study methods.
According to the physiological model ('spatial DEB', Figure 3, Appendix S7), maintenance costs are higher in winter (August) than in summer (February) due to lower temperatures that increase the demand of energy to maintain the metabolism (Kooijman, 2010).
Besides, there is less energy available in the reserve compartment to compensate for the increased maintenance costs as food availability is low in winter too (Appendix S7).
These results are strongly dependent on the assumption that metabolism performance (and therefore requested energy) follows Arrhenius laws as determined with summer acclimated individuals (Motreuil et al., 2018). For some Antarctic sea urchins, such as Sterechinus neumayeri (Meissner, 1900) it was reported a sharp metabolic switch during winter conditions. During this hypothesized non-feeding period, metabolic rates are decreasing with lower recorded oxygen consumption and slow or absent somatic growth (Brockington et al., , 2007Brockington & Peck, 2001 February and August 2017, has long been suggested in SDMs Franklin, 2009), but it is seldom achieved due to limited data availability . Conversely, ignoring the effect of seasonality in ecological niche estimation has been recently shown to reduce prediction performance (Smeraldo et al., 2018). Seasonality is a fundamental feature of environmental systems. It is particularly critical to life in temperate and high latitudes, and one key phenomenon to consider for studying both species distribution (Morelle & Lejeune, 2015;Zuckerberg et al., 2016) and metabolism (Bahlburg et al., 2021).  (Jansen et al., 2018) to estimate the proportion of organic matter that reaches the seafloor based on the knowledge of water currents. It could be also interesting to have some information about benthic detritic organic matter that the sea urchins could consume (Pascal et al., 2021). These data were, however, not available for the study area, but such models offer promising perspectives.

| Study improvements
(2) Detailed information on the link between temperatures and physiological performances is still missing, as we only have and use here the results of a survival experiment performed at different temperatures in 2018 (Motreuil et al., 2018). DEB modelling has the potential to include five Arrhenius parameters to precisely characterize the link between temperature and metabolism (Kooijman, 2010;Thomas & Bacher, 2018), but available experimental data on A. cordatus do not permit measuring them with precision. More data are still needed for our case study to reach this precision and improve the performance of the DEB model. (3) Finally, there is a lack of presence data to correctly calibrate the model, to validate it and to accurately consider model predictions as accurate likelihood of species distribution. Generating ecological models with small datasets was indeed shown to reduce modelling performances (Liu et al., 2019;Stockwell & Peterson, 2002) as it truncates predicted distribution and niche definition (El-Gabbas & Dormann, 2018;Hortal et al., 2008), and may lead to a reduction in model accuracy because the presence and background datasets would not differ markedly (Luoto et al., 2005) and constrain the evaluation process (Pearson et al., 2007) (reviewed in Guillaumot et al., 2021). Therefore, common validation approaches such as the cross-validation method (that uses a part of the dataset to train the model and another part to test it independently, Hijmans, 2012;Guillaumot et al., 2019) could not have been used for our study, which limited the power of our model evaluation.

| CON CLUS IONS
Our results suggest good performances of 'integrated Bayesian' approaches to estimate species potential ecological niche, compared with single correlative approaches or 'integrated SDM-DEB' approaches that might be biased by the subjective choice of the DEB layer used as an input into the SDM. More data are still necessary to better evaluate the model, to more accurately establish the relationship between the environmental conditions and the species physiology and to better represent the whole environment. However, this study showed the possibility to apply the method for a data-poor case study, which opens perspectives for future applications to a broad panel of natural examples, noteworthy in the context of decision making and conservation management. Indeed, better simulating the ecological range of species, including environmental changes and species physiological tolerance might constitute interesting approaches for future conservation issues on ecosystems facing global change.

ACK N OWLED G EM ENTS
We acknowledge the use of imagery from the NASA Worldview application (https://world view.earth data.nasa.gov), part of the NASA Earth Observing System Data and Information System (EOSDIS).
We are thankful to Bas Kooijman for his helpful comments to adapt the DEB model. This work was supported by a 'Fonds pour la formation à la Recherche dans l'Industrie et l'Agriculture' (FRIA) and 'Bourse Fondation de la mer' grants to C. Guillaumot Commerson for their invaluable help and support in the field.

CO N FLI C T O F I NTE R E S T
All authors disclose any potential sources of conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
R codes developed for this study are available in a public repository quoted in the material and methods section of the manuscript.

PE E R R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/ddi.13617.