Using species distribution models only may underestimate climate change impacts on future marine biodiversity Ecological Modelling

In face of global changes, projecting and mapping biodiversity changes are of critical importance to support management and conservation measures of marine ecosystems. Despite the development of a wide variety of ecosystem models capable of integrating an increasing number of ecological processes, most projections of climate-induced changes in marine biodiversity are based on species distribution models (SDMs). These correlative models present a significant advantage when the lack of knowledge on the species physiology is counterbalanced by the availability of relevant environmental variables over the species geographical range. However, correlative SDMs neglect intra- and inter-specific interactions and thereby can lead to biased pro- jections of changes in biodiversity distribution. To evaluate the influence of trophic interactions on projections of species richness and assemblage composition under climate change scenarios, we compared biodiversity pro- jections derived from an ensemble of different SDMs to projections derived from a hybrid model coupling SDMs and a multispecies trophic model in the Mediterranean Sea. Our results show that accounting for trophic in- teractions modifies projections of future biodiversity in the Mediterranean Sea. Under the RCP8.5 scenario, SDMs tended to overestimate the gains and underestimate the losses of species richness by the end of the 21st century, with marked local differences in projections, both in terms of magnitude and trend, in some biodiversity hotspots. In both SDMs and hybrid approaches, nestedness with gains in species richness was the main pattern driving dissimilarity between present and future fish and macro-invertebrate species assemblages at the Mediterranean basin scale. However, at local scale, we highlighted some differences in the relative contribution of nestedness vs replacement in driving dissimilarity. Our results call for the development of integrated modelling tools that can mechanistically consider multiple biotic and abiotic drivers to improve projections of future marine biodiversity.


Introduction
Given the ongoing growth of human population, anthropogenic pressures on ecosystems such as exploitation, pollution, habitat fragmentation and loss, and biological invasions are expected to intensify in the future, especially in a business-as-usual scenario (Bellard et al., 2012;Leadley et al., 2010;Shin et al., 2019;Steffen et al., 2018). Associated with climate change, such pressures are expected to exacerbate the deleterious impacts already observed on world's ecosystems with potentially dramatic cascading consequences for human societies (Bindoff et al., 2019;Díaz et al., 2019;Doney et al., 2012;Lotze et al., 2019;Serpetti et al., 2017;Thiault et al., 2019;. Climate-induced changes on marine ecosystems have been accelerating at such a rapid pace that they may dominate impacts resulting from other drivers of change in the coming decades (Cheng et al., 2019;IPCC, 2019;Pinsky et al., 2019). In response to climate change, species adapt, shift their geographical range or face local extinction, this causes substantial changes in species composition, abundance and biotic interactions (e.g., trophic interactions) and leads to profound reshuffles of marine ecosystem structure and functioning (Blowes et al., 2019;Bryndum-Buchholz et al., 2019;Du Pontavice et al., 2019;Pereira et al., 2010;Román-Palacios and Wiens, 2020). Fish and other marine species have already altered their geographical distributions to stay within suitable environmental conditions (Lenoir et al., 2020), and these shifts will likely continue or accelerate in the future (Baudron et al., 2020;Cheung et al.,2009Cheung et al., , 2016García Molinos et al., 2016;Pecl et al., 2017;Pinsky et al., 2020;Poloczanska et al., 2013). It is, therefore, critical to develop our capacity to project and map biodiversity changes to support management and conservation measures of marine ecosystems (Bellard et al., 2012;IPBES, 2016;Parmesan et al., 2011), a major challenge for the mitigation of current biodiversity losses and the achievement of the United Nation's Sustainable Development Goal 14 (Pinsky et al., 2018;Singh et al., 2019).
Despite the development of increasingly sophisticated models that are able to integrate an increasing number of ecological processes, a majority -64% according to the recent global assessment of the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services (IPBES) -of global and regional scale scenarios of climateinduced changes in marine biodiversity rely solely on correlative approaches (Gavish et al., 2017;IPBES, 2019). Such approaches, based on the ecological niche conceptualised by Hutchinson (1957) and commonly named Species Distribution Models (SDMs) or Environmental niche models (ENMs), are empirical statistical models relating species distribution data (occurrences or abundances) to a set of environmental descriptors to predict where a species is likely to be present in unsampled locations or time periods (see Guisan and Thuiller 2005 or Elith and Leathwick 2009 for a review). SDMs have become very popular to project past, present and future diversity distributions under various climate change scenarios (Bellard et al., 2012;Gavish et al., 2017;Guisan and Thuiller, 2005). The development of ready-to-use and generic software has greatly contributed to the widespread use of SDMs in ecology (e.g., Biomod R package (Thuiller et al., 2009(Thuiller et al., , 2016, as well as the accessibility to global occurrence datasets for an increasing number of taxa. With the growing number of studies involving the use of SDMs, however, came a growing recognition of their limitations in projecting future distributions under new environmental conditions (Araújo and Luoto, 2007;Davis et al., 1998;Dormann, 2007;Elith and Leathwick, 2009;Heikkinen et al., 2006;Hof et al., 2012;Morin and Lechowicz, 2008). In particular, and even though observed presence-absence data result from complex processes, SDMs do not explicitly consider key ecological processes such as species interactions (e.g., predation, competition, mutualism or facilitation), adaptation (e. g., phenotypic plasticity), or population dynamics which greatly influence the distribution of species at local to global scales Araújo and Luoto, 2007;Elith and Leathwick, 2009;Giannini et al., 2013;Guisan and Thuiller, 2005;Lavergne et al., 2010;Palacio and Girini, 2018;Wisz et al., 2013). De facto, by ignoring biotic interactions and most often relying on presence-absence data rather than abundance data, SDMs may under-or overestimate the number of species that can coexist at a given location, time and environmental condition, and thus lead to distorted spatial projections of future biodiversity (Engelhardt et al., 2020;Gavish et al., 2017;Gilman et al., 2010;Urban et al., 2013). While the role of biotic interactions, population dynamics and life history processes in shaping distributions and species assemblages is well recognized, very few studies and essentially all in the terrestrial realm, have attempted to investigate the effects of considering biotic interactions on projections of climate change effects on species richness and community composition (Araújo and Luoto, 2007;Hof et al., 2012;Kissling et al., 2010). So far, the inclusion of biotic interactions (e.g., using patterns of species co-occurrence (Wisz et al., 2013)) has been shown to affect SDMs' performance but we still lack an understanding of how such interactions, when explicitly considered in the modelling of species distributions, will alter projections of biodiversity distribution under climate change. In addition, when biotic interactions were considered, they were, in most cases, assumed to be static in time and space, an assumption often unverified in Nature (Davis et al., 1998;Klanderud and Totland, 2005).
The Mediterranean Sea is rare in that it is both a biodiversity and global change hotspot (Coll et al., 2012(Coll et al., , 2010Cramer et al., 2018;Giorgi, 2006;Soto-Navarro et al., 2020). Several studies have projected that fish, invertebrates and zoo-and phyto-plankton species will shift their distribution under climate change leading to local changes in species richness and composition, alterations in food web structure and decreases in functional and phylogenetic diversity Ben Rais Lasram et al. 2010;(Albouy et al., , 2015aHattab et al., 2016;Benedetti et al., 2018Benedetti et al., , 2019. These studies, however, invariably relied on environmental niche modelling that either overlooked intra-and inter-specific interactions or considered trophic links and food web structure emerging from SDM projections of co-occurrence and not from climate-driven processes driving individual species distributions and dynamics (Albouy et al., 2014b;Hattab et al., 2016).
In this study, we evaluated and quantified the influence of considering trophic interactions on changes in species richness and assemblage composition under climate change by comparing biodiversity projections (α and ß diversities) derived from an ensemble of SDMs and those from a hybrid model coupling SDMs and a multispecies trophic model in the Mediterranean Sea, under the "no mitigation" climate change scenario (RCP8.5). We used a new hybrid model, OSMOSE-MED, which allows to represent the spatial dynamics of interacting marine fishes and macroinvertebrates in the Mediterranean Sea (Moullec et al., 2019a(Moullec et al., , 2019b, including key life history processes such as growth, reproduction and mortality.

Materials and methods
We used SDMs to estimate present-day (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013) and future (2071-2100) climatic suitability (here after called species distributions) of one hundred species currently inhabiting the Mediterranean Sea (Appendix A1 in Supplementary data) under the high emission, "no mitigation policy" IPCC Representative Concentration Pathway (RCP) 8.5 scenario. The distribution maps produced by the SDMs were then used to force the spatial distribution of species within the OSMOSE-MED hybrid model (Fig. 1). As the development and calibration of SDMs and OSMOSE-MED hybrid model are fully described in Moullec et al. (2019b) and Moullec et al. (2019a), only a brief presentation of the models is given here.

Species distribution modelling
The spatial occurrence data of each species were compiled and merged from multiple databases: the Global Biodiversity Information System (GBIF), the Ocean Biogeographic Information System (OBIS), the Food and Agriculture Organization's Geonetwork portal (www.fao.org /geonetwork) and the FishMed database (Albouy et al., 2015a). Only occurrences recorded since 1975 were kept. Decadal climatologies of temperature and salinity within 40 depth layers were used to generate predictor variables used to calibrate SDMs. These data were extracted from the global World Ocean Database 2013 Version 2 at a spatial resolution of 1/4 • and were then interpolated by bilinear methods on a 1/12th (ca. 7 km in the Mediterranean Sea) grid. In order to take into account environment experienced by differences in the vertical distribution of species in the water column, six environmental metrics were derived from temperature and salinity climatologies from 1975 to 2012: average sea surface temperature and salinity (0-50 m depth), average vertical temperature and salinity (0-200 m depth), and average sea bottom temperature and salinity (bottom 50 m depth). We used temperature and salinity to calibrate SDMs (i.e., environmental niche models) as these climatic variables are typically the major factors determining large-scale patterns of species distribution and are commonly used as explanatory variables in SDMs in the Mediterranean Sea and across the world's oceans Ben Rais Lasram et al., 2010;Benedetti et al., 2019;Cheung et al., 2009;Schickele et al., 2020).
To account for the uncertainty associated with the choice of a particular correlative statistical method for climatic niche modelling, we used an ensemble forecasting approach implementing a set of eight different methods (generalized linear models, generalized additive models, classification tree analysis, boosted regression trees, random forests, multivariate adaptive regression splines, artificial neural networks and flexible discriminant analysis) embedded in the BIOMOD2 multi-model platform (Thuiller et al., 2016). Models were calibrated using data with a global coverage to capture all climate conditions encountered by the species and avoid truncated climate niche (Thuiller et al., 2004). Pseudo-absences were randomly generated in order to better characterize the environmental conditions experienced by species within their current ranges (see (Hattab et al., 2013(Hattab et al., , 2014 for methodology and Moullec et al., 2019b for more details). Models were evaluated according to the True Skill Statistic criterion (TSS; Allouche et al. 2006) with a three-fold cross-validation. First, only models with a TSS higher than 0.6 were kept for each species. Second, the consensus distribution was obtained with an ensemble forecast approach, meaning that results were weighted according to the TSS. Finally, consensus distributions were transformed into presence/absence maps by using the probability threshold that maximizes the model's TSS, meaning that the occurrence probability for grid cells under this threshold was set to zero.
SDMs calibrated under present-day conditions were then used to make projections of the spatial distribution of species in 2071-2100 based on projected future environmental predictors (Moullec et al., 2019a). Sea temperature and salinity values at different depths were acquired for the historical period  and the end of the century (2071-2100, under RCP8.5 scenario) from the regional climate system model CNRM-RCSM4 (Sevault et al., 2014) that was driven by atmosphere and ocean boundary conditions extracted from the general circulation model CNRM-CM5 (Voldoire et al., 2013). CNRM-RCSM4 showed high skill in reproducing the interannual to decadal Conceptual representation of the modelling and analytical frameworks adopted in the study. Species distribution models (SDMs) were independently calibrated for one hundred marine species and derived species binary maps (presence/absence; p/a) were used as inputs to drive species spatial distributions in the food web model OSMOSE. Future environmental data (temperature and salinity climatologies under RCP8.5 scenario) used to project future species distributions were derived from the CNRM-RCSM4 model. CNRM-RCSM4 drove the Eco3M-S biogeochemical model, its outputs (planktonic biomass) serving as potential prey field for species in OSMOSE. Changes in species richness and composition were calculated between the present period (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013) and the end of the 21st century (2071-2100) using outputs from both modelling approaches. Differences between projections were then compared.
Mediterranean variability for sea surface temperature and salinity (Sevault et al., 2014). Anomalies between the historical simulated period  and the future period (2071-2100) were calculated and added to present climate and salinity climatologies to project future species geographical distribution (Moullec et al., 2019a).
To match the spatial resolution of the OSMOSE-MED hybrid model, all species distribution layers were resampled on a regular grid of 20 × 20 km using nearest neighbour interpolation.

The osmose-med hybrid model
OSMOSE-MED, considered as a "hybrid model" in this paper, is a consistent, end-to-end modelling chain, including a general circulation model (CNRM-CM5; Voldoire et al., 2013), a regional climate model (CNRM-RCSM4; Sevault et al., 2014), a regional biogeochemistry model (Eco3M-S; Auger et al., 2011) and a multispecies trophic model (OSMOSE; Moullec et al., 2019b). Eco3M-S simulates the lower trophic part of the food wed (i.e., phyto-and zoo-plankton groups) and is driven one-way by the ocean and atmosphere outputs of CNRM-RCSM4. The food web model OSMOSE, which simulates the higher trophic level part of the food web (i.e., fish, cephalopods, crustaceans), is in turn driven by the biogeochemistry outputs of Eco3M-S (phyto-and zoo-plankton biomass). The modelling chain has been implemented under the RCP8.5 climate change scenario, allowing to consider not only potential changes in species distribution (from SDMs outputs) but also changes in planktonic productivity by the end of the 21st century in the Mediterranean Sea (Moullec et al., 2019a).
OSMOSE is a multispecies and individual-based model, spatially explicit and representing the whole life cycle of several interacting marine fish and macro-invertebrates species from eggs to adult stages Cury, 2004, 2001). Major ecological processes of the life cycle, namely growth, predation, reproduction, various sources of mortality (including fishing mortality), are modelled step by step (15 days in our study) (Appendix A2 in Supplementary data and see osmose-model.org for a full description of the OSMOSE model). OSMOSE assumes opportunistic predation based on spatio-temporal co-occurrence and size adequacy between a predator and its prey, allowing the emergence of complex trophic interactions. The model is forced by species-specific spatial distribution maps (one unique map per species throughout its life cycle) built from the SDMs approach described above (see Section 2.1). Our application covered the whole Mediterranean Sea on a regular grid of 20 × 20 km (6229 cells) and modelled the same one hundred marine species (mainly fish (85 species) but also cephalopods (5 species) and crustaceans (10 species)) than in the SDMs approach. Overall, the species represent around 95% of total declared catches in the region over the 2006-2013 period. A full description of the parameterization and calibration of OSMOSE-MED can be found in Moullec et al. (2019b).
To our knowledge, despite its complexity and associated uncertainty, the modelling chain used in the present study is the first able to project climate change and fishing effects on Mediterranean marine biodiversity in an integrated way, i.e., considering explicit and consistent changes in regional climate, ocean dynamics, nutrient cycle, plankton production, shifts in species distributions, their life cycles and their trophodynamic interactions (Moullec et al., 2019a).

Comparing projections of species richness and assemblage composition between SDMs and the osmose-med hybrid model
Following Albouy et al., (2012); Benedetti et al., (2018) and Le Marchand et al. (2020), we used species richness and the ß ratio to characterize temporal community change related to α and ß diversities.
The two diversity indices obtained within each grid cell using the SDMs approach were then compared with those obtained using the OSMOSE-MED approach, i.e., additionally constraining the trophic interactions and the whole species life cycle. α diversity was calculated by tailing species numbers for each time period (present and future) and modelling approach (SDMs and OSMOSE-MED) by stacking all the presence/absence maps (0/1) and counting the number of present species within each grid cell. To derive presence/absence maps from OSMOSE-MED outputs, we transformed abundance distribution map per species in presence/absence map, a species with a null abundance in a cell being considered as absent from this cell. Differences in species richness were computed as the difference between the future (2071-2100, with climate change) and present species richness (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013).
In a second step, we computed the ß ratio index in order to quantify the relative contribution of species replacement (i.e., turnover) vs. nestedness (i.e., species loss or gain) in temporal changes in community composition (Baselga, 2012). Disentangling the contribution of these two components of ß-diversity patterns is indeed fundamental to understanding how communities react to spatial, environmental and temporal changes (Carvalho et al., 2012). To calculate the ß ratio, we computed the pairwise Jaccard's dissimilarity measure and calculated the ratio between its two additive components: the dissimilarity due to species replacement (or turnover (ßjtu)) and the dissimilarity due to nestedness (ßjne) (see Baselga, 2010Baselga, , 2012. A ß ratio value smaller than 0.5 indicates than dissimilarity in assemblage is mainly driven by species replacement and a ß ratio value higher than 0.5 indicates that dissimilarity is mainly driven by nestedness. Furthermore, as any ß ratio value can be associated with either gains or losses in species richness (Baselga, 2010), to simultaneously analyse changes in species richness and assemblage composition, we analysed and mapped the ß ratio values and related changes in species richness Benedetti et al., 2018).
In a final step, α and ß diversity projections were compared between the SDMs and OSMOSE-MED approaches, at different spatial scales: Geographical Sub-Areas (GSAs, the spatial scale used for assessment and management of fisheries resources in the Mediterranean Sea), marine ecoregions (as listed in the Marine Strategy Framework Directive; MSFD) and the whole Mediterranean basin (Appendix A2 in Supplementary data).
All analyses were performed with R (R version 3.5.1, R Core Team 2018), using the betapart package for ß diversity analyses (Baselga and Orme, 2012).

Projecting changes in species richness
Under the RCP8.5 climate change scenario, 16 and 15 species are projected to increase their geographic range by the end of the century (2071-2100), with average gains in range area of ca. 163% and ca. 131% according to SDMs and OSMOSE-MED, respectively. On the contrary, the range of 25 (SDMs) and 27 (OSMOSE-MED) species is expected to shrink, with average losses in range area of ca. 29% (SDMs) and ca. 37% (OSMOSE-MED). amongst the loser species, 3 species are projected to become extinct (no more presence) in OSMOSE-MED while SDMs projected only a reduction of their climate suitable areas (in average -45%).
At the scale of the whole Mediterranean basin, the average differences in species richness (SR) between the present and future periods (RCP8.5 scenario) were 0.37 species per grid cell when running the SDMs only (mean SR for the present period being 28.45 species per grid cell) and -0.37 species per grid cell when running OSMOSE-MED (mean SR for the present period being 28.1 species per grid cell). With both modelling approaches, a contrasted spatial pattern emerged from the simulations between the western Mediterranean Sea, characterized by an overall decrease in SR, and the central and eastern Mediterranean Sea, characterized by an overall increase in SR. The differences between SR projections from SDMs and OSMOSE-MED were predominantly localized on the continental shelf with some exceptions, e.g., in the Levantine Sea (Fig. 2). At both marine ecoregions and Geographical Sub-Areas scales, significative differences between SR projections from SDMs and OSMOSE-MED were apparent, however under different orders of magnitude (Appendix A3 in Supplementary data). Projected changes in SR diverged the most between SDMs and OSMOSE-MED in the Tunisian Gulfs of Hammamet and Gabes, the northern Adriatic Sea, the north Aegean Sea and along the Egyptian coasts. In the northern Adriatic Sea (GSA 17), the average loss in SR was -0.61 species per grid cell in SDMs while it was -3.4 species per grid cell according to OSMOSE-MED. In the Gulfs of Hammamet and Gabes (GSA 13 and 14), SDMs projected average gains in SR of 0.34 and 1.18 species per grid cell respectively and OSMOSE-MED, on the contrary projected corresponding losses in SR of -1.85 and -1.1 species per grid cell.
At the Mediterranean Sea scale, 44.15% and 48.93% of the cells (mainly in the central and eastern Mediterranean Sea) exhibited an increase in SR by the end of the century according to SDMs and OSMOSE-MED, respectively. On the contrary, cells showing a decrease in SR represented 39.75% of the Mediterranean Sea in SDMs and 47.09% in OSMOSE-MED and were mainly located on the continental shelf of the western Mediterranean Sea and in the Adriatic Sea. Net differences in SR result from gains and losses in SR that spatially differed between OSMOSE-MED and SDMs ( Fig. 2 and Appendix A4 in Supplementary data). Using SDMs, higher gains in species richness were projected by the end of the century in the Gulfs of Hammamet and Gabes (on average 2.7 and 4.1 times more, respectively, p<0.05), in the northern Adriatic Sea (on average 1.9 times more, p<0.05), in the Levantine Sea (on average 1.2 times more, p<0.05) and to a lesser extent in the Aegean Sea (on average 1.12 times more, p<0.05). On the other hand, lower losses in SR were projected by SDMs on the continental slope of the gulf of Gabes (on average 3.7 times less, p<0.05), in the northern Adriatic Sea (on average 2.6 times less, p<0.05), the Gulf of Lions (on average 1.4 times less, p<0.05), the continental shelf of the Catalan Sea (on average 1.3 times less, p<0.05) and the Alboran Sea (on average 1.24 times less, p<0.05). In the Adriatic Sea, for example, the number of cells affected by a decrease in SR was ca. 29% lower in SDMs compared to OSMOSE-MED. In the Gulf of Gabes, ca. 41% more cells were concerned by gains in SR in SDMs compared to OSMOSE-MED. In other words, under climate change scenario, SDMs overestimated the gains and underestimated the losses in species richness compared to OSMOSE-MED.

Projecting changes in species composition
Comparing projections of simultaneous changes in species richness (α diversity) and species composition (ß diversity) from SDMs and OSMOSE-MED revealed similar changes in species composition between SDMs and OSMOSE-MED at the Mediterranean scale (mean ß Jaccard = 0.02), albeit with more than 35.6% of the cells showing a non-zero dissimilarity. At a finer spatial scale (GSAs), dissimilarities in species composition by the end of the century between the two approaches were more apparent (Fig. 3 and Fig. 4 and Appendix A5 in Supplementary data), notably in the southern Alboran Sea, the Gulf of Gabes, the Adriatic Sea and the Levantine Sea. For both modelling approaches, species nestedness (i.e., ß ratio >0.5) accompanied by gains in species richness contributed more than other patterns in driving dissimilarity between present and future assemblages (nestedness with gains in SR represented ca. 37% and ca. 34% of the Mediterranean Sea surface according to SDMs and OSMOSE-MED, respectively) ( Fig. 3 and Appendix A5 in Supplementary data). In both approaches, nestedness represented ca. 64% of the whole Mediterranean Sea but with marked local differences in projections through a west-east gradient: the western Mediterranean Sea was dominated by nestedness and a clear decrease in species richness (mean ß ratio = 0.76 and 0.77 for SDMs and OSMOSE-MED, respectively, p>0.05), the central Mediterranean Sea (Gulfs of Gabes and Sidra) showed a combination of nestedness and replacement with gains and losses in SR (mean ß ratio = 0.24 and 0.23 for SDMs and OSMOSE-MED, respectively, p>0.05) and the eastern Mediterranean Sea was characterized by nestedness associated to a net increase in species richness (mean ß ratio = 0.78 for both SDMs and OSMOSE-MED, p>0.05. Whether for nestedness or replacement, SDMs projected greater gains and lower losses in species richness compared to OSMOSE-MED. When analysing beta-diversity patterns at a finer spatial scale (i.e., at the GSAs scale), the relative contribution of nestedness and replacement did not drastically differ between SDMs and OSMOSE-MED projections (Appendix A5 in Supplementary data). However, in the Gulf of Gabes, SDMs projected that species replacement with increases in species richness was the main pattern driving dissimilarity between time periods (but with many local patches on the continental shelf where nestedness associated to gains in SR was the main driver) while OSMOSE-MED projected a decrease in SR associated with species replacement (only 10% of the cells in SDMs and three-fold less (3%) in OSMOSE-MED were associated with nestedness). Projections also markedly differed in the southern Alboran Sea (28% of the cells in SDMs and 50% in OSMOSE-MED were associated with nestedness), the northern Adriatic Sea (7% of the cells in SDMs and almost double in OSMOSE-MED (13%) were associated with nestedness), the northernmost parts of the Aegean Sea and at the mouth of the Nile.

Discussion
Most studies projecting climate change impacts on species assemblages have used environmental niche modelling approaches that overlook key biological mechanisms such as demography, dispersal, eco-evolution, and species interactions shaping the distribution and productivity of species (Bellard et al., 2012;Urban et al., 2016). Our analysis suggested that SDMs represent a relevant and parsimonious statistical tool to project the environmental suitability of species at large spatial scales. At local scales, however, projections from SDMs alone and from a model using SDMs to drive biotic (trophic) interactions markedly differed. At the local scale (geographical sub-area in our study), projections of climate-driven changes in assemblage composition and richness were both qualitatively and quantitatively different when species interactions were considered in the modelling of species distributions. When applied alone, SDMs produced more optimistic projections of climate-induced changes in species richness in the Mediterranean Sea, overestimating gains and underestimating losses of species richness compared to the model including trophodynamics. In particular, important differences were found in some areas of high conservation interest such as the northern Adriatic Sea and the Gulf of Gabes, both recognized as biodiversity hotspots in the Mediterranean Sea Coll et al., 2010;Mouillot et al., 2011). Our results agree with other studies showing that species interactions may slow climate tracking and produce more extinctions than in environmental niche models (see Urban et al., 2013). In addition, changes in species composition and single-species extinctions may modify the web of interactions within a community and may lead to cascading effects and potential coextinctions, processes not reflected in SDMs but captured by OSMOSE-MED (Bellard et al., 2012;Brook et al., 2008). . Temporal change in species composition was quantified using the ßratio index. Grid cells with a ßratio <0.5 were dominated by the replacement process while grid cells with a ßratio >0.5 were dominated by the nestedness process. White colour on the map represents a Jaccard dissimilarity index equal to zero (meaning no change in species composition between time periods), in this case ßratio was not defined. The green to blue to purple colour gradient (matching the 0, 0.5 and 1 ßratio values) was used for grid cells showing a decrease in species richness while the yellow to red to brown colour gradient was used for grid cells showing an increase in species richness. Numbers represent the proportion of the Mediterranean Sea concerned by the corresponding change in species richness and composition (top right panel). Changes in species richness and composition between OSMOSE-MED and SDMs for the future period (2071-2100) (bottom left and bottom right panels). Dissimilarity between OSMOSE-MED and SDMs projections by the end of the century was quantified using the Jaccard's dissimilarity index (ßJaccard). Colors on the maps (left panels) correspond to colors on the scatterplots (right panels). The projected extent of expansion or contraction of geographical ranges of species by the end of the century differed depending on whether species interactions were considered, with lower gains and higher losses in our model including trophodynamic interactions. Various theoretical modelling studies have suggested that changes in intensity and nature of biotic interactions such as an increase in predation or a release of competition, modulate not only the direction but also the magnitude of species range shifts (Lenoir et al., 2010;Svenning et al., 2014;Urban et al., 2013). In the OSMOSE model, the presence of a species in a given grid cell (i.e., at least one individual) is based on the one hand, by specific distribution maps (from SDMs) and, on the other hand, by its ability to maintain itself in its suitable climatic space by feeding, competing, displaying somatic growth, reproducing and resisting natural and anthropogenic pressures (i.e., fishing mortality). These added processes can cause more rapid changes in the sizes of home ranges compared to spatial changes predicted merely by SDMs that consider only the suitability of abiotic factors. For example, sea warming could lead to the arrival of new species in an area which could directly compete with and displace native species (Alexander et al., 2018). However, there is still no consensus on whether considering trophic interactions produces more pessimistic climate-induced shifts in biodiversity distribution. For example, Fernandes et al., (2013) showed that the inclusion of trophic interactions in the species-based Dynamic Bioclimate Envelope Model (DBEM) led to lower projected poleward latitudinal shifts in the distribution of marine fish, especially for pelagic species, while Engelhardt et al., (2020) reported that ignoring trophic interactions tended to overestimate the effects of climate-induced changes in species distributions.
Climate change directly impacts species physiology and spatial distribution, but also engage species in novel interactions such as competition for food, with potential feedbacks on species distribution (Araújo and Luoto, 2007;Araújo and Rozenfeld, 2014;Gallardo et al., 2016;McKnight et al., 2017;Valiente-Banuet et al., 2015;Wiens, 2016). Even if it is still an immense challenge for ecologists and modellers, mechanistically considering these interactions is essential to explore the effects of future climate change on biodiversity distribution (Jackson et al., 2009;Pellissier et al., 2013;Van der Putten et al., 2010;Wisz et al., 2013). The OSMOSE-MED model explicitly considers the whole life cycle of species. Climate change can have different and potentially more severed effects on specific life stages. For example, the recruitment (year-class success) of marine fish is normally established by processes impacting early life stages such as eggs and larvae that are planktonic (free-floating) and more sensitive than juveniles and adults to changes in abiotic factors and prone to higher rates of mortality due to starvation and/or predation (Durant et al., 2007;Llopiz et al., 2014). Furthermore, predation is modelled as a size-based and opportunistic process in OSMOSE-MED (not predefined based on diet matrices of species), making it a powerful tool to predict potential interactions of species that do not currently co-occur or interact, but may do so in the future. With all these processes represented, our OSMOSE-MED hybrid model intends to better capture the future realized niche of the species for producing more robust biodiversity projections.
Our results pointed that in both SDMs and OSMOSE-MED approaches, at the Mediterranean Sea scale, dissimilarity between presentday and future assemblages was mainly driven by nestedness with increases in species richness. In contrast to what has been suggested for endemic species by Ben Rais Lasram et al. (2010), climate change might thus increase or decrease the richness of the initial assemblage without a total turnover of assemblages. The projected patterns were spatially contrasted, with for instance, in the Gulf of Gabes, three times as many grid cells concerned by nestedness in the SDMs compared to OSMOSE-MED. In this area, while SDMs projected changes in species composition almost entirely driven by species replacement with gains in richness, OSMOSE-MED projected that changes could be mainly caused by replacement with species losses.
Depending on the spatial scale considered, our results can highly differ from previous studies in the region that focused on projections based on environmental niche modelling approaches Ben Rais Lasram et al., 2010;Benedetti et al., 2019Benedetti et al., , 2018. For instance, under the SRES A2 scenario,  found that species replacement was the dominant pattern of change projected in Mediterranean coastal fish assemblages (accounting for 98% of the continental shelf). In contrast and in agreement with our estimates, Benedetti et al., (2018) found that nestedness was the main pattern driving dissimilarity between present and future copepod assemblages. Our results also differed from Albouy et al., (2013) who projected that ca. 70% of the continental shelf could experience a decrease and ca. 27% an increase in species richness by the end of the century (47% and 44%, respectively, at the whole Mediterranean Sea scale, according to OSMOSE-MED). Compared to  and Benedetti et al., (2019), the average rate of change in species richness that we found at the Mediterranean Sea scale was also much lower, with higher losses in the western regions and lower losses in the eastern ones. This pattern was partly related to the large increase in abundance of thermophilic and/or exotic species in the eastern Mediterranean and to the decrease in abundance of some species of high commercial interest in the western part (see Moullec et al., (2019a) for an extended discussion on climate change consequences on the structure and functioning of the Mediterranean Sea by the end of the 21st century). There are several explanations for the differences in projections. First, the environmental niche models used by  were not trained at the global scale but at the Mediterranean Sea scale which can lead to truncated environmental niches and overestimation of extinction rates (Thuiller et al., 2004). Second, not the same set of species (endemic species vs. one hundred miscellaneous species in our study) and the area (continental shelf vs. whole Mediterranean Sea in our study) were modelled which could also contribute to the discrepancies between the studies. Third, previous studies used a different IPCC climate change scenario, namely the Special Report on Emission Scenario (SRES) A2 (vs the RCP8.5 scenario), and an older version of the regional circulation model CNRM-RCSM4 (IPCC, 2014; Rogelj et al., 2012). Fourth, in contrast to previous niche modelling approaches, the end-to-end model OSMOSE-MED that we developed is fully integrated and considers, for the first time in the region, major biotic and abiotic processes that drive changes in population dynamics and community structure and composition. Despite attempts to make projections more realistic and credible, several caveats persist in our approach. For example, we did not consider other drivers of fish distribution (e.g., seafloor, oxygen concentration) that may prevail at local scales to prevent life-cycle closure and population persistence (Petitgas et al., 2013) and which will most likely be impacted by future anthropogenic pressures. We also did not consider the potential arrivals of thermophilic non-native species from the Atlantic or the Red Sea which will likely influence the future biodiversity in the Mediterranean Sea (Ben Rais Lasram and Parravicini et al., 2015) and could dramatically affect future assemblages (Edelist et al., 2013). Another limit in our projections was that phenotypic plasticity and genetic evolution were not included into OSMOSE-MED to mechanistically depict the potential eco-evolutionary dynamics shaping future communities and food webs (Peck et al., 2018). Yet, these adaptive processes might allow species to evolve in a new ecological niche, thus dampening climate change impacts on marine organisms and ecosystems (Beaugrand and Kirby, 2018;Cotto et al., 2017). Although our end-to-end approach considers two major components of global change (climate change and fishing pressure), we neglected other anthropogenic pressures such as pollution or habitat destruction which can act in synergy and exacerbate the vulnerability of species, ultimately accelerating changes in species composition and richness and making the estimates from our projections conservative in terms of biodiversity loss. Finally, for practical and pragmatic reasons (Moullec et al., 2019a(Moullec et al., , 2019b and, similar to most studies conducted in the region, we did not thoroughly investigate the uncertainty associated with our projections. Although some aspects of structural uncertainty were addressed by using an ensemble estimate of SDMs, other sources of uncertainty stemming from differences amongst scenarios, global climate models, and regional physical-biogeochemical models were not included (Peck et al., 2020). This is clearly one of the future challenges that need to be addressed when tackling biodiversity projections, especially with increasingly complex models (IPBES 2016;Tang et al., 2018;Thuiller et al., 2019).

Author contributions
FM and YJS conceived and designed the research. FM developed the models and performed the analyses. SD, FG, TH, MP and YJS helped in data analyses and/or interpretation. SD and FG helped in the development of species distribution models. NB helped with programming the code of OSMOSE and the use of the HPC cluster DATARMOR. All authors contributed to the writing of the article and approved the submitted version.

Data availability statement
The raw data supporting the conclusions of this manuscript will be made available by the authors, without undue reservation, to any qualified researcher.

Declaration of Competing Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.