Evaluating Potential Distribution and Niche Divergence among Populations of the World’s Largest Living Damselfly, Megaloprepus caerulatus (Drury, 1782)

Megaloprepus caerulatus is a Neotropical species with a highly specialised niche, found from Mexico to Bolivia, primarily in mature tropical forests lower than 1500 masl. It is also the damselfly with the largest wingspan in the world. Recent studies found strong genetic isolation among populations of M. caerulatus. Further studies found genetic and morphological divergence, but ecological divergence was not tested. Here, we test for ecological divergence by evaluating niche differences among populations of M. caerulatus in Los Tuxtlas (Mexico), Corcovado (Costa Rica), Barro Colorado (Panama), and La Selva (Costa Rica). We used Ecological Niche Modelling (ENM) to compare potential distribution ranges, and we estimated the breadth and overlap of the ecological niche using equivalence and similarity tests. The potential distributions estimated with ENM were heavily fragmented and we found no geographic overlap of potential distributions among populations. However, we found geographic correspondence between populations with a close phylogenetic relationship. Even though all similarity tests were non-significant, the results of the equivalence tests suggest niche divergence between Corcovado and the other three populations, but also between Barro Colorado (Panama) and La Selva. These results show evidence of strong ecological divergence in Corcovado and Barro Colorado populations.


Introduction
The niche is integral in ecology for the study of species and as tools for conservation. Many contributions have been made on niche conceptualization, all of which relate to the basic taxonomic unit, the species. Since Grinnell [1], niche theory has been constantly developed and multiple studies at fine (i.e., physiological) and coarse (geographical) scales have been explored. Although direct evaluation of physiological constraints would be ideal to niche reconstruction, specific data are unavailable for many species in all taxonomic groups. Thus, Ecological Niche Modelling (ENM) can be used as a proxy. Tolerance limits reflect the fundamental ecological niche of a species and can be studied via their coarse-resolution associations with environments manifested across their geographic distributions [1][2][3][4]. Under certain assumptions and limitations, the correlational ENM approaches can estimate limits that approach the fundamental niche [5].
Within the fundamental niche, the existing niche is the one that corresponds to environments represented within the species' distributional area [2,6]. This existing niche can be estimated using correlative modelling approaches [5]. The ENM has been a widely used approach for estimating niches by searching for associations between species localities and their corresponding environmental conditions [5,[7][8][9][10][11][12][13][14]. Estimating abiotically suitable, potential, and occupied distribution areas has been the oldest and most widely other as other odonate species, suggesting these subspecies in fact are more likely different species. Subsequently, Fincke et al. [31] identified genetic and morphological clades among the aforementioned study populations, in addition to the following populations: La Bartola Reserve (BART) from Nicaragua, Canandé Reserve (CAN) in Ecuador (both dimorphic), El Jaguar (EJ) in Nicaragua, and Cusuco National Park (HON) in Honduras (both monomorphic). Monomorphic populations showed lower adult density, lower resource defence, fewer male-male interactions and, consequently, lower sexual selection on males, suggesting sexual selection is a diverging mechanism between monomorphic and dimorphic populations [31].
However, ecological divergence among populations or subspecies of M. caerulatus has not been tested. In this study we evaluate whether niche characteristics of Megaloprepus remain constant or instead are divergent. We analyse how intraspecific genetic-level ENMs of Megaloprepus might differ in terms of (1) potential distribution ranges, and (2) the breadth and overlap of the ecological niche between populations or subspecies.

Biological and Environmental Data
In order to have a complete database of the historical records of M. caerulatus, an exhaustive search was carried out from literature [31][32][33]35,[43][44][45][46][47][48][49] and from the GBIF database [50], which is an online portal providing access to primary biodiversity data. Even though only "research grade" records were selected, all GBIF records were verified whenever possible through photographs of the specimens to make sure they were correctly identified to species and morphotype level. Records lacking GPS coordinates but with detailed locations were georeferenced with Google Earth Pro. The resulting database had 139 records (80 from GBIF, 59 from taxonomic records) ranging from Mexico to Ecuador (Figure 1a).
Our study groups were chosen from the eight genetically and morphologically distinct populations from the subspecies M. c. caerulatus, M. c. latipennis, and M. c. subsp. nov. suggested by Feindt et al. [35] and Fincke et al. [31]. For this purpose, we analysed the 139 records in our database to find geographic (based on terrestrial ecoregions; [51]) and morphological (i.e., dimorphic and monomorphic) correspondence with the populations of the former studies ( Figure 1a). We prioritised those populations with a strong reference of genetic and morphological difference and with an adequate number of records to complete the ENM. Therefore, our study groups are four, the northernmost in Mexico and the southernmost in Panama (Figure 1b). There are two monomorphic groups, Los Tuxtlas (16 data points; Figure 1c) and Corcovado (15 data points; Figure 1d), as well as two dimorphic groups, La Selva (44 data points; Figure 1d), and Barro Colorado (9 data points; Figure 1d). For simplicity, we refer to our study groups as "populations", although the monomorphic groups Los Tuxtlas and Corcovado are different subspecies (M. c. latipennis and M. c. subsp. nov., respectively), while dimorphic groups La Selva and Barro Colorado are populations of the same subspecies, M. c. caerulatus.
We used spatial environmental data relevant to the biology of the species to produce the dataset that would represent the ecological niches [52][53][54][55]. Climatic data were obtained from WorldClim Version 1.4 [56] (http://www.worldclim.org/ (accessed on 19 May 2018)). Land use was obtained from EarthEnv [57] (http://www.earthenv.org/ (accessed on 3 October 2019)) and soil evapotranspiration from the Global High-Resolution Soil-Water Balance dataset [58]. Topographic data were obtained from GMTED2010 [59] (https://on.doi.gov/30VfqfR (accessed on 10 October 2019)) and HydroSHEDS [60] (https: //hydrosheds.org/ (accessed on 24 November 2019)). All spatial data were downloaded at a spatial resolution of 30 s (~1 km 2 ). To avoid duplication of environmental information, a correlation analysis was performed. The environmental values associated with the occurrence data (one per pixel) were used to perform the Spearman non-parametric correlation test with the "Hmisc" package [61] and a principal component analysis in R [62]. The variables that were correlated with a greater number of other variables were eliminated.
We also conducted an exploratory Jackknife test with the niche modelling algorithm of Maximum Entropy (Maxent, v3.3.3.e; [63]) in order to identify the variables that contribute the least to the construction of the model, and thus eliminate them prior to modelling the definitive ecological niche model and the similarity analyses. This resulted in the selection of 15 environmental variables: mean diurnal range, isothermality, temperature seasonality, temperature annual range, annual precipitation, precipitation of wettest month, precipitation of driest month, precipitation seasonality, precipitation of wettest quarter, precipitation of driest quarter, precipitation of warmest quarter, evergreen broadleaf tree cover, mean annual evapotranspiration, drainage direction, and flow accumulation. In fact, many of these variables have also been of importance in other odonate ENMs [11,14]. We used spatial environmental data relevant to the biology of the species to produce the dataset that would represent the ecological niches [52][53][54][55]. Climatic data were obtained from WorldClim Version 1.4 [56] (http://www.worldclim.org/ (accessed on 19 May  [35] and Fincke et al. [31] as well as all historical records for the species (black circles). This study focused on Los Tuxtlas (TUX), La Selva (SELVA), Corcovado (CNP-SIR), and Barro Colorado (BCI); (b) Occurrence data used for niche modelling and similarity tests for TUX (purple), SELVA (yellow), CNP-SIR (green), and BCI (orange); their ecological niches were projected to a greater geographic extension (black solid line polygon) which encompasses all known records for the species; (c,d) Occurrence data and models for each population, which were calibrated in M regions delimited by biological and geographic conditions, also shown as black solid line polygons.

Ecological Niche Modelling
In ENM, three major factors are considered to explain distributions of species: biotic (B), abiotic (A), and mobility (M) constraints (BAM; [2]). Biotic factors are denoted by B; however, at coarse resolutions, the biotic component is frequently diffuse and non-limiting, in contrast to how it is manifested at finer spatial resolutions [64]. On the other hand, the remaining two components have broad-scale effects. The abiotic factors, called A, represent the geographic region presenting favourable conditions, and, finally, M is the area that has historically been accessible to the species via dispersal over relevant periods [5]. Although this approach is simplified based on static approximations to the three classes of factors [65], it has proven to be a useful heuristic test. Since our hypothesis was tested on geographic extents, a coarse resolution, the component B was not considered. Evaluations of niche similarity were made in terms of whether two niches are more similar than expected given the set of environments accessible to each population across their M [66,67]. We tested the niche similarity under environmental space [68] to compare the niches of Megaloprepus populations and test for divergence.
For niche model calibration, we applied the Maximum Entropy algorithm implemented in Maxent [63], which fits a distribution of probabilities across the study area subject to the constraints of the environmental characteristics of known occurrences. Evaluation data were separated a priori, therefore no data were assigned for evaluation within Maxent. Other settings included: regularization value = 1, maximum number of points for the background = 10,000, maximum iterations = 500 and convergence threshold = 0.00001. We turned off the clamping and extrapolation options, following Owens et al. [69] to avoid artificial extrapolations of extreme values of environmental variables. To convert the output into a binary map, we used the maximum training sensitivity plus specificity threshold, because it has been shown to be an optimal method for presence-only data [70].
Only in the case of the potential distribution of SELVA, occurrences were divided into random subsets: 80% for model calibration and 20% for model evaluation. Conversely, for the rest of the Megaloprepus populations, fewer than 25 occurrence points were available, therefore models were calibrated with all data. Models were calibrated across regions posited as historically accessible to each population ( Figure 1b). M regions were delimited considering aspects of the distribution and life history of the population [65], using the limits of the surrounding ecoregions [51] and watershed boundaries and sub-basin delineation [71] as a guide (Figure 1c,d). Once the models were calibrated, they were transferred to a broader region, including the union of the M hypotheses across all of the Megaloprepus study populations (Figure 1b).
The model evaluation was performed differently according to the sample sizes of populations. The potential distribution of SELVA had a large sample size (N = 44), therefore the model was evaluated using a modification of the area under the curve of the receiver operating characteristic (partial ROC AUC ratios; [72]) using the graphical interphase Niche Toolbox [73]. This test only evaluates over the spectrum of the prediction and allows for differential weighting of the two error components (omission and commission; [72]). Thus, AUCs were limited to the proportional areas over which the models truly made predictions, and we only considered models that presented omission errors of less than 5% [72]. In the case of BCI, CNP-SIR and TUX, evaluation was accomplished via the Jackknife strategy developed for small sample-sizes by Pearson et al. [74]: significance was evaluated over n models, each excluding one locality from among the n available and evaluating the success of the model in terms of anticipating the excluded locality. The probability of these observed levels of success and failure was calculated using scripts provided by Pearson et al. [74]. This test was applied to binary models created by applying Minimum Training Presence (MTP) approaches [74].

Niche Similarity Tests
The niche similarity tests were performed in an environmental space as proposed by Broennimann et al. [68]. This method uses Schoener's D metric [75] as a measure of environmental niche overlap and includes a statistical framework to test for niche similarity as proposed by Warren et al. [66]. The value of D ranges between 0, when two populations have no overlap in the environmental space, and 1 when two populations share the same environmental space. With this method, a multivariate environment grid is created using the first two axes of a PCA that summarise all the environmental variables previously selected (PCA-env). A Gaussian kernel density is applied to estimate the occupancy of each cell (z ij ), and the D metric is calculated based on the different z ij values obtained [68]. We tested for niche equivalence to assess whether the ecological niches of a pair of Megaloprepus populations were significantly different from each other and if the two niche spaces were interchangeable. This test only assesses if the two populations are identical in their niche space by using the environmental data of their exact locations and does not consider the surrounding M space, unlike the similarity test [76]. Equivalency test compares the niche overlap values (D) of a pair of populations to a null distribution of 100 overlap values. On the other hand, niche similarity test assesses if the ecological niches of any pair of populations are less similar than expected by chance, accounting for the differences in the surrounding environmental conditions in the M regions, which are the geographic areas where the populations are distributed [66]. We determined that ecological niches in comparison were less equivalent and/or less similar if the niche overlap value of the populations being compared was significantly lower than the overlap values form the null distribution (p ≤ 0.05). This analysis was performed using the Ecospat package [77] in R. Since our goal was to test for niche divergence, we selected alternative = "lower". Thus, we present six cross-comparisons between Megaloprepus populations.

Results
The niche model performance of all populations was high and better than expected by chance. In the case of SELVA, the Partial ROC value was 1.606 (p = 0) and in the case of BCI, CNP-SIR and TUX the jackknife test (used for populations with fewer occurrence data) showed high success rates (0.77, 0.87, 0.94, respectively), as well as statistical significance (p = 0).
The environmental variables used contributed differently to the ENM of each population. For example, in the case of CNP-SIR, the variables that contributed the most were precipitation of the wettest month (45%) and mean diurnal range (10.3%), whereas precipitation seasonality contributed with 8.6%, mean annual evapotranspiration with 8.2%, isothermality with 6.3% and temperature annual range with 5.7%. For BCI, mean diurnal range also contributed significantly to the model with 57%, as well as evergreen broadleaf tree cover with 26.5%, precipitation of the driest quarter with 6.3% and flow accumulation with 5.1%. For SELVA, the two major contributors to the model were mean annual evapotranspiration (30.6%) and precipitation of the driest quarter (20.6%), although other variables that contributed significantly were mean diurnal range (15.5%), temperature annual range (12.2%) and flow accumulation (9.2). For TUX, mean diurnal range was the variable with the most significant contribution of 58.1%, others were mean annual evapotranspiration (11.7%) and precipitation seasonality (9.6%). Other variables had a contribution of <5% in each model.
The niche projection resulted in the potential distribution of the four populations studied. We found that the potential distributions of the four Megaloprepus populations do not overlap. Additionally, all potential distributions are fragmented and, in most cases, discontinuous ( Figure 2). The potential distribution of TUX is restricted to Sierra Los Tuxtlas in Veracruz at the coast of the Gulf of Mexico (Figure 2a). Interestingly, the model does not predict suitable conditions towards other regions to the north or south where nearby historical records for Megaloprepus are located (see biological and environmental data in methods; Figure 1a). The potential distribution of SELVA is mainly located in Eastern Costa Rica and North Panama (next to the Atlantic Coast), but partly extends to Colombia, relatively close to the Canandé population of Fincke et al. [31], where other historical records in Colombia are found (Figure 2b). The potential distribution of CNP-SIR is predominantly found in Western Costa Rica, except for a few small, patchy areas towards the Pacific coast of Guatemala and the Southern Mexican border (Figure 2c). The potential distribution of BCI is predominantly found in Panama but extends discontinuously towards the east of Nicaragua along the Atlantic Coast (Figure 2d). This potential distribution corresponds with La Bartola population studied by Fincke et al. [31], but not with La Selva, which is located between the two (Figure 2d).  In order to compare the niches at the environmental space, we performed pairwise tests for equivalency and similarity, therefore there are six comparisons for the four study populations. What we found can be summarised in three points. First, in all cases the PCA-env plots show low overlap between M environments (D = 0, except in Figure 3c) (Figure 3, left graphs). Second, the equivalence tests show evidence of niche divergence in four of the comparisons: TUX and CNP-SIR, SELVA and CNP-SIR, SELVA and BCI and between CNP-SIR and BCI (p ≤ 0.05; Figure 3a,d-f). It is important to keep in mind that this evidence of divergence occurs when comparing the niche with environmental data of exact locations and without considering the surrounding space, contrary to the analysis of similarity. Third, no evidence of niche divergence was found in the similarity tests, since all pairwise comparisons were non-significant (p ≥ 0.05) (Figure 3a-f).   [68]. Shaded areas in each plot show the density of the occurrences of the populations by cell. The solid contour lines illustrate 100% of the available (background) environment. The green colour represents the niche of the first population and the blue colour the niche of the second population. The red shaded areas are the niche intersections among kernel densities of occurrences. The histograms correspond to the results of equivalency (centre) and similarity tests (right) to test for niche divergence. Histograms show the observed niche overlap D between the two ranges (red lines with a diamond) and simulated niche overlaps (grey bars) on which tests of niche divergence are calculated from 100 iterations. Test significance is shown (ns, non-significant; *** p < 0.001).

Discussion
This study shows that the potential distributions of the Megaloprepus populations studied in this work are formed by highly fragmented areas of suitable habitat ranging from South Mexico to South Colombia. According to Fincke et al. [31] no Megaloprepus subspecies are known to occur sympatrically. Our results support this fact, since the potential distributions do not present any geographic overlap (Figure 2a-d). However, we found geographic correspondence between populations with a close phylogenetic relationship [31,35] despite being rather separated populations (e.g., La Selva and Canandé; Barro Colorado and La Bartola).
Our results suggest that the potential distribution of Los Tuxtlas (M. c. latipennis), which is monomorphic, is heavily isolated. There is no potential distribution nearby areas where historical records have been found, such as localities in the state of Chiapas, Mexico [48] (Figures 1a and 2a). This may be explained by the biogeographic history of the region to which the subspecies is restricted, the Sierra Los Tuxtlas. This is an isolated mountainous area on the Gulf of Mexico that was originated due to intense volcanic activity in the Miocene [78]. It is also the northernmost relict of moist tropical forest in the continent [79], in which vertical colonization [78], divergent evolutionary processes [80], and high rates of endemism have been documented [45].
The CNP-SIR potential distribution from Costa Rica (M. c. subsp. nov.), also monomorphic, is predominantly found on the Pacific coast of Costa Rica, but despite a large distance, it also reaches small areas in the north of Guatemala and the state of Chiapas (Mexico) next to the Pacific (Figure 2c). This northern range of suitable conditions for Corcovado (Figure 2c) almost coincides with several records of Megaloprepus found in Guatemala (Figure 1a). This can be partially explained because the Sierra Madre of Chiapas and the Central Chiapas Massif constitute the northern projection of the Central American mountain system [79]. These elevations yield a specific floristic composition ranging from Tropical-Subtropical moist forest to Tropical-Subtropical coniferous forest [81].
We also found interesting results for dimorphic populations (M. c. caerulatus). Despite the geographic proximity between the populations of La Selva and Barro Colorado, their potential distributions do not overlap (Figure 2b,d). All the records used in the niche modelling for both populations are located in the East Central America biogeographic province sensu [82]. The disruption of the potential distribution of La Selva may have originated from biogeographic barriers formed during the Pleistocene, which is a long period wherein landscapes had a dynamic history of dramatic changes, due to geological and climatic processes that impacted a wide variety of taxa through fragmentation and displacement of populations [83]. The potential distribution of the La Selva population shows the Napo province [84] in South America as an ideal habitat if individuals could have access (Figure 2b), which corresponds with the Napo Pleistocene Refuge proposed by Haffer [85]. However, it is important to mention that there is a debate regarding the validity of the Refugial Hypothesis in South America [86]. Therefore, further research is needed to unveil the underlying mechanisms driving the potential distribution of La Selva.
Additionally, we found significant ecological divergence between La Selva and Barro Colorado, despite being very similar populations, behaviourally and morphologically [31]. Perhaps the main difference between La Selva and Barro Colorado is their environmental conditions, particularly rainfall and seasonality. La Selva is a wet tropical forest where the larval habitat of the species seems not to dry up, whereas Barro Colorado is a tropical moist forest where larval habitats dry up for 2-3 months each year [31]. Adults from Barro Colorado population could become inactive during the dry season, which would imply an ecological adaptation via seasonal changes in behaviour. Behavioural traits have been shown to affect niche similarity or divergence in other species at inter and intraspecific level [24,29]. Even though many of the environmental variables that contributed significantly to the niche models are related to rainfall and seasonality, we cannot conclude from our data that these factors play a role in the ecological divergence found, because the PCA-env used in similarity and equivalence tests summarises all environmental variables in a multivariate environment grid using two PCs. This makes it difficult to ascertain which factors are driving ecological divergence. Alternatively, the ecological divergence between La Selva and Barro Colorado could be due to possible physiological differences in the larvae or adults, which would need to be verified in future studies.
We must mention that this study largely focuses on the abiotic factors driving niche divergence. However, Megaloprepus depends entirely on water-filled tree holes for reproduction. These holes are formed from indentations in the tree bole or buttress-particularly in tree species that are susceptible to burl formation-or in fallen trees in which flutings are filled with water [38,39]. Trees with smooth boles, such as Bursera simarouba, are unlikely to form tree holes [39]. Megaloprepus reproduces predominantly in tree holes found in gaps where other trees or large branches have fallen [87]. Moreover, large holes (>1 L) are relatively rare but support a greater number of emerging adults per season and provide more resources for the larvae, producing larger adults than small holes (<1 L) which usually only support one emerging adult [38,87]. In fact, this is the reason why territorial males primarily defend large tree holes as breeding sites [87]. Only a subset of available tree species can provide water-filled holes (especially >1 L), such as Ceiba pentandra, Dipteryx panamensis, Platypodium elegans, Ficus trigonata, and Alseis blackiana, among others [39]. Hence, tree holes are a limiting population resource for Megaloprepus [38]. Even though we included the general habitat type in the "Evergreen Broadleaf Tree Cover" variable, we did not consider the distribution of the tree species that yield the water-filled holes needed for Megaloprepus to thrive. This is a crucial biotic factor that could play a key role in their distribution and divergence, and should be accounted for in future studies.
One of the most important findings in this study is the niche differences between Megaloprepus populations in the geographic and environmental space in the equivalence test, particularly in Corcovado (M. c. subsp. nov.) and Barro Colorado populations (M. c. caerulatus) (Figure 3f). However, we must emphasize that ecological divergence is rare or very slow-occurring in species with highly specialised niches [30]. Numerous studies have found that niches generally remain constant throughout phylogeny, at least in the short-to-medium term [18,24,66,88,89]. Therefore, ecological divergence could be a result of an ancient and complex geographic event. Toussaint et al. [90] tested three different clock partitioning schemes and two different tree models, which suggest an ancient origin for the diversification of Neotropical giant damselflies in the Paleogene-Eocene (50-40 Ma). Feindt et al. [35] assume that the historical distribution of helicopter damselflies might have been in the northern portion of South America. According to morphological and phylogenetic studies, Megaloprepus is closely related to Anomisma [91] and Microstigma [90]. These taxa present an exclusive current distribution in South America with a northern limit in Colombia, without having any presence in Central America [92]. It has been suggested that the closing of the land bridge connecting south Nicaragua and Colombia was during the late Pliocene 3 Ma [93,94], although new palaeontological studies suggest this occurred [13][14][15]96]. Until the closing of this land bridge, migration of montane entomofauna between the Central American Nucleus and South America represented a very difficult undertaking and was only possible through the mountains of Talamanca [97]. To this day, fauna from Costa Rica and Panama show a higher affinity towards South American biodiversity than fauna from Mexico and other countries in Central America [97]. The mountains of Talamanca are over 1500 masl and, considering M. caerulatus is only found in areas below 1500 masl, these mountains may have constituted an important barrier for dispersal of Megaloprepus between East Central America and the Pacific region. Geographic barriers limit gene diffusion and populations become geographically isolated [30]. This could be driving divergence in the Corcovado population, located west to the mountains of Talamanca.
On the other hand, Fincke et al. [31] found the Tuxtlas population as the most ancestral, suggesting the species dispersed in a N-S direction. This would change the perception of the dispersion of Megaloprepus described thus far, which highlights the need to investigate the areas of endemism, along with the primary biogeographic homology and if possible, the relationships between the areas of secondary biogeographic homology sensu [98]. Due to the complex history of Central America, another scenario is also possible where dispersion was first S-N, then extinction occurred in intermediate fragments, and when the temperature and other environmental conditions changed in southern Mexico [97], the populations dispersed from N-S along the continuous fragments remaining. If our results of geographic correspondence between populations with a close phylogenetic relationship [31,35] are evidence of an ancient, continuous distribution, then they should be subject to discussion and testing. This is particularly necessary if the distribution eventually became fragmented as a result of natural events that occurred long ago, perhaps during the Pleistocene, which is likely given the correspondence of La Selva potential distribution with a Pleistocene refuge in South America [85,99].
In order to decipher the biogeographic history of M. caerulatus, it is paramount to identify more populations representative of South America, particularly monomorphic populations found in the northern portion of South America (M. caerulatus brevistigma; [42]), since very few records were found in this area (Figure 1a). It is also important to characterise populations found north of Los Tuxtlas [43,47], Chiapas (Mexico) and Guatemala to verify if their phylogenetic relationship is closer to Los Tuxtlas (M. c. latipennis) or Corcovado (M. c. subsp. nov.), as suggested by the potential distributions found in this study. This would shed light on the divergence of Corcovado (M. c. subsp. nov.) and provide more information to characterise monomorphic populations from Los Tuxtlas (M. c. latipennis) or from South America (M. c. brevistigma) as the most basal node.
Conversely, no difference was found in niche similarity among the populations studied. This approach not only considers the record's environmental data, but also the surrounding area, which can be heterogeneous in fragmented landscapes [100]. Megaloprepus caerulatus is a species with a highly specialised niche [35] which shows limited dispersal in open areas [39] and is susceptible to forest conversion [34]. Landscape fragmentation is a wellknown issue especially in Los Tuxtlas [39], although fragmentation rates and connectivity to other forested areas vary considerably among the populations studied [35]. Therefore, environmental heterogeneity may constitute an excluding factor, which might explain the non-significant niche similarity and significant niche equivalence found in this study. Nevertheless, the relevance of surrounding environments (i.e., M environments) in highly specialised species must be tested in order to relate niche theory to niche conservatism or divergence.
So far, we have speculated about the possible drivers of the ecological divergence found in this study, which is likely a result of complex processes that occurred long ago. However, Megaloprepus was found to be fairly sensitive to recent land use change [13], which could be due to various factors. Firstly, odonates-particularly zygopterans-with a large body size, such as Megaloprepus, are more prone to extinction [101,102]. This could be due to the long development period required to achieve such large size, which makes them more vulnerable to predators [101]. On the other hand, it may also be related to their thermal tolerance. Larger water-breathing ectotherms may be more susceptible to impaired heat tolerance by oxygen limitation [103], which could restrict the capacity of large damselflies to obtain oxygen from the environment [102].
Additionally, their highly specialised niche can make them particularly vulnerable to current anthropogenic land use change. As previously mentioned, Megaloprepus is highly dependent on mature, moist Neotropical forests and it is most common in primary forests [34]. Although they can disperse considerable distances in forest understorey, they show relatively low flight endurance in open areas (almost 1 km in open water) and their dispersal largely depends on tree-hole species [39]. Therefore, they are extremely susceptible to forest conversion [34,39]. Land use change is perhaps the most important factor constricting and/or fragmenting their distribution range, as has been shown in other odonate species [13,54,104]. In addition, climate change may be a more unpredictable threat and could exacerbate the effects of forest fragmentation [39]. This is important because we are barely beginning to understand the patterns of divergence in Megaloprepus, and while we do not know for certain how these populations will respond to land use change, these divergent populations might differ in their response, or in other words, some populations could be more prone to extinction than others. Further research is essential to identify the impact of human activities on these populations and increase conservation efforts.

Conclusions
In conclusion, we found that Megaloprepus populations show some potential ecological divergence, which is in line with previous studies that found genetic and morphological divergence in these populations. This prepares the basis for examining the specific biotic and abiotic factors limiting the distribution of these populations and driving niche divergence. Further studies are also encouraged to (1) disentangle the biogeographic history of Megaloprepus, and (2) evaluate the impact of current anthropogenic land use changes and climate change in the distribution and conservation of Megaloprepus. This will allow a better understanding of the divergence in these populations and predict their future under anthropogenic disturbance.