Environmental Drivers of Diversification and Hybridization in Neotropical Butterflies

Studying how the environment shapes current biodiversity patterns in species rich regions is a fundamental issue in biogeography, ecology, and conservation. However, in the Neotropics, the study of the forces driving species distribution and richness, is mostly based on vertebrates and plants. In this study, we used 54,392 georeferenced records for 46 species and 1,012 georeferenced records for 38 interspecific hybrids of the Neotropical Heliconius butterflies to investigate the role of the environment in shaping their distribution and richness, as well as their geographic patterns of phylogenetic diversity and phylogenetic endemism. We also evaluated whether niche similarity promotes hybridization in Heliconius. We found that these insects display five general distribution patterns mostly explained by precipitation and isothermality, and to a lesser extent, by altitude. Interestingly, altitude plays a major role as a predictor of species richness and phylogenetic diversity, while precipitation explains patterns of phylogenetic endemism. We did not find evidence supporting the role of the environment in facilitating hybridization because hybridizing species do not necessarily share the same climatic niche despite some of them having largely overlapping geographic distributions. Overall, we confirmed that, as in other organisms, high annual temperature, a constant supply of water, and spatio-topographic complexity are the main predictors of diversity in Heliconius. However, future studies at large scale need to investigate the effect of microclimate variables and ecological interactions.


INTRODUCTION
Understanding how the environment shapes species distribution and affects patterns of biological diversity is still a challenging task, especially in species rich regions, such as the Neotropics (Hawkins et al., 2003;Gotelli et al., 2009;Brown et al., 2020). To date, information on this topic is mostly based on vertebrates and plants, and suggest that the combination of high annual temperature with a constant supply of water and spatio-topographic complexity are the main predictors of species distribution, richness, and endemism (Hawkins et al., 2003;Kreft and Jetz, 2007;Qian, 2010;Vasconcelos et al., 2019). Within the Neotropics, the Amazon and the foothills of the North-eastern Andes are examples of regions that combine these conditions, and consequently, they exhibit high levels of species richness and phylogenetic diversity in monkeys, snakes, birds, amphibians, palms, and vascular plants (Kreft and Jetz, 2007;Fenker et al., 2014;Vallejos-Garrido et al., 2017;Velazco et al., 2021). Similarly, regions such as the Biogeographic Choco, Costa Rica, and the Amazon show high levels of phylogenetic endemism (e.g., Rosauer and Jetz, 2014;López-Aguirre et al., 2019;Varzinczak et al., 2020). However, these patterns have not been deeply evaluated in Neotropical invertebrates, and particularly butterflies (Pearson and Carroll, 2001;Mullen et al., 2011).
The environment, and especially climatic niche, has also been suggested to have an effect on gene flow. For example, phylogenetic discordance in multiple loci in beetles of the genus Mesocarabus seems to be the result of hybridization between species sharing the same climatic niche (Andújar et al., 2014), while in armadillos of the genus Dasypus, asymmetric gene flow appears to be facilitated by niche conservatisms at both sides of a geographic barrier (Arteaga et al., 2011). Additionally, climaticbased selection likely plays a role in maintaining mosaic hybrid zones in Quercus oaks, where climatic heterogeneity favors the co-occurrence of parental species and their hybrids (Swenson et al., 2008;Ortego et al., 2014).
Heliconius butterflies are a diverse insect group found across southern United States, Central, and South America, where they occupy divergent habitats (Jiggins, 2017). Due to the recent radiation of this butterfly genus, species pairs have different levels of reproductive isolation, which are used as proxies for different stages of speciation (Kronforst et al., 2013;Martin et al., 2013). In total, ∼25% of Heliconius species are known to hybridize in nature (Mallet et al., 1998(Mallet et al., , 2007, but the role of abiotic variables in facilitating or hampering such hybridization has been poorly studied (Mallet et al., 1990;Rosser et al., 2014).
In this study, we combined an extensive database of occurrences of species and hybrids in Heliconius as well as environmental data to investigate: (1) how the environment shapes the distribution of Heliconius at a regional scale, (2) how the environment molds species richness, phylogenetic diversity, and phylogenetic endemism in these butterflies, and (3) whether niche similarity promotes hybridization.

Species Data and Environmental Variables
We included occurrence data of 46 species of Heliconius and generated a database of the localities where these butterflies have been collected across their entire distribution range. The data were obtained from: (1) entomological collections and (2) the Heliconiinae checklist of Rosser et al. (2012). For those regions in Colombia that we identified as undersampled, we conducted field trips to improve our geographic coverage. The nomenclature of all records was updated to the most recent taxonomic checklist when needed (Lamas and Jiggins, 2017). We also included occurrence data for all interspecific hybrids documented in Heliconius. All individuals were photographed and identified based on their color pattern. We used the point-radius method to georeference specimens with missing coordinates following Wieczorek et al. (2004). Although Heliconius is widely represented in databases, such as global biodiversity information Facility (GBIF), we did not include such records to ensure the use of data that have been curated by specialists both in terms of georeference and taxonomy, or that have images of each specimen that would allow us to confirm the taxonomy.
We used the 19 climatological variables from climatologies at high resolution for the earth's land surface areas (CHELSA) at spatial resolution of 1 km (Karger et al., 2017) to characterize climatic variation across the occurrence range of Heliconius, and altitude was obtained from Jarvis et al. (2008). Collinearity between variables was avoided by estimating the Pearson correlation coefficient among all 20 variables, and the absolute value of this correlation was used to create a dissimilarity matrix (1-correlation values). We used this matrix to perform a hierarchical clustering analysis with the hclust function in R (R Core Team, 2021). We then chose one variable per cluster that had a pairwise distance <0.5. Using the selected variables, we calculated the variance inflation factor (VIF) (Dormann et al., 2013) with the HH package in R (Heiberger, 2020) and chose those variables with VIF <5 (Kubota et al., 2015).

Species Distribution Modeling and Environmental Variables Importance
First, we used R pipelines (Assis, 2020) to reduce sampling bias and spatial autocorrelation among occurrences in our species distribution models using the variables that passed the filters mentioned before. The minimum non-significant autocorrelated distances were used to prune species databases. H. nattereri and H. tristero were not modeled because they had <32 occurrence records.
Then, we generated a second database that included pseudoabsences data following Phillips et al. (2009), Soberón and Nakamura (2009), Barbet-Massin et al. (2012), and Lake et al. (2020). Because Heliconius is a very well-sampled genus we had enough information to select pseudo-absences points for each species in places where: (i) Heliconius other than the focal species have been collected, (ii) environmental conditions may not be optimal for its occurrence, and (iii) absence is not caused by dispersal limitation. Using these criteria, we defined a minimum convex polygon with a 50 km buffer area for each species and selected 10,000 pseudo-absences only in this buffer.
Then, we estimated the ensemble species distribution models (ESDMs) of Heliconius with the R package stacked species distribution models (SSDM) (Schmitt et al., 2017), equally weighting presences and pseudo-absences (prevalence weights = 0.5). Individual species distribution models (SDM) were implemented using four algorithms that optimize the use of pseudo-absences in a similar way (Barbet-Massin et al., 2012): (1) Generalized Linear Models (GLMs) (McCullagh and Nelder, 1989), (2) Generalized Boosting Models (GBMs) (Friedman et al., 2000), (3) Maximum Entropy Models (MAXENT) (Phillips et al., 2006), and (4) Generalized additive model (GAM) (Hastie and Tibshirani, 1990). Each algorithm was run 10 times. In each run, models were calibrated using 75% of the occurrence data and their accuracy was evaluated with the remaining 25%; the "holdout" method was used to ensure independence between training and evaluation sets. The data set randomly changes between runs. An ensemble model (ESDM) was obtained for each species by averaging the best SDM outputs (highest Area Under the Curve-AUC-score), and the ensemble models were evaluated with the AUC score and the Cohen's Kappa coefficient (k). Following Smith and Santos (2020), we did not model species with n < 32 or that occupy >70% of the background region (i.e., entire distribution range for the genus).
We used the relative importance values of the variables provided by SSDM to evaluate the influence of each of them within all models. The importance is estimated with a randomization process, where SSDM calculates the correlation between a prediction using all variables and a prediction where the independent variable being tested is randomly removed; this is repeated for each variable. The calculation of the relative importance is made by subtracting this correlation from one, therefore higher values are the best variables for the model (Schmitt et al., 2017).
Diversity Metrics: Species Richness, Diversity, Endemism Phylogenetic Maps, and Environmental Variables Importance Species richness, phylogenetic diversity, and phylogenetic endemism were calculated by superimposing the distribution maps of all species using the R package phyloregion (Daru et al., 2020b). In order to avoid overestimation of the diversity metrics, we created alpha hulls with the R package rangeBuilder (Davis Rabosky et al., 2016) and following (Paz et al., 2021). Briefly, we used occurrence data available for all species (54,392 georeferenced records) that had more than 10 locality points, a dynamic selection of alpha for each species, and an alpha that varied in steps of 1 (Meyer et al., 2017). We next generated a community matrix using the alpha hulls of all species with the function polys2comm in the R package phyloregion (Daru et al., 2020b).
We used the community matrix to calculate species richness by summing all species present in each cell, and also, with this matrix and the best Maximum Likelihood tree estimated with 20 nuclear and 2 mitochondrial loci for Heliconius (Kozak et al., 2015), we estimated phylogenetic diversity and phylogenetic endemism (Faith, 1992;Rosauer et al., 2009), with the functions phylogenetic diversity (PD) and phylo_endemism of the R package phyloregion (Daru et al., 2020b). To investigate whether these metrics are scale dependent, we performed the above analyses at three consecutive grain sizes (5, 10, and 20 km). We performed a linear regression model using phylogenetic diversity as response variable and species richness as predictor variable to investigate their relationship and plotted the residuals to highlight areas where these metrics are different.
We also used four machine learning algorithms to generate correlative models and then we created an ensemble prediction of each diversity metric to identify the environmental variables that best explain them (Paz et al., 2021). The algorithms used were: Random Forests (Liaw and Wiener, 2002), Neural Network (Venables and Ripley, 2002), Support Vector Machines (Karatzoglou et al., 2004), andGLM (McCullagh andNelder, 1989). The models were built with the R package caret 6.0-86 (Kuhn, 2008), and we used the varImp function to compute the weighted average of the contribution of each variable.

Evaluating the Environmental Effect in the Hybridization on Heliconius Butterflies
We estimated the Schoener's niche equivalency test (D) and Warren's niche background test (I) between pairs of hybridizing species to determine if they share environmental niches. We used the R package humboldt (Brown and Carnaval, 2019) and we followed the concept of environmental niche sensu (Phillips et al., 2006;Soberón and Nakamura, 2009), where the niche consists of the subset of conditions currently occupied and where environmental conditions at the occurrence localities constitute samples from the realized niche. The niche overlap metric Schoener's D ranges between 0 and 1, meaning no overlap and complete overlap, respectively (Rödder and Engler, 2011). The environmental overlap was visualized with a principal component analysis (PCA). We tested the significance of this metric by comparing the realized niche overlap against a null distribution of 1,000 overlaps randomly generated from the reshuffled occurrence dataset and tested whether niche background and niche equivalency were different from those expected by chance at α = 0.05 (Brown and Carnaval, 2019). This was done using the entire distribution of the entities under comparison (niche overlap test = NOT) and using only the area where they overlap (niche divergence test = NDT) (Brown and Carnaval, 2019).
From the species records we discarded 13,476 records as they could not be reliably georeferenced, thus leaving us with 54,392 records. For species modeling, these were further subject to pruning, which left a total of 13,671 records (Supplementary Table 3). There was considerable variation in the sampling effort across the phylogeny. For example, species of the erato and silvaniform clades are well-represented, whereas species from the aoede clade had the lowest number of records (Supplementary Figure 1). The variables retained and used to model species distributions and diversity metrics were: (i) minimum temperature of coldest month, (ii) altitude, (iii) precipitation of coldest quarter, (iv) isothermality, and (v) precipitation seasonality (Supplementary Figure 2). The maximum absolute pairwise correlation between minimum temperature of coldest month and precipitation of coldest quarter was 0.436. The four algorithms we implemented were accurate in predicting the distribution of species, but their combination (ensemble) was the most accurate (Supplementary Figure 3). In total, we generated 44 species distribution models for Heliconius species. These are deposited in ZENODO. 1 We found that environmental variables are better predictors of the distribution of Heliconius compared to topography. For instance, current temperature (isothermality) explains the distribution of 14 species (Figures 1A,B) and precipitation explains the distribution of 24 species (Figures 1C,D). In contrast, altitude explains the distribution of only five species ( Figure 1E). No single variable was correlated with the entire distribution of the genus (Figure 1F), but we observed some general patterns. For example, isothermality explained the distribution of widely distributed species and trans-Andean species (i.e., west of the Andes; Figures 1A,B). Also, precipitation of the coldest quarter explains the distribution of species that occur in the biogeographic Choco + Costa Rica while precipitation seasonality explains the distribution of cis-Andean species (i.e., east of the Andes) + the Pacific of Ecuador (Figures 1C,D). Altitude explains the distribution of species restricted to the eastern foothills of the Andes and highland Andean species (Figure 1E). Interestingly, we did not find a single variable that was better correlated with the distribution of H. charitonia (Supplementary Table 4).
Diversity Metrics: Species Richness, Diversity, Endemism Phylogenetic Maps, and Environmental Variables Importance We found that higher values of Heliconius species richness are concentrated in the foothills of the eastern Andes from Colombia to Ecuador, and into the Amazon basin mainly along the course of the Amazon River (Figure 2A). These results were consistent but more striking in the phylogenetic diversity maps ( Figure 2B). Also, species richness has a strong and significant effect on phylogenetic diversity (adjusted R 2 0.9887, p ≤ 2e-16; Supplementary Figure 4). Interestingly, the residuals map showed values of phylogenetic diversity below those expected from species richness in the same regions, indicating that phylogenetic diversity, although high, is underestimated (blue grids; Figure 2C). In contrast, this metric was overestimated mainly in the Central Andes, the southern Amazon in Brazil, and the northern Chaco in Bolivia (red grids; Figure 2C). The highest values of phylogenetic endemism were concentrated in: (i) the Pacific coast of Costa Rica and Panama, (ii) the central foothills of the Eastern Cordillera in Colombia, and (iii) the biogeographic Choco of Colombia ( Figure 2D). The pattern of these metrics was not scale dependent, and the results were highly congruent at 5, 10 (Supplementary Figures 5, 6, respectively), and 20 km ( Figure 2C).
The ability of the machine learning models to predict species richness, phylogenetic diversity, and phylogenetic endemism varied between algorithms (Supplementary Figure 7). The best algorithms for all diversity metrics were the ensemble model followed by random forest, while the GLM algorithm had the lowest predictive accuracy in all metrics (Supplementary  Figure 7). The best models predicted that altitude and isothermality were the most important variables for species richness and phylogenetic diversity (Figures 3A,B). In contrast, the most important variable for phylogenetic endemism was precipitation seasonality, followed by isothermality ( Figure 3C). Finally, the residuals from the spatial regression between phylogenetic diversity (response variable) and species richness (predictor variable) were explained by isothermality ( Figure 3D).

Evaluating the Environmental Effect on Hybridization in Heliconius
We found 18 pairs of hybridizing species in Heliconius. The results of the NOT and NDT tests based on Schoener's D revealed that the niches of three of these pairs (H. melpomene/H. cydno, H. melpomene/H. hecale, and H. hecalesia/H. hortense) are equivalent (Figure 4 and Table 1) and overlap climatically (D > 0.40). In contrast, 12 of these pairs did not show evidence of niche equivalency. These included both pairs that have extensive geographic overlap (such as H. ethilla and H. numata) (Supplementary Figure 8) and pairs with a narrow overlap (such as H. erato and H. himera) (Figure 5). The remaining three pairs (H. beskei/H. ethila, H. timareta/H. melpomene, and H. charitonia/H. peruvianus) showed inconclusive results (Figure 4 and Table 1). The results of these analyses were deposited in ZENODO (see text footnote 1).

DISCUSSION
We found that Heliconius butterflies display five general distribution patterns, namely: (i) wide distribution, (ii) trans-Andes, (iii) biogeographic Choco + Costa Rica, (iv) cis-Andes + Pacific of Ecuador, and (v) highland Andes. We also found that three variables (isothermality, precipitation and altitude) explain these patterns. Isothermality is a variable that quantifies how daily temperatures oscillate relative to the annual oscillations (O'Donnell and Ignizio, 2012), and its importance as one of the most explanatory variables of species distribution is not without precedent. For example, this variable explains the distribution of frugivorous bats (Chattopadhyay et al., 2019), mealybugs (Heya et al., 2020), Opiliones (Simó et al., 2014), and American monkeys (Vallejos-Garrido et al., 2017). Although all Heliconius species are strongly affected by isothermality, its effect seems to be stronger for widely distributed species and those with trans-Andean distribution. Interestingly, these species occur in regions with high and medium isothermality (>460%), that is, in regions that experience temperature changes throughout the day but keep a constant temperature throughout the year (O'Donnell and Ignizio, 2012). This suggests that these butterflies are particularly sensitive to long term changes in temperature, thus limiting their range to tropical areas. The distribution of species occurring in the biogeographic Choco of Colombia, Costa Rica, cis-Andes and the Pacific of Ecuador is also strongly limited by precipitation. Consistently, these regions have either rainforest, monsoon, or savanna climate, and they are the Neotropical regions with the highest precipitation [precipitation in the driest month (Pdry) > 60 mm]  Kozak et al. (2015), and branches that contribute the most to the phylogenetic endemism are labeled as H1-H5, both in the phylogeny and the map. All maps were plotted in grid cells of 20 km × 20 km. (Beck et al., 2018). Previous studies have suggested that cloudiness and precipitation decrease flying bout duration in butterflies and, consequently, limit their dispersal (Cormont et al., 2011). Therefore, exceptionally high levels of precipitation in such regions may act as population traps, preventing butterflies from flying over longer distances and keeping them in a single region (Rosser et al., 2014). This finding agrees with previous studies in South America, where precipitation shapes the distribution of multiple vertebrates and invertebrates (Atauchi et al., 2017;Amundrud et al., 2018;Schivo et al., 2019;de Oliveira da Conceição et al., 2020).
In addition, altitude was the best predictor for the distribution of Heliconius species that can reach elevations up to 2,600 masl, which is considerably higher than the elevational range occupied by other members of the genus (<2,200) (Rosser et al., 2012). Therefore, it is likely that these highland species have morphological or physiological modifications that allow them to expand their elevational range and occupy new niches. In fact, highland Heliconius are known to have rounder wings compared to lowland species, and this has been suggested to aid them flying dense cloud forest or compensate for the lower air pressure found at higher altitudes (Montejo-Kovacevich et al., 2019). Also, comparisons among different populations of Heliconius have revealed that highland populations are less tolerant to heat (Montejo-Kovacevich et al., 2020), which may limit their distribution range.
The foothills of the eastern Andes and the Amazon basin appeared as the regions with highest Heliconius species richness, which confirms the findings of a previous study done for the subfamily Heliconiine at a higher scale (50 km) (Rosser et al., 2012). Interestingly, both of these regions are known to present unusual concentrations of contact zones and hybrid zones (i.e., suture zones) (Dasmahapatra et al., 2010;Rosser et al., 2021), which may explain the richness they exhibit. Also, altitude, isothermality, and precipitation were the variables best correlated with this metric. This may be due to the elevational gradient found at the foothills of the eastern Andes, which offers multiple ecological niches thus favoring diversification rates (Rahbek and Graves, 2001;Jetz and Rahbek, 2002;Davies et al., 2007;Keppel et al., 2016). Additionally, there are several climate-based hypotheses that seek to explain broad-scale diversity patterns, and water and energy have emerged as crucial influencers of species richness (Silva-Flores et al., 2014). In particular, the water-energy dynamics hypothesis argues that species richness increases in places where liquid water and optimal energy conditions provide the greatest capacity for biotic dynamics FIGURE 4 | Co-occurring and hybridizing species of Heliconius. Green: species pairs with equivalent environmental niches, blue: species pairs with divergent environmental niches, and salmon: species pairs with inconclusive results. Numbers indicate the pairs of species falling into each category. (Svenning et al., 2008). The Amazon and foothills of the eastern Andes are regions with near constant hot-warm temperature throughout the year and have a permanent liquid water supply (Rosser et al., 2014;Vallejos-Garrido et al., 2017) thus ensuring an optimal water-energy dynamic. The latter translates into constant availability of plants for butterflies, including host plants for immature and pollen for adults, and continual interactions between individuals, which may be correlated with the high species richness we detected.
Similar to other studies, patterns of phylogenetic diversity were similar (although not identical) to richness (Davies Jonathan and Buckley, 2011;Fenker et al., 2014;Mendoza and Arita, 2014;Guedes et al., 2018). Interestingly, areas with highest species richness got low phylogenetic diversity (Figure 2C, blue grids), which may be a consequence of the recent increase in diversification rate in Heliconius (4.5 Ma) and the consequent cooccurrence of multiple young species in the Amazon and foothills of the eastern Andes (Rosser et al., 2012;Kozak et al., 2015). In agreement with this observation, previous research in both animals and plants have found high phylogenetic diversity in the eastern Andes of Colombia, Peru, and Ecuador (Fenker et al., 2014;Mendoza and Arita, 2014;Guedes et al., 2018;Arango et al., 2021;Velazco et al., 2021).
The highest phylogenetic endemism was found in the central eastern Andes of Colombia, and this result is possibly due to the restricted range of the species Heliconius heurippa ( Figure 2D, area H1). However, we cannot rule out this result as an overestimation since the phylogenetic tree that we used (Kozak et al., 2015) considers this taxon as a separate species and not as part of H. timareta (as recently hypothesized). If H. heurippa had been included within H. timareta, which has a wider distribution range, it is likely this result on phylogenetic endemism does not hold. Additionally, the pacific region of Costa Rica, Panama and Colombia show intermediate values of phylogenetic endemism that resulted from the presence of species that have reduced geographic range and are either longbranch species (e.g., Heliconius godmani) or species for which no close relatives are known (e.g., Heliconius hewitsoni) (Figure 2D, area H2 and H3, respectively). These regions were previously described as highly endemic phylogenetically for plants (Sandel et al., 2020), terrestrial mammals (Rosauer and Jetz, 2014), birds and amphibians (Daru et al., 2020a). Interestingly, there were several species that, although are considered as geographically endemic within Heliconius, exhibited low values of phylogenetic endemism. However, it is important to acknowledge that phylogenetic endemism is a concept based on linages rather than species, and thus, if an endemic species has a narrow range but it is closely related to a widespread species, its phylogenetic endemism will not necessarily be low (Rosauer et al., 2009). An example of this is Heliconius nattereri, an endemic species from Brazil's Atlantic Forest that, despite having a narrow distribution, is sister to the widely distributed Heliconius ethilla ( Figure 2D, area H4). Similarly, Heliconius atthis is restricted to the Ecuadorian and Peruvian Pacific, but it is sister to the widely distributed Heliconius hecale (Figure 2D, area H5). In our study we found that high precipitation and near constant hotwarm temperature throughout the year are strongly correlated with phylogenetic endemism, which agrees with studies that point a role for temperature in promoting endemism by reducing extinction rates and increasing population sizes in small areas (Jetz et al., 2004;Rosauer and Jetz, 2014;Varzinczak et al., 2020).
Our environmental niche analysis showed that hybridizing species do not necessarily share the same climatic space despite some of them having largely overlapping geographic distributions. This is the case of H. ethilla and H. numata, which frequently co-occur throughout their distribution, but there are some regions with an extreme climate, such as the Pacific coast of Colombia (a humid jungle) and the Colombian Magdalena valley (which has a marked precipitation gradient, being humid in the north and dry in the south), where H. ethilla but not H. numata occur (Supplementary Figure 9). This suggests that the former species has a broader climatic tolerance. We also detected differences in the environmental niche between pairs of hybridizing species that rarely overlap geographically, but when they do, they hybridize. For example, H. erato and H. himera occupy contrasting environmental niches in Ecuador (Jiggins et al., 1997), where H. himera lives in dry forests while H. erato inhabits wet forests of the Andes (Figure 5). Similarly, the hybridizing H. erato (H. e. venus) and H. chestertonii meet in an environmental transition zone between wet and dry forest in the Colombian Andes (Muñoz et al., 2010 ; Supplementary Figure 8).
In summary, we confirmed that, at large scales, the distribution of Heliconius, its richness, diversity, and phylogenetic endemism are mainly shaped by a combination of high annual energy (i.e., hot-warm temperature), constant water supply, and an extraordinary topographic complexity. However, species distributions are thought to result from dynamics occurring at multiple spatial scales. Therefore, including microclimate variables and ecological interactions would provide an in-depth understanding of the multiscale drivers of distribution, niche range and phylogenetic processes (Montejo-Kovacevich et al., 2020;Paz and Guarnizo, 2020). Our study confirms the richness and diversity of areas already identified in other taxa, thus strengthening the importance for their conservation as strategic hotspots of biodiversity.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://doi.org/10. 5281/zenodo.5149294.