EUROPEAN JOURNAL OF ENTOMOLOGY EUROPEAN JOURNAL OF ENTOMOLOGY Linking potential habitats of Odonata (Insecta) with changes in land use/land cover in Mexico

. Land use/land cover change (LULCC) is a major threat that affects the viability of insect populations worldwide yet our estimates of such effects are usually poor. We analysed how LULCC affected the distribution of 49 species of dragon ﬂ ies and damsel ﬂ ies in the south-central zone in Mexico during the period 2006–2012. For this, we mapped the potential species richness using ecological niche models in order to analyse predicted future changes and determined the effect of LULCC on the current and future habitats of Odonata. We also estimated current incidence of deforestation and projected its effect to 2050 using the Dinamica-EGO program. Having predicted the level of deforestation in the year 2050, we then compared current vs. expected species richness and the factors that determine it. First, roads and urban areas turned out to be the most important drivers of LULCC in our analysis. Second, deterioration occurred at all sites, but species richness remained high despite considerable habitat fragmentation. Third, there is likely to be a high species turnover rate (i.e. a high species richness, but composed of different species) even in areas where there are signi ﬁ cant changes in the vegetation. Our work illustrates both a resilience of Odonata to LULCC and provides a useful method for measuring the effects of LULCC on insects.

extent of the overlap of the images, which is expressed as a transition matrix that provides estimates of the probability of change (Cayuela, 2006;Flamenco-Sandoval et al., 2007). For estimating LULCCs Dinamica-EGO program (hereafter abbreviated as D-E) is widely used and provides spatially explicit environmental models (Soares-Filho et al., 2009). D-E is a series of algorithms known as functors, which execute the different spatial analyses (map algebra, simulations, cost maps, probability maps, etc.). D-E has a wide variety of applications (Rodrigues & Soares-Filho, 2018) particularly in the analysis of patterns in LULCC and their effects on vegetation (Soares-Filho et al., 2006;Piontekowski et al., 2012;Prieto-Amparán et al., 2019) and loss of biodiversity (Soares-Filho et al., 2002;Mas, 2017;Oliveira et al., 2019;Hemati et al., 2020;Godínez-Gómez et al., 2020).
Increased conversion and degradation of habitats due to fragmentation, climate change, agriculture and pollution has resulted in a high loss of biodiversity (e.g. Goncalvez-Souza et al., 2020). Moreover, evaluations of LULCC Linking potential habitats of Odonata (Insecta) with changes in land use/land cover in Mexico INTRODUCTION Changes in land use/land cover (LULCC) throughout the world are associated with human activities in terms of exploitation of natural resources (Turner & Meyer, 1994;Manson, 2006). The changes due to mining, agriculture and deforestation, are considered to be the main factors driving environmental degradation, for example, landscape fragmentation and loss of biodiversity at a global scale (Lambin et al., 2001).
One reliable method of measuring the degree to which terrestrial environments have been degraded is a LULCC analysis (Turner & Meyer, 1994;Lambin et al., 2001). This analysis provides information on the state of the habitats of many species and an understanding of the different effects of degradation (López, 2012;Farfán, 2015;Chang, 2017). One specifi c analytical approach in LULCC is to compare satellite images taken on different dates, which are processed and classifi ed into categories, classes, or types of land use/cover using both visual and automated methods (Mas & Correa Sandoval, 2000). This analysis reveals the ing the distributions of rare species with narrow habitat ranges (Renjana et al., 2022).
Thus, quantitative tools are needed to improve the understanding of the threats that LULCC poses to Odonata (e.g. Mendoza-Penagos et al., 2021;Pereira-Moura et al., 2021). A combination of Maxent and land use/land cover analysis is used in studies involving land use/land cover map generation (Mack & Waske, 2016;Amici et al., 2017). However, the use of Maxent in LULCC analysis for generating ecological niches is rare. Here we use D-E and Maxent to analyse the links between LULCC and patterns of Odonata richness in south-central Mexico. For this region, we previously generated distribution models for 94 species, using Maxent, which was used to construct a vulnerability index for Odonata (hereafter referred as IVO; Rodríguez-Tapia et al., 2020) and refl ects their tolerance of LULCC.
For the generation of niche models, as in Rodríguez-Tapia et al. (2020), we fi rstly cleaned the data for each species since it included redundant, repeated or clearly incorrect records. Secondly, environmental layers from Worldclim (http://www.worldclim.org/) were used to obtain environmental data associated with each record. Finally, we generated distribution models for each of the 94 species using Maxent (Phillips et al., 2006). The IVOs included the breadth of the range of habitats and ecological niche of each species. IVO ranges were between 1 and 3 in steps of 0.5. The value 1 represents least sensitive species while 3 indicates a species that is highly sensitive to environmental changes (Rodríguez-Tapia et al., 2020).
The objectives of this study were to: (a) generate a map of species richness for Odonata using the ecological niche model of Rodríguez-Tapia et al. (2020); (b) determine the change in habitats used by Odonata associated with LULCC (in this case deforestation) and its drivers (highways, rivers, water body) in south-central Mexico over the period of 2006-2010; and (c) predict the type of vegetation that will be present in 2050 and the geographic distribution of Odonata based on their tolerance of LULCC and their IVO (Rodríguez-Tapia et al., 2020).

Study area
The study area was delimited using the Freshwater Ecoregions of the World 2014 (FEOW, 2014). To estimate the extent of the potential distributions of the 94 species of Odonata selected we used the ecological niche model developed by Rodríguez-Tapia et al. (2020) (see Fig. S1). For this a database of 18,147 records was compiled using the Global Biodiversity Information Facility (GBIF, 2020). Prior to using the niche model, the data for each species was cleaned as initially it included redundant, repeated or clearly incorrect records. First, we used species for which there were at least 15 records. Second, we removed any second record of a species in the same pixel to reduce the weight of a pixel due to the repetition of presence records. To remove repeated records, we used a programme in R software (R Core Team, 2019) that detects any pixels with more than one record in each raster layer. The second process consisted of extracting the values of the environmental layers associated with each point. Then for each species the relationship between environmental variables was have mainly used geographical data on vertebrates to understand effects of LULCC on loss of biodiversity . However, other groups such as insects have been overlooked, despite the multitude of critical ecosystem functions provided by insects (Foottit & Adler, 2017;Cardoso et al., 2020) and their vulnerability that possibly makes them more threatened than other animal taxa (Wagner et al., 2021).
Odonata; with two suborders: Zygoptera, or damselfl ies and Anisoptera, or dragonfl ies are an emerging model group for determining the effect of habitat loss on invertebrate biodiversity because they can indicate large-scale changes in biodiversity and degradation of environmental conditions (Bried & Mazzacano, 2010;Gómez-Tolosa et al., 2020). Odonata are predatory insects, with aquatic larvae and terrestrial adults. The fact that larval and adult Odonata occupy aquatic and terrestrial ecosystems implies that they can provide valuable information of alterations in the conditions in these ecosystems (Butler & de Maynadier, 2008). Their complex life cycle and often high population densities make them an important regulatory component in freshwater ecosystems (Knight et al., 2005). However, many populations and species have disappeared (Clausnitzer et al., 2009) as freshwater habitats are among the most threatened ecosystems (Dijkstra et al., 2014;Reid et al., 2019). Recently, several researchers have used aquatic insects as bioindicators, especially in the Neotropics (Monteiro-Junior et al., 2013). Despite this, there is still a lack of information, especially on Neotropical Odonata (Gómez-Tolosa et al., 2021). In this regard, there are 356 species in Mexico (Cuevas-Yánez et al., 2017).
The most common and homogeneous stressor associated with this collapse is LULCC (Clausnitzer et al., 2009). Among other effects, LULCC results in the modifi cation and/or extraction of water from larval habitats or alterations to the terrestrial habitats of adults due to agricultural, road construction and other human activities. These activities affect local ambient temperature, water quality, available resting sites, and/or food availability for the larvae and adults (e.g. Monteiro-Junior et al., 2013;Chávez et al., 2015;Jakob & Poulin, 2016).
In addition to LULCC analysis, models of the distributions of potential habitats of species are important for a wide variety of applications in conservation biology and ecology (Ferrier, 2002), for example, for dragonfl ies and damselfl ies (Boys & Adam, 2021). Several models exist for predicting the distribution of potential habitats of species and each method has its own advantages and disadvantages (Guisan & Zimmermann, 2000). In particular, the maximum entropy (Maxent) model is a species distribution model based on a machine learning response (Phillips et al., 2004;Baldwin, 2009). Maxent predicts the potential distributions of species. This method has several advantages: it only requires species' presence data (even for relatively small datasets) and environmental information or variables (Elith et al., 2011;Rammou et al., 2022). Also, Maxent performs better than other methods when predict-correlated and any redundant correlations discarded. We generated distribution models for each of the species using the program Maxent version 3.4.0 (Phillips et al., 2006), with a pixel size of ≈ 1 km 2 and then reclassifi ed each model with a threshold value of 0.6 for binary values as reported in the literature (Anderson, 2003;Pearson et al., 2006;Velóz, 2009). The models were revised using ArcMap 10.2.1 software (ESRI Inc., Redlands CA, USA). The probability values of the models of the potential distributions of the 49 species ranged between 0 and 1. This range corresponds to probabilities generated by the potential distribution models that were lowest and highest respectively. Finally, we selected probability values from 0.6 to the maximum value for each species in order to compile map for each species. Then, we constructed a richness map in order to identify areas of varying richness, using species with different IVO values (Table 1, Fig. 1). This analysed the effect of LULCC on different species of Odonata with respect to their tolerance. This richness map was then compared with the deforestation map to assess the risk that each species faces due to the deterioration of their habitats. This risk is thus based on each species' distribution.
The fi rst step in creating the deforestation maps mentioned above was to develop a model for deforestation based on LULCC scenarios and variables that promote change, which is either due to human activity (i.e., roads, cities), or physical factors, which when combined can be used to infer or predict future landscape conditions. To develop the deforestation model, we used two digital land use/land cover (LULC) maps at a 1 : 250,000 scale To generate the LULCC-related maps, the vegetation series IV and V were classifi ed and grouped into three categories defi ned by the ecological affi nities of each species of Odonata. The habitat groupings were: (1) Poor or unfavourable, (2) Acceptable and (3) Optimal (see Table S1 and Fig. S2). This classifi cation allowed us to: (a) simplify the different vegetation classes and provide a better interpretation of habitat changes; and (b) facilitate operations carried out by the D-E (for the description of this software, see below), as this software yields better results when using fewer classes of vegetation (Prieto-Amparán et al., 2019). Related to the above, D-E software makes use of variables related to landscape features such as roads, rivers and vegetation, which were analysed using maps to estimate their possible infl uence on the distribution of Odonata.
According to the literature  and expert knowledge we selected the following six variables as possible drivers transitions in LULCC and the subsequent changes in the distribution of Odonata: distances in meters to nearest paved road, water body, perennial river, urban area, location in Natural Protected Areas (categorical values) and nearest distance to Natural Protected Areas (see Table S2 and Fig. 3). For each of these six variables, we transformed six layers into a raster format, with a pixel size of ≈ 1 km 2 the required input of D-E. To analyse the results of the change matrix, we transformed the map information into a matrix format (raster). The transition matrix describes the deforestation in a system over a discrete time interval (t 0t 1 ). These matrices of transition determine the extent of change, which is the percentage of land that underwent a change from one state to another (which is an attribute of LULCC).

Extent of change and transition matrix
The fi rst step in the development of the deforestation model was to produce an annual and global Markov matrix (Table 2A). The annual matrix describes yearly deforestation and the global matrix changes in deforestation over the whole period based on discrete time changes in the response/state variables. On the matrix, the value of any of the six variables (i.e., the drivers of LULCC) for a given period of time (t 1 ) must correspond to the sum of all the values from previous time periods (t 0 ), as this sum indicates the total change. These matrices can be presented as a single-pass matrix (for a given time period, for example, 10 years) or as multi-pass matrix (with a pass for each discrete interval of time, for example, one year per pass) (Soares Filho et al., 2015). The transition matrix, which compares changes or transitions between classes of vegetation over a period of time (Table  2), determines the percentage change. This is the percentage of land that changed to a different state. Based on the analysis of geographic information for each LULC map from the D-E program, the intersections of the land covers in the two LULC series corresponding to t 0 and t 1 (Series IV and Series V) (INEGI, 2012(INEGI, , 2015 were identifi ed. Using these intersections, a transition matrix was developed for the period 2007-2011 in order to identify the changes, using the equation proposed by the FAO (1996) and used by López et al. (2019), to characterize the transitions between the Optimal, Acceptable and Poor categories defi ned above (Table 2B): where r = the change, A 1 = area at date 1, A 2 = area at date 2, y = number of years between the two dates. Positive r values indicate that there was an increase in area and negative values a loss of area.

Modelling and spatial simulation using D-E
For modelling deforestation, we used the methods proposed by Mas & Flamenco Sandoval (2011) (see Fig. S4). The deforestation scenarios generated indicated the loss of optimal habitats for Odonata. The modelling was based on the analyses of past changes, which enabled an evaluation of changes in different types of LULCC and the spatial relationship between the location of these changes and the explanatory variables that infl uence the spatial distribution of these changes. Based on this analysis, areas that were most prone to change (change probability maps) were identifi ed and maps generated that indicate future trends in LULCC.

Model calibration: Calculation of the annual matrix and weight of evidence
The calibration of the deforestation model adjusts the model parameters (Rykiel, 1996). This involves providing information to enable the model to determine the incidence of change, the types of transitions and their most likely location at the geographical scale. For this, we overlapped the t 0 and t 1 maps in order to generate a map of LULCC and a change matrix that indicated the areas with each type of transition during the period of observation (four years). This change matrix can be transformed into a Markov matrix of probability of a change, which indicates the yearly probability of the occurrence of each transition (Soares-Filho et al., 2002).
A comparison of the superimposed maps t 0 and t 1 was used to identify areas that are more likely to change in terms of vegetation and these locations were indicated on another map. These areas of change were compared with each variable that promotes change as before (Fig. 3). This revealed a relationship between the potential for change and variables that promote change by different means, resulting in a weight of evidence coeffi cient (WoE). Using WoE coeffi cients and binary maps it is possible to determine local spatial correlations between several explanatory variables, based on the spatial associations among maps (Goodacre et al., 1993). Using a Bayesian model, the maps are combined to produce a map of the probability of change (Bonham-Carter, 1994;Soares-Filho et al., 2004). This method is implemented in D-E and used to determine transition probabilities. These probabilities indicate areas that are most likely to change, i.e. pix- Fig. 1. Map of the area studied with predicted species richness based on ecological niche models. els that will change from one state to another as there is a high probability that they will change from one type of vegetation to another as a function of the different combinations of probability change that the D-E software is capable of producing (Soares-Filho et al., 2002. Finally, an analysis was carried out to remove highly correlated explanatory variables and verify their spatial independence using Cramer's test (Bonham-Carter, 1994;Soares-Filho et al., 2009). The WoE indicates the correlation among variables and their ranges, which range from 0 to 1, where 0 indicates an independent and 1 a dependent variable. Values above 0.45 indicate strong correlations between variables, which were then removed from the analysis (López et al., 2019). The reason for this is that highly correlated variables have the same effect on LULCC, thus including two strongly correlated variables overemphasize the effect of those variables on the probability calculation.

Model simulation: LULCC
Using the WoE and maps of the explanatory variables indicated above, we generated a probability map (Fig. 3) for each of the three transitions between Optimal, Acceptable and Poor habitats in the LULCC, as explained above. The model for simulating the LULCC used the probability map and functors PATCHER and EXPANDER (Soares-Filho et al., 2002).
This generated a map, called the "real observed map" (here, the map of the three categories, Optimal, Acceptable and Poor) of the Series VI map (Class "Poor " is the last map in Fig. S2). This "real observed map" was compared with the map that was simulated to determine the model's capacity to predict changes in classes of vegetation and estimate how similar they are to each other (Fig.  2). Using the probability map described above, prospective maps for the three defi ned classes were generated. To obtain more realistic results, D-E uses two functors to reproduce the spatial patterns in change. The fi rst functor is PATCHER, which generates changes in patches (for example, zones that are poor habitats for Odonata in zones of optimal habitat), while the second functor, EXPANDER makes changes to existing LULC classes (for example, zones of poor habitat extending into nearby or neighbouring zones of optimal habitat).

Model validation
The map based on the previous simulation was validated using the 'real observed map'. This validation used the 'fuzzy similarity index ' (FSI;Hagen, 2003) modifi ed and implemented in D-E by Soares-Filho et al. (2009) and a constant decline function with multiple windows. This generates comparisons of window sizes of one pixel (1000 × 1000 m, i.e. 1 km 2 ) and 1, 3, 5, 7, 9, 11, 13, 15, 17, etc, up to the desired maximum size of 7 × 7 pixels, beyond 60% similarity (see Table S3). The concept of minimum similarity refers to the comparison that D-E makes between real and simulated maps, while maximum similarity is the comparison between the simulated and real map. The validation consists of comparing the spatial coincidence between the simulated and observed LULCC maps, accounting for the spatial coincidence at several pixel sizes (Eastman et al., 2005;Pérez-Vega et al., 2012). In the FSI, the representation of a pixel is infl uenced by itself and its neighbours (Ximenes et al., 2011;Yanai et al., 2012;Chadid et al., 2015), focusing on the areas of change while considering the pixels that surround them (Mas & Flamenco Sandoval, 2011). In other words, the FSI verifi es the concordance between the simulated and observed sets of LULC data by determining the number of cells that coincide within increasing window sizes in a neighbourhood (Costanza, 1989;Soares-Filho et al., 2009).

Analysis of the transition matrix
Of the six possible transitions (Table 2), only those from Optimal to Po or (3 to 1), Acceptable to Poor (2 to 1) and Optimal to Acceptable (3 to 2) that indicate habitat deterioration was considered. In the case of the annual transition, the percentage change from Optimal to Acceptable to Poor was only 0.47%. However, the transition from Optimal to Acceptable occurred in 9.6% of the pixels, which indicates a trend towards site deterioration (Table 2).
In the case of the global transition matrix, the trends were very similar, but not of the same extent. The transitions from Optimal and Acceptable to Poor occurred in just 0.12% of the pixels and the transition from Optimal to Acceptable in 2.5%, which indicates a less marked trend (Table 2).

Simulation: Analysis of the weight of evidence and map of the probability of deforestation
For the three categories of habitat (Optimal, Acceptable and Poor) there was a correlation of 0.587 between NPA and the distance to NPA, so the former was not included in order to ensure that all the remaining variables were spatially independent.

Validation of the deforestation model
The similarity of the values above 50% for the maps compared is considered acceptable for the validation of the deforestation models (Piontekowski et al., 2012). In our case, this percentage was recorded for the 7 × 7 pixel (7 km 2 ) with a maximum similarity value of 60% (Fig. S5). This means that the model is able to correctly assign sites of LULCC within a radius of at least 7 km 2 .

Trends in transitions and deforestation
The results of the deforestation trajectories from D-E indicated that the most common change in terms of area was that from Acceptable to Poor and occurred over an area of 260,644 km 2 . This change had probability values of between 0.4 to 0.9 and central Mexico was the area most strongly affected by these changes. The predictions of the level of probability for each cell analysed are based on cu-mulative values, for the cell in question. A cell with a value of 100 is the most and that of one with values close to 0 the least adequate. The probability of a transition from Optimal to Acceptable was similar and occurred over an area of 60,706 km 2 . Finally, the transition from Optimal to Poor occurred only over an area 37,950 km 2 (Fig. 3). Importantly, a substantial but fragmented area was categorised as Poor. That is, these are areas whose neighbouring pixels were not previously categorised as Poor. In addition, the following drivers of change related to the loss of vegetation cover were identifi ed as closeness to roads (highways) and urban areas. The drivers that favoured the continuance or recovery of vegetation cover were mainly the distance to rivers, followed by distances to water bodies and Natural Protected Areas. The simulated maps were very similar to real maps because the areas that remained unchanged were similar. In other words, areas that did not change at a certain scale dominate the landscape and are indicated by zero changes on the map. Thus, the comparison between the real map and the simulated map was based on the spatial fi t between the areas that changed using fuzzy focus.
Prediction for the year 2050 was based on a comparison of the initial vegetation map of series IV with that of series V. Based on this comparison, prediction of the trends in deforestation was based on modelling the environmental dynamics of both series using D-E. This simulation of the predicted deforestation enables one to infer that the areas that suffer a deterioration in vegetation will be less suitable for Odonata.

Integration of the simulated map for 2050 and the results of the general species distribution model
One of the most notable characteristics of the map simulated for 2050 was the patchiness of deforested areas (Fig.   4). This may be due to the roadway variable (highways), which generated this patchy effect.
The intersections between the predicted map and the predicted area of distribution of the species of Odonata (Fig. 5) revealed an affected area of approximately 40% of zone C. This affects areas with different ranges in species richness from areas with 3 species (green) to those with 49 species (blue) (Fig. 5). The second most affected zone was zone A, where there were also slight effects on all species richness intervals due to changes in vegetation. Finally, zones B and D were the least severely effected in the area and had the highest species richness, with around 5%. However, in these zones, areas with species richness between 3 to 20 showed a 30% reduction due to deforestation. This indicates that the effects of deforestation in these two areas with fairly high species richness are in category 2 -Acceptable, which indicates that there will be no transitions to unfavourable types of vegetation.

DISCUSSION
This prediction of the future for Odonata is novel and based on combining ecological niche and LULCC modelling. The deforestation model indicates that in general for all transitions the changes will result in a moderate deterioration in habitats.

Analysis of the transition matrix
The deforestation model based on LULCC and driving variables indicated that, of the three transitions, the change to the category Poor was less than 1%, which is not a serious deterioration. However, the trend in site deterioration, in terms of Optimal to Acceptable was nearly 10%, which confi rms a previously reported trend for Mexico (Mas & Flamenco Sandoval, 2011;López et al., 2019).

Simulation: Analysis of weight of evidence and map of probability of deforestation
The deforestation model correctly assigned sites in terms of vegetation cover, which allowed us to predict the situation in the year 2050. The trends predicted by the deforestation model based on D-E revealed deforested areas that were not close to previously deforested or Poor sites. In other words, the deforestation model predicted deforested sites that were far from those previously recorded. We believe that this is the result of using driving variables, such as roads and urban areas, which coincides with the fi ndings of Farfán (2015), Espinoza (2016) and Viégras et al. (2020). For example, Farfán (2015 reports that the deterioration of natural vegetation is associated with variables such as communication routes, which have an important effect on the dynamics of LULCC. On the other hand, Espinoza (2016) reports that the constant changes in the type of agricultural vegetation is associated with an increase in communication routes and that areas closest to the communication routes are the most affected in terms of loss of forest cover. Finally, Viégras et al. (2020) indicate that roads promote urban growth and change in vegetation. However, when there are conservation areas, these prevent the appearance of more roads and urbanization.

Trends in transitions and deforestation
Trends in transitions indicate that changes in LULCC, such as the transformation to agricultural land or urbanization have the greatest effect on freshwater ecosystems (e.g. Li et al., 2021). Thus, Odonata and other aquatic insects cannot be effectively protected without stopping the modifi cation and pollution of their aquatic and terrestrial habitats. Although Odonata are relatively resilient to environmental change , other aquatic insects are not, for example: Ephemeroptera, Plecoptera, Trichoptera, Heteroptera and Chironomidae (Albert et al., 2020;Brasil et al., 2020;Cardoso et al., 2020). Thus, there will inevitably be a reduction in insect abundance and/or richness. Although we provide evidence of environmental deterioration in Mexico, this occurs and will continue to occur in many other regions due to habitat deterioration and deforestation.

Integration of the predicted map for 2050 and the results of the species distribution model
The integration of the fi nal map with the predicted level of deforestation in 2050 and the species richness of Odonata with different IVO values indicate that the four most species-rich zones will be affected by deforestation. Thus, the south-eastern region of Mexico will be severely affected by these changes. However, it is important to point out that many other sites within the study area will be severely affected by deforestation. Our results also indicate that of the 29 species with IVO values of 1.5 or 1.0, many are not present in drastically disturbed areas. These results allow us to infer that at least half of the species of Odonata will be affected by habitat deterioration. Thus, the same Odonata species richness is present in areas that have undergone land use change, but the local species composition has changed considerably (higher beta diversity), which is similar to that reported by . In addition, however, the consequences of habitat modifi cation for each species is unknown. These effects need further study since there are species that are highly plastic and resilient, which allows them to remain at these sites (García-García et al., 2017. Both Rocha-Ortega et al. (2019) and the present study involved the same species in the same areas. On the other hand, there are other species for which even slight changes can have strong negative effects, for example, the damselfl y Hetaerina americana in central Mexico, for which the abundance and physiological condition of adults has decreased due to pollution of rivers, but they are still present (Córdoba-Aguilar & Rocha-Ortega, 2019).
This study clearly demonstrates the diffi culty of developing predictive models for regions subject to dramatic changes. There have been great changes in agriculture and a marked increase in deforestation recently due to the constant pressure to change the use of land for mainly for urban and agricultural use. In spite of this, tools like D-E enable the development of multi-temporal predictive models. One can produce suitable models, which integrate expert knowledge (i.e., they are "knowledge driven") and in which the variables can be weighted or modifi ed based on the knowledge of the study group. In our case, the "knowledge driven" consisted of our experience, relevant literature and asking experts which variables are the most important for the development and survival of Odonata. Thus, D-E is an excellent tool for identifying the drivers and rates of change, and in particular the prediction of future scenarios.
Another important contribution of this study is the development of a deforestation model that can be used to predict the situation in 2050 and produce niche models for 49 species of Odonata, which are presented in a map of species richness. Integration of both maps enabled us to locate sites where the habitats of Odonata were very poor. In this sense, there are no other studies that combine these two methods and predict future changes in land use and/or deforestation. In any case, it is important to replicate this analysis using other groups, such as aquatic invertebrates and vertebrates.

CONCLUSIONS
We estimate that zone A will be the area where the effects on geographic distribution and land use will be most marked (Fig. 5). Nevertheless, zone D will be the most affected in terms of the habitat deterioration of almost half of the species of Odonata. In terms the IVO values most species will be affected, especially those with values of 1.5-1.0 despite their tolerance of perturbation.
The effects on the species of Odonata varied in the different zones shown in Fig. 5. Zone A is close to an important Nature Reserve (RB-Sierra Gorda). However, in the central-southern part of this zone there are extensive urban and agricultural areas with a highly fragmented landscape and high risk of an increasing change in land use. Zone B is located within another important Nature Reserve (RB-Barranca de Metztitlán) with a large number of ravines, canyons and an extensive forest. The soil in this zone does not differ very much, however there is a signifi cant pressure from agriculture at the margins of the reserve. Zone C is located within one of the most highly diverse areas in Mexico, in the Sierra Madre del Sur. This area is mountainous, with a considerable number of rivers. Finally, zone D is located in one of the wettest regions in Mexico, with an area of 2,500 km 2 and abundance water bodies. Over the last two decades this zone was subjected to both deforestation and use of land for agriculture, keeping only a small area of natural vegetation as a protected nature reserve.