2.1 Species selection

In order to decide which species were going to be evaluated, we firstly needed to have a list of the species that were distributed in Mexican islands. As a first approach, the Mexican government provided us a list of 259 terrestrial vertebrate species (^{CONABIO, 2016}). Then, we made an additional review of species in reports, theses, books, and governmental programs; and we increased the species list up to more than double the number of species that was firstly given to us by the Comisión Nacional para el Conocimiento y Uso de la Biodiversidad (CONABIO) (Table I, Appendix I). Finally, the geographic distribution of the listed species (Appendix I) was double-checked with polygons of known distributions from the red list of the International Union for Conservation of Nature (^{IUCN, 2017}) and with taxonomic authority catalogs, in order to confirm their current taxonomic status and distribution on Mexican islands (^{CONABIO, 2011}, ²⁰¹³; ^{Navarro and Gordillo, 2006}; ^{Ramírez, 1999}).

Once we prepared the species list, we looked up for the species geographic records at CONABIO (2017) and the Global Biodiversity Information Facility (^{GBIF, 2017}). The first criterion for choosing species to be modeled was the number of unique geographic records. A minimum amount of records to create reliable models is still under discussion, but some experts suggest at least five (^{Pearson et al., 2007}). We decided to use a minimum threshold of 10 unique records in order to have enough data for creating and validating the model. The second criterion was to reduce uncertainty and error sources coming from scarce information and conceptual complexity, and although marine mammals and birds were listed in the occurrence data set, we decided not to consider them in the ecological niche modeling. As far as we know, there are no future marine layers for Mexico and as a result it is not possible to create ecological niche models (ENM) in the ocean. In terms of birds, most of islands’ species with over 10 unique records are migratory and consequently they need to be modeled in a particular way.

The groups of vertebrates modeled were mammals, reptiles and amphibians. Only 14 species of amphibians had more than 10 unique geographic records, so we decided to model all of them. Regarding mammals and reptiles, once we separated the species with 10 unique records we randomly selected 20 species from each taxonomic group. The random sample had the goal of giving us the representativeness over the group integrating species with different current risk category, distribution (large vs. restricted) (see Appendix II), and other biological characteristics such as population status and dispersal ability (even though this information is not available for most species). The aim was getting an insight of how the groups will be impacted in general by climate change conditions. The list of evaluated species is given in Appendix II.

2.2 Climate change scenarios

The species ecological niche was modeled under a current climatic scenario (1961-2000) and two-time series (2015-2039 and 2075-2099) under two different representative concentration pathways (RCPs 4.5 and 8.5 Wm^-2) (^{IPCC, 2014}) at a 1 km² resolution, using regionalized climatology created with a weighted ensemble (Reliability Ensemble Averaging, REA) of 21 general circulation models (^{Fernández-Eguiarte et al., 2015}). We modeled two contrasting RCPs scenarios to get an idea of the different outcome possibilities. The worst-case emission scenario is the RCP 8.5 Wm^-2, and it helps having an insight of the consequences under the greatest changes expected in climate. The moderate scenario, RCP 4.5 Wm^-2, gave us an insight of what would happen in terms of species geographic distribution if some positive actions take place to reduce future impacts.

To run ENMs we used 19 bioclimatic variables that represent annual tendencies and temporal variation of climatic values that can be related to more biological meaningful results (^{Nix, 1986}). We generated these variables with the package “dismo” (^{Hijmans et al., 2014}) of the open software R (^{R Core Team, 2014}) and monthly values of precipitation, and maximum and minimum temperatures. The geographic area in which our species were projected was the entire Mexican territory (M). The geographic area or M in a BAM diagram represents the reachable geographic area for the species, while B are the biotic factors of the species niche and A the abiotic factors (^{Soberón and Nakamura, 2009}).

To select the geographic area M we should consider which BAM diagram is more plausible for our species. The main assumption for species distribution models using climatic change scenarios is that abiotic variables will be the main species constraints, and any reachable area for our target species cannot be occupied due to a lack of favorable environments. It would be a BAM diagram in which A is a subset of M and B is disregarded (^{Saupe et al., 2012}). We considered the Mexican territory as the potential area to be occupied in a free dispersion scenario for all species (M), and climatic conditions (A) the principal species constraints to reach any region of Mexico (A is a ubset of M ). We decided to include the entire Mexican territory, because even though we are modeling island species, most of them are not endemic to islands but are also distributed in the continent. Additionally, using a single M for all evaluated species allowed us a fair comparison of loss or gain area under climate change conditions. In order to reduce collinearity between the climatic variables, we performed a Pearson correlation test and selected nine variables that were not related to each other (determination coefficient value threshold R ² = 0.8) in the geographic area where models were projected (M) (^{Soberón and Nakamura, 2009}).

The following variables were used to create the model: Bio1: mean annual temperature; Bio2: mean diurnal range; Bio4: temporal variability; Bio6: minimum temperature of the coldest month; Bio7: annual range temperature; Bio12: annual precipitation; Bio13: precipitation of the wettest quarter; Bio18: precipitation of the hottest quarter, and Bio19: precipitation of the coldest quarter. This test was performed with the function cor in the R-based package (^{R Core Team, 2014}).

2.3 Modeling with BIOMOD

The BIOMOD package (^{Thuiller et al., 2009}) is a platform for ENM that works as a multimodeling system where several algorithms can be used simultaneously and helps creating different types of model ensembles (^{Thuiller et al., 2009}). One of the uncertainties in climate change assessments is that projections of species distributions can vary between algorithms (^{Pearson et al., 2006}). Ensembles are a way to manage and reduce uncertainty (^{Araújo and New, 2007}). In this study we selected two algorithms: maximum entropy (MAXENT [^{Phillips et al., 2006}]) and GAM (general additive models) that have shown a good performance in comparison to other methods (^{Elith et al., 2006}). We used the default parameterization, creating randomly 10 000 pseudo-absences and 10 replicates per algorithm. We split species records and used 70% to create the model and 30% to validate it. As a model assembly rule, we chose the total consensus, which was weighted by the sum of probabilities of the model through the true skill statistic value (TSS) ≤ 0.7, meaning that every model used for the ensemble had at least that validation value. The threshold for the binary maps (presence/absence) was also de maximization of the TSS. We chose the TSS test because it has been shown that it is the best statistical option to interpret binary maps dealing with ecological phenomena (^{Allouche et al., 2006}). We also used the BIOMOD range change function to estimate the proportion and relative number of lost, gained and stable pixels for each time period and RCPs.

Only four out of 54 species evaluated had less than 20 records, one mammal and three reptiles (see Appendix II); however, for species with less than 20 unique records a jackknife analysis was carried out. The jackknife method helps validating the predictability of models constructed with less than 20 data (^{Pearson et al., 2007}). For species having more than 700 records, we used the SpThin packages (^{Aiello-Lammens et al., 2015}). This library helps eliminating a specific amount of records to reduce bias resulting from oversampling and keeps the useful information in terms of unique records. It works by calculating the minimum geographic distance between two points and randomly eliminating the ones which distance is equal or shorter than this minimum distance, where thin.par = 10 km.

3.1 Biodiversity in the islands of Mexico

As we described above, we had to compile the information about vertebrates registered in Mexican islands. To our knowledge, this is the most complete list of Mexican vertebrate species with an insular distribution. Our list is made up of 87 mammals, 228 birds, 117 reptiles and 20 amphibians. From this list, that includes all terrestrial vertebrates (mammals, birds, reptiles and amphibians), there are 10 species not recognized by ^{GBIF.org (2017)} and other six that do not have any geographic record (Appendix I). Most bird species that have been registered in Mexican islands are migratory (Appendix I), which represents a more complex system to model (^{Pérez-Moreno et al., 2016}); however, we encourage future assessments with a specific type of analysis.

3.2 General tendencies in climate change vulnerability

For each evaluated species we reported unique georeference data, the percentage of area loss in the future with respect to the area that is climatically suitable under current conditions, the percentage of new area that will have climatically suitable conditions under specific climate change scenarios in the future, and the net change, which is the difference between the future suitable area and the area loss (Appendix II). This information is presented for both RCPs 4.5 and 8.5 Wm^-2 and for both short- and long term time periods. If we assume that species will not have any kind of restriction to disperse and migrate (unlimited dispersion scenario), then it is possible to take into account the net change. However, for precaution in terms of conservation we cannot forget about the limited dispersion scenario, in which species are restricted by biological/ecological characteristics and anthropogenic barriers.

In general, all evaluated species from all groups presented climatically suitable areas that may be lost in both RCPs and time frames (Appendix II). Still, increases in suitability areas were in general higher than losses, meaning that in the future there will be more climatically suitable areas for most species. Significant statistical differences were observed between RCPs in terms of area loss, gain and net change at both time frames (Wilcoxon test, P ≤ 0.05). Significant differences were only found in the loss and gain of climatically suitable areas between scenarios 4.5 and 8.5 Wm^-2 in the long term (2075-2099) and in losses, gains and net change between times (4.5 Wm^-2 [2015-2039] vs. 4.5 Wm^-2 [2075-2099], and 8.5 Wm^-2 [2015-2039] vs. 8.5 Wm^-2 [2075-2099]).

To identify vulnerability between the studied groups, we created boxplots of the percentage of climatically suitable loss (Fig. 1) under a limited dispersion scenario. In this figure it is possible to distinguish that in both extremes analyzed, the “best case scenario” (4.5 Wm^-2) in the short term and “the worst-case emission scenario” in the long term (8.5 Wm^-2), reptiles are the species with greatest losses followed by amphibians and mammals. In the long term, the variability in area loss increases up to almost 100% of the area that is climatically suitable under current climatic conditions (Fig. 1).

Fig. 1 Percentage of area loss in all evaluated groups in terms of species richness (number of species in a geographic area) for all taxonomic groups. The greatest species distributional shifts are shown in central and northern Mexico. Mammals’ richest areas under current climatic conditions can be found in the Baja California peninsula. With time, the richest areas decrease and besides shifting towards the north, they show an eastward trend in colonization. For reptiles, the richest areas were found mainly in the coast of Baja California and Sonora, but under future scenarios the areas with higher number of species decrease and shift eastward toward the center of Mexico. A similar pattern can be found for amphibians (Appendix III). 

3.2 Top ten vulnerability list

We ranked species and selected the top 10 with greatest losses under a limited dispersion scenario and the top 10 that will be more affected even when taking into account future suitable regions-i.e., unlimited dispersion scenario-(Table I). Under a limited dispersion scenario, reptiles are the group that has more species in the top 10 with greatest losses, followed by amphibians. Mammals are almost not present in the top 10 list (Table Ia, c). When net change is analyzed, reptiles continue to be a vulnerable group followed by amphibians and mammals. The current suitable area loss in these groups is so great that it cannot be compensated by future new areas with suitable conditions.

3.3 Vulnerability inside the three vertebrate groups

4.1 General tendencies in climate change vulnerability

This study is a good first approach to imrpove our knowledge on species diversity in Mexican islands and getting an insight of their vulnerability under climate change conditions. The compilation of species found in these islands allows us to understand that greater efforts are needed to increase our knowledge on which species are present in our territory. If this basic information is incomplete, it is very difficult to channel conservation efforts. Several species of different groups of vertebrates are present in Mexican islands but lack the minimum occurrence records (Appendix I). An increase in sampling efforts is needed for species mentioned in our appendix as “without records”, which are 16 reptile species and two birds (Appendix I).

In terms of vulnerability to climate change, species that will suffer in the future great losses of their current climatically suitable area do not necessarily present a negative net change. Species might be gaining new potential distribution areas bigger than area losses. However, the fact that new climatically suitable areas will exist in the future does not mean that species will be able to reach them. There are factors such as dispersion ability and ecological or anthropogenic barriers that might keep species from relocating to the new suitable regions (^{Pearson and Dawson, 2003}; ^{Parmesan 2006}). Consequently, it is important to keep in mind the question to be answered in order to decide which result is needed: the unlimited or the limited dispersion scenario. We think that in terms of conservation it will be important to use both results. Both scenarios can channel our efforts regarding conservation plans. For example, if ecological and anthropogenic barriers were obvious, then assuming that species will be able to reach new suitable areas would be unrealistic. Those species could be considered vulnerable even when its net change is positive. However, translocation could be an alternative. For species without new potential distribution areas, preservation efforts must take place in their current distribution area. In general, our results project positive net changes, meaning that new climatically suitable areas were greater than area losses for most species evaluated. We could also find significant differences in area losses between scenarios RCP 4.5 and RCP 8.5 Wm^-2 in the long term, meaning that if humanity chooses to follow an ecological way of life and reduces emissions, negative impacts in the distribution (and possible extinctions) of species will be significantly reduced. On the other hand, the richest areas for all evaluated taxonomic groups are expected to shift and change under future scenarios. In other words, hotspots of species richness will probably change and consequently current hotspots should not be the only tools used to concentrate conservation efforts.

In general, reptiles are the species with greatest losses and lower net change values (Fig. 1 and Table I). This result was very surprising, because amphibians as a group have been reported to be very vulnerable to environmental and climate changes (^{Blaustein et al., 2003}; ^{Foden et al., 2013a}). However, it is important to keep in mind that even when their net change or losses in the future are not greater than those of reptiles, most species in this group have very restricted distributions and other biological characteristics (such as the strong dependence to specific timing and amount of precipitation in order to fulfill their life history requirements) that make them sensitive to climate change (^{Wake, 2007}).

4.2 Vulnerability in the three vertebrate groups