Population genomics reveals differences in genetic structure between two endemic arboreal rodent species in threatened cloud forest habitat

Genomic tools are now commonly used to assess the genetic diversity and genetic structure of species and populations, and they provide the ability to describe and address the negative effects of population declines and fragmentation. However, such studies are lacking for arboreal mammals despite their contribution to various ecosystem services, especially in uncommon and critically endangered ecosystems such as cloud forests. The aim of this work was to evaluate and compare the genetic diversity and population structure of two endemic arboreal mice from Mexican cloud forests that are associated with areas with different levels of impacts from human activities. We performed genotyping-by-sequencing in 47 Habromys schmidlyi and 17 Reithrodontomys wagneri individuals to evaluate genetic diversity and differentiation. In both species, the genetic diversity was low compared to other cricetid species, and we observed different population structure patterns, potentially linked to the different ecological associations. We detected two genetic groups in H. schmidlyi, that is a territorial species present in areas of low incline, while a single genetic group was found in R. wagneri, which forms family groups in areas with steep slopes. Overall, these results highlight how species’ genetic diversity can be differentially impacted depending on differential ecological associations within the same ecosystem. This information is essential for the development of the adequate conservation and management of these species.


Introduction
It is well known that the current Anthropocene period is associated with continued biodiversity loss at all ecological levels, from genetic diversity to ecosystem integrity and functionality (Ceballos et al. 2010;Otto 2018;Turvey and Crees 2019). The Convention on Biological Diversity (CBD), in their 2020 global framework, stated a need to minimize biodiversity loss in three aspects: genetic, species, and ecosystem diversity (Hoban et al. 2020). The CBD also highlighted the need to focus on the state of conservation and environmental risks in places where threatened species occur, mentioning the desirability of information on genetic diversity -important for conservation decision-making at both national and international levels (Darling et al. 2017;Breed et al. 2019;Hoban et al. 2020). Genomic tools offer a powerful approach to assess the status of species and populations, as these tools can reveal ecological and demographic processes that shape genetic structure across the entire distribution of taxa, which is useful for conservation practitioners when dealing with threats to populations, species, and even ecosystems (Allendorf et al. 2010;Ekblom and Galindo 2011;Hoban et al. 2013;Hunter et al. 2018;Garrison et al. 2021). More importantly, these genomic tools can be used in non-model organisms with no prior genomic information or with low probability of detection or sampling -reducing the biases that may arise from small sample sizes, frequent in endangered species -since genome markers such as SNPs are not sensitive to these sampling effects (Ekblom and Galindo 2011;Hunter et al. 2018;Breed et al. 2019;Garrison et al. 2021). In this context, a set of four essential biodiversity variables have been proposed for genetic composition monitoring: (i) genetic diversity, (ii) genetic differentiation between populations, (iii) inbreeding, and (iv) effective population size (Hoban et al. 2022).
In terms of the health of ecosystems, it is also considered important to improve biodiversity protection (Hoban et al. 2020). Cloud forests are considered a relict ecosystem, which is threatened by climate change, as well as by habitat and land-use changes. They now represent less than 0.5% of Mexican territory, in isolated patches of less than 8 000 ha, and is the result of natural fragmentation of the ecosystem -due to the emergence of the Trans-Mexican Volcanic Belt, which started early-to mid-Miocene, with the last stratovolcanoes arising during the last 1.5 Myr (Ferrari et al. 2012) -but also by the replacement of half of its original distribution by human-induced landscapes ( Fig. 1) (Challenger 1998;Luna-Vega et al. 1999;Ruiz-Jiménez et al. 2012). Consequently, the increasing patchiness of this ecosystem has affected the distribution of species by creating smaller and more isolated populations and reducing home-ranges due to the scarcity of resources and increased interspecific competition, particularly in species with low dispersal ability and probability of detection (Castañeda-Rico et al. 2011;Leonardi et al. 2012;Marines-Macías et al. 2018).
In cloud forest, arboreal rodents (together with other arboreal mammals) play a key-role in the ecosystem as seed dispersers and pollinators (Romo-Vázquez et al. 2005;Vega et al. 2017;Marines-Macías et al. 2018), and these species can also help to enhance community diversity by exploiting canopy resources, leaving ground-level resources for non-arboreal animals (Macarthur et al. 1966;August 1983). However, these arboreal species are facing strong threats due to their extreme habitat specialization, as well as low dispersal ability and detectability during sampling . Highly specialist species are more prone to extinction due to habitat loss and fragmentation in human-modified landscapes, limiting migration and gene flow between populations and increasing inbreeding and between-population genetic divergence, which promote the reduction of effective population size and loss of genetic variability (Reed and Frankham 2003;Banks-Leite et al. 2014;Millán-Aguilar et al. 2016).
One of the largest patches of cloud forest in Mexico, and the second in priority of conservation according to national authorities, is the region of the "Upper Basin of the Amacuzac River" (CONABIO 2010;Ruiz-Jiménez et al. 2012). Within this region, three arboreal rodent species have been recorded, the Schmidlyi's crested-tailed deer mouse (Habromys schmidlyi), the Wagner's harvest mouse (Reithrodontomys wagneri), and the Goldman's diminutive woodrat (Nelsonia goldmani) (Pardiñas et al. 2017;Marines-Macías et al. 2018). Among these, it has been reported that H. schmidlyi and R. wagneri are sensitive to the scarcity of suitable substrates for nesting, so they depend on the local availability of arboreal substrates for the construction of their nests (Marines-Macías et al. 2018). In this way, these and other arboreal mice are useful indicator species of ecosystem change (Clavel et al. 2011;Linnell et al. 2018). Therefore, evaluating their genomic diversity could be beneficial for future conservation and management plans.
Even though the arboreal H. schmidlyi and R. wagneri are sympatric in some Mexican cloud forests, they occupy and use resources in different ways. Habromys schmidlyi shows strong association with cloud forest tree species in areas of low incline (e.g., Cleyera velutina and Carpinus caroliniana), together with Quercus sp. in areas surrounding the cloud forest patches (Marines-Macías et al. 2018 . pdf], even though conservation status has not been formally evaluated for R. wagneri.
All these ecological and conservation characteristics make H. schmidlyi and R. wagneri interesting systems to explore how the natural fragmentation of cloud forest ecosystems translates into differences in genetic diversity and population structure. For these, we expect that H. schmidlyi will show higher levels of genetic structure because of its preference for isolated areas of low incline, whereas R. wagneri will show a weaker genetic structure due to the connectivity among steep slopes. According to this, the aim of this study was to evaluate the genomic diversity and population structure in these two ecologically distinctive arboreal rodent species living in sympatry. This study also aims to contribute to the development of conservation guidelines for the protection of cloud forest mammalian biodiversity, and our results could serve as a basis for further studies focused on the management and conservation plans for different species associated with the cloud forest, but also with other fragmented ecosystems.

Sampling and DNA extraction
During 2012-2013 extensive fieldwork was conducted at the "Cerro del Huixteco" National State Park, within the Sierra de Taxco, Guerrero, Mexico (18° 36′ N 99° 36′ W), and the "Picacho de Oro y Plata" State Park, Zacualpan, Estado de México, Mexico, and surrounding areas (18° 43′ N, 99° 46′ W) (Fig. 1), where the study species are sympatric. Sampling at both locations was performed using 120 Sherman live-traps, with 60 traps placed at the ground level and the remaining 60 traps on tree trunks and branches, following the transect set on the ground. The arboreal trapping was done using climbing equipment to reach the highest possible point (20 m) [for further sampling information, see Marines-Macías et al. (2018)]. All trapped individuals were identified in situ using external characteristics, and ear clip samples were obtained using a medical punch previously disinfected with 10% chlorine and 96% ethanol. Tissues were fixed individually in 96% ethanol and then stored at − 70 °C until processing. DNA extractions were performed with the Mag-Bind® Blood & Tissue DNA OMEGA bio-tech kit (Omega Bio-tek, Inc., Norcross, US), following the manufacturer's protocols with some modifications in incubation temperatures and time to increase DNA yield (Supplementary 1).

Genotyping-by-sequencing (GBS)
DNA was quantified with a Qubit 2.0 Fluorometer (dsDNA Broad Range Assay) for a total of 54 H. schmidlyi and 21 R. wagneri, and the samples normalized for standard concentrations required for high-throughput sequencing, using an average of 40 μl sample volumes and concentrations between 20 and 50 ng/μl. Samples were then sent to the Genomic Diversity Facility (GDF) at Cornell University for GBS library construction (Elshire et al. 2011) and sequencing services. Briefly, DNA was sheared using the PstI enzyme for genome complexity reduction, including one blank as a negative control. Individual samples were barcoded and then combined in a single pool, which was amplified through PCR reaction, purified, and quantified for sequencing on the Illumina Hi-seq 2000/2500 sequencer (Illumina Inc., USA).
Due to lack of a reference genome for both species, SNP calling was carried out with the UNEAK pipeline for each species independently (Lu et al. 2013), defining a minimum call of five reads, maximum locus and individual missing data of 20% and 25%, respectively, following the methodology of Barbosa et al. (2017). Postprocessing included further filtering in order to avoid repetitive and paralogous loci (heterozygosity above 75%) following the pipeline of White et al. (2013).
Finally, a principal coordinates analysis (PCoA) using the ade4 package (Dray and Dufour 2007) was performed individually for each rodent species to test for the genetic variation of the SNPs obtained and their distribution among populations in a model-free framework.

Genetic diversity
We used the R package diveRsity (Keenan et al. 2013) to evaluate the fit to the Hardy-Weinberg equilibrium within each sampling location using Fisher's exact test, by calculating the rate of fixation or inbreeding (F) applying a Bonferroni correction for multiple comparisons (Rice 1989).
The SNP variation and its distribution among populations based on a model-free framework were tested using the same R package, estimating the rarefied allelic richness and private alleles, observed (Ho) and expected (He) heterozygosity, allowing 5% missing data and using 10,000 global iterations.

Genetic structure
To examine the genetic structure of each species, the software STRU CTU RE (Pritchard et al. 2000) and STRU CTU RE HAR-VESTER (Earl and vonHoldt 2012) were used to determine the most probable number of populations by calculating the ΔK value with the Evanno method (Evanno et al. 2005). STRU CTU-RE analyses were conducted testing from 1 to 5 populations (K), with the admixture model and correlated allele frequencies. We used an initial burn-in of 200,000, followed by 800,000 Markov Chain Monte Carlo iterations with 5 independent runs for each K value.
The geographic structure by species was also tested in Geneland (Guillot et al. 2012), with this analysis incorporating the maximum distance that an individual can move as an uncertainty value, together with the geographic position of each sample -information that was obtained according to home-range areas and movement distances in Marines-Macías et al. (2018). For this, we followed the recommendations of Guillot et al. (2012) for processing the genomic data, so we created a random subset of both datasets: four subsets for H. schmidlyi and three for R. wagneri. In each subset, the occurrence of different populations from 1 to 5 was tested, using 1 million generations with a thinning of 100 and a burn-in of 1,000. The convergence of the results over five independent runs and all subsets was determined, and the run with the highest likelihood value was chosen.
Additionally, hierarchical genetic structure was evaluated by analysis of molecular variance (AMOVA) in Arlequin 3.5.2.2 (Excoffier and Lischer 2010). To perform this analysis, both sampling sites were tested as the main groups, while the population identity of each individual was based on closeness or correspondence to overlapping home-ranges areas, according to results for each species obtained by Marines-Macías et al. (2018), which considers that these overlapping-home-ranges represent subpopulations. These analyses were performed using pairwise distance matrices (Fst) with 10,000 permutations to test the significance of our results.
To evaluate the distribution of the individual genetic distances of each species, Nei's distances were calculated between individuals in the adegenet R package (Jombart and Ahmed 2011). Additionally, the effective population size (Ne) per population by species were estimated in the LDNe 1.31 software (Waples and Do 2008), with a random mating model, a minor allele frequency of 5% and Jackknife with parametric corrections for independence of loci.
To test for isolation-by-distance, a Mantel test was performed correlating the geographic and genetic distances. The genetic distance matrices were calculated between each sampling site (the same Fst pairwise comparisons of populations as used in the AMOVA analysis) in GenAlEx v6.503. Both Mantel tests were performed using 9,999 permutations, and as a measure of population distance, Fst / (1-Fst) was used, following the recommendations of Rousset (1997).
Finally, the effective number of migrants for both species (Nm) was calculated by the private alleles method proposed by Barton and Slatkin (1986) with a correction for population size, using Genepop 4.4 (Rousset 2008).

Results
Of the 54 individuals of H. schmidlyi, a total of 40,665 SNPs were obtained with the UNEAK pipeline, of which 4,006 were within the defined thresholds of missing data. Similarly, of the 21 individuals of R. wagneri, a total of 60,137 SNPs were obtained, of which 2,895 fulfilled the same filtering thresholds.
According to the PCoA exploratory analysis, a clear differentiation between the two populations of H. schmidlyi was observed, indicating two well-defined genetically distinct groups. These groups corresponded to the Zacualpan and the Taxco populations. In contrast, no genetic differentiation was found between the two sampling sites for R. wagneri (Fig. 2).

Genetic diversity
The rate of fixation or inbreeding (F) for the 4,006 loci obtained from H. schmidlyi, were significantly different from zero in 925 loci (23%), indicating an excess of homozygotes. However, the paired value of the fixation index after Bonferroni correction (F = 0.01, P = 0.0001) indicates that both the Zacualpan and Taxco populations were in Hardy-Weinberg equilibrium, so we include all loci in further analyses (Supplementary 2). Similarly, from the 2,895 loci obtained from R. wagneri, 212 (7%) showed an excess of observed homozygotes, and the paired value of the fixation index after Bonferroni correction (F = 0.0083, P = 0.0002) showed that these populations were also in Hardy-Weinberg equilibrium (Supplementary 2).
In general, the Ho/He values were similar for both species when the populations of Taxco and Zacualpan were compared (Table 1). Furthermore, for H. schmidlyi, a rarefied mean allelic richness of 1.86 and a rarefied private allelic richness of 0.08 were obtained in the Taxco population, compared with 1.83 and 0.13 obtained in Zacualpan. On the other hand, for R. wagneri, a rarefied mean allelic richness of 1.73 and a rarefied private allelic richness of 0.17 were obtained in the Taxco population, with 1.77 and 0.08 for the Zacualpan population (Supplementary 3).

Genetic structure
The ΔK values of STRU CTU RE analyses indicated that the most likely value of K was 2 for both species (Fig. 3). Considering these K values, the STRU CTU RE analysis showed a clear division of two geographic clusters for H. schmidlyi ( Fig. 3A; Table 2). However, for R. wagneri, there was no clear evidence of two geographic clusters ( Fig. 3C; Table 2), since only two individuals from the Zacualpan population showed high values of belonging to a different genetic population; this suggests a lack of genetic subdivision for this species. Additionally, when the Geneland analysis are considered, there was evidence of genetic structure in H. schmidlyi, inferring the same two genetic clusters corresponding to the two sampled localities, Taxco and Zacualpan (Fig. 3B). Meanwhile, for R. wagneri, the Geneland analysis revealed a single genetic cluster (Fig. 3D).
The AMOVA analyses showed that the highest variation was within populations (i.e., the a priori identity of population membership according to the criterion of overlapping home-ranges [see Marines-Macías et al. (2018)]), and among groups (i.e., Zacualpan vs Taxco) for both species. In the case of R. wagneri, the variation within populations was higher (93.69%) compared with H. schmidlyi (85.31%), but the variation among groups was lower (4.03%) compared to H. schmidlyi (11.15%) ( Table 3). The AMOVA revealed significant variation for all F-statistics for H. schmidlyi, but only among populations (Fct = 0.04; P = 0.007) for R. wagneri (Table 3).
Nei's distance values show a low differentiation for both species. However, when observing the distribution of the distances in R. wagneri, a unimodal distribution is observed, which could indicate less variation between the comparisons, unlike H. schmidlyi in which there is a bimodal distribution, which suggests a greater variation between comparisons (Fig. 4). The effective population sizes were low for both species in both populations, with an estimate of 20.8 individuals in Taxco and 10.2 in Zacualpan for H. schmidlyi and 16.9 individuals in Taxco and 9.3 in Zacualpan for R. wagneri (Table 1).
Only for H. schmidlyi was geographic distance positively correlated with genetic distance (R 2 = 0.679, P = 0.011), but in the case of R. wagneri the correlation was close to zero and not statistically significant (R 2 = 0.288, P = 0.169). Supporting these results of isolation-by-distance, the number of migrants after size correction between Taxco and Zacualpan was lower (0.221 ± 0.15 migrants) for H. schmidlyi than for R. wagneri (0.462 ± 0.175).

Genetic diversity
The genetic diversity shown by both species seems to be lower than other endemic Mexican rodents such as Oryzomys couesi cozumelae (He = 0.690-0.864; Vega et al.  Espindola et al. 2014). It seems also low compared with other mammal species distributed in fragmented areas and endangered habitats such as those in the Australian fragmented Eucalyptus forests, Antechinus agilis (He = 0.844-0.86; Banks et al. 2005) and Petauroides volans (He = 0.59-0.77; Taylor et al. 2007). Since not all the above He values were obtained using SNPs, future analyses using homologous loci between species are needed to corroborate if genetic diversity of H. schmidlyi and R. wagneri are truly lower. Given the threatened status of the Mexican cloud forests, anthropogenic disturbances in the habitat could also negatively affect the connectivity between the patches, promoting the decrease in the genomic diversity of these two arboreal rodents, as has been documented in other animals inhabiting protected areas (Mohammadi et al. 2018;Fetene et al. 2019;Garrison et al. 2021). Furthermore, it has been documented that life history and natural demographic characteristics of populations could play an important role in influencing low levels of heterozygosity (Garrison et al. 2021). In this sense, the low levels of heterozygosity shown by other microendemic mammalian species that inhabit fragmented habitats have been attributed to the formation of family groups (Frère et al. 2010;Castellanos-Morales et al. 2014;Espindola et al. 2014). This specific demographic characteristic has been found in R. wagneri in both populations, where family groups were found using radiotracking information, although in the case of H. schmidlyi this characteristic was only documented in Zacualpan (Marines-Macías et al. 2018).

Genetic differentiation
Given the small number of samples for R. wagneri, the results must be taken with great caution; however, they can be interpreted as a trend. Taking into account the fact that the two sampling locations are separated by 22.7 km, which can be considered enough to expect high genetic differentiation in rodent species (Trizio et al. 2005;Bogdanov et al. 2009;Vargas et al. 2012;Castellanos-Morales et al. 2014), a genetic difference was detected between the H. schmidlyi populations, but was not detected in R. wagneri. While the clustering analyzes (Geneland and STRU CTU RE) showed a "highly" structured individual identity for H. schmidlyi, the AMOVA and Nei's distances suggest that this genetic structure is supported by low levels of genetic differentiation. This low genetic differentiation within populations of H. schmidlyi could be expected due to the recent divergence of this species from its sister taxon H. simulatus, between 1.63 and 2.33 Mya (Leon-Paniagua et al. 2007), and the fact that other Neotominae needed substantial time to reach high levels of genetic differentiation (Steppan and Schenk 2017).
Considering the habitat selection previously reported for H. schmidlyi, this species prefers areas of low incline (Marines-Macías et al. 2018), and the presence of steep slopes that divide the populations of Zacualpan and Taxco could limit the number of migrants, as shown in this study, favoring a greater genetic structure between populations. Interestingly, geographic features also shape the genetic structure of another co-generic species, H. simulatus, in which isolated populations act as genetic islands (Castañeda-Rico et al. 2011). We additionally detected an isolation-bydistance pattern in H. schmidlyi that could also be limiting the number of migrants between these two populations. So, the geography of the landscape and the distance could be promoting the subdivision of H. schmidlyi. Contrary to what was observed in H. schmidlyi, R. wagneri prefers trees located in steep slopes, habitat that could facilitate connectivity and gene flow along the Sierra de Taxco, where this species was sampled. This habitat selection could explain the lack of genetic structuring in R. wagneri, the higher number of migrants compared with H. schmidlyi, and the lack of isolation-by-distance observed. In summary, these results suggest that the steep slopes, which allow gene flow in R. wagneri, represent a local geographic barrier for H. schmidlyi, which prefers areas of low incline.
Additionally, different behaviors have been reported between the two populations of H. schmidlyi. Territorial behavior has been observed in Taxco: individuals have separate home ranges and males tend to expand that cruising area in the reproductive season, whereas in Zacualpan a gregarious social structure has been reported (Burt 1943;Coleman and Downs 2010;Marines-Macías et al. 2018). This supports the distinction of the two H. schmidlyi populations that are genetically, behaviorally, and ecologically different and also agree with the results found in R. wagneri: in this species, a gregarious social structure was reported in both Taxco and Zacualpan, in which females and males share the same areas (Burt 1943;Coleman and Downs 2010;Marines-Macías et al. 2018); showing a continuity and unity between sampled localities.
When considering estimated effective population sizes, these species show low numbers compared to other rodent species, such as Reithrodontomys raviventris (Ne = 15-131;Statham et al. 2022) or Cricetus cricetus (Ne = ~ 6000-7000; Melosik et al. 2017). However, it is important to note that rodent generations overlap due to their short life cycles, violating the sampling assumption of non-overlapping generations. As such, our Ne estimates could be lower than reality, since it is very likely that several age classes were sampled (England et al. 2006;Marandel et al. 2019;Garrison et al. 2021). Therefore, the Ne estimates reported in this study should be used with caution, and we strongly recommend further age-class cohort samplings to elucidate the actual trends over time, following the suggestions of England et al. (2006) and Marandel et al. (2019).

Implications for conservation
Habromys schmidlyi and R. wagneri are under special protection by international agencies or the Mexican National law. Both authorities agree that the main concerns are the low number of individuals and the declining trend of the populations (DOF 2010;Álvarez-Castañeda et al. 2018). We observed that H. schmidlyi was affected by local geographic features, so the continuity and connection of forest patches that this species inhabits needs to increase to avoid negative processes that could decrease genetic diversity. We believe our results highlight the importance of both H. schmidlyi populations as unique evolutionarily significant units, since they represent unique sources that should be conserved to protect the species, as well as contributing to overall conservation and management plans (Moritz 1995;Castañeda-Rico et al. 2011;Eizaguirre and Baltazar-Soares 2014). In the case of R. wagneri, the two sampling localities analyzed represent the same genetic unit for this taxon, and even though the number of migrants is low, the conservation of the natural habitats located in steep slopes is also important to maintain the genetic connectivity and avoid the negative effects of genetic drift. The conservation and management of these two species must be focused on maintaining the health and continuity of the cloud forest, which also will play an important role in the conservation and the continuity of the ecological processes in it (Mitrovski et al. 2007;Eizaguirre and Baltazar-Soares 2014;Marines-Macías et al. 2018). Furthermore, future monitoring is necessary to evaluate the conservation risk of both species in order to protect their genetic diversity (Eizaguirre and Baltazar-Soares 2014).
In this work we analyzed all known populations of H. schmidlyi and compared them with sympatric populations of R. wagneri. We showed how the complex topography in the region, which includes multiple interconnected steep slopes that isolate multiple areas of low incline, promotes different evolutionary responses in Neotominae rodents. Although both species belong to the Reithrodontomyini tribe, their contrasting social interactions and habitat selection seems to lead to a differential impact on the connectivity of their populations. This affects their genetic diversity and genetic structure and highlights that the idiosyncrasy of each species determines their evolutionary responses. Likewise, this study showcases a pervasive phenomenon across microendemic species with low population density, low detection probability, and associated with endangered ecosystems, such as Peromyscus melanurus, Reithrodontomys brevirostris, and Nelsonia goldmani . Hence, future research should also address the entire distribution of R. wagneri and be especially concerned with the conservation of other microendemic species and the habitats that they inhabit.
Finally, the information shown here is of great importance to generate management plans and biodiversity conservation strategies not only for the species themselves, but also to be used as a proxy for evaluating the health of the ecosystem, since these arboreal rodents have key roles in the maintenance of cloud forest as seed dispersers, pollinators, and ecosystem regenerators (Schupp 1993;Williams et al. 2000;Vander-Wall et al. 2005).

Supplementary Information
The online version contains supplementary material available at https:// doi. org/ 10. 1007/ s13364-022-00667-x. Funding This work was supported by grants of Consejo Nacional de Ciencia y Tecnología (CONACyT) 239482 and Programa de Apoyo a Proyectos de Investigación e Innovación Tecnológica (UNAM-PAPIIT) 216713. There are no financial or commercial conflicts of interest. Pablo Colunga-Salas and Tania Marines-Macías were master's students at the time of study, based at the Posgrado en Ciencias Biológicas, UNAM, and were supported by a fellowship from CONACyT and Apoyo a los estudios de Posgrado (UNAM-PAEP). Giovani Hernández-Canchola was supported by the Programa de Becas Posdoctorales de la Dirección General Asuntos del Personal (DGAPA, UNAM), the US National Science Foundation (DEB-1754393), and CONACyT.

Conflict of interest The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.