One hundred years of climate change in Mexico

RESEARCH ARTICLE One hundred years of climate change in Mexico Angela P. Cuervo-RobayoID1,2, Carolina Ureta3, Miguel A. Go´ mez-AlboresID4, Anny K. Meneses-MosqueraID1, Oswaldo Te´llez-Valde´ s5, Enrique Martı´nez-Meyer1* 1 Departamento de Zoologı´a, Instituto de Biologı´a, Universidad Nacional Auto´ noma de Me´ xico, Ciudad de Me´xico, Me´ xico, 2 Comisio´ n Nacional para el Conocimiento y Uso de la Biodiversidad (Conabio), Ciudad de Me´xico, Me´ xico, 3 Ca´tedras-Departamento de Ciencias Atmosfe´ricas, Centro de Ciencias de la Atmo´ sfera, Universidad Nacional Auto´noma de Me´xico, Ciudad de Me´xico, Me´ xico, 4 Instituto Interamericano de Tecnologı´a y Ciencias del Agua, Universidad Auto´ noma del Estado de Me´ xico, Toluca, Estado de Me´xico, Me´xico, 5 Facultad de Estudios Superiores Iztacala, Unidad de Biotecnologı´a y Prototipos, Laboratorio de Recursos Naturales, Universidad Nacional Auto´ noma de Me´xico, Tlalnepantla, Me´ xico * emm@ib.unam.mx Abstract Spatial assessments of historical climate change provide information that can be used by scientists to analyze climate variation over time and evaluate, for example, its effects on bio- diversity, in order to focus their research and conservation efforts. Despite the fact that there are global climatic databases available at high spatial resolution, they represent a short tem- poral window that impedes evaluating historical changes of climate and their impacts on bio- diversity. To fill this gap, we developed climate gridded surfaces for Mexico for three periods

that cover most of the 20th and early 21st centuries: t1-1940 (1910–1949), t2-1970 (1950– 1979) and t3-2000 (1980–2009), and used these interpolated surfaces to describe how cli- mate has changed over time, both countrywide and in its 19 biogeographic provinces. Results from our characterization of climate change indicate that the mean annual tempera- ture has increased by nearly 0.2˚C on average across the whole country from t2-1970 to t3- 2000. However, changes have not been spatially uniform: Nearctic provinces in the north have suffered higher temperature increases than southern tropical regions. Central and southern provinces cooled at the beginning of the 20th century but warmed consistently since the 1970s. Precipitation increased between t1-1940 and t2-1970 across the country, more notably in the northern provinces, and it decreased between t2-1970 and t3-2000 in most of the country. Results on the historical climate conditions in Mexico may be useful for climate change analyses for both environmental and social sciences. Nonetheless, our cli- matology was based on information from climate stations for which 9.4–36.2% presented inhomogeneities over time probably owing to non-climatic factors, and climate station den- sity changed over time. Therefore, the estimated changes observed in our analysis need to be interpreted cautiously. PLOS ONE PLOS ONE | https://doi.org/10.1371/journal.pone.0209808 July 16, 2020 1 / 19 a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 OPEN ACCESS Citation: Cuervo-Robayo AP, Ureta C, Go´mez-

Albores MA, Meneses-Mosquera AK, Te´llez-Valde´s O, Martı´nez-Meyer E (2020) One hundred years of climate change in Mexico. PLoS ONE 15(7): e0209808. https://doi.org/10.1371/journal. pone.0209808 Editor: Lucas C.R. Silva, University of Oregon, UNITED STATES Received: December 7, 2018 Accepted: June 19, 2020 Published: July 16, 2020 Copyright: © 2020 Cuervo-Robayo et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Data Availability Statement: The National Commission for Biodiversity (CONABIO) is the repository for the data. We have provided all the information to CONABIO. We have included the URLs to access the data in Table 1. Funding: The Mexican Consejo Nacional de Ciencia y Tecnologı´a (CONACyT) provided a Ph.D. scholarship to A.P.C-R (No. 217650). This work was funded by PAPIIT-DGAPA-UNAM to O.T-V. (Grant number IN-225010).

ings that impede evaluating the historical climate change and their impacts on biodiversity; for example, they span only a single period of time in the late 20th century. Furthermore, the qual- ity and reliability of global databases are compromised in many regions because the data from weather stations used for producing global interpolations generally come from open, public databases, whereas in some places–such as Mexico–there is more information available, but needs to be requested. Also, there is more time and effort devoted to detecting and correcting errors in databases if people are focused on a particular region than worldwide; for example, reconciling location (coordinates) with elevation of stations has been more thoroughly done at regional levels [14,15]. Until now, gridded climate surfaces for Mexico are only readily avail- able for single time slices [15–17], so they are not suitable for evaluating temporal climatic changes. Regardless, time series information at the level of climate stations is available for cli- mate change analysis [18,19], although some regions across the country lack adequate spatial coverage. One option to overcome this limitation is to analyze data from climatic surfaces at different periods from the recent past and the present in an attempt to understand where, when and to what degree the climate has been changing, in an effort to evaluate its environmental and social implications [2,20–28]. For instance, combining data of climate surfaces with species

inventories can be useful for detecting drivers of change in biological communities [28,29] and responses of biodiversity to climate change [30]. Historic climate surfaces can be used to iden- tify areas that have been under the effect of climate change for some time to analyze the response and resilience of natural and social subsystems [31]. Furthermore, coupling this information with future projections can help to inform decision-makers [32–34]. This approach has been successful in Australia, where historical climate data facilitated conserva- tionists to detect recent range shifts in bird communities that may be informative for future biotic responses [35]. To date, most studies evaluating the impacts of climate change on biodiversity have focused, first, on future projections rather than under a retrospective approach, and second, on individ- ual species and ecosystems [1], and not at other organizational levels [36]. A more inclusive approach would be to assess climate change impacts in biogeographic units, which represent PLOS ONE Historical climate change in Mexico PLOS ONE | https://doi.org/10.1371/journal.pone.0209808 July 16, 2020 2 / 19 Competing interests: No authors have competing interests.

natural spatial areas that integrate physiographic, evolutionary and ecological features of biodi- versity [37]. However, the impacts of climate change at this level have yet to be performed [32]. To fill this gap, we developed average historical climate surfaces for Mexico, for three peri- ods covering most of the 20th and early 21th centuries: t1-1940 (1910–1949), t2-1970 (1950– 1979) and t3-2000 (1980–2009). We used these surfaces to describe historical climate change countrywide and in the 19 biogeographic provinces of Mexico. Materials and methods Climate data We analyzed monthly minimum and maximum temperature and accumulated rainfall data gathered from weather stations from the National Meteorological Office that were previously organized by Cuervo-Robayo [17], to derive monthly mean climate surfaces for three periods: t1-1940 (1910–1949), t2-1970 (1950–1979) and t3-2000 (1980–2009). We selected these periods based on previous global [38,39] and regional [19] climate change analysis and we also consid- ered the number of stations available for each period [40]. For t1-1940, we used a 40-year period, instead of 30 years as in the other time slices, due to the limited number of stations available. Because interpolation estimations are normally poor toward the edges of a region, the analysis should expand beyond the limits to the target region. Consequently, we included weather stations from southern portions of the United States, northern Belize and Guatemala. Data from the U.S. were gathered from the United States Historical Climatology Network

(USHC: http://cdiac.ornl.gov/epubs/ndp/ushcn/access.html), and data from Belize and Guate- mala were gathered through the FAOCLIM 2.0 software (http://www.fao.org/nr/climpag/pub/ en1102_en.asp) and the National Climatic Data Center (http://www.ncdc.noaa.gov data were organized and averaged into their corresponding periods using the “Structuration” extension from the Integrated Water Resources Management tool implemented in Idrisi Selva software, which is freely available at http://idrisi.uaemex.mx [41]. For the period t1-1940, weather sta- tions were located mainly in central Mexico, in the states of Guanajuato, Jalisco, Michoaca´n, and Me´xico. For the other two time periods, weather stations distributed throughout the coun- try, mainly at or near densely populated areas (Fig 1). We used ANUSPLIN software version 4.37 [42] to generate countrywide continuous cli- mate surfaces for the continental area (i.e., excluding oceans). This program uses the thin- plate smoothing splines interpolation technique, which integrates climatic and topographic data for making spatial climatic estimations, thus performing better than other interpolation methods [43,44]. Nonetheless, methods that incorporate expert knowledge, e.g. PRISM [45], sometimes perform better to capture features, such as rain shadow, which can be advantageous in complex geographical regions. Finally, we only used climate data from stations that operated for more than 10 years (S1) in at least one variable (temperature or precipitation). We used a

second-order spline with three independent variables (latitude, longitude and elevation) [17]. The value of the smoothing parameter was determined by minimizing a measure of the predic- tive error of the fitted surface given by the generalized cross-validation (GCV). For precipita- tion, we used a square-root transformation that reduces positively skewed values and ignores all negative values in precipitation data [42,46]. It also applies more smoothing to large precipi- tation values and less smoothing to small precipitation data values [47]. Given that we had more than 2000 stations for t2-1970 and t3-2000, we used SPLINB, as recommended by Hutch- inson [46], and used SELNOT to select a set of knots to reduce the complexity of the fitted spline. For t1-1940 we used SPLINA, which uses all available stations [46,47]. PLOS ONE Historical climate change in Mexico PLOS ONE | https://doi.org/10.1371/journal.pone.0209808 July 16, 2020 3 / 19

ANUSPLIN produces a list of the stations with the largest data residuals that denote fitting errors. With this list, Cuervo-Robayo et al. [17] corrected the geographic position for 100 erro- neous stations by using online gazetteers and Google Earth; however, for the two later periods, we eliminated about 20 stations that maintained high residual values (< 1) during the fitting process, without increasing error in the diagnostic statistics. This residual error may have resulted from the lack of a homogenization process in the stations’ data. We assessed the accu- racy of the fitted surfaces by examining ANUSPLIN diagnostic measures [42]. The signal indi- cates the degrees of freedom associated with the surfaces, which reflects the complexity of the surface and varies between a small positive integer and the number of stations used to generate the surface [46,48]. Hutchinson & Gessler [46] suggested that the signal should not be greater than about half the number of data points. Models with a signal below this threshold tend to be more robust and reliable in regions where data are scarce [49]. We also examined the root mean square error (RTMSE) and the square root of the GCV (RTGCV). RTMSE is an optimis- tic assessment of predictive error because some interpolation error is removed, therefore the true error is somewhere between the RTMSE and the RTGCV error [42]. Gridded climate Fig 1. Location of weather stations used to generate climate surfaces of precipitation, maximum and minimum temperature for three periods: a. t1-1940 (1910–1949), b.

t2-1970 (1950–1979) and c. t3-2000 (1980–2009). Shaded areas represent mountainous regions above 1500 masl. Stations used were not corrected for inhomogeneities and densities vary between periods and climate variables. For t1-1940 we used 803 stations for precipitation and 500 for temperature, for t2-1970, we used 3411 and 3670 and in t3-2000, 3870 and 4200 for precipitation and temperature, respectively. https://doi.org/10.1371/journal.pone.0209808.g001 PLOS ONE Historical climate change in Mexico PLOS ONE | https://doi.org/10.1371/journal.pone.0209808 July 16, 2020 4 / 19

surface (i.e., gridpoints) of land areas were generated with the function LAPGRD, using a 30-arc second of spatial resolution digital elevation model (https://lta.cr.usgs.gov/GTOPO30). As an additional output, we derived 19 bioclimatic variables for each period with the dismo library [50] of R software [51]. These variables included annual, quarterly and monthly sum- maries of temperature and precipitation that represent more biologically meaningful variables than the original climate surfaces and have been widely used in several studies of climate change impacts on species and ecosystems [12]. For t1-1940, we used 803 stations for precipita- tion and 500 for minimum and maximum temperatures. For t2-1970, the number of stations was 3411 and 3670 for precipitation and temperature, respectively. For t3-2000, the number of stations was 3870 and 4200 for precipitation and temperature, respectively (Fig 1). Across the country, 47% of the weather stations are located in mountainous regions above 1500 masl; the remaining 63% are located in the lowlands. Data used for the analysis are available upon request, and the original dataset is available at the National Meteorological Office (http://smn. cna.gob.mx/es/, S1 File). Monthly climate surfaces and bioclimatic variables for each period are freely available (see URLs in Table 1). In order to evaluate differences between periods, we performed a discriminant analysis on deseasoned climatic surfaces using the Statistica software version 10 [52]. Deseasoning is rec-

ommended to eliminate seasonality from the analysis that may obscure temporal changes. To do so, we calculated standardized anomalies (z-scores) for each month of the three periods using the Deseason Panel of the Earth Trends Modeler (ETM) of Idrisi Selva [53]. We also evaluated if the difference in the number of weather stations between t2-1970 and t3-2000 had an effect on the climate surfaces; we did not compare t1-1940 against the other two periods because station instrumentation and especially density was quite different for the first period from that of the last two periods [40]. In order to compare t2-1970 and t3-2000, we created a Table 1. URLs for downloading monthly gridded surfaces of maximum temperature, minimum temperature and precipitation, and 19 bioclimatic parameters. Metadata can be accessed through the metadata tab. Variable Period url Monthly maximum temperature (MaxT) t1-1940 (1910– 1949) http://geoportal.conabio.gob.mx/metadatos/doc/html/ tm191049gw.html t2-1970 (1950– 1979) http://geoportal.conabio.gob.mx/metadatos/doc/html/ tm195079gw.html t3-2000 (1980– 2009) http://geoportal.conabio.gob.mx/metadatos/doc/html/ tm198009gw.html Monthly minimum temperature (MinT) t1-1940 (1910– 1949) http://geoportal.conabio.gob.mx/metadatos/doc/html/ tmi191049gw.html t2-1970 (1950– 1979) http://geoportal.conabio.gob.mx/metadatos/doc/html/ tmi195079gw.html t3-2000 (1980– 2009) http://geoportal.conabio.gob.mx/metadatos/doc/html/ tmi198009gw.html Monthly precipitation (Ppt) t1-1940 (1910– 1949) http://geoportal.conabio.gob.mx/metadatos/doc/html/ p19101949gw.html

t2-1970 (1950– 1979) http://geoportal.conabio.gob.mx/metadatos/doc/html/ p19501979gw.html t3-2000 (1980– 2009) http://geoportal.conabio.gob.mx/metadatos/doc/html/ p19802009gw.html 19 Bioclimatic variables (Bio) t1-1940 (1910– 1949) http://geoportal.conabio.gob.mx/metadatos/doc/html/ b19101949gw.html t2-1970 (1950– 1979) http://geoportal.conabio.gob.mx/metadatos/doc/html/ b19501979gw.html t3-2000 (1980– 2009) http://geoportal.conabio.gob.mx/metadatos/doc/html/ b19802009gw.html https://doi.org/10.1371/journal.pone.0209808.t001 PLOS ONE Historical climate change in Mexico PLOS ONE | https://doi.org/10.1371/journal.pone.0209808 July 16, 2020 5 / 19

set of climate surfaces in which we used the same number of weather stations for both periods. We found that differences between t2-1970 and t3-2000 were not due to the number of stations, thus we included all available stations for the second and third periods, mainly because a larger number of stations improve model signal [17,54]. Biogeographic provinces of Mexico Mexico is recognized as a megadiverse country [55]. Its extraordinary biodiversity is due to its high environmental heterogeneity and the confluence of the Nearctic and Neotropical biogeo- graphic regions connected by a transition zone (Fig 2) [56]. The Nearctic region encompasses the arid subtropical areas of the north and includes the provinces of: California, Baja Califor- nia, Del Cabo, Sonorense, Altiplano Norte, Altiplano Sur, and Tamaulipeca. The Neotropical region covers the humid and sub-humid tropical areas of the south, including the provinces: Costa del Pacı´fico, Golfo de Me´xico, Depresio´n del Balsas, Oaxaca, Altos de Chiapas, Soco- nusco, Yucata´n, and Pete´n. The Transition Zone includes the mountainous provinces in the central portion of the country: Sierra Madre Oriental, Sierra Madre Occidental, Eje Volca´nico, and Sierra Madre del Sur [56,57]. Fig 2. Biogeographic regions and provinces of Mexico: Nearctic and Neotropical regions, and the transition zone. Reprinted from [57] under a CC BY license, with permission from [CONABIO], original copyright [1997]. https://doi.org/10.1371/journal.pone.0209808.g002 PLOS ONE Historical climate change in Mexico

PLOS ONE | https://doi.org/10.1371/journal.pone.0209808 July 16, 2020 6 / 19

A map of the biogeographic provinces of Mexico (Fig 2) was obtained from the Comisión Nacional para el Conocimiento y Uso de la Biodiversidad [57]. This regionalization was based on the distributional pattern of four taxonomic groups (vascular plants, amphibians, reptiles, and mammals) and the main morpho-tectonic features of Mexico [56]. Each unit represents a relatively homogeneous area with high levels of endemic species sharing similar historical, physiographic, climatic, edaphic, and vegetation features [57]. Climatic changes between time periods We estimated annual and seasonal changes of maximum temperature, minimum temperature and precipitation, at three scales: Mexico, the three biogeographic regions (i.e., Nearctic, Neo- tropical and the Transition Zone) and the 19 biogeographic provinces, using as source data the monthly ANUSPLIN interpolated climate surfaces (grid points). For estimating changes from t1-1940 to t2-1970 and from t2-1970 to t3-2000 at the country level and the three regions, first, we calculated descriptive statistics using the grid point values of the raster surfaces to summarize temperature and precipitation changes. We also estimated extreme warming and cooling at the 90th and 10th percentile distribution of the grid point val- ues of temperature changes and the median fractional increases and decreases in precipitation. Then, for the 19 provinces, we calculated seasonal changes from t1-1940 to t2-1970 and from t2-1970 to t3-2000 using the average values of the grid points of the climate surfaces calculated

with the software Earth Trends Modeler (ETM) of Idrisi Selva [52]. Seasons were defined as: (i) December, January, February (DJF); (ii) March, April, May (MAM); (iii) June, July, August (JJA); and (iv) September, October, November (SON). Changes in precipitation were esti- mated in millimeters (mm) and percentage (%), and changes in temperature in Celsius degrees (˚C). To visualize a monotonic upward or downward change through time in the interpolated cli- mate surfaces of annual precipitation (Bio 12), maximum temperature of the warmest month (Bio 05) and minimum temperature of the coldest month (Bio 06), we produced a map show- ing the consistency of the sign between periods as follows: at each grid point, we obtained the difference between t2-1970 to t1-1940 and between t3-2000 to t2-1970, if the sign was consistent between these two transitions we classified them either as Positive or Negative, otherwise, we classified them as Not consistent. Finally, we performed a Pettitt Test for Change-point Detection on annual averages time- series of all individual weather stations used for creating the climate surfaces across the coun- try [58]. The Pettitt test is a non-parametric analysis for detecting inhomogeneities in a time- series and identifying the possible year of breakpoint [59]. Inhomogeneities in climatic time- series may occur owing to climatic factors, or due to instrumental, observer, or positional changes in the stations, and they may distort the signal of climatic change [60]. We used the pettitt.test function [61] implement in the R software [51]. Results

Climate surfaces Diagnostic measures obtained from ANUSPLIN indicated that model fit of splines varied across historical time slices. The average ratio of the signal to the number of data points was <0.5 for monthly temperatures and precipitation, mostly for t2-1970 and t3-2000 (S1 Table in S1 File). For t1-1940, only a small number of weather stations (<900) had available precipita- tion data. Thus, we eliminated only those stations with higher residual error (>1) to maintain the largest number of stations (>800). The average signal for precipitation was above the per- mitted threshold (values greater than about half the number of data points), indicating that the PLOS ONE Historical climate change in Mexico PLOS ONE | https://doi.org/10.1371/journal.pone.0209808 July 16, 2020 7 / 19

number of stations is insufficient to depict the complex spatial patterns of mean precipitation [47]; however, for some months (i.e., from January to May) the ratio signal was acceptable. Therefore, the precipitation surface for period t1-1940 must be used with caution. The monthly average RTMSE for temperatures was below 0.6˚C and below 10 mm for precipitation of all three-climate periods. We found differences between the three time periods in the deseasoned climatic surfaces of monthly precipitation (Wilks’ Lambda = 0.295, Fdf = 24,178478 = 6249.079, P < 0.001), maxi- mum temperature (Wilks’ Lambda = 0.419, Fdf = 24,178478 = 4050.582, P < 0.001) and minimum temperature (Wilks’ Lambda = 0.398, Fdf = 24,178478 = 4346.174, P < 0.001) averaged across all Mexico. The Pettitt test showed that the number of weather stations holding inhomogeneities varied among the three climatic variables and biogeographic provinces. For maximum temperature, 899 stations out of 4152 (21.7%) presented inhomogeneities, ranging from 10.6–28.2% among provinces (S2 Table in S1 File). These numbers were higher for minimum temperature, for which 36.2% of 5079 presented inhomogeneities. For some provinces inhomogeneities occurred in more than 40% of the stations (S3 Table in S1 File). Inhomogeneities were lower for precipitation (9.4% of 5239 weather stations, ranging from 4–19% among provinces) (S4 Table in S1 File). The breakpoint years–which refers to the estimated year in which an abrupt transition occurred in the time series–were more frequent during the second half of t2-1970

and the first half of t3-2000 in all three climatic variables and across all provinces (S1, S2 and S3 Figs in S1 File). Changes between time periods The magnitude of change in annual mean temperature (Bio 01) and precipitation (Bio 12) dif- fered between the first and second time periods and between the second and third time peri- ods. At the country level, the mean annual temperature decreased, on average, 0.32˚C from t1- 1940 to t2-1970, and increased 0.22˚C from t2-1970 to t3-2000 (Table 2). This tendency was also consistent across regions, but the Neotropics experienced the greatest change between t1- 1940 to t2-1970, (-0.64˚C), while the Nearctic did so between t2-1970 and t3-2000 (-0.3˚C; Table 2). The largest increases (>1.5˚C) from t2-1970 to t3-2000 were found in the desert mountains of northern Mexico and Baja California (S4 Fig in S1 File). Additionally, extreme cooling and warming increased in magnitude since t2-1970 (Table 2). Across Mexico, maximum and minimum temperatures have increased since t2-1970 in most provinces and for almost all seasons, even when at the beginning of the century maxi- mum temperature decreased mainly in the southern provinces of the country (Fig 3). Between t1-1940 and t2-1970, Baja California and Yucata´n provinces suffered the greatest changes, with mean annual temperature decreases ranging from almost 2˚C to almost 4˚C in the summer (Fig 3). Between t2-1970 and t3-2000, there was a tendency of temperature increase in the northern provinces with no consistent change observed in the southern ones (Figs 3 and 5).

Within seasons, the pattern was similar to the changes observed in the mean annual tempera- ture (Table 2), as we observed an increase in temperatures throughout periods, particularly for the second half of the 20th century, although the changes in minimum temperature were not consistent between periods. This happened mainly in the biogeographic provinces of Depre- sio´n del Balsas, Oaxaca, Sierra Madre del Sur, Oaxaca, and Socunusco. Moreover, the increase observed in the minimum temperature in winter (DJF) and spring (MAM) between t2-1970 to t3-2000 was much larger in the northern provinces than elsewhere in the country (Fig 3). In general, across Mexico, the total area where positive changes in precipitation occurred during the century is larger than the area with negative changes. However, the magnitude of PLOS ONE Historical climate change in Mexico PLOS ONE | https://doi.org/10.1371/journal.pone.0209808 July 16, 2020 8 / 19

change in mean precipitation is smaller between t2-1970–t3-2000 than between t1-1940–t2- 1970 (Table 3). Mean changes in annual precipitation since t2-1970 ranged between 4% increase for the Nearctic to -5% for the Neotropical region. Despite the change in mean precip- itation in the Nearctic was positive, median conditions between t2-1970–t3-2000 have decreased. Changes in precipitation were higher in the Transition Zone and the Neotropical region than in the Nearctic region, with a reduction of precipitation of 3% to 8%, respectively (Table 3). Precipitation has generally decreased in the tropical provinces of the south and in the mountainous provinces of the Transition Zone, whereas it has increased in the Nearctic prov- inces for both temporal changes (Fig 4). In several southern biogeographic provinces, the mag- nitude of precipitation decreased in t2-1970–t3-2000, mainly during spring (MAM) and autumn (SON; Fig 4). Mean precipitation change was higher from t1-1940 to t2-1970 in five biogeographic provinces (Sierra Madre del Sur, Costa del Pacı´fico, Oaxaca, Soconusco, and Golfo de Me´xico) compared to the mean precipitation change from t2-1970 to t3-2000. In the remaining provinces, the precipitation change was positive in t1-1940–t2-1970, but magnitude decreased in t1-1940–t2-1970 (Fig 4). Monotonic upward or downward change in Bio 12, Bio 05 and Bio 06 occurred over 50% of the interpolated climate surfaces (Fig 5). A consistent reduction in maximum temperature (Bio 05) occurred in 13% of the country, whereas 12% of the country experienced an increase

in this parameter. Minimum temperature (Bio 06) decreased monotonically in 26% of the country and consistently increased in 12% of the country. For precipitation (Bio 12), consis- tent increases and reductions were observed in 18% and 5% of the country, respectively, mainly in the Altiplano Norte, Sierra Madre Occidental, Sonorense, Baja California, and Pete´n provinces. Positive changes in Bio 06 were particularly consistent in the mountainous regions Table 2. Average temperature changes from t1-1940 to t2-1970 and from t2-1970 to t3-2000. Annual mean temperature changes (˚C) Region Statistics t1-1940� to t2-1970 t2-1970 to t3-2000 Mexico Mean -0.32 0.22 Median -0.28 0.30 10th percentile -1.10 -0.46 90th percentile 0.45 0.69 Nearctic Mean -0.43 0.30 Median -0.33 0.28 10th percentile -1.34 -0.26 90th percentile 0.49 0.89 Transition Zone Mean -0.39 0.03 Median -0.42 0.06 10th percentile -1.33 -0.45 90th percentile 0.56 0.49 Neotropical Mean -0.64 0.06 Median -0.70 0.08 10th percentile -1.45 -0.43 90th percentile 0.26 0.48 Summary statistics were computed for Mexico and for three biogeographic regions. To calculate the statistics, we used the grid points (at 30-arc second spatial resolution) that fell within a particular province. Stations used were not corrected for inhomogeneities and densities vary between periods, especially between t1-1940 to t2-1970. �This period is longer than t2-1970 and t3-2000. https://doi.org/10.1371/journal.pone.0209808.t002 PLOS ONE Historical climate change in Mexico

PLOS ONE | https://doi.org/10.1371/journal.pone.0209808 July 16, 2020 9 / 19

of the country, at the highest peaks of the Eje Neovolca´nico province: Popocate´petl, Iztaccı´- huatl and Nevado de Toluca volcanoes. Negative changes for the minimum temperature occurred mainly in Yucata´n, Los Altos de Chiapas and Pete´n provinces (Fig 5). Finally, we identified a consistent increase in Bio 05 in the California, Baja California, Sonorense, Del Cabo, and Altiplano Norte provinces throughout the century (Fig 5). Discussion Climate change is currently the environmental issue of greatest concern worldwide, thus extensive research efforts are underway. Surprisingly, there are very few formal evaluations of historical climate change effects on biodiversity [21]. For this purpose, historical climate inter- polations have proven useful to analyze climate variation during the recent past and its effects on biodiversity and agrobiodiversity [21,62]. Here, we developed historic climate interpolated Fig 3. Seasonal change of maximum (♦) and minimum (bars) temperature from t1-1940 to t2-1970 (Δt1-t2) and from t2-1970 to t3-2000 (Δt2-t3) in the biogeographic provinces of Mexico. Negative values indicate a decrease in temperature from the previous period and positive values indicate an increase. DB: Depresio´n del Balsas, SMOr: Sierra Madre Oriental, DC: Del Cabo, C: California, Sse: Sonorense, SMOcc: Sierra Madre Occidental; AN: Altiplano Norte, T: Tamaulipeca, BC: Baja California, CP: Costa del Pacı´fico, AS: Altiplano Sur, EN: Eje Neovolca´nico, GM: Golfo de Me´xico, O: Oaxaca, Y: Yucata´n, P: Pete´n, AC: Altos de Chiapas, Ssc: Soconusco, and

SMS: Sierra Madre del Sur. Stations used were not corrected for inhomogeneities and densities vary between periods, especially between t1-1940 to t2-1970. Biogeographic provinces reprinted from [57] under a CC BY license, with permission from [CONABIO], original copyright [1997]. https://doi.org/10.1371/journal.pone.0209808.g003 PLOS ONE Historical climate change in Mexico PLOS ONE | https://doi.org/10.1371/journal.pone.0209808 July 16, 2020 10 / 19

surfaces to evaluate to what extent Mexico and its biogeographic regions and provinces have been exposed to climate change during the 20th and early 21st centuries. Significance for ecological studies Mexico has been projected to maintain medium to high climate stability [36], but, as expected, historically the magnitudes of change have not been geographically uniform across Mexico. Understanding the history of climate change at the regional level can provide useful informa- tion for analyzing how climate change impacts biodiversity [62]. Our results show that, on average, minimum and maximum temperatures decreased and precipitation increased coun- trywide from t1-1940 to t2-1970, whereas from t2-1970 to t3-2000, average minimum and maxi- mum temperatures rose and precipitation decreased. Similar results were obtained for temperature by Pavia et al. [19] and Englehart and Douglas [18] from analysis of weather sta- tions’ data. Furthermore, these authors found that warming has not been consistent Fig 4. Seasonal precipitation changes in % (♦) and millimeters (bars) from t1-1940 to t2-1970 (Δt1-t2) and from t2-1970 to t3-2000 (Δt2-t3) in the biogeographic provinces of Mexico. Negative values indicate a decrease in precipitation from the previous period and positive values indicate an increase. DB: Depresio´n del Balsas, SMOr: Sierra Madre Oriental, DC: Del Cabo, C: California, Sse: Sonorense, SMOcc: Sierra Madre Occidental; AN: Altiplano Norte, T: Tamaulipeca, BC: Baja California,

CP: Costa del Pacı´fico, AS: Altiplano Sur, EN: Eje Neovolca´nico, GM: Golfo de Me´xico, O: Oaxaca, Y: Yucata´n, P: Pete´n, AC: Altos de Chiapas, Ssc: Soconusco, and SMS: Sierra Madre del Sur. Stations used were not corrected for inhomogeneities and densities vary between periods, especially between t1-1940 to t2-1970. Biogeographic provinces reprinted from [57] under a CC BY license, with permission from [CONABIO], original copyright [1997]. https://doi.org/10.1371/journal.pone.0209808.g004 PLOS ONE Historical climate change in Mexico PLOS ONE | https://doi.org/10.1371/journal.pone.0209808 July 16, 2020 11 / 19

throughout the country during the 20th century, as we also observed. They reported that cool- ing occurred mainly in the central and southern provinces of the country in the mid-century (1940–1969). In turn, we observed a decrease of the maximum temperature at the Sierra Madre Oriental and Soconusco provinces from t1-1940 to t2-1970, and an increase since t2- 1970 in the Sierra Madre Oriental, Golfo de Me´xico, Altos de Chiapas, and Yucata´n provinces. Our analyses demonstrate that the Neotropical region (Costa del Pacı´fico, Golfo de Me´xico, Depresio´n del Balsas, Oaxaca, Altos de Chiapas, Soconusco, Yucata´n, and Pete´n provinces) has exhibited a more pronounced decline in precipitation, probably as a consequence of an increase in the frequency and intensity of El Niño events in the last two decades [63]. Con- versely, precipitation and temperature [19] have shown an increase in some provinces of the north, namely, California, Sonorense, Altliplano Norte, and Tamaulipeca. However, precipita- tion is counterbalanced by evapotranspiration, and the combined effect of an increase in pre- cipitation and temperature can cause a larger vapor pressure deficit and evaporation, reducing Fig 5. Climate change in Mexico from t1-1910 to t3-2000. Negative and positive areas represent downwards and upwards changes, respectively, in annual precipitation (Bio 12), maximum temperature of the warmest month (Bio 05) and minimum temperature of the coldest month (Bio 06). Blue areas indicate negative signs (decrease)

for the difference between t1-1940–t2-1970 and t2-1970–t3-2000. Red areas indicate positive signs (increase) for the difference between t1-1940–t2-1970 and t2-1970–t3- 2000. White areas indicate that the sign between t1-1940–t2-1970 and t2-1970–t3-2000 were inconsistent. Biogeographic provinces reprinted from [57] under a CC BY license, with permission from [CONABIO], original copyright [1997]. https://doi.org/10.1371/journal.pone.0209808.g005 PLOS ONE Historical climate change in Mexico PLOS ONE | https://doi.org/10.1371/journal.pone.0209808 July 16, 2020 12 / 19

water availability [62], as has been observed in Tamaulipeca, Baja California and Sonorense provinces [64,65]. Arid and semiarid provinces are highly dependent on water availability, reg- ulating net ecosystem productivity [62,66] and agriculture [36,64]. Human population and agriculture have increased in these provinces in the last three decades, exceeding water avail- ability in some areas [63,67]. There are large spatial variations in the projected changes of temperature and precipitation in most of the Altiplano Sur and Eje Neovolca´nico provinces, although within the latter we identified an increase in minimum temperature in several high-elevation volcanoes. The rapid retreat of glaciers during the 20th century confirms this observation [68]. The confluence of flora and fauna from the Nearctic and Neotropical regions makes the Transition Zone particu- larly rich and unique [60]. Many species are endemic to this region and occupy naturally nar- row distributional ranges. This condition coupled with the fact that the Transition Zone concentrates the highest human density and largest urban nuclei of the country, its biodiver- sity is highly vulnerable to the synergistic effects of multiple stressors, including climate change [10]. Furthermore, survival of many species depends on their capacity to keep pace with cli- mate; however, their ability to respond to climatic changes via range shifts are limited due to human-induced obstacles, making them highly vulnerable to the current warming event [36,62].

There are several factors other than global warming that also influence climatic variations through time, for instance, the Pacific Decadal Oscillation, the Atlantic Multidecadal Oscilla- tion and the El Niño Southern Oscillation. At a regional scale, the direct attribution of anthro- pogenic climate change to the observed patterns is complicated; nonetheless, our results are in line with those obtained by the IPCC, showing increase in temperature since mid-20th century, Table 3. Average precipitation changes from t1-1940 to t2-1970 and from t2-1970 to t3-2000. Anual Precipitation (mm) Region Statistics t1-1940� to t2-1970 t2-1970 to t3-2000 Mexico Mean 123 (16) -13 (-2) Median 48 (8) 1 (0) 10th percentile 28 (11) 23 (8) 90th percentile 26 (2) 164 (10) Nearctic Mean 3 (1) 18 (4) Median 35 (10) 18 (5) 10th percentile 3 (1) 18 (13) 90th percentile 41 (7) -1 (0) Transition zone Mean 141 (16) -12 (-1) Median 121 (15) -13 (-2) 10th percentile 40 (8) 10 (2) 90th percentile 262 (21) -36 (-3) Neotropical region Mean 289 (21) -62 (-5) Median 122 (10) -8 (-1) 10th percentile 139 (16) -28 (-4) 90th percentile 740 (33) -160 (-8) Summary statistics are reported as millimeters and percentages (in parenthesis) and were computed for Mexico and for three biogeographic regions. To calculate the statistics, we used the 30-arc gridpoints that fell within a particular province. Stations used were not corrected for inhomogeneities and densities vary between periods, especially between t1-1940 to t2-1970. �This period is longer than t2-1970 and t3-2000.

https://doi.org/10.1371/journal.pone.0209808.t003 PLOS ONE Historical climate change in Mexico PLOS ONE | https://doi.org/10.1371/journal.pone.0209808 July 16, 2020 13 / 19

with greater changes in the north of the country [9]. We also identified a general increase in precipitation throughout the century, although its magnitude decreased between t2-1970 and t3-2000. This change might be related to the fact that a negative Atlantic Multidecadal Oscilla- tion and a Pacific Decadal Oscillation were observed during the last period [9]. Pervasive climatic change leads to changes in species composition, as has been observed for birds in Mexico [28]. Understanding recent climate and ecological changes are useful to focus research efforts and conservation actions; therefore, the climate changes detected in this study may be useful to establish priorities for conservation. For instance, the Transition Zone includes the Sierra Madre Occidental, Sierra Madre Oriental and Eje Neovolca´nico provinces, which harbor diverse types of temperate vegetation, such as oak, pine and cloud forests; all of them, but particularly the latter, are highly vulnerable to climate change [68,69]. The Sierra Madre Oriental harbors the richest coniferous forests in the world [70]. Our results show that this province has been warming since mid-20th century and is projected to be affected by cli- mate change in the future [69,71], thus conservation efforts should be prioritized in this region [70]. Changes in climate are now occurring simultaneously with other environmental disrup- tions. It will not be possible to fully understand biodiversity responses to climate change with- out considering the interactions with other components of environmental change [72]. Species

responses will depend on their exposure and sensitivity to other human-induced pressures, their inherent capacity to adapt to new conditions, the magnitude and speed of environmental changes, and time lags in their responses [1,10,62]. This work aims to serve as a baseline for improving our knowledge regarding historic climate change in Mexico in a spatially-explicit fashion, to identify specific areas where climate change is occurring and to identify its direc- tion and magnitude. In sum, climate change is occurring unevenly in Mexico. Provinces in the Nearctic region showed a higher and more consistent warming over time than provinces from the Transition Zone and Neotropical region. Precipitation has also generally increased more consistently in northern provinces over the whole period, whereas it tended to decrease in most of the Neo- tropical region and Transition Zone between t2-1970 and t3-2000. Nonetheless, it is important to bear in mind that the first period (t1-1940) holds higher uncertainty than the other two given the lower number of weather stations used to produce the climatic surfaces. Finally, this information can be used to improve regional projections of future climate impacts on biodi- versity, which would provide scientists and authorities with more reliable data and informa- tion for making better decisions in the face of climate change. Limitation of the analysis We have generated three sets of climatic surfaces for Mexico with error estimates based on withheld data across the country [42] comparable to errors obtained in two other climate sur-

face datasets developed for the country for different time frames [15,17]. However, it is impor- tant to mention the main shortcomings in this type of data and analysis for a better use of them. First, it has been demonstrated that long-term station records that are used to interpo- late climate surfaces generally suffer from incompleteness and inconsistencies through time caused by equipment failure, replacement or station relocation [40,73], limiting their use for calculating multi-decadal climate trends [72,74]. This is the case for the data that we used for generating these climatic surfaces. Our results of the Pettitt test indicate that a considerable proportion of the stations presented inhomogeneities along the time-series, especially for min- imum temperature (see S1 File). Furthermore, the breakpoint years occurred mostly during the second half of t2-1970 and the first half of t3-2000 (S1, S2, S3 Figs in S1 File). To know how PLOS ONE Historical climate change in Mexico PLOS ONE | https://doi.org/10.1371/journal.pone.0209808 July 16, 2020 14 / 19

much of these inhomogeneities are due to climatic or non-climatic factors requires further investigation. Regardless, it is likely that non-climatic variations may have influenced to some extent the observed changes between these periods. Consequently, these non-climatic varia- tions also increase uncertainty in ecological studies [75]; hence, additional sources, such as remote sensing-based climate products can be used to reduce uncertainty, e.g., [76]. Second, climate surfaces for t1-1940 were interpolated for a longer period but with fewer weather stations, thus are more uncertain than the climate surfaces for the other two periods [17,40]. We advise caution when using the surfaces from this period. Finally, there is a geo- graphical bias in the location of weather stations in the three periods, since density of stations was higher in or nearby human settlements and lower in remote locations, some of which can be important for biodiversity analyses [74,77]. This geographical bias has implications in the spatial structure of error and uncertainty, particularly in areas with low station density [75], where scarcity of stations impedes to capture the effect of terrain complexity in climate vari- ability [78]. Considering all these limitations in the data used for the construction of the cli- matic surfaces, the use of them and the interpretation of the results should be made with caution. Supporting information S1 File. ANUSPLIN statistics and additional spatial distribution of climate data. This file

contains the http of a kmz with the distribution of weather stations across Mexico and ANUS- PLIN statistics. We also included results of the Pettitt test for change-point detection in precip- itation, maximum and minimum temperature stations, and three figures that describe the frequency distribution of the detected breakpoints years. Finally, we share a map showing extremes of warming. (DOCX) Acknowledgments We thank the Servicio Meteorolo´gico Nacional (SMN) for providing data from weather sta- tions. We also thank the Comisio´n Nacional para el Conocimiento y Uso de la Biodiversidad (CONABIO) for serve as the repository of the climate surfaces developed in this study. Caro- lina Ureta thanks Catedras-Conacyt for its support. Comments of three anonymous referees are gratefully acknowledged. Author Contributions Conceptualization: Angela P. Cuervo-Robayo. Data curation: Angela P. Cuervo-Robayo, Miguel A. Go´mez-Albores, Anny K. Meneses- Mosquera. Formal analysis: Angela P. Cuervo-Robayo, Enrique Martı´nez-Meyer. Funding acquisition: Angela P. Cuervo-Robayo, Oswaldo Te´llez-Valde´s, Enrique Martı´nez- Meyer. Investigation: Angela P. Cuervo-Robayo, Carolina Ureta, Miguel A. Go´mez-Albores, Oswaldo Te´llez-Valde´s, Enrique Martı´nez-Meyer. Methodology: Angela P. Cuervo-Robayo, Miguel A. Go´mez-Albores, Anny K. Meneses-Mos- quera, Enrique Martı´nez-Meyer. Project administration: Angela P. Cuervo-Robayo, Enrique Martı´nez-Meyer. PLOS ONE Historical climate change in Mexico PLOS ONE | https://doi.org/10.1371/journal.pone.0209808 July 16, 2020 15 / 19

Resources: Angela P. Cuervo-Robayo, Oswaldo Te´llez-Valde´s, Enrique Martı´nez-Meyer. Software: Oswaldo Te´llez-Valde´s. Supervision: Angela P. Cuervo-Robayo, Oswaldo Te´llez-Valde´s, Enrique Martı´nez-Meyer. Validation: Miguel A. Go´mez-Albores, Enrique Martı´nez-Meyer. Visualization: Angela P. Cuervo-Robayo. Writing – original draft: Angela P. Cuervo-Robayo, Carolina Ureta, Enrique Martı´nez- Meyer. Writing – review & editing: Angela P. Cuervo-Robayo, Carolina Ureta, Enrique Martı´nez- Meyer. References 1. Sala OE, Chapin FS, Armesto JJ, Berlow E, Bloomfield J, et al. (2000) Global biodiversity scenarios for the year 2100. Science 287: 1770–1774. https://doi.org/10.1126/science.287.5459.1770 PMID: 10710299 2. IPCC (2014) Climate change 2014: Impacts, adaptation and vulnerability. Summary for policy makers. https://doi.org/10.1016/j.renene.2009.11.012 3. Thuiller W (2004) Patterns and uncertainties of species’ range shifts under climate change. Global change biology 10: 2020–2027. 4. Thuiller W, Albert C, Arau´jo MB, Berry PM, Cabeza M, et al. (2008) Predicting global change impacts on plant species´ distributions: Future challenges. Perspectives in Plant Ecology, Evolution and Systemat- ics: 1–16. 5. Williams J, Jackson S, Kutzbach J (2007) Projected distributions of novel and disappearing climates by 2100 AD. PNAS 104: 5738–5742. https://doi.org/10.1073/pnas.0606292104 PMID: 17389402 6. Ureta C, Martorell C, Hortal J, Fornoni J (2012) Assessing extinction risks under the combined effects of