Introduction
The relationship between biodiversity and climate has become a concern for the scientific community, due to the negative effects of climate change on vegetation (Intergovernmental Panel on Climate Change [IPCC], 2014). Increases in temperature and changes in precipitation affect species in various ways, and these changes would modify their distribution, diversity, and abundance in the future (Burgmer, Hillerbrand, & Pfenninger, 2007; Saenz-Romero, Rehfeldt, Ortega-Rodriguez, Marin-Togo, & Madrigal-Sanchez, 2015).
Climate change scenarios for Mexico in the 21st century predict an increase in average annual temperature of 0.5 to 2 °C by 2030, 2.3 °C by 2060 and up to 3.7 °C by 2090 (Sáenz-Romero et al., 2010). These scenarios predict an increase in the area occupied by dry forests of 7.4 % from 2025 (Villers-Ruiz & Trejo-Vázquez, 2000), while the tropical zones of the southeast such as the Selva el Ocote Biosphere Reserve would show an increase in average annual temperature between 0.41 and 0.83 °C and a decrease in precipitation of 35 to 71 mm by 2030 (Manzanilla-Quiñones & Aguirre-Calderón, 2017).
The Meliaceae family is a group of plants of great importance in the tropics, of which red cedar (Cedrela odorata L.) and mahogany (Swietenia macrophylla King) stand out as valuable woods (Mendizábal-Hernández, Alba-Landa, & Suárez-Dorantes, 2009; Pennington & Sarukhán, 2005; Secretaría de Medio Ambiente y Recursos Naturales [SEMARNAT], 2016). Cedrela odorata is naturally distributed in an altitudinal range from 0 to 1 200 m from Central America in Belize and Panama to South America in Venezuela, Colombia and part of the Amazon in Brazil and Peru (Salazar, Soihet, & Méndez, 2000). In Mexico, the species is located in southern Tamaulipas, southeast of San Luis Potosi, Sinaloa, Guerrero, the Central Depression, the coasts of Chiapas and the Yucatan Peninsula (Gómez, Monterroso, & Tinoco, 2007; Pennington & Sarukhán, 2005; Romo-Lozano, Vargas-Hernández, López-Upton, & Ávila-Ángulo, 2017).
Timber production in Mexico has fluctuated from 6.7 to 7 million cubic meters in logs (m3r) annually, where the main producing states are Durango (35.1 %), Chihuahua (13.2 %), Veracruz (7.8 %), Michoacán (6.7 %) and Oaxaca (5.9 %), which contributed 68.6 % of total production. Of this production, only 0.5 % corresponds to precious woods such as mahogany and red cedar, with Campeche (18 373 m3r), Chiapas (836 m3r), Jalisco (2 784 m3r), Veracruz (4 517 m3r) and Quintana Roo (4 087 m3r) being the states with the highest timber production and which generated revenues of MXN 919 406 up to MXN 18 373 000 (SEMARNAT, 2016, 2017).
Currently, after mahogany, red cedar is the most important timber species in the forest industry in Mexico (Pennington & Sarukhán, 2005; SEMARNAT, 2016). Red cedar has been a species of high economic value, since its wood is considered fine and aesthetically precious, for the manufacture of fine luxury furniture, veneer, and cabinetmaking. Due to the excellent quality of the wood, red cedar has been selectively cut in tropical forests for multiple uses, which has resulted in the larger diameter and straight shaft trees being scarce in today's natural populations.
Due to poor forest management of C. odorata, the species has fragmented its habitat and reduced its natural populations significantly (Rodríguez et al., 2003). This situation has caused the red cedar to be included in the Official Mexican Standard NOM-059 as a species subject to special protection (SEMARNAT, 2010) and in Appendix III of the Convention on International Trade in Endangered Species of Wild Fauna and Flora (CITES, 2011) so from February 2011, in order to make any kind of use, a series of legal procedures must be followed. These involve obtaining permits in which population studies, sampling, and specifications for restoration, repopulation, and reintroduction programs, projects, or activities must be attached (Ruiz-Jiménez, De los Santos-Posadas, Parraguirre-Lezama, & Saavedra-Millán, 2018).
Researchers have used Geographic Information Systems (GIS) and ecological niche modeling algorithms to generate habitat quality suitability maps to identify geographic areas with potential for restoration or conservation activities for species of ecological and economic interest (Garza-López et al., 2016; Manzanilla et al., 2019). The use of programs for ecological niche modeling such as GARP (genetic algorithm), BIOCLIM (climate envelope algorithm) and MaxEnt (maximum entropy algorithm) have allowed the prediction of species distribution areas with good statistical probability (Franklin, 2010; Peterson, 2011; Phillips & Dudík, 2008). These algorithms allow the identification of areas with favorable biotic and abiotic conditions, both for the conservation of a taxon and for the establishment of areas suitable for species with high ecological and economic value such as C. odorata (Morales, 2012; Perosa et al., 2014; Soberón, Osorio-Olvera, & Peterson, 2017). The delimitation of these areas is based on the environmental zoning carried out on the basis of an exploratory analysis of eco-geographical variables (latitude, longitude, altitude, slope, climate and soil) (Castellanos-Acuña et al., 2018).
In order to have elements to make decisions for the management and conservation of C. odorata in the Yucatan Peninsula, the objective of this work was to delimit the potential current and future distribution (2030), as well as to identify and propose areas with environmental and morphological characteristics suitable for the conservation and production of seeds.
Materials and methods
The study area is located in the Yucatan Peninsula, situated in southeast Mexico between 22° 31' 43'' and 17° 48' 51'' LN and 92° 20' 11'' and 86° 42' 36'' LO; it covers the states of Campeche, Quintana Roo and Yucatan, and has an area of 139 840 km2 (Figure 1). The relief is flat with an average altitude of 50 m and only in the center-south there are elevations of up to 350 m (Instituto Nacional de Estadística y Geografía [INEGI], 2015). The dominant climate is Aw1(x') warm sub-humid with an average annual temperature of more than 22 °C and with annual accumulated precipitation intervals of 700 to 1 600 mm (Cuervo-Robayo et al., 2014; García 1998).
Data collection
Two hundred twenty-seven records of C. odorata were obtained from the databases of the National Forest and Soil Inventory (INFyS) 2004-2009 (CONAFOR, 2009), the book "Árboles tropicales de México" (Pennington & Sarukhán, 2005) and the Global Biodiversity Information Facility (GBIF, 2017) (Figure 1). The databases were purified leaving a record per cell of 1 km2 (30 arc seconds), their correct position was checked eliminating duplicated records, badly georeferenced, as well as those located in urban areas and water bodies. These procedures were carried out on the Niche ToolBox platform of the Comisión Nacional para el Conocimiento y Uso de la Biodiversidad (CONABIO) (Osorio-Olvera, Vijay, Narayani, Soberón, & Falconi, 2016) and in the ArcMap program version 10.3 (Environmental Scientific Research Institute [ESRI], 2014). After the debugging process, 216 valid records were obtained for the calibration of the models. For a better evaluation of the data, 50 models with different modeling criteria were generated (Table 1).
Models | Source | Data |
---|---|---|
10 | National Forest and Soil Inventory (INFyS) | 166 |
10 | Global Biodiversity Information Facility (GBIF) and Árboles tropicales de México | 50 |
10 | GBIF, Árboles tropicales de México and INFyS | 216 |
10 | GBIF, Árboles tropicales de México and INFyS 25 % | 91 |
10 | GBIF, Árboles tropicales de México and INFyS 50 % | 133 |
Current and future bioclimatic variables
The 19 bioclimatic variables (http://visualcrow.com/descarga-mexico.html) generated for Mexico for the period 1910-2009 were downloaded (Cuervo-Robayo et al., 2014). For the evaluation of the potential distribution in the future (2030), the bioclimatic variables of the general circulation models GFDL_CM3 and MIROC_ESM (http://www.ccafs-climate.org/data_spatial_downscaling/) projected to the near future (2030), which simulate radioactive forcing trajectories, in greenhouse gas emissions, similar to the current ones (RCP 4.5 W·m-2) were used. These bioclimatic variables have a resolution of one km2 (30 arc seconds); in addition, the slope (%), orientation (°) and altitude (m) variables obtained from the digital elevation model with a resolution of 90 m in GeoTIFF format were included. The type of vegetation was obtained from the series VI "Land use and vegetation" of INEGI (2016) with a scale of 250 m. The edaphology (soil type) of the series I "Edaphological vector data set", with a scale of 250 m, was obtained from CONABIO (1995). These last two variables were downloaded in vector format (shapefile). All variables were standardized and transformed into ASCII format at a scale of 1 km2 (30 arc seconds) with the help of ArcMap version 10.3 (ESRI, 2014).
Delimitation of the area M
The variables were adjusted to the size of the spatial modeling area M (Yucatan Peninsula), which has been described as the space where a species is or is assumed to be, based on biological knowledge and its ability to disperse (Martínez-Méndez, Aguirre-Planter, Eguiarte, & Jaramillo-Correa, 2016; Soberón & Peterson, 2005). From area M, environmental variables were cut and adjusted to the same 1 km2 pixel size (30 arc seconds) with the help of ArcMap version 10.3 (ESRI, 2014).
Modeling of current and future potential distribution
The models were generated in the MaxEnt algorithm version 3.3.3k, which uses an exploratory technique and identifies sites with similar values according to the principle of maximum entropy (Phillips, Anderson, & Schaphire, 2006). This algorithm is one of the most efficient for modeling species distribution using presence data only (Elith et al., 2011). In order to increase the reliability of the MaxEnt results, 50 distribution models were generated and tested under 10 modeling criteria (Table 2) using different threshold application rules and a maximum number of replications of 1 000, with a convergence limit of 0.00001 (Manzanilla et al., 2019; Phillips et al., 2006). This made it possible to determine which model best fits and explains the current distribution of the species.
Model | Internal replication | Model | Internal replication | Application Rule Threshold | Replicas |
---|---|---|---|---|---|
M1 | Bootstrap | M6 | Cross validation | Equal training sensitivity and specificity | 1 000 |
M2 | Bootstrap | M7 | Cross validation | Maximum training sensitivity plus specificity | 1 000 |
M3 | Bootstrap | M8 | Cross validation | Equal test sensitivity and specificity | 1 000 |
M4 | Bootstrap | M9 | Cross validation | Maximum test sensitivity plus specificity | 1 000 |
M5 | Bootstrap | M10 | Cross validation | Without threshold rule application | 500 |
The models were calibrated using 75 % of the species presence data in csv format, selected at random to train the models, and the remaining 25 % for the validation test (Alba-Sánchez et al., 2010; Martínez-Méndez et al., 2016; Phillips et al., 2006). The relative contribution (%) of each variable in the models was estimated using the jackknife test (Shcheglovitova & Anderson, 2013). The output was of the logistic type, which indicates an index of similarity (probabilistic) of suitable conditions for the species with values from 0 to 1, where values close to 1 indicate suitable conditions for the development and presence of the species (Coitiño, Montenegro, Fallabrino, González, & Hernández, 2013; Phillips et al., 2006).
To generate the future model (2030), the parameters of the current distribution model with the best fit and statistical performance were transferred to the MaxEnt program (Morrone & Escalante, 2016).
Model validation
The models were evaluated by means of the Area Under the Curve (AUC) test, which is obtained directly from the Receiver Operating Characteristic (ROC) technique (Coitiño et al., 2013; Phillips et al., 2006); however, this validation has been questioned because it does not consider true absence data (Peterson, Papes, & Soberón, 2008). Therefore, it was necessary to perform a more robust analysis using a partial ROC test in the Tool for Partial-ROC version 1.0 program (Narayani, 2008), which was designed to counteract the shortcomings of the AUC (Peterson et al., 2008).
Partial ROC analyses were performed on 50 % of the records with 95 % reliability, with 1 000 iterations by bootstrapping and fixing a 5 % error of omission. The test generates values from 1 to 2; the closer it is to 2, the more reliable the model is (Martinez-Mendez et al., 2016; Peterson et al., 2008). Also, a Z test was performed, where the calculated Z value of the AUC proportions of partial ROC must be greater than the Z value of tables (Z > 2.3 = 99 %), which provided greater support to the models (P < 0.01) (Martínez-Méndez et al., 2016; Monterrubio-Rico et al., 2016).
Potential areas of distribution
The values of the logistic output of the model with the best statistical fit were reclassified into three categories of habitat quality with equal intervals (low, medium and high), through the reclass module of the ArcMap 10.3 program (ESRI, 2014). The value of the high habitat quality category was used to transform the models from continuous to binary (fit and unfit). Based on the reclassification, the area occupied by the potential distribution areas of C. odorata was estimated for the current and future periods.
Conservation and seed production areas
The conservation and seed production areas were selected based on the identification and delimitation of the areas with the highest habitat quality in the current and future periods. An additional criterion for determining locations with the possibility of serving as a climatic refuge for the conservation of the species was the location of natural protected areas (ANPs) (SEMARNAT & Comisión Nacional de Áreas Naturales Protegidas [CONANP], 2017), whose boundaries were used and combined with the areas of distribution determined for both periods.
The selection of suitable areas for seed production was based on the evaluation of stand quality (Huang, Titus, & Wiens, 1992), which reflects the productive capacity of the sites through the dominant height of the trees and the diameter-height ratio. Based on the proposals of Huang et al. (1992) and Yuancai and Parresol (2001), it has been observed that sites with larger diameter trees tend to have a greater dominance in total height. The suitable areas for seed production were identified through the INFyS maximum diameter data in raster format, from the data obtained from the 400 m2 rectangular sampling clusters, designed for forests (CONAFOR, 2009). This file was transformed to vector format in the ArcMap 10.3 program (ESRI, 2014); subsequently, the polygons of the current and future distribution were used to identify the areas where the species is found. Finally, the areas for seed production were chosen based on the reclassification of the values of the diameter category ≥50 cm; according to Huang et al. (1992) and Yuancai and Parresol (2001), this value would correspond to a tree of 20 m in height, which would be in the reproductive stage.
Results and discussion
The results of the current AUC models recorded values of 0.893 for training data and 0.805 for validation, while for the future models, GFDL_CM3 and MIROC_ESM, the values were 0.962 to 0.957 for training data and 0.731 to 0.751 for validation data. The models generated for the current and future periods had a good predictive performance as a function of the AUC, better than a randomly generated model. According to Araújo and Guisan (2006), Peterson (2011), and Miranda, Geada, and Sotolongo (2016), AUC values close to or equal to 0.50 show equal or worse performance than a randomly generated model, while 0.60-0.70 are considered bad, 0.71-0.80 regular, 0.81-0.90 good, and >0.90 are considered excellent.
The models with the best statistical performance were the 10 models generated from the 50 records of the GBIF and "Árboles tropicales de México " databases. Of these, the one that presented the best fit and statistical performance was model 7 (M7VC); this result was corroborated by the partial ROC and Z tests, which were statistically significant (P < 0.01) (Table 3).
Period | Models | Average radii partial ROC | Standard deviation | Z Test |
---|---|---|---|---|
Current | M1B | 1.447 | 0.106 | P < 0.01 |
M2B | 1.550 | 0.110 | P < 0.01 | |
M3B | 1.552 | 0.112 | P < 0.01 | |
M4B | 1.549 | 0.113 | P < 0.01 | |
M5B | 1.531 | 0.112 | P < 0.01 | |
M6CV | 1.533 | 0.112 | P < 0.01 | |
M7CV | 1.554 | 0.106 | P < 0.01 | |
M8CV | 1.546 | 0.112 | P < 0.01 | |
M9CV | 1.538 | 0.111 | P < 0.01 | |
M10CV | 1.533 | 0.114 | P < 0.01 | |
Future | GFDL_CM3 | 1.111 | 0.119 | P < 0.01 |
MIROC_ESM | 1.111 | 0.087 | P < 0.01 |
Relevant variables
The jackknife test indicated that the variables with the greatest contribution in the generation of the current distribution models of C. odorata were vegetation with 34.7 %, followed by the precipitation of the wettest month with 14.6 %, edaphology with 8.8 %, average temperature of the coldest quarter with 8.6 % and slope with 7 %. In the analysis of the potential future distribution of the GFDL_CM3 and MIROC_ESM models, the most relevant environmental variables were vegetation with 84 % and 83.7 %, mean temperature of the coldest quarter with 6.5 % and 9 %, and edaphology with 2.8 % and 2.5 %, respectively (Figure 2). These variables determine the current and future (2030) ecological niche of the species in the Yucatan Peninsula.
According to the contribution of the environmental variables, C. odorata requires Leptosol type soils with precipitation in the wettest month of 216 mm and oscillation in the average temperature of the coldest quarter of 22 to 24 °C. The species has a preference for ecosystems with secondary tree and shrub vegetation of medium deciduous, sub-deciduous, evergreen and sub-perennial forests. These results are consistent with the environmental requirements of a tropical forest species, which have been described in the literature, providing a degree of confidence that the species is found in areas congruent with sites with habitat conditions similar to those described by Rodríguez (2003), and Pennington and Sarukhán (2005). The results of the environmental variables already mentioned are similar to those reported by Gómez et al. (2007) and Garza-López et al. (2016), who mention that tropical species such as red cedar and mahogany tend to develop better during the dry seasons (no longer than four months), in sites with average annual temperatures of 24 to 28 °C, maximum temperature of the hottest month of 24 to 32 °C and minimum temperature of the coldest month of 11 to 22 °C. Andrade and Solis (2004) point out that C. odorata is both a secondary and a primary species, because it is usually a pioneer in areas of secondary vegetation in medium and high forests, and also as a frequent element in the upper stratum of mature forests. Salazar et al. (2000) mention that C. odorata requires annual rainfall of 1 200 to 2 500 mm and an altitudinal gradient of 0 to 1 500 m for its optimal development and growth, conditions that are located within the densest forested areas towards the southeast of the Yucatan peninsula. The most relevant variables for the future (2030), in both evaluated models, were vegetation, soil and average temperature of the coldest quarter, which remained as the main determining variables in the distribution of the species. This indicates that changes or alterations in the vegetation associated with C. odorata would have a great impact on the species' future distribution areas.
Current and future potential distribution
The results indicate a current potential distribution of C. odorata of 405 065 ha, equivalent to 2.9 % of the study area, where this area meets the highest environmental conditions of habitat quality, suitable for the species. The average habitat quality represents an area of 19.2 % and the low habitat quality represents 77.9 %, which represents marginal areas of almost no suitability. The results of this classification allowed us to generate a map of the current habitat quality available for C. odorata in the Yucatan Peninsula (Figure 3).
On the other hand, the results of the future models GFDL_CM3 and MIROC_ESM, projected to 2030, estimate a potential suitable area of 277 037 and 226 138 ha, equivalent to 1.6 and 2 % of the study area, respectively. These results represent a significant reduction between 31 and 44.8 % of the suitable habitat for the species with respect to its estimated current potential distribution (Figure 4).
Although climate change scenarios predict a significant reduction in the species' ecological niche, they also indicate that natural populations of C. odorata would find a climatic refuge into the denser parts of the mid and upper forests, located in the southeast of Campeche State, specifically in the Calakmul Biosphere Reserve. The results on ecological niche decline are consistent with the projections of Estrada-Contreras, Equihua, Laborde, Meyer, and Sánchez-Velásquez (2016), who modeled the current and future distribution of C. odorata under climate change scenarios in Mexico, and reported a decline in its range of 3.8 % under the CGCM3 B1 scenario by 2030 from its current distribution. Garza-López et al. (2016) estimated a 60 % reduction in S. macrophylla climate habitat in the Yucatan Peninsula by 2030. The results of these authors and those obtained in the present study suggest a medium-high vulnerability for the tropical species analyzed in southeastern Mexico. The combination of high temperatures and decreased precipitation would be the main cause for the reduction of suitable conditions for C. odorata in much of the peninsula; such a combination would make the species vulnerable to forest fires and attack by pests and diseases.
Villers-Ruiz et al. (2000) point out that the tropical forests of southeast Mexico could be converted into dry forests. In turn, Manzanilla-Quiñones and Aguirre-Calderón (2017) predict an increase in average annual temperature of 0.41 to 0.83 °C and a decrease in annual precipitation between 35 and 71 mm in tropical areas of southeastern Mexico by 2030, so it is very likely that similar conditions will occur in the Yucatan Peninsula for the same year.
Due to global warming, modifications in the distribution of species are already taking place, with significant changes in the taxa that inhabit mountainous regions, where they tend to move altitudinally in order to subsist (Lenoir, Gegout, Marquet, de Ruffray, & Brisse, 2008; Manzanilla-Quiñones, Aguirre-Calderón, Jiménez-Pérez, Treviño-Garza, & Yerena-Yamallel, 2019; Parolo & Rossi, 2008; Telwala, Brook, Manish, & Pandit, 2013). For this reason, it is a priority to carry out and analyze the simulations of climate change scenarios to determine which populations are at risk of disappearing due to the effects of climate change.
Conservation and seed production areas
Future models indicate that C. odorata will conserve a more suitable area of 88 603 to 64 342 ha, within the Calakmul, Los Petenes and Ría Celestún Biosphere Reserves (Figure 5). These reserves are located mainly in the municipalities of Hopelchén, Champotón and Calkiní in the state of Campeche, and in the municipality of Halachó in the state of Yucatán, where environmental conditions would serve as climatic refuges for the species in the near future (2030). Regarding the Calakmul Biosphere Reserve, results are similar to those reported by Garza-López et al. (2016, 2018) for S. macrophylla and Lysiloma latisiliquum L. Benth. These results may be useful in selecting areas for migration or assisted recolonization for conservation and restoration purposes for C. odorata in the Yucatan Peninsula, as suggested by Saenz-Romero, Rehfeldt, Crookston, Duval, and Beaulieu (2009) and the Secretariat of the Convention on Biological Diversity (CDB, 2009); furthermore, proper management of these areas will directly benefit the protection of the species.
The suitable areas for seed production were mainly located in the Calakmul and Los Petenes Biosphere Reserve, in the state of Campeche, with maximum diameters between 50 and 80 cm and heights of 10 to 30 m within an area of 341 481 ha (Figure 6). In these areas there are red cedar trees with morphological and phenotypical conditions suitable for seed production and collection and with which it is feasible to establish in situ conservation and restoration activities. Navarro-Cerrillo, Clavero, Vidaña, Quero, and Duque-Lazo (2016) mention that defining germplasm-producing areas would help conservation and safeguard the genetic material of the species of interest.
Conclusions
The ecological niche models provided valuable information on the areas with the highest habitat quality suitability of Cedrela odorata in the current and future (2030) distribution in the Yucatan peninsula. The models helped to define areas of distribution with phenotypically desirable specimens for seed production and suitable for the conservation of the species in the natural protected areas of Calakmul, Los Petenes and Ría Celestún. The most relevant variable was vegetation, which means that the slightest change in this variable would severely affect the distribution of the species in the future. Niche models of C. odorata can be used to support restoration and conservation programs for the species, as well as to propose management strategies for the establishment of seed-producing zones and in-situ protection of the species' genetic pool. This type of study is extremely relevant in the conservation of species in the face of climate change in Mexico.