Predicting the potential distribution of the endemic seabird Pelecanus thagus in the Humboldt Current Large Marine Ecosystem under different climate change scenarios

Species records

Pelican nesting and occurrence data were compiled from four main sources: the Neotropical Waterbird Census (https://lac.wetlands.org/), eBird (https://ebird.org/), the Global Biodiversity Information Facility (https://www.gbif.org/), and the literature. The geo-coordinates for each data point were referenced from the information in the literature or through the use of coordinates in Google Earth. We excluded duplicate or unclear locations and verified the accuracy of the data. We found a total of 4,818 georeferenced data points referring to pelican sightings (in resting place, nesting sites, coves, beaches, etc.), encompassing its entire geographic distribution from 2000 to 2015. Of these records, a subsampling was performed at a distance of 15 km (cell size), obtaining a total of 264 records, with which the modeling was performed. This subsampling were conducted in R, version 3.0.2 (R Development Core Team, 2013). The breeding distribution of the pelican was modeling with information for 34 nesting sites (Vinueza, Sornoza & Yáñez Muñoz, 2015; Zavalaga, 2015; Cursach et al., 2018).

Environmental variables

The environmental variables used to characterize the current distribution (and breeding) of the pelican were selected based on climate and oceanography. The climate variables used in this study were downloaded from the EcoClimate database (http://www.ecoclimate.org) (Lima-Ribeiro et al., 2015). These variables were represented by maximum, minimum, and mean values of monthly, quarterly, and annual temperatures, and the precipitation values recorded between 1950 and 2000. These parameters provided a combination of means, extremes, and seasonal differences in variables known to influence the distribution of species (Root et al., 2003). With the species distribution modeling toolbox extension implemented in ArcGIS, all bioclimate variables that showed a correlation higher than 0.7 were eliminated (Brown, 2014). Finally, six climate variables were selected: annual mean temperature, mean daytime temperature range, isothermality, seasonality in temperature, annual precipitation, seasonality in precipitation. The oceanographic variables used were sea surface temperature (SST) and marine net primary productivity (mg C m⁻² day⁻¹), as they are considered the main descriptors of the spatial distribution of seabirds (Quillfeldt et al., 2015; Ingenloff, 2017). These variables were obtained from the National Oceanic and Atmospheric Administration (NOAA, http://www.ngdc.noaa.gov/). For the analyses, we used mean values per climate season for a period of nine years (2004 to 2013), totaling eight oceanographic variables. All environmental variables used in this study were interpolated by the kriging method, with a uniform resolution of 0.5° × 0.5° using the QGIS 3.2.0 software (Lima-Ribeiro et al., 2015; Varela, Lima-Ribeiro & Terribile, 2015).

To evaluate the effects of the different climate change scenarios on the spatial distribution of pelicans, we did not include the oceanographic variables. The future climate scenarios corresponded to those proposed by the Intergovernmental Panel on Climate Change (IPCC, 2014). These scenarios were obtained from the ecoClimate website (http://ecoclimate.org/), which contains climate models available for different temporal intervals. To do this, we used the model developed by the Community Climate System Model version 4 of the National Center for Atmospheric Research (Gent et al., 2011). This is due to the good results for the South-East Pacific (Larson, Pegion & Kirtman, 2018; Zheng et al., 2018).

The projections for the six preselected variables and the projected minimum and maximum trajectories of the concentrations of greenhouse gases were obtained. That is 2.6 and 8.5 rcp (representative concentration pathways), respectively. These values indicate increases in the heat absorbed by the planet Earth due to the concentration of greenhouse gases up to 2100, in each trajectory and expressed in watts per square meter. Thus, 2.6 rcp is the moderate projection for the scenario with the least climate change, whereas 8.5 rcp is a more pessimistic projection and represents a severe scenario with the greatest climate change (Taylor, Stouffer & Meehl, 2012).

Modeling of the potential geographic distribution

The MaxEnt software (MaxEnt version 3.3.3k, http://www.cs.princeton.edu/ schapire/maxent/) has been frequently used for species distribution models under current and future climate scenarios (Phillips & Dudík, 2008). We used MaxEnt to model the geographic distribution of the pelican, including under two previously described climate change scenarios (Elith et al., 2006; Taylor, Stouffer & Meehl, 2012). The model was elaborated by MaxEnt auto-features (5,000 iterations). Logistic output was used for all analyses. The quality of the model was evaluated using the area under the curve (AUC) and the continuous Boyce index (Hirzel et al., 2006). AUC values can vary from 0 to 1, where a value greater than 0.9 is considered an indicator of “good” discrimination skills (Peterson et al., 2011). Values of the Boyce index vary between −1 and 1, where positive values indicate a model with predictions that are consistent with the distribution of observed presences in the evaluation dataset (Boyce, 2002). Both analyses were conducted in R using the “biomod2” package (R Development Core Team, 2013).

For each distribution model, a 30-fold cross-validation was used, with a data proportion of 25% for training and 75% for evaluation. The most important environmental variables were identified by estimating the relative contribution (%) to the model (Phillips, Anderson & Schapire, 2006). Jackknife test was used to evaluate the importance of the environmental variables for predictive modeling (Almalki et al., 2015).

Model yield for potential distribution

The model of presence with the best fit showed a gain of 3.04 and a Boyce Index of 0.99. Also, an AUC_training of 0.98 and an AUC_evaluation of 0.98 and a standard deviation of 0.004. While the modeling of breeding distribution showed a gain of 2.24 and a Boyce Index of 0.98, with an AUC_training of 0.98 and an AUC_evaluation of 0.98 and a standard deviation of 0.003. The AUC values were relatively similar, so the models used are appropriate for predicting the presence and breeding distribution of the species. AUC_evaluation 0.98 indicates that the pelican has a wide geographic distribution and breeding in relation to the area corresponding to the environmental data. The model predicts that the potential geographic distribution of the pelican reaches an approximate surface area of 466,836 km², latitudinally distributed from southern Ecuador (2°13′09″S) to southern Chile (46°59′07″S). Over this extensive marine–coastal surface, the probability of occurrence for this species varied between a 0.16 (minimum) and 0.84 (maximum) (Table 1). Areas with the highest probabilities of occurrence for the pelican are represented with intense red colors in Figs. 1A and 2A. These areas are mainly distributed from northern Peru to central Chile.

Importance of environmental variables

Among the six climatic variables and eight oceanographic variables, the mean daytime temperature range (Bio2) and the summer marine primary productivity, contributed the most to the current and potential distribution of the pelican (Table 2). These two factors explained 78.47% of the modeled distribution. The mean daytime temperature responded to the probability of the presence of the pelican, with a high probability of finding the species in areas where the mean daytime temperature ranges between 6 and 8 °C. In turn, the summer marine primary productivity also influenced the probability of the presence of the pelican, with a greater probability of finding the species during the summer season in areas with high primary productivity. The other factors such as, spring marine primary productivity, isothermality, and seasonality in temperature, contributed 9.24%, 3.23%, and 1.74%, respectively, to the modeled distribution. Therefore, thermal and primary productivity conditions are more important than other variables for mapping pelican distribution (Table 2).

The modeling of breeding distribution showed that the mean daytime temperature range contributed with 91.5% to the model, while the summer marine primary productivity contributed with 8.5%.

Potential geographic distribution of the pelican as a function of climate change

Based on the six climatic variables selected in the study, the model predicts that the projected pelican distribution currently attains an area of 596,753 km² (Table 3). This area is larger than that initially projected (466,836 km²), where the oceanographic variables were integrated. Regarding the projections of climate change for 2100, under the moderate scenario of 2.6 rcp a slight decrease (−0.68%) in pelican spatial distribution is predicted (Table 3). Under the severe scenario of 8.5 rcp, a slight increase (4.51%) in pelican spatial distribution is predicted (Table 3).

The projected habitable surface area under climate change of 2.6 rcp does not presents a major change with respect to the current geographic distribution of the pelican (Table 3). Under the severe scenario, the model predicts that the pelican habitable surface will vary depending on geographic area (Fig. 1). For example, in northern Chile its habitable surface would decrease, whereas in central and southern Chile it would increase over time (Fig. 2). The projected habitable surface area and the probabilities of occurrence for the pelican are spatially schematized in Fig. 2.

For the case of the modeling of breeding distribution, an area of 435,640 km² is projected (Table 4). Regarding the projections of climate change for 2100, under the moderate scenario of 2.6 rcp an increase (8.77%) in pelican breeding distribution is predicted (Table 4). Under the severe scenario of 8.5 rcp, a decrease (−19.30%) in pelican breeding distribution is predicted (Table 4). Under the severe scenario, the model predicts a decrease of occurrence probability of nesting sites of the pelican in northern Ecuador and north-central Chile (Fig. 2).