Charlène Guillaumot
email: charleneguillaumot21@gmail.com

Alexis Martin

Marc Eléaume

Thomas Saucède

Methodological clues for improving species distribution models in data-poor areas: example of sub-Antarctic benthic species on the Kerguelen Plateau

Keywords: Species distribution modeling, model performance, historical datasets, Kerguelen Plateau, presence-only data

Methods for improving species distribution models in

data-poor areas: example of sub-Antarctic benthic species on the Kerguelen Plateau.

INTRODUCTION

Today, species distribution models (SDM) constitute essential tools for conservation biologists to understand species distribution patterns and their underpinning drivers (see [START_REF] Guillera-Arroita | Is my species distribution model fit for purpose? Matching data and models to applications[END_REF] for a review), assess the combined effects of environmental changes and direct human pressures (i.e. economic activities including tourism) on natural habitats [START_REF] Gutt | Correlative and dynamic species distribution modelling for ecological predictions in the Antarctic: a cross-disciplinary concept[END_REF], define conservation priorities (Vierod et al. 2014[START_REF] Greathead | Environmental requirements for three sea pen species: relevance to distribution and conservation[END_REF], and develop relevant management plans [START_REF] Reiss | Benthos distribution modelling and its relevance for marine ecosystem management[END_REF][START_REF] Koubbi | Ecoregionalisation of the Kerguelen and Crozet islands oceanic zone Part I: Introduction and Kerguelen oceanic zone[END_REF]. SDM allow to interpolate the known distribution of single species, assemblages or communities [START_REF] Ferrier | Spatial modelling of biodiversity at the community level[END_REF] to little-accessed or under-sampled areas [START_REF] Reiss | Species distribution modelling of marine benthos: a North Sea case study[END_REF][START_REF] Robinson | Pushing the limits in marine species distribution modelling: lessons from the land present challenges and opportunities[END_REF]) and help improve our knowledge of the distribution of rare species [START_REF] Mccune | Species distribution models predict rare species occurrences despite significant effects of landscape context[END_REF].

In regions subject to fast environmental changes and significant anthropogenic activities, SDM can provide useful tools for conservation purposes [START_REF] Guisan | Predicting species distributions for conservation decisions[END_REF][START_REF] Reiss | Benthos distribution modelling and its relevance for marine ecosystem management[END_REF]). However, modeling species distribution over vast and remote areas is challenging and questions the relevance of the method compared to more traditional and qualitative approaches [START_REF] Koubbi | Ecoregionalisation of the Kerguelen and Crozet islands oceanic zone Part I: Introduction and Kerguelen oceanic zone[END_REF]). In such regions, our knowledge of species distribution usually relies on historical and heterogeneous presence-only datasets, which concentrate many gaps and can induce methodological biases altering the level of SDM performance [START_REF] Loiselle | Predicting species distributions from herbarium collections: does climate bias in collection sampling influence model outcomes?[END_REF][START_REF] Costa | Sampling bias and the use of ecological niche modeling in conservation planning: a field evaluation in a biodiversity hotspot[END_REF][START_REF] Newbold | Applications and limitations of museum data for conservation and ecology, with particular attention to species distribution models[END_REF]. The use of historical data in SDM has been widely discussed [START_REF] Reutter | Modelling habitat suitability using museum collections: an example with three sympatric Apodemus species from the Alps[END_REF][START_REF] Hortal | Limitations of biodiversity databases: case study on seed-plant diversity in Tenerife, Canary Islands[END_REF], 2008), for instance with regards to the spatial and temporal heterogeneities induced by the practice of different sampling strategies. Limitations to SDM performance are mainly due to uncertainties in data location and detection [START_REF] Costa | Sampling bias and the use of ecological niche modeling in conservation planning: a field evaluation in a biodiversity hotspot[END_REF][START_REF] Naimi | Where is positional uncertainty a problem for species distribution modelling?[END_REF], Tessarolo et al. 2014), to over-estimations of habitat suitability in intensively sampled areas [START_REF] Guillera-Arroita | Is my species distribution model fit for purpose? Matching data and models to applications[END_REF], and to artefacts in niche descriptions [START_REF] Hortal | Historical bias in biodiversity inventories affects the observed environmental niche of the species[END_REF]. The lack of available data from remote areas also constitutes a limitation to SDM, which are restricted to presence-only data, and are regarded as less reliable and less efficient than presence-absence and abundance-based models [START_REF] Brotons | Presence-absence versus presenceonly modelling methods for predicting bird habitat suitability[END_REF]). Over the past few years, many methodological developments in SDM procedures have been produced to correct for such biases [START_REF] Dormann | Effects of incorporating spatial autocorrelation into the analysis of species distribution data[END_REF][START_REF] Phillips | Sample selection bias and presence-only distribution models: implications for background and pseudo-absence data[END_REF][START_REF] Barbet-Massin | Selecting pseudo-absences for species distribution models: how, where and how many? Method[END_REF]) but no single procedure emerged [START_REF] Qiao | No silver bullets in correlative ecological niche modelling: insights from testing among many potential algorithms for niche estimation[END_REF] and few practical solutions have been proposed to deal with poor and heterogeneous datasets.

Our knowledge of Southern Ocean species distribution remains patchy [START_REF] Koubbi | Ecoregionalisation of the Kerguelen and Crozet islands oceanic zone Part I: Introduction and Kerguelen oceanic zone[END_REF]. Therefore, the growing interest of marine biologists and biogeographers for the region has led to the conception of collaborative projects compiling past and present marine biodiversity data in information networks like the SCAR-Marine Biodiversity Information Network (SCAR-MarBIN) [START_REF] Griffiths | Quantifying Antarctic marine biodiversity: The SCAR-MarBIN data portal[END_REF], the Biogeographic Atlas of the Southern Ocean [START_REF] Broyer | Biogeographic Atlas of the Southern Ocean[END_REF]) and other open access databases [START_REF] Danis | Connecting biodiversity data during the IPY: the path towards e-polar science[END_REF][START_REF] Gutt | Antarctic macrobenthic communities: A compilation of circumpolar information[END_REF], Van de Putte et al. 2014). However, running species distribution models in the region still requires a significant effort of data compilation [START_REF] Guillaumot | Echinoids of the Kerguelen Plateau: occurrence data and environmental setting for past, present, and future species distribution modelling[END_REF] to complement the existing open access data sources, and check for data quality. In addition, modeling Southern Ocean species distribution poses auxiliary problems due to the paucity of data and model performances that can vary with ecological niche width [START_REF] Qiao | No silver bullets in correlative ecological niche modelling: insights from testing among many potential algorithms for niche estimation[END_REF]. Recent works have developed methodologies to adapt SDM to rare species and poorly-sampled areas but none was tested for the Southern Ocean (Pokharel et al. 2016[START_REF] Phillips | Applying species distribution modelling to a data poor, pelagic fish complex: the ocean sunfishes[END_REF].

In this work, we analysed the reliability of modeling procedures with regards to the heterogeneous nature of data available and the gaps in our knowledge of species distribution.

We compiled echinoid presence-only data collected from several ancient and recent oceanographic campaigns carried out on the Kerguelen Plateau (sub-Antarctic region) for one and a half century. The distribution of four echinoid species with contrasting ecological niches was modeled, and the reliability and the performance of modeling procedures were tested. We propose methodological clues to correct for spatial and temporal biases and assess the sensitivity of modeling procedures to species ecological niche width. This is the first methodological approach to correct for potential biases in SDM in the Southern Ocean. Our objective is to offer useful perspectives for future modeling works along with a practical and transferable protocol to test for the reliability and performance of modeling procedures.

MATERIAL AND METHODS

Biological data

Species occurrence data were taken from [START_REF] Guillaumot | Echinoids of the Kerguelen Plateau: occurrence data and environmental setting for past, present, and future species distribution modelling[END_REF][START_REF] Phillips | Applying species distribution modelling to a data poor, pelagic fish complex: the ocean sunfishes[END_REF].

The dataset includes presence-only data of echinoid species collected during scientific cruises carried out on the Kerguelen Plateau (63°/81°E; -46°/-56°S) since 1872 (Fig. 1). Scientific objectives, dates, sampling efforts, gears, and surveyed areas have differed between cruises, leading to spatial and temporal heterogeneities [START_REF] Guillaumot | Echinoids of the Kerguelen Plateau: occurrence data and environmental setting for past, present, and future species distribution modelling[END_REF]). In the dataset, four echinoid species with contrasting ecological preferences and a high number of presence-only records were selected. Species include two sediment feeders of the family Schizasteridae, one shallow water species, Abatus cordatus, and a deeper one, Brisaster antarcticus, one carnivorous/detritivorous and eurybathic species of the family Cidaridae, Ctenocidaris nutrix, and one omnivorous and eurybathic species of Echinidae, Sterechinus diadema [START_REF] David | Antarctic Echinoidea -Synopses of the Antarctic benthos[END_REF]) (Fig. 1). A. cordatus is a coastal species that is endemic to the Kerguelen Plateau, B.

antarcticus is known in the Kerguelen and Crozet archipelagoes and has broader environmental preferences than A. cordatus (Fig. 1). C. nutrix and S. diadema are widespread in the Southern Ocean and have contrasting environmental preferences (Fig. 1).

Environmental descriptors

Environmental descriptors were taken from [START_REF] Guillaumot | Echinoids of the Kerguelen Plateau: occurrence data and environmental setting for past, present, and future species distribution modelling[END_REF]. The dataset covers the geographic extent of the Kerguelen Plateau (63°/81°E and -46°/-56°S) and compiles environmental data for six decades included in . Environmental data are available at a grid cell resolution of 10km precision. Environmental layers include no data pixels, particularly in seafloor related descriptors. Data were not interpolated to avoid the potential biases due to interpolation procedures.

Collinearity between descriptors can alter modeling performances [START_REF] Phillips | Maximum entropy modeling of species geographic distributions[END_REF] because collinear data may (1) inflate standard errors, (2) induce the violation of residual independency during model validation and (3) generate noise that can be interpreted as a link between descriptors [START_REF] Dormann | Collinearity: a review of methods to deal with it and a simulation study evaluating their performance[END_REF]). To reduce the collinearity effect, we computed the Variance Inflation Factor (VIF) and the Spearman correlation coefficient (rs) between all available descriptors from [START_REF] Guillaumot | Echinoids of the Kerguelen Plateau: occurrence data and environmental setting for past, present, and future species distribution modelling[END_REF]. VIF analysis was performed using a stepwise procedure, using the vifstep function proposed in the 'usdm' R package [START_REF] Naimi | Where is positional uncertainty a problem for species distribution modelling?[END_REF]). Descriptor pairs with high VIF and rs values were omitted based on the commonly used thresholds of VIF < 5 and rs < 0.85 [START_REF] Phillips | Applying species distribution modelling to a data poor, pelagic fish complex: the ocean sunfishes[END_REF][START_REF] Dormann | Collinearity: a review of methods to deal with it and a simulation study evaluating their performance[END_REF][START_REF] Duque-Lazo | Transferability of species distribution models: The case of Phytophthora cinnamomi in Southwest Spain and Southwest Australia[END_REF]. Environmental descriptors finally selected to model species distribution are displayed in Table 1.

Environmental changes were tested between 1955 and 2012. The comparison of pixel values between periods was generated using a Wilcoxon signed-rank test with the Bonferroni correction.

Analytical procedures

The flow chart of Figure 2 details the analytical procedure used in the present work.

Model selection

Due to the growing interest of ecologists for species distribution modeling, a large range of modeling techniques is now available [START_REF] Reiss | Species distribution modelling of marine benthos: a North Sea case study[END_REF][START_REF] Guillera-Arroita | Is my species distribution model fit for purpose? Matching data and models to applications[END_REF][START_REF] Qiao | No silver bullets in correlative ecological niche modelling: insights from testing among many potential algorithms for niche estimation[END_REF]. Running the most appropriate model involves selecting the best modeling technique for the data under analysis and also involves considering the scientific objectives to be addressed [START_REF] Reiss | Species distribution modelling of marine benthos: a North Sea case study[END_REF][START_REF] Qiao | No silver bullets in correlative ecological niche modelling: insights from testing among many potential algorithms for niche estimation[END_REF].

Here we compared several modeling techniques using the 'biomod2' R3.3.0 library (Thuiller et al. 2016) and we tested the performance of these approaches with regards to the chronological addition of new data and to the transferability of models between areas. Several models were generated with an increasing number of occurrence data (Fig. S1). The best modeling techniques were then compared with each other using a non-random crossvalidation procedure (Fig. S2, Wengen and Olden 2012) in order to determine the approach with the best accuracy in transferability performances (Randin et al. 2006, Wengen and[START_REF] Wenger | Assessing transferability of ecological models: an underappreciated aspect of statistical validation[END_REF].

Results show high performance and stability values for Random Forest (RF) and Boosted Regression Trees (BRT) in our case study (Appendix 1). However, BRT performed better in transferability in comparison with RF [START_REF] Heikkinen | Does the interpolation accuracy of species distribution models come at the expense of transferability?[END_REF], and previous works showed that RF does not deal correctly with missing values and patchy datasets [START_REF] Breiman | Random forests[END_REF][START_REF] Barbet-Massin | Selecting pseudo-absences for species distribution models: how, where and how many? Method[END_REF][START_REF] Qiao | No silver bullets in correlative ecological niche modelling: insights from testing among many potential algorithms for niche estimation[END_REF], see Table S1 for a review). Therefore, BRT was chosen in the present work to generate the analyses.

BRT calibration was realised using the 'gbm' R package [START_REF] Ridgeway | gbm: Generalized Boosted Regression Models[END_REF][START_REF] Elith | A working guide to boosted regression trees[END_REF].

The three main parameters (learning rate lr, tree complexity tc and bag fraction bg) were selected using the method developed by [START_REF] Elith | A working guide to boosted regression trees[END_REF] to figure out the combination of values minimizing the predicted deviance of the models [START_REF] Elith | Boosted Regression Trees for ecological modeling[END_REF]. The parameters were finally set at respectively lr= 0.0001, tc=2 and bf= 0.75.

Following [START_REF] Barbet-Massin | Selecting pseudo-absences for species distribution models: how, where and how many? Method[END_REF], we sampled the same number of background data as the number of presence data available for computing BRT models. Considering the low number of presence data available, 100 model replicates (i.e. background sampling) were generated for each analysis. Finally, to correct for data aggregation in space, presence duplicates were removed when present in a same 10km resolution pixel.

Model performance was assessed by measuring AUC score values (Area Under the Curve of the Receiver Operating Curve) of each model replicate using the 'dismo' R library [START_REF] Hijmans | dismo: Species Distribution Modeling. R package version 1[END_REF]. AUC expresses the relationship between model sensitivity and the commission error (1-specificity). The sensitivity corresponds to the number of « presence » pixels correctly predicted as present and the specificity the number of « absence » pixels correctly predicted as absent [START_REF] Fielding | A review of methods for the assessment of prediction errors in conservation presence/absence models[END_REF]. The use of the AUC to evaluate SDM performance has been repeatedly discussed [START_REF] Lobo | AUC: a misleading measure of the performance of predictive distribution models[END_REF], Peterson et al. 2008) but the AUC remains the most appropriate metric for presence-background models as values stay stable with low-prevalence datasets and are not sensitive to threshold effects [START_REF] Hand | Measuring classifier performance: a coherent alternative to the area under the ROC curve[END_REF][START_REF] Proosdij | Minimum required number of specimen records to develop accurate species distribution models[END_REF]. Following the recommendation of Jimenez-Valverde et al. (2012), we used the AUC to estimate the robustness of models but not for direct comparisons between models that were generated for different species, on different studied areas and with different training samples.

Correcting for sampling bias

Data collected during the various scientific cruises led over the Kerguelen Plateau for the last 145 years present conspicuous spatial heterogeneities. The resulting bias can generate an unequal number of records in the different sectors of the study area and heterogeneous patterns in record distribution. Such heterogeneities might increase the risk of over-estimating the contribution of environmental conditions to the models in the most sampled areas [START_REF] Araújo | Five (or so) challenges for species distribution modelling[END_REF].

The effect of spatial heterogeneities on the quality of distribution models was tested using a null model approach. A first null model, null model #1, was generated by sampling presence data at random within the total set of sites that were visited during the different campaigns, whether echinoid specimens were collected at these sites or not (Fig. 3). Because absence data are not available, this approach allows us to assess the weight of sampling bias in the models.

If a sampling bias is significant, null model #1 is expected to produce distribution maps with higher suitability values in the most sampled areas [START_REF] Merckx | Null models reveal preferential sampling, spatial autocorrelation and overfitting in habitat suitability modelling[END_REF].

A second null model, null model #2, was built by simulating presence data sampled at random over the entire studied area. Null model #2 is expected to produce distribution maps of equal suitability over the entire study area. If sampling is spatially biased, we expect that null model #1 deviates from null model #2 [START_REF] Raes | A null-model for significance testing of presence-only species distribution models[END_REF].

The two null models were generated for the four selected species. The number of presenceonly data used in the models was contained between the number of data collected until the MD04 campaign and until the PROTEKER campaign, between 1974 and 2015, which corresponds to periods of high sampling effort (Fig. 1). In each null model, 100 replicates were produced. Time-averaged environmental descriptors were used for the analysis.

To correct for sampling bias when null models #1 and #2 significantly differ between each other, we used the methodology proposed by [START_REF] Phillips | Sample selection bias and presence-only distribution models: implications for background and pseudo-absence data[END_REF], which has been shown to improve modeling performances [START_REF] Phillips | Sample selection bias and presence-only distribution models: implications for background and pseudo-absence data[END_REF][START_REF] Aguirre-Gutiérrez | Fit-forpurpose: species distribution model performance depends on evaluation criteria-Dutch hoverflies as a case study[END_REF]. A grid layer was built using a kernel density estimation (KDE) to represent sampling spatial bias.

The layer was calculated from the map of visited sites. The estimated proportion of presenceonly data present in each pixel was determined using the 'kde2d' function of 'MASS' R package (Venables and Ripley 2002). Background data were sampled according to the weighting scheme of the KDE layer, to reduce discrepancies between presence-only records and background data [START_REF] Phillips | Sample selection bias and presence-only distribution models: implications for background and pseudo-absence data[END_REF][START_REF] Barbet-Massin | Selecting pseudo-absences for species distribution models: how, where and how many? Method[END_REF]. In order to test for the efficiency of model correction based on the KDE, the Pearson r correlation was computed between pixel values of the KDE layer (a proxy for the sampling effort) and the predicted probabilities of models, before and after the KDE correction.

Spatial heterogeneities in data collection can also generate spatial autocorrelation (SAC)

between presence records, which can violate model calibration assumptions, and affect model accuracy with wrong parameter estimations [START_REF] Segurado | Consequences of spatial autocorrelation for niche-based models[END_REF][START_REF] Dormann | Effects of incorporating spatial autocorrelation into the analysis of species distribution data[END_REF][START_REF] Crase | A new method for dealing with residual spatial autocorrelation in species distribution models[END_REF]). Several approaches have been developed to take into account SAC in SDM [START_REF] Crase | A new method for dealing with residual spatial autocorrelation in species distribution models[END_REF] for a review). They consist in including an additional term in models (the autocovariate), which represents the influence of neighboring records on modeling predictions.

The significance of SAC was tested using the Moran I autocorrelation index computed on model residuals [START_REF] Luoto | Uncertainty of bioclimate envelope models based on the geographical distribution of species[END_REF][START_REF] Crase | A new method for dealing with residual spatial autocorrelation in species distribution models[END_REF] for both original and corrected models. Models were built using time-averaged environmental descriptors .

Testing for the effect of the chronological addition of new records on model performance

Our dataset compiles presence-only data collected during various scientific cruises and with distinct sampling protocols, which may alter the performance of the models (Fig. 1). To test for model reliability, we separately analysed (1) the influence of the chronological addition of presence records, (2) the influence of data number alone and (3) the influence of sampling patterns (distribution of data in space). The analyses were performed for A. cordatus, C. nutrix and S. diadema. Not enough data were available for B. antarcticus. We used timeaveraged environmental descriptors to generate the models.

To test for the potential effect of the chronological addition of new data on model performance, we followed the protocol proposed by [START_REF] Aguiar | Effect of chronological addition of records to species distribution maps: The case of Tonatia saurophila maresi (Chiroptera, Phyllostomidae) in South America[END_REF]. The dataset was split into distinct subsets corresponding to main periods of sampling effort (1975: including Marion Dufresne campaigns, 1993: including ANARE campaigns, 2010: including POKER II campaign, 2015: including PROTEKER campaigns). New presence data were progressively added to the models, following the chronological collection of new records. The influence of the chronological addition of data was assessed by measuring the correlation between models using the Schoener's D statistic. The Schoener's D is a correlation metric adapted to the study of niche similarities (Warren et al. 2008, Rödder and[START_REF] Rödder | Quantitative metrics of overlaps in Grinnellian niches: advances and possible drawbacks[END_REF]. It evaluates the similarity of pixel values between two distribution grids. A D value of 0 means that the two maps are perfectly different, a D value of 1 means that maps are perfectly similar. Values were computed using the niche.overlap function of the 'ENMeval' R package [START_REF] Muscarella | ENMeval: An R package for conducting spatially independent evaluations and estimating optimal model complexity for ecological niche models[END_REF]).

The significance of correlations was tested following a null model protocol, using 100 replicates, pairwise-compared using the Schoener's D statistic [START_REF] Raes | A null-model for significance testing of presence-only species distribution models[END_REF][START_REF] Warren | Environmental niche equivalency versus conservatism: quantitative approaches to niche evolution[END_REF][START_REF] Ficetola | From introduction to the establishment of alien species: bioclimatic differences between presence and reproduction localities in the slider turtle[END_REF].

The distinct effect of data addition and sampling patterns were tested separately. To test for the effect of data addition alone, models were built by sampling an increasing number of presence data at random in the total area for A. cordatus [n=54, 76, 95], C. nutrix [n=46, 54, 106, 114] and S. diadema [n=54, 66, 98]. These thresholds correspond to the number of presence-only data used in the chronological addition analysis.

Finally, to test for the effect of sampling patterns, different models were produced by sampling presence data at random either within a subset of real data collected along transects (MD03 campaign) or within a subset of real data collected at random (POKER, PROTEKER campaigns). All models were compared between each other.

Testing for the effect of temporal variations on model performance

To test for the effect of environmental shifts on models, different distribution models were generated using distinct environmental descriptors for four periods ([1955-1964]; [1965][1966][1967][1968][1969][1970][1971][1972][1973][1974]; [1975][1976][1977][1978][1979][1980][1981][1982][1983][1984][1985][1986][1987][1988][1989][1990][1991][1992][1993][1994]; [2005][2006][2007][2008][2009][2010][2011][2012]) and the complete set of presence data available.

Similarities between models were measured using the Schoener's D statistic.

RESULTS

Environmental shifts

Mean sea surface temperature and amplitude, mean seafloor temperature and amplitude, mean sea surface salinity and amplitude were all tested significantly different between all the studied decades (p < 0.001). Seafloor temperature amplitude only was not proved significantly different between decades [2005][2006][2007][2008][2009][2010][2011][2012] and [1955][1956][1957][1958][1959][1960][1961][1962][1963][1964]. These results indicate that significant environmental shifts occurred during the studied time period and may induce important variations in models as the dataset extends over 145 years.

Spatial bias

Null model #1 predicts higher suitability values in areas with most intense sampling effort corresponding to the northern part of the Kerguelen Plateau and the vicinity of the Kerguelen archipelago (Fig. 3A). In contrast, null model #2 predicts equally medium suitability values over the entire Kerguelen Plateau because presence data were sampled at random in the area (Fig. 3B). The difference between null models #1 and #2 was tested significant for the four species (Fig. 3) showing that sampling bias has a significant impact on model outputs, which will over estimate environment suitability in areas with the highest number of sampling sites if no correction is applied.

Correlation between visited areas and predicted probability distribution decreases in models built with the KDE-correction compared to non-corrected models (Table 2), showing that the correction is efficient to reduce the influence of sampling bias on modeling performances. However, the correction was proved less efficient in models of the coastal and narrow niche species A. cordatus for which correlation values after the KDE correction remain high (r=0.44) (Table 2). Spatial autocorrelation (SAC) was tested significant for non-corrected models (Moran index, Imin=0.05, Imax=0.16) but values were not significant in corrected models (Imin=0.04, Imax=0.06), except for A. cordatus (Table S2, Figure S3). This shows that the KDE procedure also corrected for SAC for three of the four studied species.

Chronological addition of new records

The different models built with a chronological addition of new data show high AUC values for C. nutrix and A. cordatus (0.814±0.018 < AUCC.nutrix < 0.883±0.024 and 0.908±0.023 < AUCA.cordatus < 0.909±0.018 respectively) demonstrating the relevance of all models (Fig. 4, Fig. S4). For these two species, Schoener's D correlation values are high (𝐷 ̅ A.cordatus= 0.978±0.023, 𝐷 ̅ C.nutrix=0.968±0.020) and were tested significant, showing that models are similar to each other.

In contrast, models generated for S. diadema significantly differ between each other with lower Schoener's D statistics (𝐷 ̅ S.diadema = 0.932±0.036) (Fig. S5). Therefore, the chronological addition of new data has contrasting impacts on model outputs according to the studied species, which may be explained by a various sensitivity of models to data addition and to sampling patterns.

Data addition and sampling patterns

Comparison of models produced with an increasing number of data presents high and significant Schoener's D values (minimum=0.979 ± 0.031 for S. diadema, maximum=0.985 ± 0.020 for C. nutrix), showing that model outputs do not vary significantly with increasing data in our case study (Table 3).

To test for the influence of sampling patterns, models built using subsets with contrasting distribution patterns (radial versus random patterns) were compared. Schoener's D statistics measured between these two types of models present low values. This suggests a significant influence of sampling patterns on model outputs (Table 3).

Environmental change and model performance

The different models generated with contrasting environmental descriptors are highly similar as shown by high Schoener's D and low standard deviation values (𝐷 ̅ = 0.981±0.005). This proves that environmental shifts have no significant impact on model outputs. In addition, the respective contributions of environmental descriptors to models do not vary significantly between periods for the four species. However, A. cordatus seems to be less impacted by environmental shifts in comparison with other species (Fig. 5).

Finally, the contribution of time-averaged environmental descriptors over the total studied period tends to differ from contributions computed for each decade separately (Fig. 5).

Final species distribution models

Sampling bias analyses and model corrections show that reliable distribution models can be built for C. nutrix only. It is the only dataset in which spatial and temporal heterogeneities do not impact prediction performances significantly. A final, reliable model was produced for C. nutrix over the Kerguelen Plateau (Fig. 6).

DISCUSSION

Data scarcity and heterogeneity

First research surveys of the Kerguelen Plateau date back to the oceanographic campaign of the HMS Challenger in 1872. One and a half century later, our knowledge of benthic species distribution on the Kerguelen Plateau has significantly increased but remains patchy [START_REF] Koubbi | Ecoregionalisation of the Kerguelen and Crozet islands oceanic zone Part I: Introduction and Kerguelen oceanic zone[END_REF]. As in most parts of the Southern Ocean, modeling species distribution on the Kerguelen Plateau faces significant limitations due to data gaps and heterogeneities [START_REF] Guillaumot | Echinoids of the Kerguelen Plateau: occurrence data and environmental setting for past, present, and future species distribution modelling[END_REF]. Such limitations seriously question the relevance of modeling procedures, which are required by environmental managers for conservation purposes [START_REF] Féral | PROTEKER: implementation of a submarine observatory at the Kerguelen Islands (Southern Ocean) -Underwater Technol[END_REF][START_REF] Koubbi | Ecoregionalisation of the Kerguelen and Crozet islands oceanic zone Part I: Introduction and Kerguelen oceanic zone[END_REF]). In the present work, we follow a step by step protocol to assess, quantify, and correct the potential effects of data scarcity and heterogeneity on models, a critical issue when considering the growing interest for modeling approaches in Antarctic and Sub-Antarctic regions [START_REF] Gutt | Correlative and dynamic species distribution modelling for ecological predictions in the Antarctic: a cross-disciplinary concept[END_REF]. Our results demonstrate that such approaches can prove feasible and reliable in certain case studies, when data quality and sampling bias can be tested and corrected.

Coping with spatial and temporal bias in presence-only datasets

Spatial bias and spatial autocorrelation (SAC)

Building SDM for remote and little-accessed regions often implies the use of spatially biased datasets conditioned by sampling caveats. Because parts of these regions that are the most easily accessed aggregate most of the available presence data, more weight is given to the most sampled sites and model performance is reduced [START_REF] Phillips | Sample selection bias and presence-only distribution models: implications for background and pseudo-absence data[END_REF]). In the present work, a significant difference was measured between the two null models (that were generated by selecting presence data either from visited stations only or at random over the total investigated area), highlighting the strong heterogeneity of sampling effort with more data collected in the northern part of the Kerguelen Plateau and in coastal, shallow areas.

The significant spatial autocorrelation (SAC) values that were computed from model residuals also reveal the impact of sampling bias. The significance of SAC on uncorrected model residuals can be partly explained by the relative accumulation and high density of presence data in shallow areas of the Kerguelen Plateau where species presence probability is overpredicted. One could argue that SAC analysis does not apply to SDM as species presence proximities must be considered in the environmental niche space, not in the geography.

However, in the present study, the difference between null models constitutes an operational evidence of the impact of sample clumping on model outputs, which is also revealed by significant SAC values.

To correct for sampling bias, we used a background-based correction method [START_REF] Phillips | Sample selection bias and presence-only distribution models: implications for background and pseudo-absence data[END_REF]) that was already used in former studies based on presence-only and limited datasets [START_REF] Mateo | Profile or group discriminative techniques? Generating reliable species distribution models using pseudo-absences and targetgroup absences from natural history collections -Divers[END_REF], Pokharel et al. 2016[START_REF] Phillips | Applying species distribution modelling to a data poor, pelagic fish complex: the ocean sunfishes[END_REF]). These methods allow to reduce the effect of sample spatial bias on modeling performance by weighting background records according to sampling patterns. In the present study, the correction was proved efficient to correct both for the influence of the uneven sampling effort on predicted distributions (Table 2) and for SAC on all SDM except for models of A. cordatus. A. cordatus is a coastal, shallow marine species that was mainly sampled in the northern part of the Kerguelen Plateau. Species presence records are strongly conditioned by the location of most important sampling efforts. This is in line with previous studies that already highlighted the difficulties of modeling the distribution of narrow niche species with low prevalence distribution (i.e, corresponding to the proportion of the area where presence records are located) [START_REF] Barbet-Massin | Selecting pseudo-absences for species distribution models: how, where and how many? Method[END_REF][START_REF] Qiao | No silver bullets in correlative ecological niche modelling: insights from testing among many potential algorithms for niche estimation[END_REF]. In small presence-only datasets, the methodologies used to correct for spatial bias are not as efficient for narrow niche species as for broader niche species. Reducing the extent of distribution modeling of narrow niche species to the boundaries of their environmental limits could prove a good alternative.

Influence of record addition

The chronological addition of new data has a limited impact on certain model outputs as demonstrated by high similarities between the chronological models generated for A. cordatus and C. nutrix. In contrast, chronological models of S. diadema significantly differ between each other. The detailed analysis of data increment proved that the increasing number of presences has no impact on modeling performance, which is not in line with previous works (Stockwell andPeterson 2002, Wisz et al. 2008). However, these results can be altered by our incomplete knowledge of species full distribution due to the limited number of data available and to sampling bias [START_REF] Hernandez | The effect of sample size and species characteristics on performance of different species distribution modeling methods[END_REF], Bean et al. 2012). In models of S. diadema, differences between the chronological models are due to contrasted spatial patterns between datasets (transect versus random patterns).

Historical data and environmental change

Significant environmental shifts were measured for the descriptors analysed between 1955 and 2012 over the Kerguelen Plateau (i.e. mean sea surface temperature and amplitude, mean surface salinity and amplitude). However, for all species, distribution models built for each decade are highly similar between each other. These results confirm that temporal heterogeneities in datasets do not necessarily impact the robustness of models, because species preferences for their environment may be wider than the magnitude of changes in time. Working with both present and historical data to improve the completeness of occurrence records proved reliable when assuming that species niche and distribution have not significantly changed during the studied time period.

Between the five decades, the respective contributions of temperature and salinity to the models did not vary over the range of within-decade variation for B. antarcticus, C. nutrix and S. diadema. Variations between decades are more marked in models produced for A. cordatus. This near-shore species is found in shallow waters of Kerguelen and Heard islands, where environmental descriptors include many no data pixels [START_REF] Guillaumot | Echinoids of the Kerguelen Plateau: occurrence data and environmental setting for past, present, and future species distribution modelling[END_REF].

Consequently, the varying contributions of temperature and salinity to the models of A. cordatus between decades cannot be attributed with certainty to the effect of environmental change but to modeling limitations.

Sea surface temperature and salinity amplitudes have significant contributions to the models, contributing more than averaged parameters (i.e. A. cordatus and B. antarcticus, Fig. 5). This is in line with the results of [START_REF] Bradie | A quantitative synthesis of the importance of variables used in MaxEnt species distribution models[END_REF] who tested for the contribution of several environmental descriptors on a wide panel of taxa. They showed the importance of including seasonal means and extremes in models to further depict species distribution, considering their stronger relationships with species niche width and ecological traits (i.e. growth and survival, see [START_REF] Franklin | Mapping species distributions: spatial inference and prediction[END_REF].

Using time-averaged descriptors over the entire period could have been considered the best approach to produce representative models, independent from short-term environmental variations. Unexpectedly, our results show that for all species, contributions of time-averaged descriptors to the models are much more different than all differences between decadal descriptors (Fig. 5). This shows that using time-averaged descriptors for long time periods does not necessarily improve model reliability in comparison with using descriptors averaged for a shorter time period. This also shows the importance of the descriptor selection in modeling procedures, a critical issue for improving model performance as already stressed in previous studies [START_REF] Bradie | A quantitative synthesis of the importance of variables used in MaxEnt species distribution models[END_REF]. This is particularly relevant for certain regions of the Southern Ocean like the Western Antarctic Peninsula that has experienced among the most significant environmental changes in the world ocean during the last decades (Turner 2015).

Influence of species niche width in modeling performances

Among the four studied species, A. cordatus has the narrowest ecological niche and most restricted distribution in the vicinity of coastal areas of the Kerguelen and Heard archipelagoes. Such limited geographic and environmental distributions compared to the total extent of the studied area implies that similar environmental conditions prevail in geographically close occurrence sites. This induces a strong SAC pattern that explains the difficulties to correct for spatial bias in comparison to other species models. Moreover, the limited environmental variability between coastal sampling sites of the different oceanographic surveys can also explain the absence of data addition effect on modeling performances for A. cordatus.

In contrast, C. nutrix and S. diadema have wider ecological niches than A. cordatus (Fig. 1).

For these two species, record data are more widely distributed and show contrasting sampling patterns (i.e. transect-like versus random patterns) that have been shown to influence modeling performance in S. diadema only (Table 3). This can be explained by the higher number of presence records available for C. nutrix (n=114 and n=98 for C. nutrix and S. diadema respectively) that allowed a more complete survey of C. nutrix distribution. Finally, C. nutrix dataset only presents a quality and number of occurrence records that fulfill all methodological requirements to produce reliable distribution models.

Considering species niche width to cope with spatial and temporal bias in SDM is important, as already shown by Tessarolo et al. (2014) who studied the influence of survey designs on the performance of distribution models for endemic species with narrow ecological niches.

They concluded that survey designs have a low impact on models in comparison with the effect of niche width, data number, and type of modeling technique used. However, they did not generate any analysis of species with broad ecological niche as a comparison. Our results are also in line with other modeling studies in which distribution models of species with broad niche were the least stable [START_REF] Reiss | Species distribution modelling of marine benthos: a North Sea case study[END_REF][START_REF] Qiao | No silver bullets in correlative ecological niche modelling: insights from testing among many potential algorithms for niche estimation[END_REF][START_REF] Guo | Uncertainty in ensemble modelling of largescale species distribution: Effects from species characteristics and model techniques[END_REF][START_REF] Ranc | Performance tradeoffs in targetgroup bias correction for species distribution models[END_REF].

Conclusion

The use of SDMs has gained in importance during the last decades. They can provide complementary information for environmental managers. Modeling results can help interpolate species distribution, identify the potential drivers of species distribution and predict the potential effects of environmental changes on habitat suitability. However, modeling species distribution over vast and remote marine areas like the Southern Ocean using poor and heterogeneous datasets remains challenging and improvement of biological and environmental datasets is still required.

In the present study, we showed that reliable species distribution models can be produced in such areas as long as the number and quality of data allow testing and correcting for the effects of biases. Using historical data requires proper environmental descriptors for modeling the effect of environmental changes on species distribution. Using time-averaged predictors over long time periods can generate unrealistic models.

Model selection is also crucial at this stage and the statistical performance of models is not the only criteria to be considered. Modeling procedures must be chosen with regards to the scientific issues that are addressed. Two procedures (BRT and RF) performed best in our case study, but one of them (BRT) was proved to be more relevant because it better dealt with transferability and data patchiness.

Modeling species distribution in data-poor areas poses the practical problem of the minimum number of presence-only data required to run reliable models, although this is not the only and most critical issue. Beforehand, the number of occurrence records must be high enough for testing model robustness and reliability. In regions with limited access, sampling effort may be heterogeneous, which influences model performance. We showed that sampling bias can be corrected, but the efficiency of correction depends on species niche width, narrow niche species models being more troublesome to correct. In our study, A. cordatus is a species limited to coastal shallow areas, which implies a strong correlation between species occurrence and sampling patterns. Restricting the model to a more reduced area could allow to correct for spatial bias and improve modeling performances.

There is also a crucial need for improving the quality of datasets [START_REF] Kennicutt | Six priorities for Antarctic science[END_REF]) and running more accurate models to better tackle conservation issues [START_REF] Guisan | Predicting species distributions for conservation decisions[END_REF][START_REF] Rodríguez | The application of predictive modelling of species distribution to biodiversity conservation[END_REF]. For the time being, producing uncertainty maps can be an alternative [START_REF] Rocchini | Accounting for uncertainty when mapping species distributions: the need for maps of ignorance[END_REF], Tessarolo et al. 2014) and can provide additional information to environmental managers and stakeholders [START_REF] Addison | Practical solutions for making models indispensable in conservation decision-making[END_REF][START_REF] Guisan | Predicting species distributions for conservation decisions[END_REF].

Model reliability and performance also rely on the interaction between dataset completeness and species intrinsic ecological properties. Hence, we showed that the type and width of ecological niches were important to consider, distribution of narrow niche species being easier to model and less sensitive to incomplete datasets [START_REF] Guo | Uncertainty in ensemble modelling of largescale species distribution: Effects from species characteristics and model techniques[END_REF][START_REF] Ranc | Performance tradeoffs in targetgroup bias correction for species distribution models[END_REF]). However, narrow niches usually imply that species are distributed over small areas for which distribution models will be highly sensitive to extrapolations.

Our protocol shows that reliable SDMs can be produced when enough data are available and dataset bias can be tested and corrected. In the present study, one SDM only (C. nutrix) could be corrected for spatial and temporal heterogeneities to generate reliable distribution predictions. However, our results stress the need to consider methodological issues when modeling species distribution based on poor and spatially biased datasets. They should contribute to bring new insights and enhance modeling performance in future studies. Table 1: List of environmental descriptors selected for SDM. Asterisks (*) indicate that environmental layers are available for the following time periods: , [1955][1956][1957][1958][1959][1960][1961][1962][1963][1964], [1965-1974], [1975-1994], [2005][2006][2007][2008][2009][2010][2011][2012]. Minimum and maximum values are indicated for . Spatial resolution of layers: 10km resolution grid-cell pixels.

Figures legend list

Figure 1 :

 1 Figure 1: (A) Map showing occurrence data of the four studied echinoid species over the

Figure 2 .

 2 Figure 2. Tests and procedures carried out in the present work. Arrows indicate the stepwise

Figure 3 :

 3 Figure 3: (A) Null model #1 and (B) Null model #2 for the different species under study.

Figure 4 :

 4 Figure 4: First row: Distribution models of Ctenocidaris nutrix species with increasing

Figure 5 :

 5 Figure 5: Mean contributions of environmental descriptors to the models with standard

Figure 6 :Figure 1 :

 61 Figure 6: Species distribution model generated for Ctenocidaris nutrix using all presence-

Figure 2 .

 2 Figure 2. Tests and procedures carried out in the present work. Arrows indicate the stepwise

Figure 3 :

 3 Figure 3: (A) Null model #1 and (B) Null model #2 for the different species under study.

Figure 4 :

 4 Figure 4: First row: Distribution models of Ctenocidaris nutrix species with increasing

Figure 5 :

 5 Figure 5: Mean contributions of environmental descriptors to the models with standard

Figure 6 :

 6 Figure 6: Species distribution model generated for Ctenocidaris nutrix using all presence-

Table 2 :

 2 r Pearson correlation of pixel values between the KDE layer and the predicted probability of each species model. Statistic probabilities are all < 0.05.

		Before KDE correction After KDE correction
	Abatus cordatus	0.72	0.44
	Brisaster antarcticus	0.60	-0.17
	Ctenocidaris nutrix	0.80	0.11
	Sterechinus diadema	0.61	0.20

Table 3 :

 3 Influence of data addition and sampling patterns on models for Abatus cordatus, Ctenocidaris nutrix and Sterechinus diadema. Column 1: mean Schoener's D and associated p-value computed between models (100 replicates) produced respectively with {n=54, 76, 95}, {n=46, 54, 106, 114} and {n=54, 66, 98} occurrences randomly sampled from the total dataset. Column 2: mean Schoener's D and associated p-value computed between models (100 replicates) produced with subsets contrasting in data distribution patterns (transect versus random sampling).

	Species	occurrence number	sampling pattern
		Mean Dobs	Mean p-value	Dobs	p-value
	Abatus cordatus	0.981 ± 0.025	<0.05	-
	Ctenocidaris nutrix	0.985 ± 0.020	<0.05	0.941 ± 0.030	0.147
	Sterechinus diadema	0.979 ± 0.031	<0.05	0.842 ± 0.040	0.941

Acknowledgements

This work is a contribution to the IPEV program PROTEKER funded by the French polar institute (IPEV program n°1044) and contribution no. 21 to the vERSO project (www.versoproject.be) funded by the Belgian Science Policy Office (BELSPO, contract n°BR/132/A1/vERSO).

Smith WH and Sandwell DT (1997) Global sea floor topography from satellite altimetry and ship depth soundings. Science. 277: 1956-1962. Stockwell DR and Peterson AT (2002)