Ensemble ecological niche modelling
Presence-pseudoabsence ensemble ecological niche models were generated for each species based on the current climatic data using the BIOMOD framework (Thuiller et al., 2009), implemented in the biomod2package (Thuiller et al., 2021) in R. The ensemble included four algorithms known to be robust and perform well across a range of distribution scales (Meller et al., 2014): multivariate adaptive regression splines (MARS), artificial neural networks (ANN), random forests (RF), and maximum entropy (MAXENT, implemented usingmaxnet ; Phillips et al., 2017). This diversity in computational models provides low inter-correlation in model components, which is ideal for higher predictive performance in ensemble ENMs (Elith, 2019; Valavi et al., 2022). The models were calibrated according to the default parameters and parameterization processes in biomod2 .
To validate model results we used five-fold cross-validation, dividing the occurrence data into five subsets, and testing each subset against a model calibrated on the remaining four subsets. Hold-out validation, where a small subset of the data is used only for evaluation, while the remaining data are used for calibration-validation, was unsuitable due to the small sample size of occurrences in several species. Model performance was evaluated using the area under the receiver operating characteristic (ROC) curve (AUC) and the true skill statistic (TSS). Variable contributions were calculated by performing five permutations for each algorithm, each involving removal of the focal variable from the model to calculate the difference in performance. We report the average permutation importance.
The individual models with a final TSS value > 0.7 were then combined into ensemble consensus models for each species which we used with future climate data to predict the future climatically suitable areas. Predicted suitability was reclassified into binary suitable-unsuitable predictions using the threshold at which TSS was maximised (values above this threshold represent climatic suitability). These binary maps were used to characterise changes in climatically suitable areas from current to future climatic scenarios using four metrics: the relative change in climatically suitable areas (% change), the percentage of the current area retained into the future, and the distance and azimuth (angle of direction) between centroids of current and future suitable areas. Change and retention were calculated respectively as:
\begin{equation} Change\ =\ \frac{Current\ suitable\ area\ -\ Future\ suitable\ area}{\text{Current\ suitable\ area}}\ *\ 100\nonumber \\ \end{equation}\begin{equation} Retention\ =\ \frac{Current\ suitable\ area\ \land\ Future\ suitable\ area}{\text{Current\ suitable\ area}}\ *\ 100\nonumber \\ \end{equation}
To estimate distance and azimuth we first defined distinct fragments within suitable areas. Fragments represent connected suitable cells from the binary maps, estimated by first applying morphological dilation using a buffer of radius 2.5 km, followed by a morphological erosion using a negative buffer of the same radius. As a result of this process, suitable areas within 2.5 km of each other were connected to reflect the assumption that bats can move easily within this distance. The resulting number of fragments varied among species depending on the configuration of the suitable area. For all fragments we then identified centroids and the shortest straight paths between each centroid in the current scenario and its nearest neighbour in the future scenario, assuming zero movement costs. The length and azimuth of that path was used to estimate distances and angles of direction. For each species we calculated a mean distance and mean azimuth. Azimuths were classified into cardinal and intercardinal directions based on angle, assuming North to be 0 and 360. All centroid-based calculations were performed using thegeosphere package (Hijmans, 2022) in R.
Summary models of climatically suitable areas were generated for all species together by creating a map of species richness for each time period. All binary models for each time period were summed together into a model that ranged from 0 to 110, representing the number of species for which a cell is climatically suitable according to their binary models. This model was then rescaled to a 0 to 1 scale for consistent comparison and converted into a binary map using a threshold of 0.3, such that in the final binary models for current and future, the positive cells represented suitability hotspots - areas that were climatically suitable for ≥ 30% (more than 33) of the study species. Certainty in projected models was calculated as a continuous model by first averaging individual replicates within a species to assess model agreement, and then averaging these species models into a final model ranging from 0 to 1, which can be interpreted as high certainty of climatic unsuitability through uncertainty to high certainty of climatic suitability.