Study area

The study area comprises the entire territory of Mexico, between 14° N and 32° N, 86° W and 118° W. Mexico shares land borders with the US (to the north) and Guatemala and Belize (to the south). The total area of Mexico is approximately 1 900 000 km² and is characterized by a great diversity of landscapes and climates. Approximately 85 % of the country (except the Yucatán Peninsula and the coastal plains of the east and northwest) is formed by mountain ranges, plateaus and numerous valleys.

Desmodus rotundus occurrence data

To generate the present potential distribution map for D. rotundus, the study area was delimited using the known distribution of the species in Mexico.

We collected D. rotundus occurrence data from the National Mammal Collection of the Institute of Biology and the Museum of Natural History of the National University of Mexico (UNAM), the Mexican Commission for the Knowledge and Use of Biodiversity (Conabio), and the Global Biodiversity Information Facility.²⁶ Data collected for each record included locality (geographic coordinates in decimal degrees), height above sea level and source of information. We used only records from 1980 to 2013, avoiding old and imprecise data.

Bovine paralytic rabies and population data

We used animal health weekly reports from the National Agro-Alimentary Health, Safety and Quality Service (Senasica),²⁷ to develop a database of cases of bovine paralytic rabies from 2007 to 2015. We also used information about the national pastoral range-livestock population from the Agricultural, Livestock and Forestry Census 2007.²⁸ This data will be used to identify areas with more cattle and potential riesk areas. The information in both data sources was provided at the municipal level.

Environmental data

For this study, a total of 20 environmental variables were considered including 19 bioclimatic variables from temperature that reflect seasonal and annual trends of temperature and precipitation (Table 1). These variables were obtained from the WorldClim database ²⁹ for the 1950-2000 period. We used data for elevation above sea level, obtained from the Digital Elevation Model generated by the National Institute of Statistics, Geography and Informatics³⁰; and used ArcGIS 10.2³¹ to perform spatial analyses at a resolution of 30 arc-seconds (~ 1km²).

Table 1 Contributions of environmental variables used in the model of Desmodus rotundus in Mexico. 

Note: Bold font indicates variables that contribute ≥ 1 % to the model.

Source of data: WorldClim.²⁹

We used two Global Circulation Models (GCMs) to determine the likely distribution patterns of D. rotundus in Mexico for the decades of 2050 (average 2041-2060) and 2070 (average 2061-2080) at a spatial resolution of 30 seconds: Geophysical Fluid Dynamics Laboratory Model, version 3 (GFDL-CM3) and Hadley Centre Global Environment Model-ES version 2 (HadGEM2).²⁹ These climate change predictive models were proposed in the Fifth Assessment Report of the Intergovernmental Panel on Climate Change (IPCC) and the Coupled Model Intercomparison Project, Phase 5.³² Both GCMs considered two greenhouse gas concentration trajectories or Representative Concentration Pathways (RCP); 4.5 and 8.5. The former assumes annual emissions will decline near the year 2050 and a conservative increase in population growth, while the latter assumes CO₂ and methane emissions will continue to rise throughout the twenty-first century and a rapid increase in the human population.³²

Ecological niche model

We used the correlative maximum entropy based model or MaxEnt program ^33,34 to assess the potential distribution of D. rotundus in Mexico. MaxEnt is a robust algorithm that provides informative results and requires presence-only data.^20,35-37 MaxEnt is a correlative model based on an adjustment Bayesian procedure under the principle of maximum entropy that is used to estimate probability distributions (i.e., more uniform distributions) based on records of the presence of a species. It is subject to restrictions given by the environmental information or other variables such as topography.^33,34 To build the model, we divided the total number of records into two sets, a calibration data set (random selection of 75 % of the data as training data) and a set of test data (25 %) in 100 replicates in 500 iterations with different random partitions (bootstrap method). To test the performance of the model and its discriminatory capacity, we used the AUC_test (area under the ROC curve based on the testing data), where values close to one show a good prediction and values close to 0.5 indicate that the model behaves randomly.^33,38 It is important to note that although some studies have identified the limitations of the AUC_test as a performance measure of models,^39,40 this test is often used as a single threshold-independent measure for model performance and is a highly effective measure of the performance of ordinal score models.^36,41 The statistical significance of the model was evaluated using a binomial test which determines whethe the models are better than randomly obtained models.^33,38 Likewise, a Jackknife test was conducted to measure and determine which of these implemented variables were the most important for the presence of the species.³³ Subsequently, we generated models of current and future distribution based on the probability threshold that minimizes the error of omission (≤ 10 % training and validation) and the fraction of the predicted area to avoid overprediction.^20,38,42 We used a multivariate environmental similarity surface (MESS) to assess the level of similarity of the environmental conditions of the predicted areas between the two GCMs employed.³⁶ The results allowed us to distinguish areas of model extrapolation from those with similar environmental conditions at any site in the world.

An additional component in developing niche models is the mobility area (M). This area is frequently much larger than the actual distribution of the species in question.⁴³ Therefore, the predictions represent hypotheses about conditions that are similar to those where the species has been observed. These conditions are likely to be among the existing fundamental niche and the niche occupied by the species.²⁰ Because Desmodus rotundus has a wide dispersal capacity and information about the natural history is available, we considered the entire country of Mexico as our M in the models. Currently the species is not present in the northern states of Mexico such as Baja California, Sonora, Chihuahua and Coahuila.^13,15 However, these areas represent a potential opportunity to be colonized, if they present suitable environmental conditions in the future.

Areas at risk of bovine paralytic rabies

Predicted areas were identified as high risk due to probability estimates of the potential occurrence of D. rotundus in the study area. We superimposed maps of current and future D. rotundus distributions with database records of bovine paralytic rabies from Senasica²⁷ to identify municipalities with potential risk of disease. We then input the layer of estimated number of free-ranging cattle from the Agricultural, Livestock and Forestry Census 2007,²⁸ in each area to assess the potential extent of cattle susceptibility to the disease.

Bovine paralytic rabies distribution

Bovine paralytic rabies remains a threat for livestock production in Latin America and Mexico where the vampire bat, Desmodus rotundus, is the main vector.^4,11 From 2007 to 2015, 119 389 cases of bovine paralytic rabies were reported, of which 1 872 were confirmed, distributed in 477 municipalities in 25 Mexican states. States with the highest number of confirmed cases were Veracruz (342 confirmed cases, 25 923 reports), San Luis Potosi (211 confirmed cases, 2 624 reports), Yucatán (190 confirmed cases, 15 783 reports), Tabasco (181 confirmed cases, 16 132 reports), and Chiapas (127 confirmed cases, 13 181 reports). States with the lowest number of reported cases were the State of Mexico, Guanajuato and Baja California Sur. Municipalities with the greatest number of positive cases of bovine paralytic rabies in the same period were: Tierra Blanca (75) and Papantla (61) in Veracruz, Huimanguillo in Tabasco, Mapastepec (60) in Chiapas, and Candelaria (60) in Campeche.

The annual number of confirmed cases fluctuated between 138 and 288, recorded three maxima: 288 cases in 2009, 275 in 2010 and, 243 in 2012 (Figure 1). When comparing seasons, we observed that most cases were reported during the dry season; a total of 1064 cases were reported between January and May, with two peaks in March and May (200 and 234 cases, respectively; Figure 2). Note that in this study we used correlative maxima entropy models to develop maps of the current distribution of suitable Desmodus rotundus habitats in Mexico. We then used climate change predictions to determine how changes in bioclimatic factors influence potential distribution of D. rotundus, the principal vector of rabies virus within pastoral range-livestock production systems, to identify potential bovine paralytic rabies risk areas.¹²

Figure 1 Annual number of cases of bovine paralytic rabies in Mexico, 2007-2015. 

Figure 2 Number of monthly cases of bovine paralytic rabies in Mexico, 2007-2015. 

Current distribution of Desmodus rotundus

The current potential distribution of D. rotundus we generated from the 649 confirmed records for the period 1980-2013, performed well and had good adjustment (AUC_test value of 0.95 ± 0.029) and a low rate of omission (P < 0.0002). The most important variables associated with the presence of D. rotundus (Table 1) were temperature seasonality (49.5 %), followed by minimum temperature of coldest month (16.6 %) and annual temperature range (10 %). The areas identified with a greater likelihood of the presence of D. rotundus maintain a temperature difference of less than 4 °C between the colder and warmer quarters, and both variables had low values, indicating low variation in temperature throughout the year.

Model prediction results for the current scenario indicate habitat suitability for Desmodus rotundus covers an area of 834 965 km² (Figure 3), and the areas with the most appropriate conditions (probability > 0.49) are located along the coastal areas: The Pacific coast, from southern Sonora to the southwestern corner of the country, and the Gulf of Mexico coast, from central Tamaulipas to the Isthmus of Tehuantepec, Chiapas and the Yucatán Peninsula (Figure 3). Areas with a high probability of suitable conditions were primarily in the states of San Luis Potosi and Veracruz and in the southeastern region of the country. There were areas in the northern regions of the country, such as the states of Tamaulipas, Sonora and Sinaloa, in which D. rotundus are present in tropical and subtropical environments with warm temperatures that do not fluctuate more than 2 °C annually. However, the transformation of natural environments to grasslands for the development of the cattle industry, and the micro-climatic changes fostered by these disturbances have favoured the expansion of D. rotundus into other regions.^12,15,44

Figure 3 Current areas in Mexico identified as suitable for Desmodus rotundus using MaxEnt model. Blue points show localities of D. rotundus occurrence from 1980 to 2013. 

Future distribution of D. rotundus

Predicted risk maps for the 2050 and 2070 climate scenarios show a reduction in potential distribution of Desmodus rotundus as previously identified in the current scenario. The best predictive model of D. rotundus distribution under each of the two GCMs, GFDL-CM3 and HadGEM2, had AUC_test values of 0.97 ± 0.02, indicating a high predictive power for both models. Under a conservative scenario of RCP 4.5, the distribution models under both GCMs show a marked reduction of total area compared to the current values: a reduction of 20 % (170 317 km²) for the decade of 2050 and 28 % (233 989 km²) for the decade of 2070 (Figures 4, 5). A similar pattern was observed when the more development-intensive RCP 8.5 scenario was used: a reduction of 27 % (229 486 km²) for the decade of 2050 and 34 % (283 471 km²) for the decade of 2070 (Figures 4, 5). Comparison of the species distribution models revealed that with respect to the current geography, GFDL-CM3 GCM predicted reductions of 15 % and 33 % for both periods of time, whereas HadGEM2 predicted a reduction of 9 % and 12 % for 2050 and 2070, respectively.

Figure 4 Modelled suitability for future distribution of Desmodus rotundus according to Global Climate Model GFDL-CM3 for two time periods (2050 and 2070), and two Representative Concentration Pathways (RCP 4.5 and 8.5). Left-hand column shows suitability values, with blue indicating more suitable conditions. 

Figure 5 Modelled suitability for future distribution of Desmodus rotundus according to Global Climate Model HadGEM2 for two time periods (2050 and 2070), and two Representative Concentration Pathways (RCP 4.5 and 8.5). Left-hand column shows suitability values, with blue indicating more suitable conditions. 

The MESS analysis, which showed broad areas of similar environmental conditions under current and predicted future climatic conditions, also identified suitable areas where D. rotundus is currently absent, including central Mexico, Coahuila, Baja California, western Sonora, and northern Chihuahua. Under all IPCC climate change scenarios analyzed, D. rotundus expanded its range into the northwest of the country, occupying new areas in the states of Sonora, Sinaloa, Zacatecas, Aguascalientes, and Guanajuato, the northeast regions of Mexico in the states of Coahuila, Nuevo León and San Luis Potosi, and the Central Mexican Plateau, principally the states of Mexico, Puebla and Hidalgo.^8,12 In contrast, the reduction in distribution occurred primarily in the southeastern portion of the country, notably in the Gulf Coastal Plain in southern Veracruz, Tabasco, southwestern Campeche, and the Yucatán Peninsula (Figures 4, 5).

It is important to note that while there may be some shift in the geographic patterns of the distribution of D. rotundus under future climate change scenarios, these changes are not significant on a nation-wide scale, and the species will continue to be located primarily in tropical and subtropical regions of the country.^12,15 However, an increase in the suitable area for other vectors of the rabies disease and other biotic factors driving the dynamics of the disease, such as shelter and food availability and biological corridors, may contribute to an increased risk of the disease.^15,44,45

Our results indicate significant changes in the pattern of distribution and incidence of bovine paralytic rabies in Mexico. The highest risk for bovine paralytic rabies occurs in the dry season (April and May). This pattern is consistent with reports of a spike in the number of rabies cases in Brazil, Colombia and Costa Rica during the dry season.^46,47 We identified temperature as a principal determinant in the distribution of D. rotundus; it affects the spread of the rabies disease among tropical and semi-tropical environments and, the likelihood of its spread into more temperate areas that may experience increased temperature. We estimate that more than 600 000 head of cattle in Mexico could be exposed to bovine paralytic rabies virus by 2050-2070, according to the current number of free-ranging cattle from the Agricultural, Livestock and Forestry Census 2007.²⁸ Nevertheless, an increase of suitable areas for other vectors of the rabies disease and other biotic factors driving the dynamics of the disease, such as shelter and food availability and biological corridors, may contribute to the increased risk of the disease.^8,44,48 This number is likely an underestimation because the cattle population is likely to exceed the values reported by the 2007 census that were used in the model and because of the ongoing and projected conversion of natural ecosystems to grazing land to accommodate livestock production.