Mario M Tagliari
email: mario.tagliari@posgrad.ufsc.br|ghislain.vieilledent@cirad.fr

Pascal Danthu

Jean-Michel Leong

Pock Tsy

Cyrille Cornu

Jonathan Roger Michel

Henri Lenoir

Vítor Carvalho-Rocha

Ghislain Vieilledent

Mario Muniz Tagliari

Jean-Michel Leong

Jonathan Lenoir

M Muniz

J.-M Leong

Fofifa

Not all species will migrate poleward as the climate warms: The case of the seven baobab species in Madagascar

Keywords: Baobabs, climatic anomaly, climate change, elevation, extinction risk, latitude, Madagascar, temperature seasonality, species distribution models, species range shift

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Climate change has already modified the spatial distribution of tropical biodiversity [START_REF] Chen | Elevation increases in moth assemblages over 42 years on a tropical mountain[END_REF][START_REF] Fadrique | Widespread but heterogeneous responses of Andean forests to climate change[END_REF][START_REF] Feeley | Where are the tropical plants? A call for better inclusion of tropical plants in studies investigating and predicting the effects of climate change[END_REF]. Increasing temperatures, anomalous precipitation regimes [START_REF] Anderson-Teixeira | Altered dynamics of forest recovery under a changing climate[END_REF] and more frequent and severe extreme events (e.g. heatwaves, droughts and wildfires; [START_REF] Garcia | Multiple dimensions of climate change and their implications for biodiversity[END_REF] all pose significant challenges to biodiversity by pushing species towards the limits of their climatic tolerances [START_REF] Rodríguez-Castañeda | The world and its shades of green: a meta-analysis on trophic cascades across temperature and precipitation gradients[END_REF]. Temperature has been the main variable considered when studying biotic responses to climate change for several reasons. First, it is an easy to measure variable. The first reliable thermometers have been used to measure air temperature since the beginning of the 18th century [START_REF] Fahrenheit | VIII. Experimenta & observationes de congelatione aquae in vacuo factae a D. G. Fahrenheit[END_REF]. Second, temperature generally decreases with elevation and latitude and has been historically used to define habitat types on Earth (FAO, 2010;[START_REF] Von Humboldt | Des lignes isothermes et la distribution de la chaleur sur le globe[END_REF][START_REF] Holdridge | Determination of world plant formations from simple climatic data[END_REF]. Third, temperature is a known determinant of species biology and distribution [START_REF] Sentinella | Tropical plants do not have narrower temperature tolerances, but are more at risk from warming because they are close to their upper thermal limits[END_REF][START_REF] Tewksbury | Ecology: Putting the heat on tropical animals[END_REF]. Fourth, temperature is strongly correlated to CO2 concentration in the atmosphere and is expected to significantly increase at the global scale in the future, in association with increasing CO2 emissions due to human activities [START_REF] Ipcc | Climate change 2014: Synthesis Report[END_REF]. As a consequence, it is commonly accepted that species will move towards higher elevations (i.e. upslope) and latitudes (i.e. poleward) to track shifting isotherms as the climate warms [START_REF] Colwell | Global warming, elevational range shifts, and lowland biotic attrition in the wet tropics[END_REF][START_REF] Lenoir | A significant upward shift in plant species optimum elevation during the 20th century[END_REF][START_REF] Lenoir | Climate-related range shifts -a global multidimensional synthesis and new research directions[END_REF].

However, it has also been observed that species may go "against the flow" under the effect of climate change: towards lower elevations (i.e. downslope) and lower latitudes (i.e. equatorward) to find suitable climate conditions [START_REF] Lenoir | Going against the flow: Potential mechanisms for unexpected downslope range shifts in a warming climate[END_REF]. These unexpected directional range shifts may involve several potential determinants, e.g. indirect biotic response due to the combined effect of both climate warming and land-use change; changes in interspecific interactions such as competition release; sensitivity to other environmental gradients not conforming with upslope and poleward range shifts; physiological or evolutionary adaptations; and random shifts due to stochastic ecological processes [START_REF] Crimmins | Changes in climatic water balance drive downhill shifts in plant species' optimum elevations[END_REF][START_REF] Lenoir | Going against the flow: Potential mechanisms for unexpected downslope range shifts in a warming climate[END_REF][START_REF] Pinsky | Marine Taxa Track Local Climate Velocities[END_REF]). Yet, studies reporting species range shifts in response to anthropogenic climate change usually focus on two geographical dimensions solely -latitude and elevation [START_REF] Lenoir | Climate-related range shifts -a global multidimensional synthesis and new research directions[END_REF] -and one single climatic dimension, namely mean annual temperature. Hence, these studies disregard other relevant climatic predictors such as changes in precipitation regime, water balance, or temperature seasonality, which may force species to shift downward in elevation [START_REF] Crimmins | Changes in climatic water balance drive downhill shifts in plant species' optimum elevations[END_REF][START_REF] Lenoir | Going against the flow: Potential mechanisms for unexpected downslope range shifts in a warming climate[END_REF]. Given that, it is extremely important to account for additional climatic variables different from the mean annual temperature and consider other potential ecological processes that could explain species range shifts in multiple directions.

The most commonly employed tools to predict current and future distribution of species under climate change from a set of observations and climatic predictors are correlative species distribution models (SDMs) [START_REF] Elith | Do they? How do they? Why do they differ? on finding reasons for differing performances of species distribution models[END_REF][START_REF] Foden | Climate change vulnerability assessment of species[END_REF][START_REF] Guisan | Predicting species distributions for conservation decisions[END_REF][START_REF] Porfirio | Improving the Use of Species Distribution Models in Conservation Planning and Management under Climate Change[END_REF]. The main outputs of these SDMs are maps of species potential distributions in the present and future. Nowadays, many easy-to-use softwares (the JAVA Maxent; Phillips, Anderson, & Schapire, 2006) or libraries (e.g. "sdm" and "biomod2"; [START_REF] Naimi | sdm: a reproducible and extensible R platform for species distribution modelling[END_REF][START_REF] Thuiller | BIOMOD -A platform for ensemble forecasting of species distributions[END_REF] have been made available to easily derive such redistribution maps. Comparing present with future species distribution maps, one can assess species vulnerability to climate change by looking at species range shift, contraction, or expansion [START_REF] Vieilledent | Vulnerability of baobab species to climate change and effectiveness of the protected area network in Madagascar: Towards new conservation priorities[END_REF].

Surprisingly, conservation studies that employ correlative SDMs to assess species vulnerability to climate change mostly do not disentangle the respective effect of predictor variables in explaining species range shift [START_REF] Fourcade | Paintings predict the distribution of species, or the challenge of selecting environmental predictors and evaluation statistics[END_REF]. Thus, they fail at explaining the possible underlying mechanisms behind such changes. For example, in an article studying the vulnerability of three baobab species to climate change in Madagascar, [START_REF] Vieilledent | Vulnerability of baobab species to climate change and effectiveness of the protected area network in Madagascar: Towards new conservation priorities[END_REF] showed that Adansonia suarezensis H. Perrier, and A. perrieri Capuron, will likely become extinct by 2085. However, the study does not analyze the respective role of each variable in explaining species range shift and does not suggest potential mechanisms that could explain the species extinction. In addition, a recent study assessing the vulnerability of Madagascar endemic baobabs to future climate change also failed to explore the underlying mechanisms behind expected species range shifts [START_REF] Wan | Modeling impacts of climate change on the potential distribution of six endemic baobab species in Madagascar[END_REF].

To help fill this gap, we investigate here the specific role of a set of climatic variables in explaining shifts in species distribution associated with climate change. We chose the seven emblematic baobab species that can be found in Madagascar (Adansonia L. genus; Malvaceae family) for this purpose. Six of the seven species are endemic to Madagascar (Adansonia grandidieri Baill., A. madagascariensis Baill., A. perrieri, A. rubrostipa Jum. and Perr., A. suarezensis, and A. za Baill), while the remaining species also occurs in the African continent: A. digitata L. [START_REF] Wickens | The Baobabs: Pachycauls of Africa, Madagascar and Australia[END_REF]. Each of the seven baobab species are located in different regions of Madagascar, being adapted to different climates and could potentially have different responses to climate change. We gathered an extensive and unprecedented occurrence dataset obtained from more than 15 years of field prospection and photo-interpretation of highresolution satellite images. We used an ensemble modelling approach to model the climatic niche and to predict the distribution of these seven species. Using SDMs and climatic projections, we assessed whether each species could experience range shift, contraction or expansion, and in which direction. Looking at the relative climatic variable importance and future climatic anomalies, we disentangled the role of each climatic variable in explaining species range shift. Based on species' extinction risk, we made suggestions to update the current baobab species conservation status. Finally, we attempted to generalize our results to other species in the tropics that should experience similar climatic anomalies in the future.

Material and methods

Presence and pseudo-absence data

We used photo-interpretation of very high-resolution QuickBird (61 cm resolution: most of the cases, especially for baobab identification from the crown or by their projected shadow) and Spot5 (2.5 m resolution: only in few specific cases, such as dense and homogeneous forests) satellite images available on Google Earth (http://www.google.com/earth/index.html; see [START_REF] Yu | Google Earth as a virtual globe tool for Earth science applications at the global scale: progress and perspectives[END_REF] to locate A. grandidieri and A. suarezensis individual trees in Madagascar. To validate occurrence data from photo-interpretation, ground-truth verifications were conducted identifying baobabs trees by the basis of their crown size, shape and color during flowering period (see [START_REF] Vieilledent | Vulnerability of baobab species to climate change and effectiveness of the protected area network in Madagascar: Towards new conservation priorities[END_REF]. Groundtruth verification was conducted during the flowering period to facilitate species identification and validation of species occurrence data (see [START_REF] Vieilledent | Vulnerability of baobab species to climate change and effectiveness of the protected area network in Madagascar: Towards new conservation priorities[END_REF] for further details). For the other five Malagasy baobab species (A. digitata, A. madagascariensis, A. perrieri, A. rubrostipa, and A. za) we used an extensive presence only data-set available thanks to prospective fieldwork (2000 to 2015) from the Cirad Madagascar team. During fieldwork, baobab trees were identified at the species level and georeferenced with a GPS to generate a unique occurrence data-set for all Malagasy baobab species.

Our raw data-set contains 137,285 occurrence records encompassing all seven Malagasy baobab species. First, we removed all points with coordinates outside Madagascar (only for A. digitata because occurrence records were also collected in Comoro islands). Then, for each of the seven species separately, we created a grid with 1-km² cell resolution covering the Madagascar territory and identified all cells that had at least one occurrence record for the focal baobab species. Finally, we removed all cells and respective presence observation data with incomplete bioclimatic information. For instance, the initial set of 1,686 occurrence records available for A. suarezensis was reduced to a total of 170-pixel units of 1-km 2 each (Table S1 for all baobabs species). Our observation sample size was sufficient to perform SDMs because the recommended minimum sample size (see [START_REF] Van Proosdij | Minimum required number of specimen records to develop accurate species distribution models[END_REF] for narrow-ranged species (as for A. perrieri -21 1-km 2 grid cells) is as low as 3 while it is recommended to have at least 13 occurrence points for widespread species (as for A. grandidieri or A. za -3,772 and 460 1-km 2 grid cells, respectively). We randomly sampled 10,000 pseudo-absences (i.e. virtual absence data which are drawn to be representative of the environmental variability in the study-area; [START_REF] Barbet-Massin | Selecting pseudo-absences for species distribution models: how, where and how many? Methods in Ecology and Evolution[END_REF]) across all Madagascar for each species to constitute a presence/pseudo-absence data-set. By using pseudo-absences we used both presence and pseudo-absence information to predict species' habitat suitability and distribution, optimizing spatial and environmental discrimination [START_REF] Senay | Novel Three-Step Pseudo-Absence Selection Technique for Improved Species Distribution Modelling[END_REF]. Consequently, we aimed to have a good representativity of the climate variability in Madagascar and to be able to compute a relative probability of presence across the country.

Bioclimatic data

We used current (~1950-2000) and future (2055 and 2085) climatic data at 30 arc-seconds resolution (about 1 km at the equator) over the entire spatial extent of Madagascar. This data is freely available on MadaClim (https://madaclim.cirad.fr/). The MadaClim website provides climatic data for Madagascar obtained from the WorldClim (http://worldclim.org/bioclim/) and CGIAR-CCAFS climate data portal (http://www.ccafs-climate.org/). We selected four bioclimatic variables [START_REF] Hijmans | Very high resolution interpolated climate surfaces for global land areas[END_REF] to model species distribution which were weakly correlated among each other and easy to interpret with regard to baobab species distribution. Three of them were previously selected via a principal component analysis among all the 19 WorldClim bioclimatic variables (following [START_REF] Vieilledent | Vulnerability of baobab species to climate change and effectiveness of the protected area network in Madagascar: Towards new conservation priorities[END_REF]: (1) mean annual precipitation -prec (mm.y -1 ); (2) mean annual temperature -tmean (°C); (3) temperature seasonality -tseas (sd x 100 ⁰C). In addition, we included a synthetic variable reflecting (4) climatic water deficit -cwd (mm). The cwd variable was computed from monthly precipitation (mcprec) and potential evapotranspiration (pet) using the following formula (Equation 1):

𝑐𝑤𝑑 = -∑ 𝑖 𝑚𝑖𝑛 (𝑚𝑐𝑝𝑟𝑒𝑐 𝑖 -𝑝𝑒𝑡 𝑖 , 0) (1)
Potential evapotranspiration is defined as the evaporation amount that would occur if a sufficient water source was available. We used the Thornthwaite formula [START_REF] Thornthwaite | An approach toward a rational classification of climate[END_REF] to compute the monthly potential evapotranspiration. The four selected bioclimatic variables are widely used (i) to define biomes globally, known as Holdridge Life Zones System [START_REF] Holdridge | Determination of world plant formations from simple climatic data[END_REF] and (ii) as proxies for other bioclimatic variables. For instance, the mean annual temperature (tmean) is a proxy for solar radiation and temperature stress [START_REF] Haigh | The Sun and the Earth's Climate[END_REF]. Additionally, the mean annual temperature may indicate potential losses of plant productivity [START_REF] Hatfield | Temperature extremes: Effect on plant growth and development[END_REF]. The temperature seasonality (tseas) can be interpreted as a proxy for the growing season [START_REF] Hatfield | Temperature extremes: Effect on plant growth and development[END_REF] while the annual precipitation (prec) is a proxy for potential water availability [START_REF] Amissah | Rainfall seasonality and drought performance shape the distribution of tropical tree species in Ghana[END_REF]. Finally, the climatic water deficit (cwd) can be indicative of water stress and drought periods [START_REF] Fayolle | Patterns of tree species composition across tropical African forests[END_REF][START_REF] Stephenson | Climatic control of vegetation distribution -the role of the water balance[END_REF].

For future climate data (2055 and 2085) we selected three different global circulation models (GCMs) from the World Climate Research Programme (CMIP5) (i.e., NorESM1-M, GISS-E2-R, and HadGEM2-ES) under two representative concentration pathways (RCPs: carbon dioxide emission scenarios) (i.e., RCP 4.5 and 8.5). The RCP 8.5 scenario is characterized by high concentration and increasing CO2 gas levels emissions [START_REF] Riahi | Scenarios of long-term socio-economic and environmental development under climate stabilization[END_REF][START_REF] Van Vuuren | The representative concentration pathways: An overview[END_REF] and can be considered the most likely emission scenario in the absence of effective mitigation policies regarding CO2 emissions, whereas RCP 4.5 is known as the "mitigation scenario" because of projected reduction of CO2 gas levels emissions [START_REF] Van Vuuren | The representative concentration pathways: An overview[END_REF]. A recent discussion has been brought in the literature affirming that RCP 8.5 is a problematic scenario for near-term (2030-2050) emissions and indicate that RCP 4.5 is more likely than RCP 8.5 [START_REF] Hausfather | RCP 8.5 is a problematic scenario for near-term emissions[END_REF]. Despite this recent discussion, we projected our main results under RCP 8.5 because: (i) we projected for long-term climate change (i.e. 2085) where projections presented by RCP 8.5 in 2100 are more probable than RCP 4.5 (Schwalm, Glendo & Duffy, 2020a); (ii) historical cumulative CO2 emissions from 2005 to 2020 are more in accordance with RCP 8.5 than RCP 4.5 (Schwalm, Glendo & Duffy, 2020)b; (iii) RCP 4.5 underestimate biotic feedbacks (e.g. changes in soil dynamics, forest fires frequency and severity, permafrost thaw) which accelerates warming, further supporting RCP 8.5 (Schwalm, Glendo & Duffy, 2020b); (iv) in our study we used RCP 8.5 for the sake of risk assessment and not to compare RCPs effectiveness, despite RCP 4.5 projected temperature by 2100 is 1.7-3.2 °C, whereas for RCP 8.5 is 3.2-5.4 °C [START_REF] Fuss | Betting on negative emissions[END_REF]). As a consequence, the outputs of RCP 4.5 in our study are presented in the supplementary material.

Species distribution modeling: statistical algorithms, model performance and importance of bioclimatic variables

We selected four statistical algorithms to model the bioclimatic niche and distribution of the seven studied baobab species: generalized linear models (GLMs); generalized additive models (GAMs); random forests (RF); and Maxent. Algorithms selection included standard regression models, such as the parametric GLM and the non-parametric GAM, classification tree (RF), and maximum entropy approach (Maxent). We aimed to quantify output uncertainty and generate a gradient from robustness (GLM and GAM) to complex algorithms, i.e. RF and Maxent [START_REF] Elith | Do they? How do they? Why do they differ? on finding reasons for differing performances of species distribution models[END_REF]. The uncertainty quantification of predictive modeling follows the premise of the ensemble modelling approach [START_REF] Araújo | Ensemble forecasting of species distributions[END_REF], which enables a consensus identification among all forecasts and the exploration of the full breadth of intermodal variability [START_REF] Kujala | Conservation Planning with Uncertain Climate Change Projections[END_REF].

As we used two regression models (GLM and GAM) and two machine learning approaches (RF and Maxent) to fit SDMs, the inclusion of 10,000 pseudo-absence points (background points for Maxent algorithm) is advised for better SDM outputs to obtain more accurate results [START_REF] Barbet-Massin | Selecting pseudo-absences for species distribution models: how, where and how many? Methods in Ecology and Evolution[END_REF]. We randomly split our presence/pseudo-absence data-set using 70% for model calibration (training data-subset) and 30% for model validation (testing data-subset) to evaluate the predictive performances of our SDMs [START_REF] Hijmans | Cross-validation of species distribution models: removing spatial sorting bias and calibration with a null model[END_REF]. We repeated the cross-validation procedure five times.

Model performance in predicting species presence-absence was estimated using four different and complementary metrics: Area Under the Receiving Operating Characteristics Curve (AUC); True Skills Statistics (TSS); Sensitivity (Sen); and Specificity (Spe; [START_REF] Liu | Measuring and comparing the accuracy of species distribution models with presence-absence data[END_REF]. We thus calculated the mean value of AUC and TSS metrics across the five testing data-subsets obtained from the crossvalidation procedure for each selected algorithm. By doing this we were able to describe the modelling performance in predicting species presence-absence. We also computed AUC and TSS metrics across the full data-set.

The AUC computes the model probability to rank a randomly chosen presence site instead of a randomly absent site [START_REF] Liu | Measuring and comparing the accuracy of species distribution models with presence-absence data[END_REF][START_REF] Pearce | Evaluating the predictive performance of habitat models developed using logistic regression[END_REF] and is commonly used as an accuracy index for SDMs using ensemble modelling approaches [START_REF] Hao | A review of evidence about use and performance of species distribution modelling ensembles like BIOMOD[END_REF]. It is a threshold-independent index, and it is also independent to prevalence [START_REF] Allouche | Assessing the accuracy of species distribution models: Prevalence, kappa and the true skill statistic (TSS)[END_REF], which is the proportion of samples representing species presence [START_REF] Mcpherson | The effects of species' range sizes on the accuracy of distribution models: Ecological phenomenon or statistical artefact[END_REF]. If AUC values are >= 0.9, the model is commonly considered as highly accurate [START_REF] Thuiller | BIOMOD -A platform for ensemble forecasting of species distributions[END_REF]. The TSS metric is a threshold-dependent index [START_REF] Liu | Measuring and comparing the accuracy of species distribution models with presence-absence data[END_REF] and is computed with a probability threshold maximizing its values. TSS values range from -1 to 1, and accurate models (correctly predicting both presences and absences) lead to values close to one [START_REF] Thuiller | BIOMOD -A platform for ensemble forecasting of species distributions[END_REF].

The TSS index is equal to Sensitivity + Specificity -1. Sensitivity is the probability of correctly predicting a presence while specificity is the probability of correctly predicting an absence [START_REF] Liu | Measuring and comparing the accuracy of species distribution models with presence-absence data[END_REF]. As well as the AUC index, TSS is not sensitive to prevalence [START_REF] Allouche | Assessing the accuracy of species distribution models: Prevalence, kappa and the true skill statistic (TSS)[END_REF][START_REF] Lawson | Prevalence, thresholds and the performance of presence-absence models[END_REF], so we used both accuracy indexes to evaluate SDM outputs for rare (such as A. perrieri) or abundant (such as A. grandidieri) baobab species.

To evaluate the performance of the ensemble model based on committee averaging we used three previously mentioned metrics: TSS, Sen, and Spe, following [START_REF] Araújo | Ensemble forecasting of species distributions[END_REF]. We previously defined an evaluation threshold using the accuracy index TSS (i.e. minimum score of 0.6 or 60%) in order to: (i) remove "bad algorithms/models"; (ii) build our ensemble model; (iii) test and evaluate the ensemble model forecasting capability (i.e. predicting species presence-absence); and (iv) make the binary transformation for the committee averaging computation [START_REF] Thuiller | BIOMOD -A platform for ensemble forecasting of species distributions[END_REF].

For each statistical algorithm, we calculated the relative variable importance among the four studied bioclimatic variables selected for the SDMs. The computation principle follows the one used for RFs, where one bioclimatic variable is shuffled over the full data-set or the testing data-subset. To compute variable importance (I), the model prediction is calculated in this shuffled data-set, and a correlation (Pearson's correlation) is computed between baseline predictions (pred_ref) and the shuffled predictions (pred_shuffled; see Equation 2).

𝐼 = 1 -𝑐𝑜𝑟 (𝑝𝑟𝑒 𝑟𝑒𝑓 , 𝑝𝑟𝑒𝑑 𝑠ℎ𝑢𝑓𝑓𝑙𝑒𝑑 ) (2)
We thus generated a rank according to the variable importance over the four statistical algorithms for each species. The rank was defined by calculating the mean obtained from 6 model runs (5 runs from the testing data-subset and 1 run from the full data-set). The rank with higher mean values indicates which bioclimatic variable is more important to explain the species distribution. We used the Biomod2 R package [START_REF] Thuiller | BIOMOD -A platform for ensemble forecasting of species distributions[END_REF]) to generate the SDMs.

Current species distribution, climatic niche and elevational range

For each species and each modelling algorithm, predicted probabilities of occurrence during the current period were binary transformed (0 for species absence and 1 for species presence) using the probability threshold maximizing TSS. Then, the current species distribution area (SDAp in km²) was defined as the set of 1-km 2 pixels where two out of the four modelling algorithms predicted the presence of the focal species. When only one algorithm out of the four predicted a presence for a given species, it was considered as uncertain. The species was considered absent when none of the four algorithms predicted a presence.

To characterize each species bioclimatic niche, we randomly sampled 1,000 points in the current species distribution area and computed the density (i.e. frequency), mean values, and 95% quantiles for each of the four studied bioclimatic variables, as well as for elevation. Elevation data in Madagascar was obtained from the SRTM (Shuttle Radar Topography Mission) 90 m Digital Elevation Data available from the CGIAR-CSI GeoPortal. Elevation data was aggregated at 1 km resolution and is also available for download on the MadaClim website.

Ensemble forecasting and future species distribution

To predict species distribution area in the future, we used an ensemble forecasting approach (Araujo & New 2007). For each of the two RCPs separately (RCP 4.5 and RCP 8.5), we combined climatic projections obtained from the three different GCMs (NorESM1-M, GISS-E2-R, and HadGEM2-ES) and across the four modelling algorithms (GLM, GAM, Random Forest, and Maxent). We thus obtained, for each species under each RCP and for each year (2055 and 2085), 12 maps of the future probability of presence. Again, the probability of presence was converted into binary data (0 for species absence and 1 for species presence) using the same probability threshold which maximizes TSS during current conditions. Species distribution area in the future (SDAf) was defined as the set of 1-km 2 pixels where most projections (6 out of 12) predicted the presence of the focal species. When less than 6 models out of the 12 predicted a presence, the species presence in the future was considered uncertain. The species was considered absent in the future when none of the 12 projections predicted a presence.

When predicting future range maps under each RCP for 2055 and 2085, we also considered two contrasting dispersal hypotheses, with the reality likely to fall in between. The full-dispersal hypothesis considers the possibility for all baobab species to colonize new climatically favorable sites outside the current species distribution range. The zero-dispersal hypothesis considers the impossibility for baobab species to naturally colonize new climatically favorable sites outside the current species distribution range. This can be due to unsuitable conditions (other than climate, such as land-use) outside the current species distribution area, or to that baobab species may not be able to disperse seeds due to geographical barriers or in the absence of animal dispersers (see [START_REF] Vieilledent | Vulnerability of baobab species to climate change and effectiveness of the protected area network in Madagascar: Towards new conservation priorities[END_REF].

Species range shift and vulnerability to climate change

To evaluate the effect of climate change and the vulnerability of the seven Malagasy baobab species, we calculated the percentage of area change between the future (SDAf) and present (SDAp) distribution range (in km²). To compute SDAf, SDAp, and mean elevational shifts, we extracted all presence points indicated as 'presence' by the ensemble modelling and calculated the changes for each of the investigated future scenarios (mean, 95% quantiles interval, and percentage of area change in km²). Focusing on the year 2085, under RCP 8.5 for the full-dispersal and zero-dispersal scenarios, we suggested updates for baobabs conservation management strategies given their future distribution and vulnerability to climate change according to the International Union for Conservation of Nature Red List (IUCN, 2012a). Finally, we calculated potential latitudinal and elevational species range shifts by extracting 1,000 random points inside each species distribution range projected for the present and the future (2055 and 2085), under both RCPs 4.5 and 8.5 scenarios and for the full-dispersal hypothesis.

Results

Range contraction and vulnerability of baobab species to climate change

For all seven baobab species we obtained high True Skill Statistics (TSS) values for both the model cross-validation (Table S2) and the ensemble model (Table S3). For the ensemble model, the TSS was >=0.83 for all the species except for A. za, for which the TSS was equal to 0.67 due to a relatively lower specificity (Table S3). High TSS values indicate that the models can then be confidently used to predict the vulnerability of species to climate change. Four baobab species are expected to be highly vulnerable to climate change under RCP 8.5, whereas for RCP 4.5 three species are expected to be highly vulnerable. Adansonia perrieri and A. suarezensis might experience a complete range loss modelled by 2085 (under RCPs 4.5 and 8.5), and could face extinction, under both the full and zerodispersal hypothesis (Figures 1,S2; Tables 1,S4). Still, under RCP 8.5, A. madagascariensis and A. rubrostipa could experience a contraction in modelled suitable range by 2085, and thus might be threatened by extinction, under both the full (≥ -32% area) and zero-dispersal hypothesis (≥ -71% area).

The three other baobab species modelled under RCP 8.5 (A. za, A. grandidieri and A. digitata), are expected to be resilient to climate change (Figure S1). These species might experience a small range contraction (down to -3%) under the zero-dispersal hypothesis and a strong range expansion (from +118% for A. za up to +300% for A. grandidieri and A. digitata) when modelled under the full-dispersal hypothesis (Table 1). Under RCP 4.5 and zero-dispersal hypothesis A. madagascariensis is predicted to lose 8% of its current distribution area, indicating that the species might not be threatened (Figure S2).

The projections under RCP 4.5 and full-dispersal hypothesis also indicated a strong range expansion from +78% for A. za up to +268% for A. grandidieri and A. digitata (Figure S2; Table S4). Among these three species, only A. digitata is expected to contract its modelled suitable range by 2085 under zerodispersal hypothesis: -20% (Figure S2). Both A. grandidieri and A. za might not be affected according to this predictive scenario.

Baobab species range shift in latitude and elevation

Among the four vulnerable baobab species, A. suarezensis and A. perrieri are expected to contract their range equatorward in 2085 under RCP 4.5 and for the full-dispersal hypothesis (both species are expected to go extinct under RCP 8.5 whatever the dispersal hypothesis). Adansonia madagascariensis is also expected to move equatorward while A. rubrostipa is expected to move poleward under RCP 8.5 and full-dispersal hypothesis (Figures 1,2). For the three resilient species (A. za, A. grandidieri and A. digitata), their distributional ranges should expand equatorward (A. za), poleward (A. digitata), or in both directions for A. grandidieri under RCP 8.5 and the full-dispersal hypothesis (Figure S1).

Regarding shifts in elevation, all baobab species are expected to shift upward in the future (from +100 m for A. za to +351 m for A. rubrostipa; Figures 2 andS3) when considering the full-dispersal hypothesis and scenario RCP 8.5 (except for A. perrieri and A. suarezensis where we considered RCP 4.5). When considering RCP 8.5 and the zero-dispersal hypothesis all threatened baobab species are expected to shrink their distribution upward (Figure 2; Table 1). Under RCP 4.5 for 2085, baobabs are expected to remain at current elevational gradient or to move upwards until 2085 (Figure S4). The potential redistribution of these species in 2085 and RCP 4.5 indicates that the A. digitata and A. rubrostipa might move poleward, whereas the other species are expected to remain in similar latitudes in the future (Figure S4).

Climatic gradients and future climatic anomalies in Madagascar

We identified four main climatic gradients in Madagascar (Figure 3) which are important to subsequently interpret species distribution on the island. A first gradient shows a decrease of the mean annual temperature with elevation. A second one shows an increase in temperature seasonality associated with higher latitude (lower seasonality at the North, toward the Equator, higher seasonality at the South, toward the South pole). A third gradient shows a decrease in annual precipitation from East to West associated with dominant Eastern winds and orographic precipitation (higher precipitation in the East). Finally, a fourth gradient shows an increase of the climatic water deficit from East to West due to the combining effects of precipitation and temperature. Computation of future climatic anomalies in 2085 under RCP 8.5 shows a general increase in the mean annual temperature over the whole Madagascar (from +2.5 to +4.0 °C), with a stronger increase in the inner-land than on the coast. Temperature seasonality should also generally increase over the whole Madagascar, especially at the North of the island where the temperature seasonality anomaly should reach up to +3 °C. Precipitation should generally decrease over the island (from 0 to -300 mm/yr) with a stronger decrease in the North-East. Associated with the general increase in temperature and decrease in precipitation, the climatic water deficit should generally increase (from 0 to +1500 mm/yr) over the island and decrease should be stronger in the Western part of Madagascar (> 500 mm/yr).

Importance of each bioclimatic variable in explaining species redistribution

For A. madagascariensis, A. perrieri, and A. suarezensis, three out of the four potentially threatened baobab species in 2085 under RCP 8.5, the most important variable for explaining species distribution and thus redistribution was temperature seasonality (Figure 4 and Table 2). These three species are currently distributed in the Northern part of Madagascar (Figure 1) and might experience, by 2085 under RCP 8.5, a strong increase in temperature seasonality (from +1,13 to +1,48 °C) inside their current distribution range (Figures 1 and4, Table S5). For the fourth threatened species (A. rubrostipa), the most important variable was climatic water deficit (Table 2). This species could experience, by 2085 under RCP 8.5, a strong increase in climatic water deficit inside its current distribution range (+870 mm/yr; Figures 1, 3 and4). The second most important variable for the four threatened species were either annual mean temperature (A. madagascariensis), mean annual precipitation (A. rubrostipa and A. suarezensis), or climatic water deficit (A. perrieri).

For A. grandidieri, and A. za, two out of the three non-threatened baobab species under RCP 8.5 in 2085, the most important variable for explaining these species distribution and redistribution was mean annual precipitation (Figure 4 and Table 2). Mean annual precipitation should not significantly change (from -31 to -70 mm/yr) inside the current distribution range of A. grandidieri and A. za (Figures 3,4, and Table S5), which are two species already adapted to dry climate (precipitation < 1000 mm/yr, Figure S5; Table S5). The second most important variable for these two species was annual mean temperature which is expected to significantly increase inside the two species' distribution range in the future (+3.5 °C, Figure 3; Table S5), although not affecting their distribution (Figure 1). For A. digitata, the third non-threatened species, the two most important variables were temperature seasonality (+0,87 °C) and the climatic water deficit (+680 mm/yr). These two variables are expected to increase inside the species distribution range in the future (Figures 3,4 and Table S5), although not affecting its distribution.

Discussion

Vulnerability of baobab species to climate change and conservation status

We showed that four out of the seven Malagasy baobab species are expected to experience a strong range contraction under the effect of climate change (> 70% for year 2085 under RCP 8.5 and the zerodispersal hypothesis) and could be strongly threatened with extinction according to our predictive scenarios. These four species are: A. madagascariensis, A. perrieri, A. suarezensis, and A. rubrostipa. Among these four species, A. perrieri and A. suarezensis could face a complete loss of their habitat by 2085 due in particular to an increase in the temperature seasonality in the future. The three other Malagasy baobab species, A. grandidieri, A. za, and A. digitata did not indicate any significant range contraction when modelled under climate change scenarios, except for A. digitata which might reduce its modelled distribution by 20% in 2085 under RCP 4.5 (Figure S2). For A. grandidieri and A. za, this resilience can be easily explained. A. za is a generalist species that can be found in a large range of climatic conditions and has a large distribution over Madagascar (Figures S1,S5). For A. grandideri, while it has a much narrower climatic niche than A. za and can be considered as a specialist species, it is already adapted to hot and dry climates (Figures 4,S1 and S5). Conversely, we did not find a simple explanation for the resilience of A. digitata under RCP 8.5 in 2085. Temperature seasonality and climatic water deficit (which are expected to increase in the future, Figure 3) were the most important climatic variables in explaining the distribution of this species (Figure 4). We hypothesize that the potential combination of the four climatic variables (which are rather close in terms of importance, Table 2) determines the large suitable habitat for the species in the future and its predicted resilience to climate change (Figure S1).

These results are in line with those obtained by [START_REF] Vieilledent | Vulnerability of baobab species to climate change and effectiveness of the protected area network in Madagascar: Towards new conservation priorities[END_REF] who have previously demonstrated, with different statistical models and IPCC climate scenarios, that both A. perrieri and A. suarezensis should become extinct by 2085 due to climate change, and that A. grandidieri should not be vulnerable to climate change. In a recent study, [START_REF] Wan | Modeling impacts of climate change on the potential distribution of six endemic baobab species in Madagascar[END_REF] found contradictory results to ours: a resilience of A. perrieri and A. suarezensis, and a vulnerability of A. za to climate change. However, these conclusions were taken by using a much smaller dataset (245 occurrence points distributed among the 6 endemic Malagasy baobab species) than ours (4830 in total, see Table S1). In addition, their statistical approach was very limited in comparison to ours. [START_REF] Wan | Modeling impacts of climate change on the potential distribution of six endemic baobab species in Madagascar[END_REF] used only one algorithm (Maxent), while we used four algorithms for both the ensemble modelling and the forecasting approach. In particular, the vulnerability of A. za to climate change found by [START_REF] Wan | Modeling impacts of climate change on the potential distribution of six endemic baobab species in Madagascar[END_REF] seems to be in contradiction with its known biology and distribution as a generalist species (Figures S1,S2, S5 and Table S5).

In light of these results, we recommend updating the IUCN conservation status for the four threatened baobab species (Table 1) based on the risk assessment under RCP 8.5. We base our recommendations on the IUCN Red List Categories and Criteria version 3.1 (IUCN, 2012b). We recommend updating A. madagascariensis from "Near Threatened" to "Endangered A3c" (population reduction >= 50% in the future). For A. suarezensis, we recommend updating the conservation status from "Endangered B1ab + B2ab" to "Critically Endangered B1ab + B2ab" (complete extent of occurrence loss by 2085). For A. perrieri, we recommend updating the conservation status from "Critically Endangered C2a(i)" to "Critically Endangered C2a(i) + A3c" (few mature individuals and potential extinction in the long-term). Finally, for A. rubrostipa, we recommend updating the conservation status from "Least Concern" to "Endangered A3c" (> 85% habitat loss in 2085). Climate change is not the only threat to Malagasy baobab species. The seven baobab species are also severely threatened by habitat loss associated with the rampant deforestation in Madagascar [START_REF] Vieilledent | Combining global tree cover loss data with historical national forest cover maps to look at six decades of deforestation and forest fragmentation in Madagascar[END_REF] which prevents species from regenerating.

Species range shifts in latitude and elevation under climate change

We have shown in our study that temperature seasonality was the most important variable at explaining species distribution for three out of the four threatened Malagasy baobab species. These three species (A. madagascariensis, A. perrieri, and A. suarezensis) are currently distributed in the Northern part of Madagascar, close to the Equator line, where the seasonality is lower. These three baobab species might experience, by 2085 under RCP 8.5, a strong increase in temperature seasonality. This strong increase in temperature seasonality is expected to be general to all the Northern region of Madagascar. To track the change in temperature seasonality, these three species might move equatorward, where the temperature seasonality is lower. Specifically for A. madagascariensis, mean annual temperature was the second most important variable explaining its distribution, and could also influence the species redistribution equatorward, thanks to suitable habitats in Northeastern Madagascar in 2085. Several studies on climate change have considered that the general trend for species, under the effect of climate change, will be to shift their distribution upward or poleward to escape from the increasing mean temperature globally [START_REF] Chen | Rapid range shifts of species associated with high levels of climate warming[END_REF][START_REF] Lenoir | A significant upward shift in plant species optimum elevation during the 20th century[END_REF][START_REF] Parmesan | A globally coherent fingerprint of climate change[END_REF][START_REF] Pecl | Biodiversity redistribution under climate change: Impacts on ecosystems and human wellbeing[END_REF][START_REF] Vanderwal | Focus on poleward shifts in species' distribution underestimates the fingerprint of climate change[END_REF]. Using emblematic Malagasy baobab species as an example, we demonstrate that this might not always be the case. Depending on both the bioclimatic variables that preferentially determine their distribution (the temperature seasonality in our study) and the future climatic anomalies (increase in temperature seasonality in our study), some species are expected to move in the opposite direction under the effect of climate change, i.e. equatorward.

In a review article, [START_REF] Lenoir | Going against the flow: Potential mechanisms for unexpected downslope range shifts in a warming climate[END_REF] have examined the potential mechanisms that could push species to go "against the flow" under climate change. They underlined the importance of additional ecological processes, in addition to climate change, to explain observed downslope range shifts in a warming climate. These mechanisms encompass biotic interactions (release of the competition associated with species range shift under climate change) and land-use change (new suitable habitats available downward). In our study, we show that it is not necessary to invoke other processes than climate change to explain shifts in species distribution in opposite directions. Studying the past range shift of 464 Australian bird species, [START_REF] Vanderwal | Focus on poleward shifts in species' distribution underestimates the fingerprint of climate change[END_REF], showed that complex interactions between temperature, precipitation, and species-specific tolerances could result in multi-directional distribution shifts, including equatorward. In our study, we illustrate one simple climatic mechanism, based on the change in temperature seasonality, by which species can shift their distribution equatorward.

We have also shown that it is not contradictory for a species to move both equatorward and upward under climate change, as is the case for A. perrieri, A. suarezensis, and A. madagascariensis (Figure 2). This counter-intuitive range shift has already been reported for sub-mountainous forest plant species in France which have shifted their distribution both southward (i.e. equatorward in France) and upward [START_REF] Kuhn | Early signs of range disjunction of submountainous plant species: An unexplored consequence of future and contemporary climate changes[END_REF]). The explanation lies in the presence of mountainous areas towards the south of the species' current distribution areas. In Madagascar, the explanation is different. We have shown that temperature seasonality is much more correlated to latitude than to elevation and that conversely, mean annual temperature is much more correlated to elevation than to latitude (Figure S6). Consequently, a species can shift its distribution both equatorward and upward to track changes in both temperature seasonality and mean annual temperature, respectively. Moreover, assuming that a species moves towards the equator to track changes in temperature seasonality, it might be that the lands towards the equator are located at higher elevation, thus leading to an upward shift of the species.

Vulnerability of tropical species to change in temperature seasonality

Our findings could have strong implications regarding species response to climate change in the tropics. In tropical regions, species are adapted to low temperature seasonality (Hua, 2016;[START_REF] Janzen | Why mountain passes are higher in the tropics[END_REF][START_REF] Pacifici | Species' traits influenced their response to recent climate change[END_REF][START_REF] Sheldon | The influence of temperature seasonality on elevational range size across latitude: a test using[END_REF]. Because sunlight duration (~12 hours a day) and solar incidence do not change significantly throughout the year, the temperature seasonality in tropical regions close to the Equator is narrower in comparison with subtropical or temperate regions (Figure 5a). For instance, temperature seasonality influences plant species biology and traits as it determines the length of the growing season and their phenology, such as the date of foliation, flowering and fruiting [START_REF] Pacifici | Species' traits influenced their response to recent climate change[END_REF][START_REF] Wright | Phenological Responses do Seasonality in Tropical Forest Plants[END_REF].

Here we have shown that an increase in temperature seasonality could force species to shift their distribution equatorward. Looking at the projected change in temperature seasonality in 2085 under RCP 8.5, a general increase in temperature seasonality across the tropics is expected (up to +10°C for the standard deviation of the monthly temperatures), with a particularly marked change in the Amazon region (Figure 5b). As for A. perrieri and A. suarezensis in Madagascar, the species redistribution equatorward to track changes in temperature seasonality might be impeded by several geographic and climatic barriers. This might also happen through several tropical lands located on islands in the Caribbean, Indian Ocean, or Southeast Asia for example. On these islands, the absence of land equatorward could act as a geographical barrier for species moving equatorward due to climate change. In addition, much of the tropical natural areas have been degraded, largely because of deforestation [START_REF] Hansen | The fate of tropical forest fragments[END_REF]. This could also prevent species from finding suitable habitats equatorward in tropical forests. Finally, species already distributed at the equator will not be able to move toward areas with lower seasonality, in analogy with species already at the top of the mountain which cannot shift their distribution upward to track temperature changes.

To conclude, our study shows that not all species should migrate poleward or upward as the climate warms, which reinforces the results of previous studies [START_REF] Kuhn | Early signs of range disjunction of submountainous plant species: An unexplored consequence of future and contemporary climate changes[END_REF][START_REF] Vanderwal | Focus on poleward shifts in species' distribution underestimates the fingerprint of climate change[END_REF]. Depending on the relative importance of the bioclimatic variables at explaining species distribution and regional climatic anomalies, shifts in species distribution can be multi-directional, including equatorward. We also underline the importance of the increase in temperature seasonality in the tropics that could potentially drag to extinction a large number of species adapted to low seasonality, among which three out of the seven emblematic baobab species of Madagascar. predicting the species presence. For the distribution maps, the species is assumed to be present (green areas) when a majority of models predicts a presence (votes >= 2 in the present, and >= 6 in the future). The species is considered absent (grey areas) when no model (votes = 0), or a minority of models (votes < 2 in the present, and Decrease in precipitation is expected to be stronger in the East (between -150 and -300 mm.y -1 ). and a strong increase in the climatic water deficit. indicates the number of models (0-12) predicting the species presence. For the distribution maps, the species is assumed to be present (green areas) when a majority of models predicts a presence (votes >= 2 in the present, and >= 6 in the future). The species is considered absent (grey areas) when no model (votes = 0), or a minority of models (votes < 2 in the present, and < 6 in the future), predicts a presence. Maps for A. perrieri and A.

Declarations

Figures

suarezensis, two species distributed at the extreme North of Madagascar, have been zoomed in (black squares).

Figure 1 .

 1 Figure 1. Species range contraction under climate change for the four threatened baobab species under RCP 8.5. The four species are A. madagascariensis, A. perrieri, A. rubrostipa, and A. suarezensis (one species per row).

  Occurrence points over Madagascar elevation map (elevation in m); (b, f, j, n) Current predicted species distribution. Legend indicates the number of models (0-4) predicting the species presence; (c, g, k, o) Projected species distribution in 2085 under scenario RCP 8.5 and the full-dispersal hypothesis. Legend indicates the number of models (0-12) predicting the species presence; (d, h, l, p) Projected species distribution in 2085 under scenario RCP 8.5 and the zero-dispersal hypothesis. Legend indicates the number of models (0-12)

Figure 2 .Figure 3 .

 23 Figure 2. Change in elevation and latitude for the most threatened baobab species. We randomly sampled 1000 points inside the species predicted occurrence area in the present and in the future (due to the extremely reduced distribution area for A. perrieri and A. suarezensis, we only sampled 416 and 105 points, respectively for year 2085, and another 15 points for A. suarezensis in 2055). For A. madagascariensis and A. rubrostipa we considered the scenario RCP 8.5 and the full-dispersal hypothesis. Under RCP 8.5, both A. perrieri and A. suarezensis became extinct in 2055 and 2085. As a consequence, we used RCP 4.5 to show change in elevation and latitude for these two species. (a, b) A. madagascariensis could migrate to higher elevation under climate change scenarios for 2055 and 2085 and also change its latitudinal range to lower latitudes (i.e. equatorward) in 2085. (e, f) Adansonia rubrostipa is expected to move to higher elevations and shift its latitudinal range to higher latitudes, i.e. poleward. (c, g) Both A. perrieri and A. suarezensis could shift their range to more elevated areas. (d) Adansonia perrieri might shift its range towards lower latitudes, i.e. equatorward. (h) Adansonia suarezensis is expected to retain its niche at lower latitudes

Figure 4 .

 4 Figure 4. Comparison of current (blue) and future (green) bioclimatic conditions experienced by each species within the current extent of their respective geographical ranges. We selected only the two most important climatic variables determining species distribution. Left column (panels a, c, e, g, i, k, m) shows the first most important variables; Right column (panels b, d, f, h, k, l, n) shows the second most important variables. Horizontal axis represents one of the four bioclimatic variables: mean annual temperature (temp -°C x 10), temperature seasonality (°C, standard deviation [sd] x 100), mean annual precipitation (mm/y), and climatic water deficit (mm). Vertical axis shows the distribution of values for that bioclimatic variable. Density plots show current (blue density plots) and future (2085, RCP 8.5, mean of 3 GCMs; green density plots). Dark-green shaded areas show current and future overlapping values. Bioclimatic envelope (current and future) of each variable was calculated within the current species distribution area (1000 random points extracted within occurrence areas indicated by the Ensemble approach). Vertical lines: dashed lines represent 95% bioclimatic envelope of future (green) and current (blue) variables within each species current distribution area; solid lines represent the computed mean value for current (blue) and future (green) bioclimatic data.

Figure 5 .

 5 Figure 5. Map of temperature seasonality and future anomaly in the tropics. (a) Map of the temperature seasonality at 10' resolution across the tropics. Black dashed line represents the equator. Intertropical regions in South America, Africa, Southeast Asia, and Oceania have similar low temperature seasonality values and temperature seasonality is decreasing from the poles to the equator. (b) Temperature seasonality anomaly obtained while comparing future temperature seasonality in 2085 under RCP 8.5 with current temperature seasonality. Future temperature seasonality in 2085 under RCP 8.5 is the mean of three GCMs (NorESM1-M,GISS-E2-R, and HadGEM2-ES). Temperature seasonality will increase throughout most of the tropics, with a particular strong increase in South America. (c) Future temperature seasonality in 2085 under RCP 8.5. (mean of the three selected GCMs). Despite changes in temperature seasonality in the future, the gradient of temperature seasonality, with a lower temperature seasonality at the equator, will be conserved in the future.

Figure S2 .

 S2 Figure S2. Species range contraction under climate change (RCP 4.5) for all baobab species. The seven species are A. digitata, A. grandidieri, A. madagascariensis (A. mada.), A. perrieri, A. rubrostipa, A. suarezensis, and A. za (one species per row). (a, e, i, m, q, u, y) Occurrence points over Madagascar elevation map (elevation in m); (b, f, j, n, r, v, z) Current predicted species distribution. Legend indicates the number of models (0-4) predicting the species presence;(c, g, k, o, s, w, a') Projected species distribution in 2085 under scenario RCP 4.5 and the fulldispersal hypothesis. Legend indicates the number of models (0-12) predicting the species presence; (d, h, l, p, t, x, b') Projected species distribution in 2085 under scenario RCP 4.5 and the zero-dispersal hypothesis. Legend

Figure S3 .

 S3 Figure S3. Change in elevation and latitude for the climate-resilient baobab species. We randomly sampled 1000 points inside the species predicted occurrence area in the present and in the future under RCP 8.5 scenario and full-dispersal hypothesis. (a, b) For A. digitata, the species should maintain stable its average elevational and latitudinal gradient from the current projection until 2085. (c, d) Adansonia grandidieri might move to higher elevations and shift its latitudinal range to lower latitudes. (e, f) Adansonia za could slightly shift its elevational gradient to more elevated areas in 2085 and higher latitudes, i.e. equatorward, from the present to 2085.

Figure S5 .

 S5 Figure S5. Bioclimatic niche of the seven baobabs species found in Madagascar. We computed species bioclimatic niche by randomly sampling 1000 points within each species distribution area for the four climatic variables selected in this study: annual mean temperature, temperature seasonality, mean annual precipitation, and climatic water deficit. Dashed lines (red) represent the climatic conditions through Madagascar. Baobab niche breadth indicates difference between species' climatic niche, with climatic specialist ones (A. digitata and A. madagascariensis -adapted to warmer temperatures -; A. grandidieri -adapted to higher seasonality in temperature and drier climate -; A. perrieri and A. suarezensis -adapted to lower seasonality in temperature), and generalist/ubiquitous species (A. rubrostipa and A. za), with a wider climatic niche which encompass almost every other baobab species bioclimatic niche.

< 6 in the future), predicts a presence. Maps for A. perrieri and A. suarezensis, two species distributed at the

extreme North of Madagascar, have been zoomed in (black squares).

Acknowledgments -Authors would like to thank all the people who have contributed to baobab occurrence data collection in Madagascar since the 2000s, in particular Wilfried Ramahafaly, Lucien Rasoanaivoson, Marcelo Saïdoune Maharitra, and Emilson Rakotoarisoa. MMT thanks to Juliano Bogoni and Maurício Cantor for their useful insights during the writing of the manuscript. MMT also thanks Emanuel Sinic and Sayuri Kyota for their tremendous support between 2019 and 2020.

Funding -This work was funded by the FRB -FFEM (Fondation pour la Recherche sur la Biodiversité -Fond Français pour l'Environnement Mondial) through the BioSceneMada project (project agreement AAP-SCEN-2013 I). MMT and VCR received their PhD grants from the "Coordenação de Aperfeiçoamento de Pessoal de Nível Superior'' -Brazil (CAPES) -Finance Code 001. MMT received additional fundings from MAPA's Private Fund -Number 8052-7. MMT also received a fellowship from CIRAD ("Action Incitative: Soutien aux doctorants du Sud'') for a three-month visiting period at UMR AMAP in Montpellier in 2019.

Data availability -Data (including the CIRAD baobab occurrence database in Madagascar) and code used to obtain the results of the present study are available on GitHub (https://github.com/ghislainv/baobabs_mada). The GitHub repository has been permanently archived in the CIRAD Dataverse (https://doi.org/10.18167/DVN1/LIALRR).

Conflicts of interest -The authors declare they have no conflicts of interest. Table 1. Baobabs' vulnerability to climate change and elevational range shift in 2085 under scenario RCP 8.5.

Tables

We calculated the species distribution area (km²) in the present (SDAp) and future (SDAf) to describe the change in the species distribution area (Change SDAp SDAf %) according to two dispersal hypotheses (full and zerodispersal). See IUCN Red List categories and criteria, version 3.1, second edition | IUCN Library System for criteria explanation. *IUCN criteria suggested for Madagascar; ** IUCN status defined by populational aspects for the referred species once our models did not predict vulnerability for this species.

Baobab species IUCN status

SDAp (km²) Legend indicates the number of models (0-12) predicting the species presence. For the distribution maps, the species is assumed to be present (green areas) when a majority of models predicts a presence (votes >= 2 in the present, and >= 6 in the future). The species is considered absent (grey areas) when no model (votes = 0), or a minority of models (votes < 2 in the present, and < 6 in the future), predicts a presence. Table S1. Presence data-set used for the species distribution modeling (SDM). We computed the number of presence points with GPS coordinates for each species. Considering a 1-km² spatial grid covering Madagascar, we identified the cells including at least one presence point for each species independently.

Supplementary Tables

Doing so we obtained a presence data-set at 1-km² resolution for each species. We removed cells with incomplete associated bioclimatic data from the presence data-set.

Species