High-resolution topographic variables accurately predict the distribution of rare plant species for conservation area selection in a narrow-endemism hotspot in New Caledonia

Species distribution models (SDMs) represent a widely acknowledged tool to identify priority areas on the basis of occurrence data and environmental factors. However, high levels of topographical and climatic micro-variation are a hindrance to reliably modelling the distribution of narrow-endemic species when based on classic occurrence and climate datasets. Here, we used high-resolution environmental variables and occurrence data obtained from dedicated field studies to produce accurate SDMs at a local scale. We modelled the potential current distribution of 23 of the 25 rarest species from Mount Kaala, a hotspot of narrow-endemism in New Caledonia, using occurrence data from two recent sampling campaigns, and eight high-resolution (10 m and 30 m) environmental predictors in a Species Distribution Modelling framework. After a first sampling operation, we surveyed six additional areas containing, overall, 13 of the 20 species modelled at this stage, to validate our projections where the highest species richness levels were predicted. The ability of our method to define conservation areas was largely validated with an average 84% of predicted species found in the validation areas, and additional data collected enabling us to model three more species. We therefore identified the areas of highest conservation value for the whole of Mount Kaala. Our results support the ability of SDMs based on presence-only data such as MaxEnt to predict areas of high conservation value using fine-resolution environmental layers and field-collected occurrence data in the context of small and heterogeneous systems such as tropical islands.


Introduction
More than half of plant extinction events over the last 500 years occurred in oceanic islands, which are home to more than 30% of critically endangered plants (Tershy et al. 2015). In this regard, islands are priority territories for the development of plant conservation programs (McGeoch et al. 2016). In an era of global biodiversity crisis (Ceballos et al. 2015), the world's biodiversity hotspots have been defined according both to the level of richness and endemism of their flora and the threats they are facing, nearly half of them being island territories (Myers et al. 2000). As stressed by Cañadas et al. (2014), some hotspots are composed with micro-or nano-hotspots due to high levels of narrow-endemism and should be the focus of conservation actions. The presence of such narrow-hotspots was highlighted in New Caledonia by Wulff et al. (2013) and was recently identified as in need of urgent conservation planning (Gâteblé et al. 2018). In these areas of high conservation value, management planning often suffers from a lack of high-resolution background data on precipitation and temperature regimes, resulting in a poor understanding of biogeographical patterns and species ecology (Heywood 2011). Knowledge about the ecology of rare species is even more limited in the case of narrow endemic species (NES) occurring in inaccessible locations (Kier et al. 2009). Pressures from global changes are particularly concerning in islands (Weigelt et al. 2016) and NES are particularly threatened (Wulff et al. 2013;Cañadas et al. 2014). There is a need to characterise the ecological niche of rare insular plants and their potential response to environmental changes.
Species distribution models (SDMs) are now widely used to produce transferable projections of the potential distribution of known species, which helps identify important areas for conservation among other goals (Guisan et al. 2013;Leroy et al. 2014). In this regard, Heywood (2011) proposed the use of this tool to guide the protection of island flora in data-poor territories. However, modelling the distribution of rare plant species is challenging in several aspects. First, the resolution of the most commonly used bioclimatic data (Booth et al. 2014) such as the Worldclim Database (Fick and Hijmans 2017) or Chelsa climate data (Karger et al. 2017) are too coarse to properly represent highly heterogeneous habitats. The most up to date WorldClim dataset is available at a 1 km resolution (World-Clim2-Fick and Hijmans 2017), while it is known that rare plant distributions can be locally driven by topo-climatic variations at a 10-100 m scale (Franklin et al. 2013). The resolution of these global datasets is therefore insufficient to properly model species distributions and guide decision-making tools at the local scale (Kier et al. 2009). Secondly, occurrence data of rare or narrow-ranging plant species are difficult to obtain from existing databases. Generalist datasets and herbarium collections are the most common sources of such data (Loiselle et al. 2003;Gallardo and Aldridge 2013). However, they often lack occurrence data on rare species, and herbarium collections may be geographically biased (Ter Steege et al. 2011). Ad hoc collection of input presence data may thus form a better basis for high resolution modelling of rare plant species in narrow ranges. Finally, the selection of absence data is a key point in SDM projections, hence calling for a precise description of the survey effort underlying each occurrence dataset (Phillips 2008). Poorly defined absence data may indeed result in overpredictions in surveyed areas and produce "precise answers to the wrong questions" (Yackulic et al. 2013).
While SDM studies have advanced toward more efficient algorithms at large scales, there have been relatively few applications to the projection of NES distributions. Recent studies are answering some of these challenges by integrating microvariation predictors into their projections (Tomlinson et al. 2019). However, high-resolution background data 1 3 are not always available. When such background data are lacking, a combination of dedicated field collection of presence data and post hoc field validation of projections should produce reliable guidelines for in situ conservation purposes.
The New Caledonia archipelago is one of the world's biodiversity hotspots of the highest conservation value (Myers et al. 2000;Wulff et al. 2013). Almost 1% of the world's flora is endemic to this southwest Pacific archipelago (Munzinger et al. 2020). In this region, about 20% of the indigenous flora is considered as narrow-endemic, i.e., restricted to one, two or three locations (Wulff et al. 2013). This high endemism is the consequence of a unique geological history, resulting in a high topographic and climatic heterogeneity. Ultramafic soils are very frequent in the south of Grande Terre (the main island) and sparsely distributed northwards along the west coast in isolated and steep mountains (Bonvallot et al. 2012). In terms of plant ecology, these mountains are considered as ultramafic islands within a non-ultramafic matrix (Isnard et al. 2016). The original composition of ultramafic soil, particularly the abundant presence of nickel ore, offers substantial economic value to the mining industry that is widely developed in New Caledonia.
Mining represents a major threat to biodiversity through vegetation removal, topographic alteration, and persisting pollution (Lefcort et al. 2010;Palmer et al. 2010;Duran et al. 2013). In New Caledonia, the number of mining tenements issued rose by about a hundred between 1997 and 2017, resulting in more than 45% of all ultramafic areas currently being under mining concessions . Furthermore, no protected areas have been implemented on ultramafic substrates in the northern half of the main island , which emphasises the need for urgent conservation programs to identify and protect its endemic-rich flora. In 2019, the conservation status of one third of the New Caledonian flora was assessed by the local Red List Authority (RLA-NC). Out of a thousand species, 43% were classified as threatened, mainly by open-cast mining, bushfires or invasive alien species, and 11% were considered "data-deficient" (RLA-NC 2019). This lack of information needs to be addressed urgently for the design of efficient management strategies (Collen et al. 2008;IUCN 2012;IUCN/SSC 2013).
In this study, we coupled empirical field data collection and species distribution modelling to produce high-resolution "conservation priority" maps and promote the conservation of narrow-endemic taxa in one of New-Caledonia's hotspots of narrow-endemism threatened by mining activity. We predicted the distribution of 23 species among the rarest taxa of Mt. Kaala, using presence data from two recent sampling campaigns, and eight available habitat variables (10 m and 30 m resolution). To evaluate the ability of our projections to define areas for conservation, we carried out six additional surveys in areas where the highest species richness was predicted. Our results are discussed in light of the local and regional strategies for plant conservation. The study provides an integrative framework that should be considered by environment managers in hotspots of narrow-endemism.

Study site
The New Caledonia archipelago is located in the Coral Sea, 1500 km east of Australia and 2000 km north-northwest of New Zealand. It encompasses the main island of Grande Terre, the Loyalty Islands, and a number of smaller islands and islets. Grande Terre is 400 km long and 50-70 km wide.

3
Our study site is defined by the geological entity known as Mount Kaala, located in northwest Grande Terre (164°23′E, 20°36′S), with an altitude of 1070 m (Fig. 1). The site is dominated by ultramafic substrates with a serpentinite outcrop at the bottom, and is isolated in a large volcano-sedimentary coastal plain. The vegetation of Mount Kaala is a mosaic of varyingly degraded shrublands called "maquis miniers" (ca. 61% of the study area based on Barrière et al. 2007), dense humid forest (ca. 19%) and degraded secondary herbaceous patches (ca. 16%) on the slopes. The summit plateau, exploited for ore extraction, is composed of bare soils (ca. 4% of the study area) and secondary herbaceous patches.

Studied species
The species list we used in this study was determined in a previous study (Table 1, Lannuzel, unpublished.), and is based on a comprehensive screening of all herbarium data  1 3 available worldwide (mainly Paris-P, Nouméa-NOU, Missouri Botanical Garden-MO) for the Mount Kaala area. We filtered the resulting species list using the IUCN criteria and methodology (IUCN 2012). This process highlighted 23 taxa assessed as either Critically Endangered or Endangered in 2017 by the local Red List Authority (RLA-NC). We also kept two non-evaluated taxa because we assumed their extinction risk was very high. We identified four specimens collected in the early stages of the study as new species and the publication of their formal description is pending. We assume the species list we used is up to date for this area. However, given that new species are continuously described in New Caledonia (Gâteblé et al. 2018), our study provides recommendations for conservation on the basis of current knowledge.

Environmental variables
We used six environmental predictors in this study, all based on topography and remote sensing. Because of the high resolution needed, the use of global bioclimatic data was not appropriate. Recently, Tomlinson et al. (2019) showed the relevance of using edaphic proxies versus downscaled climatic variables in the case of narrow-range endemic species. As a consequence of the lack of reliable edaphic predictors in New Caledonia, we used proxies representing the local topo-climatic variations (Pouteau et al. 2012), as well as a geological layer, as proxies for edaphic conditions. We computed four abiotic layers all derived from a Digital Elevation Model (DEM) freely available at a 10 m resolution (accuracy ± 2 m) (Gouvernement de la Nouvelle-Calédonie 2017), following the TopoTo Raster methods (Juffroy 2012). The first proxy is the DEM per se, because elevation is linked to temperature at the local scale (Bonvallot et al. 2012). Slope steepness (hereafter referred to as "slope") affects the velocity of both surface and subsurface flow, and hence soil water content, erosion potential and soil formation (Pouteau et al. 2012;Wilson and Gallant 2000). Slope exposure (hereafter "exposure") was used as a proxy of both insolation and windwardness, which can influence the drying rate of vegetation (Wilson and Gallant 2000). We used the topographic wetness index (TWI) as a proxy for soil moisture, and therefore water availability for plants (Western et al. 1999;Radula et al. 2018). In New Caledonia, TWI was recently proven to affect species assemblages at local scales (Blanchard et al. 2019). Slope and exposure were calculated using the QGIS 2.18.7 "terrain analysis" tool. TWI was computed with the SAGA GIS dedicated algorithm. It is based on slope and catchment area, the latter also computed using SAGA GIS. We also used a geological polygon layer (1/50,000 scale) with substrate type categories (Maurizot et al. 2005) rasterized at a 10 m resolution. The only biotic predictor used was the Normalised Difference Vegetation Index (NDVI), a continuous index representing vegetation productivity (Pettorelli et al. 2005). Recent studies showed a correlation between tree species richness and single NDVI scene taken in the wet season in New Caledonia (Pouteau et al. 2018). We therefore interpreted it as a proxy for habitat suitability because many taxa are restricted to specific vegetation types (Isnard et al. 2016). NDVI was calculated from 30 m-resolution Landsat 8 images from 24th March 2017 (USGS 2017), corresponding to one of the only scenes without clouds in the study area during the wet season. We accounted for collinearity between environmental predictors using Pearson correlation coefficients (

Sampling effort
We acquired a first dataset during the five field studies in 2017 and 2018. We sought the most cost-effective method for data collection considering the steepness of the mountain and limited accessibility. Field surveys were designed to equally sample all altitudinal (DEM-100 m classes) and vegetation classes (NDVI-0.1 index value classes). We used strip transects, which is recommended in the case of highly clustered multi-species assemblages (Melville and Welsh 2001;Ogutu et al. 2006). Because strip transect accuracy depends on width, we assumed that detection was possible within a 5 m wide strip on either side of the observer. We chose a 5 m limit as the maximum distance at which we confidently detected a shrub species in dense shrublands, and small herbaceous species in herbaceous maquis. Each GPS mark was then located at the centre of a 10 m wide square-cell within which all species of interest were recorded.
A total of twenty-one strip transects were carried out on Mount Kaala (Fig. 1a). Transect locations were recorded using a GPS device (Garmin GPSMap64; spatial accuracy ± 3 m). All data were transmitted to the local Red List Authority and are available on the www.endem ia.nc website.

Modelling
To broaden the appeal of our approach to the widest user community, we chose a modelling algorithm that is reproducible without needing high computational power and which can be undertaken by any stakeholder with moderate background knowledge of SDMs. We used MaxEnt (v.3.3.3; Phillips et al. 2006), a user-friendly and well documented software. This algorithm is widely acknowledged to be efficient for rare species and performs well with small datasets (Phillips et al. 2004(Phillips et al. , 2006Pearson et al. 2007;Wisz et al. 2008). The convergence threshold was set at 0.00001 with a maximum number of iterations of 500. We chose the logistic output format with suitability values ranging from 0 (unsuitable) to 1 (optimal; Phillips and Dudík 2008). We ran five cross-validation replicates for each species. When the number of occurrences was lower than five, the software automatically reduced this number. For each replicate, it randomly splits the dataset into two parts: 80% for training, 20% for testing. The remaining parameters were kept as default settings.

Sampling bias correction
The sampling methodology inevitably produced spatial and ecological auto-correlation of occurrences, which could represent a major bias (Kramer-Schadt et al. 2013). One expected consequence of this bias is to make predictions biased toward over-sampled areas (Phillips et al. 2009). As stated by Phillips et al. (2009), presence-only modelling is fairly robust as long as background points have the same sample selection bias as the presence points (Yackulic et al. 2013). However, this background selection method was hard to implement where the sampling effort is not precisely reported (Phillips et al. 2009;Kramer-Schadt et al. 2013;Fourcade et al. 2014). In our case, the dedicated sampling inherently provided sampled and unsampled areas. Moreover, this method enabled the use of all presence locations, as opposed to split or spatial filtering methods (Fourcade et al. 2014). It is therefore preferable for rare species with fewer occurrence 1 3 data. We added a 5 m buffer to our GPS tracks to generate 10-m wide polygons corresponding to our sampling design. The resulting polygon layer was rasterised to obtain a 10-m resolution raster layer with a value of 1 in sampled pixels, and 0 in unsampled pixels. This layer was implemented in Maxent as a "Bias file". The software then selected background points within the surveyed area, while results were projected on the overall study area.

Model performance
Selection of appropriate validation metrics such as AUC and TSS is a matter of debate in presence-pseudo-absence modelling, especially for rare species when prevalence is low (Leroy et al. 2018). To reduce misleading results, Leroy et al. (2018) advocated for the use of similarity tests in such cases. To evaluate model performance independently from prevalence, we used the continuous Boyce index (B cont(0,1) ; Hirzel et al. 2006) which avoids arbitrary selection of class cutting thresholds by using a sliding window, and consequently might be more efficient for small sample tests. This indicator was computed for each replicate using the "ecospat" package (Di Cola et al. 2017) in R (R Core Team 2018). The Boyce index varies from − 1 to 1. A value of 1 corresponds to a perfect match between observed and predicted occurrences, while values under 0 are no better than random. Every acceptable prediction (B cont(0,1) > 0.5) layer was transformed into a binary presence/absence map using the maximum sensitivity/specificity threshold (Jiménez-Valverde and Lobo 2007; Cao et al. 2013). Considering the debate on validation metrics for SDM, we also computed the true skill statistics (TSS), because many authors consider it as prevalence-independent, but did not refer to it for further analyses. Model performance is considered as poor if TSS value is lower than 0.4, moderate between 0.4 and 0.6, and good above 0.6 (Beauregard and de Blois 2014). We used the maximum sensitivity/specificity threshold to calculate TSS (Allouche et al. 2006;Williams et al. 2009). We also provide AUC (that is automatically calculated by MaxEnt) in Appendix Table 5. We present average predicted values for each taxon.

Field validation
We produced a general map by stacking every specific binary map (Fig. 2, and Appendix Fig. 4). We obtained a raster layer giving the expected number of species present within each pixel. Six areas labelled from "a" to "f" (Fig. 2) were chosen for survey within the areas where pixel values fell within the 4th upper quartile (Rhoden et al. 2017). These areas were surveyed, during a sampling trip on 6th and 7th February 2018. We then computed a post-hoc test, i.e. the Positive Predictive Power index (PPP; Fois et al. 2015) to assess the efficiency of our SDMs. PPP represents the ratio between true positives and the sum of predicted positives. Here, we considered a true positive when a predicted species was actually observed in the area considered and a false positive when a species was predicted but not observed.
All data collected during this second field campaign were added to the total dataset to produce the final prediction maps. B cont(0,1) were re-calculated accordingly.

Plant conservation priority scoring
A lot of literature is available to determine important areas for conservation (Blasi et al. 2011;Balaguru et al. 2006;Darbyshire et al. 2017). We adapted the different criteria to match available data and specific local characteristics, such as narrow-endemism. Darbyshire et al. (2017) identified three important criteria while considering the establishment of important plant areas (IPAs): threatened species, botanical richness, and threatened habitat. The threatened species criterion is relatively straightforward if IUCN assessments are locally available, which was the case here. Other criteria are harder to implement in a data-poor context. When no exhaustive inventories are available, botanical richness can be assessed by the density of important/valuable species (Darbyshire et al. 2017). We used the predicted presence of species in our list, weighted by the IUCN status, to inform both these criteria. The IUCN (2016) recommended considering the irreplaceability of a species or population. In New Caledonia, NES that are restricted to one locality are irreplaceable, because, if they disappear from Mount Kaala in our example, they go extinct. We thus considered every predicted presence of such species as of irreplaceable value in our index. Threatened habitats are not known clearly in New Caledonia. In the absence of such a reference source, scores of the "threatened habitat" criterion were computed using two components. The first was the hierarchy of the given vegetation type within the post-fire ecological succession based on McCoy et al. (1999), from bare soil (lower level) to dense humid forest (upper level). The second was an area index representing the relative scarcity of each vegetation type on the mountain. The vegetation map used at this step was Fig. 2 Stack map of every acceptable binary occurrence prediction of 20 rare plant taxa on Mt Kaala (New Caledonia) obtained during the 1st field study, which we used to define validation areas ("a"-"f"). The location of transect surveys are also indicated. On the right, examples of validation areas with transect surveys and studied species' occurrences produced by Barrière et al. (2007). This results in the following formula where i is the vegetation type: We then proposed the list of criteria shown in Table 2. The final layer was calculated to determine an index of conservation priority I, computed with the following formula: where A, B, and C are the priority conservation criteria mentioned in Table 2 and j is the taxon considered. Criteria A and B were calculated for each taxon, while criterion C was calculated once for the whole study area. The I index therefore varies between 0 and 3, where 0 corresponds to an area with the lowest conservation priority, and 3 represents the highest level of priority.

Pre-validation modelling
At the end of the first data collection phase, the mean number of occurrences was 33 ± 11 SE per species. Three species were detected only once, and a maximum of 259 occurrences was reached for Zieria chevalieri. All elevation and NDVI classes were included in the surveys. Relative to the whole massif area, transects represented about 0.5% of Mount Kaala. A total of 1610, pixels (10 × 10 m) were sampled, for a total of 800 "presence" pixels. For elevation classes, 0.68 ± 0.19% of areas were surveyed with an over-representation of the highest classes (700 m and higher). For NDVI, 0.36 ± 0.07% of areas were surveyed with a slight over-representation of the 0.1-0.2 classes. An under-representation on higher classes (0.4 and more) corresponded to dense humid forests. Species with one or two occurrences (four species) were discarded for the modelling procedure. Distributions of twenty-one species were projected, among which 17 had a B cont(0,1) ranging from 0.58 to 1 and were thus considered reliable.The four models with a B cont(0,1) < 0.5 were not considered for the validation step.
Niche overlap values ranged from 0 to 8 species predicted per pixel.

Field validation
We conducted surveys in six areas projected to host more than five species. We selected two areas at low altitude, one in the middle altitudinal range, and three at high altitude. Mean PPP was 0.80 ± 0.02 SE, with no PPP under 0.7 (Table 3). Surveyed areas contained 15 of 17 predicted taxa, and no unpredicted species were observed.
Surveys represented 5% of each selected area on average, and 460 additional surveyed pixels (46,000 m 2 -28% of the initial dataset). Additional occurrence data were then added to the occurrence datasets. This allowed modelling of the distribution of one more species, and improved the validation metric of four others. Final results included 22 prediction maps, with B cont(0,1) values ranging from 0.51 to 1. Each predictor's contribution varied among species. DEM was the predominant predictor with 51% contribution on average, but some species' projections were mostly driven by NDVI (e.g. Gynochthodes truncata), exposure (e.g. Marsdenia kaalaensis) or slope (e.g. Myrsine novocaledoniae subsp kaalaensis). Soil contribution represented 7% on average, but emerged as a strong predictor for some species (e.g. Olax hypoleuca-28%, Gynochthodes truncata-26%). TWI was the poorest predictor with an average 2% contribution, and a maximum contribution of 6%. Detailed statistics and results are given in Appendix Table 5. Detailed response curves are also given in Appendix Fig. 5. For most species, the marginal response curves of each predictor showed no significant differences with the response obtained with models created using only the corresponding variable. In a few specific cases only (C. kaalaensis, C. velutinum, P. deciduiramus, Psychotria sp. nov., T. minutiflora), exposition and slope show different curve shapes. The final specific maps, once stacked, showed a spatial niche overlap (spatial congruence) of 71% with at least two species predicted per pixel. The number of predicted species per pixel ranged from 0 to 12 (mean = 2.77). The 4th upper quartile of the total number of species per pixel (i.e., more than 9 species predicted) represented 0.8% of the study area.

Priority areas for rare plant conservation
The I index varied within the study area and was rather low on average (median I = 0.60). Pixels with I > 1.5 represented 6.3% of the total study area and were mostly located on the slopes, at intermediate altitudinal range inside and outside active mining concessions (Fig. 3). On the western slopes, it represented some areas at the top end of the three main valleys, where vegetation is less impacted. On the southern and eastern slopes, valleys are less deep, and areas of high I index were revealed in every little thalweg protected from mining erosion and bushfires.

Discussion
In a context of highly heterogeneous habitats, our results showed that the occurrence of narrow-endemic species can be accurately predicted using topographic, geological and vegetation predictors. The Boyce validation metric seemed efficient in assessing the quality of predictions, even with only four presence data (P. pinifolia, B cont(0,1) = 0.97, verified during validation). Our field validation surveys confirmed 80% of projected occurrences, which supports the reliability of our models (Jiménez-Valverde 2012; Somodi et al. 2017;Leroy et al. 2018). Hence, we were able to use the projected maps to develop an evidencebased spatially explicit index for local environment managers. Data sampling was designed to offset the lack of presence data in international databases. However, the production of presence-only data sets depends on three parameters: (i) sampling probability, (ii) occurrence probability, and (iii) detection probability. Together, they produce a bias often occurring in presence-only modelling studies (Yackulic et al. 2013). Our bias correction method maximised sampling-probability using transect tracks as bias files, without removing any occurrence data, which is a key point when modelling Fig. 3 Distribution of priority areas for conservation (computed with the I index) at Mt Kaala (New Caledonia). High I values represent areas with a high number of rare plant species, weighted according to IUCN status, narrow-endemism status and habitat importance. Active mining concessions are circled in black the distribution of rare species (Fourcade et al. 2014). This method, here applied to transect sampling, also addressed the detection probability issue by using the 5 m-buffer added to the tracks. Finally, we measured a first likelihood of occurrence, and used it for further field surveys which enabled us to add new occurrence data for the final models. This was only possible because we started from dedicated surveys, allowing a high confidence in the definition of sampled/unsampled pixel maps. Controlling for these three parameters is obviously more difficult on SDMs at large scales when the sampling effort underlying occurrence data is often unknown (Phillips et al. 2009;Kramer-Schadt et al. 2013;Fourcade et al. 2014).
The choice of environmental predictors is a challenging issue in fine-scale SDMs (Nezer et al. 2016). Hijmans et al. (2005) showed the inaccuracy of global climate products such as Wordclim in the case of remote oceanic islands (including New Caledonia), because of the scarcity of local meteorological station networks (Turner et al. 2003;Pouteau et al. 2015). These products are thus unlikely to produce reliable predictions of very rare plants' suitable habitat, often influenced by 10-100 m scale variations (Franklin et al. 2013). The proxies we used here, mostly derived from the DEM layer, are the most relevant and locally available at high resolution for plant ecology (Pouteau et al. 2012;Wilson et al. 2013). However, three predictors were computed from the DEM which increases the risk of multicollinearity issues. We tested this issue using the Spearman correlation index, but also analysed the response curves of each predictor and species (Appendix Fig. 5) to confirm the absence of multicollinearity. For each species and predictors, two curves are computed. The first one represents the response of the model to the predictor while all other predictors are set at their mean value in the sample. The second represents the response to the predictor in a model built only with this predictor. The comparison between both curves shows whether a predictor is independent (both curves are similar) or correlated (both curves show different shapes- Elith et al. 2011;Phillips 2017). Almost all results show independence between predictors, except for some species, for which exposition and slope had slightly different shapes. In these cases, the influence of a predictor on a species' distribution must be read on the response curve of the predictor alone as it is more accurate if predictors are correlated, even slightly. The level of collinearity between variables was relatively low (maximum Spearman index = − 0.48 and 0.39) and was low across species. Therefore, we believe potential collinearity between our environmental predictors may occur at the margin of species distributions, but may not alter our conclusions.
DEM was the most influential predictor overall. However, the contribution of each predictor varied with species. Gynochthodes truncata for example, one of the most commonly recorded target species, was found all along the altitudinal gradient, and had a distribution influenced mostly by NDVI, while Phyllanthus pterocladus was equally influenced by DEM (34%) and NDVI (29%). The potential distribution of most NES restricted to Mount Kaala was influenced primarily by DEM and high altitude. This could be explained by a stronger island-like isolation of high-altitude areas on ultramafic mountains like Mount Kaala in New Caledonia. Further studies on other ultramafic mountains would give more insight on narrow-endemism patterns, and contribute to understanding the drivers of such restricted distributions in the New Caledonian flora.
DEM per se and exposure are considered as climate proxies for temperature, and windwardness/insolation respectively. Slope, along with TWI and soil, influence soil formation and water content. In New Caledonia, L'Huillier et al. (2010) showed a causal relationship between elevation, slope, and soil characteristics on ultramafic substrates. Brown hypermagnesian soils such as colluvial formations (Isnard et al. 2016) cover lower altitudes, whereas ferralsols with high rock content influenced by slope and proximity to the ridges (Isnard et al. 2016) occur at middle and higher altitudes. The combination of soil substrate, elevation, and slope can be considered as reliable proxies of soil characteristic diversity. Our models were based on edaphic and climate proxies, which are known to be more relevant than strictly climatic ones for distinguishing subtle variations in plant distributions (Illan et al. 2010;Beauregard and de Blois 2014;Meineri and Hylander 2017).
Because of the uniqueness of New Caledonian flora, with high narrow-endemism rates and topo-climatic heterogeneity (Jaffré 1993;Isnard et al. 2016), we adapted already-existing criteria for the selection of conservation areas (Balaguru et al. 2006;Blasi et al. 2011;IUCN 2016;Darbyshire et al. 2017). The criteria we used considered both habitat suitability and vegetation conservation value. Narrow-endemism represents a major conservation stake on mined ultramafic mountains (Wulff et al. 2013).Considering the potential distribution of NES in our irreplaceability index highlighted new areas that would not have emerged from vegetation and species density-based indices. Our scoring therefore accounts for both the local flora specificity and the lack of comprehensive survey data. In the same way, since the country still has no red list of ecosystems, we evaluated the relative rarity of each vegetation type and its position in the vegetation succession in order to assess conservation values. We believe that this index constitutes a strong basis for further implementation of conservation programs in Mount Kaala. Our framework could even be replicated in similar ultramafic mountain contexts, provided that the reliability of irreplaceability indices could be controlled. One may be surprised by the relatively low I index value of low altitude areas, despite the proven high presence of IUCN threatened species. This is likely to be a result of the absence of NES at low altitude, and the relative rarity of forests on these areas mainly covered by shrublands on brown hypermagnesian soils. Consequently, the lower value attributed to early-stage vegetation compared to nearly climax forests tends to decrease the global conservation value of low altitude slopes. However, recent works emphasised the sensitivity of these rich vegetation types that are highly threatened by bushfires . Further developments of vegetation conservation value indices may give more weight to these habitats harbouring high levels of IUCN threatened taxa. Here, we consider the higher I index values in higher elevations relevant with regard to the threats due to mining activity. Ibanez et al. (2019), following previous works Wulff et al. 2013), also stressed the need to consider the high narrow-endemic species diversity on the island. The local network of protected areas, composed of few relatively large areas , does not sit well with the distribution of highly restricted taxa. A network of micro-scale conservation areas, such as the ones identified in this study, could enhance the conservation status of flora across the territory overall. Our results will be transmitted to mining operators on Mount Kaala to promote the protection of areas of high conservation value, and to local environmental authorities as an operational tool for conservation planning.
The high PPP values in every checked area corroborated the statistical results and enabled us to find additional presence data points, consolidating the dataset for final analyses and unravelling additional occurrences of rare species. MaxEnt-guided surveys confirmed the relevance of our models for onsite conservation planning. However, the extent of these surveys was restricted due to time/funding limitations, and we focused on areas with the highest density of rare species. In so doing, we validated the areas of greater conservation interest on the basis of NES diversity with a cost-effective approach. This method cannot therefore be considered as a statistical validation of models stricto sensu (Rhoden et al. 2017;Hertzog et al. 2014), as we did not validate the predicted areas of NES absence on Mount Kaala. A further improvement could be a validation sampling protocol in which every quartile, including predicted absences, of projected species density is equally surveyed (Rhoden et al. 2017;Hertzog et al. 2014). On the other hand, a confirmation of "low conservation" value areas may legitimate further degradations and may not appear to be an efficient resource investment to funding institutions. We validated entire areas containing several species, but didn't check the pixel per pixel validity of the model. Such validation would certainly lead to poorer PPP values. Another solution would have been the use of coarser resolution (e.g. 30 m). However, such resolution would have required additional sampling efforts to balance the lower detection probability if we aimed at being exhaustive in bigger pixels, significantly reducing the cost-effectiveness of our approach. We thus consider our predictions as a fair trade-off between limitations and conservation objectives, i.e. spatial prioritisation rather than species-specific distribution modelling. This prioritisation relied on fewer environmental predictors than in other areas with more fine-grained data available and must be used with caution. The partial field validation process appeared here as an acceptable counterpart to the lack of ecological data.
Island floras are characterised by high levels of endemism and sensitivity often associated with greater knowledge gaps (Caujapé-Castells et al. 2010), resulting in a lack of reliable data for local-scale biodiversity management (Cayuela et al. 2009;Collen et al. 2008). Plus, in rugged terrain areas (Whittaker and Fernandez-Palacios 2007), plant distributions are driven by very fine scale topo-climatic variations. Hence, biodiversity conservation needs to be assessed at a scale that matches these variations. Tomlinson et al. (2019) showed the relevance of coupled edaphic and topographic predictors in modelling shortrange endemic species distributions through high resolution SDMs. In their case, the effect of topography was supposed to come from the correlation between soil characteristics and elevation, with certain soil types occurring only at high altitude. As shown by L'Huillier et al. (2010), New Caledonian ultramafic edaphic diversity is partly linked to elevation, slope, and position on the slope as well, justifying our choice of predictors. A detailed mapping of fine soil characteristics would be of great help to improve our understanding of the discrete drivers of plant species' distribution in New Caledonia, but the cost-benefit ratio of such data acquisition would need to be assessed. Considering the local context, and the urgent need for conservation planning, we assumed that the risk of using limited environmental proxies was lower than the benefits from the acquisition of such predictions of conservation value. Our maps may be used as guidance for the determination of new occurrences for the considered rare species, and our index as a spatial guide for the prioritisation of management actions.
The global conservation effort suffers from a lack of funding, as has been highlighted for tropical islands (Caujapé-Castells et al. 2010). One aim here was to propose a reproducible framework while trying to use the best practices in biodiversity assessment using SDMs (Araujo et al. 2019). We used only freely available data covering all of New Caledonia, and the computational power required was that of a common office computer. Hence, we believe that our study addresses the difficulties encountered by biodiversity managers in areas with high topo-climatic diversity and low data availability such as tropical islands (Koch et al. 2017;El-Gabbas et al. 2020), and would be easily transferable to decisionmakers (Guisan et al. 2013). In addition, both dedicated surveys and field validation of our results gave weight to our model projections despite the constraints encountered. Both the distribution maps and the spatially explicit index we produced represent an opportunity for local stakeholders to develop their conservation plans. Our study promotes further collaboration efforts between researchers, managers and landowners to incorporate our results into existing management strategies and reproduce similar assessments in other rich and understudied areas. Finally, we encourage the use of user-friendly modelling tools in urgent conservation contexts to stimulate further initiatives and help rare plants step into the limelight.

3
Acknowledgements We thank Boris Leroy for useful discussions and wise advice. The authors are grateful to every consulted herbarium for giving access to their data. Roy Benyon is also acknowledged for his precious advice on the English translation. The project was funded by the Société des Mines de la Tontouta, which operates a mine on the top of the mountain, in order to improve conservation efforts within their operational area, and in the surroundings. We thank the Institut Agronomique néo-Calédonien (IAC) for funding JB, MT and BF, and providing field logistical support. We thank the Société Le Nickel (SLN) for giving access to the southern parts of the mountain, and the Kaala-Gomen municipal authorities for giving