Genomic differentiation and niche divergence in the Hetaerina americana (Odonata) cryptic species complex

The evolution of reproductive barriers, that is, the speciation process, implies the limitation of gene flow between populations. Different patterns of genomic differentiation throughout the speciation continuum may provide insights into the causal evolutionary forces of species divergence. In this study, we analysed a cryptic species complex of the genus Hetaerina (Odonata). This complex includes H. americana and H. calverti; however, in H. americana two highly differentiated genetic groups have been previously detected, which, we hypothesize, may correspond to different species with low morphological variation. We obtained single nucleotide polymorphism (SNP) data for 90 individuals belonging to the different taxa in the complex and carried out differentiation tests to identify genetic isolation. The results from STRUCTURE and discriminant analysis of principal components (DAPC), based on almost 5000 SNPs, confirmed the presence of three highly differentiated taxa. Also, we found FST values above 0.5 in pairwise comparisons, which indicates a considerable degree of genetic isolation among the suggested species. We also found low climatic niche overlap among all taxa, suggesting that each group occurs at specific conditions of temperature, precipitation and elevation. We propose that H. americana comprises two cryptic species, which may be reproductively isolated by ecological barriers related to niche divergence, since the morphological variation is minimal and, therefore, mechanical barriers are probably less effective compared to other related species such as H. calverti.

recorded (Sánchez-Guillén et al., 2014a).Reproductive isolation in this group usually results from sexual and mechanical barriers (Sánchez-Guillén et al., 2012;Svensson et al., 2006) such as morphological traits associated with the mating recognition system (e.g., wing and body coloration and the shape of caudal appendages).
These barriers usually represent about 90% of the total reproductive isolation, and their evolution is associated with sexual selection driven by female choice (Sánchez-Guillén et al., 2014b).For this reason, and in addition to their high degrees of niche conservatism, the speciation process in zygopterans is primarily driven by sexual selection (Wellenreuther et al., 2012;Wellenreuther & Sánchez-Guillén, 2016).Nevertheless, some cases of speciation have also been described (Damm et al., 2010;Jordan et al., 2003;McPeek et al., 1996), where factors such as the presence of predators, temperature and waterbody type drive the evolution of reproductive isolation through ecological barriers.
Studies describing genomic differentiation and local adaptation patterns in Odonata are scarce (Dudaniec et al., 2018), which is fundamental to understand the speciation process in species complexes with a low degree of morphological differentiation.Here, we analyse the genomic differentiation patterns of a cryptic species complex of the genus Hetaerina intending to contribute to the understanding of reproductive isolation evolution in odonates.
Hetaerina (Calopterygidae) is a mainly Neotropical genus of damselflies, commonly known as rubyspots, comprising 39 species (Paulson et al., 2023).Rubyspots are characterized by territorial behaviour and high degrees of behavioural interference when two or more species occur in sympatry (Drury et al., 2015).The most emblematic member of the genus is H. americana s. l. (the American Rubyspot), which has a wide distribution from Nicaragua to Canada inhabiting a great variety of ecosystems such as temperate forests, tropical dry forests, and desert scrubs.However, Vega-Sánchez et al. (2019) reported the existence of three differentiated population groups within this taxon, based on nuclear genetic and morphological data and suggested that these groups may correspond to different species.The first group is distributed from the United States to the north of Mexico (H.americana North, hereafter), a second group is mainly distributed in the center and south of Mexico (H.americana South, hereafter), and a third group occurs from northern Mexico to Guatemala.This last group was later described as a new species, H. calverti (Vega-Sánchez et al., 2020).
Additional analyses indicated that H. calverti is reproductively isolated from the other two taxa by mechanical barriers related to the shape of caudal appendages (Vega-Sánchez et al., 2022).However, morphological variation in these structures is less conspicuous between H. americana North and South, suggesting that the three entities represent different stages of the speciation continuum and that other isolation mechanisms besides mechanical barriers could be present.Interestingly, H. americana North and South have been found in sympatry with H. calverti but there is no record of sympatric localities for H. americana North and South.Also, evidence of purifying selection on the mitochondrial gene cytochrome oxidase I was found, possibly related to the gradient of temperature found along the whole distribution of the three taxa (Vega-Sánchez et al., 2019) which suggests that environmental variation also affects the genetic differentiation patterns in this species complex, especially between H. americana North and South (Vega-Sánchez et al., 2019).
Based on this, we hypothesized that, if the taxa present effective reproductive isolation, high genomic differentiation values will be found throughout the whole genome and that, in the case of H. americana North and South, niche variation throughout the distribution will contribute to the reproductive isolation more than mechanical barriers.
To test this, in this study, we implemented population genomic analyses to unravel the evolutionary history of the American Rubyspot species complex and address the following questions: (1) How many species can be recognized based on the genomic data?(2) Is there evidence of a contribution of environmental variation to the genetic differentiation patterns?(3) Are neutral processes the main driving forces involved in the speciation of this group?

| Sampling and sequencing
We analysed a subsample of 95 individuals from Vega-Sánchez et al. (2019).Twenty-five individuals belonged to H. americana North, 39 to H. americana South, and 31 to H. calverti.These individuals were sampled from different populations (8-10) for each of the three taxa (Figure 1a).
We extracted genomic DNA from the thoracic muscle of each individual using the Pure Link Genomic DNA mini kit (Invitrogen) and followed the manufacturer's protocol.The DNA was quantified using a Qubit fluorometer system (High Sensitivity DNA kit, Life Technologies), and the samples were standardized to a final concentration of 40-50 ng/μL of DNA.Library preparation and sequencing were performed using the 3RAD method (Bayona-Vásquez et al., 2019) with the MspI, BamHI and ClaI enzymes, at the University of Georgia.Demultiplexed reads can be found in the sequence read archive (SRA) database from NCBI under BioProject number PRJNA1003955.

| SNP calling
We used STACKS ver.2.54 (Rochette et al., 2019) to identify SNPs.and M (distance allowed between stacks) for de novo assembly were chosen based on the tests of a subset of samples (r80 method) as suggested by Paris et al. (Paris et al., 2017).The final STACKS parameters (m, n, M) used in each set and the total number of SNPs obtained are presented in Table S1.
In the populations module of STACKS, the parameters were constant for every dataset: --r 0.5, --R 0.7, --p 2, --min-maf 0.05, --max-obs-het 0.6 and, to avoid linked loci, we used --write-singlesnp.The parameter -r is the minimum percentage of individuals in a population required to process a locus for that population, −R is the minimum percentage of individuals across populations required to process a locus, −p is the minimum number of populations a locus must be present in to process a locus; --min-maf specifies a minimum minor allele frequency required to process a nucleotide site at a locus and --max-obs-het specifies a maximum observed heterozygosity required to process a nucleotide site at a locus.

| Genomic diversity and differentiation analyses
We calculated the heterozygosity, nucleotide diversity and the inbreeding coefficient for datasets 1-4 with the snpR v.1.2.7.1 package in R (Hemstrom & Jones, 2023), to evaluate the genetic diversity in each group and populations within the three taxa.For the genetic structure, we performed AMOVA analyses with 10,000 permutations in Arlequin v.3.5.2.2 (Excoffier & Lischer, 2010)   and absolute (D XY ) genomic divergence for every locus between pairs of taxa and between populations within taxa.To correct for missing data, for these analyses, we used the variant plus invariant sites VCF files.These analyses were performed in pixy v.1.2.6.(Korunes & Samuk, 2021).
We also estimated the number of genetic clusters through Bayesian analyses in STRUCTURE v.2.3.4 (Pritchard et al., 2000).These analyses were performed for dataset 1 to explore the genetic structure among the three taxa.Also, to analyse the structure within each group we performed STRUCTURE analyses for datasets 2, 3 and 4 to analyse the structure within each group.Finally, to evaluate the differentiation between H. americana North and South we performed this analysis with dataset 7. The program was set to run with values of assumed genetic clusters (K) from 1 to 4, 9 or 10 (depending on the number of populations in each data set) with a burn-in of 10,000 steps and chain length of 100,000 using the admixture model with correlated allele frequencies, without prior population information.Ten independent runs were performed for each value of K.The most likely number of genetic clusters was determined using the method of Evanno (ΔK) (Evanno et al., 2005), on the STRUCTURE HARVESTER website (Earl & vonHoldt, 2012).
We performed principal component analyses (PCA) in PLINK v.1.90(Chang et al., 2015) to visualize genetic structure for the same datasets as for STRUCTURE analyses.Moreover, we carried out a discriminant analysis of principal components (DAPC) (Jombart et al., 2010) in adegenet v.2.1.4 in R (Jombart, 2008) for dataset 1 to test for the presence of three taxa; this method maximizes the variance among groups while minimizing the variation within groups.
We tested for K = 3 as the optimal group number and to estimate the number of principal components (PCs) to be retained, we used the cross-validation method.
Finally, we also estimated the phylogenetic relationships between and within the taxa (datasets 1-4).We used the concatenated consensus sequence per locus (that includes all variant and invariant sites) per population (--phylip-all-var in populations module), resulting in a single sequence per population.We estimated the mutation model with ModelFinder (Kalyaanamoorthy et al., 2017) and the phylogenetic tree was estimated under a maximum likelihood approximation with IQ-TREE (Nguyen et al., 2015).We also performed 10,000 bootstrap alignments using the ultrafast bootstrap method (Hoang et al., 2018).All these analyses were run in the IQ-TREE server (Trifinopoulos et al., 2016).The trees were visualized and edited with the ggtree v3.5.0.901 package (Yu et al., 2017).

| Ecological niche modelling
To determine niche differences among the taxa within the Hetaerina americana complex, we performed an ecological niche model (ENM) for each group and tested for niche divergence among them.For these analyses, we only used the records of populations sampled by Vega-Sánchez et al. (2019, 2020), because of the uncertainty of a correct species assignment in other databases due to the cryptic nature of this complex.In total, we used 21 occurrences for H. americana North, 14 for H. americana South and 55 for H. calverti.
To perform the ENMs, we used 19 variables from WorldClim 2.1 with a resolution of 30 arc seconds (Fick & Hijmans, 2017).Because the choice of the calibration area has an important influence on model training and validation (Barve et al., 2011), we hypothesized that the area of accessibility (M) for H. calverti and H. americana North includes the biogeographic provinces within the Nearctic and Neotropical regions (Escalante et al., 2021;Morrone et al., 2022) where the presence of these taxa has been recorded (Table S2).In the case of H. americana South, the calibration area encompassed only the biogeographic provinces in Mexico (Morrone et al., 2017) where the taxon occurs.Then, the 19 variables were cropped using the M of each taxon as a mask.
To reduce redundancy among climate variables, one from each pair with a high Pearson correlation value (>0.85) was discarded, preferentially retaining those with a higher contribution to the ENM according to a preliminary Jackknife analysis using all variables, implemented in MAXENT 3.4.1 (Phillips & Dudík, 2008).The final set of variables used for each group is presented as supplementary material (Table S3).To conduct the ENM, we used the maximum entropy algorithm implemented in MAXENT.The R package kuenm (Cobos et al., 2019) was employed to explore 90 possible combinations of three feature class combinations (linear = l, quadratic = q, product = p) and 15 values of regularization multipliers (from 0.1 to 6).From all models obtained, those with better evaluation based on the Akaike Information Criterion, the area under the curve (AUC) of the Receiver Operating Characteristic (ROC) and the omission rate, were selected.Once the best parametrization was defined, MAXENT was run with these characteristics to generate the final models.The rasters obtained were thresholded based on 10th percentile values to obtain binary maps for each taxon.To identify areas where taxa could be potentially co-distributed, the three potential distribution models were mapped together.
We used three approaches to test for differences in climatic niches among taxa: niche overlap, niche identity and niche similarity tests.The niche overlap test makes multiple pairwise comparisons between the models of two taxa to estimate their observed similarity (D).We measured niche overlap using a PCA-environment approach (Broennimann et al., 2012;Warren et al., 2008).The most important and uncorrelated variables in the calibration area of the three groups according to the Jackknife analysis were selected to perform the PCA (AMT, TS, MxTWaM, MTWeQ, MTCQ, PCQ and Elev; Table S3).Then, the first two principal components were used to define the environmental space based on a Kernel density function and measures of climatic niche overlap among taxa were estimated using the Schoener's D and Hellinger's I indexes, ranging from 0 (no overlap) to 1 (complete overlap) (Warren et al., 2008).
The niche identity (equivalency) test evaluates the hypothesis that two groups are statistically significantly different (divergence) or equivalent (conservatism).This test first measures the overlap index of random pseudo-replicates in the environmental space of the study area.Then, a null distribution is created based on those replicates and compared with the observed value between the two taxa in a one-tailed test (Broennimann et al., 2012;Warren et al., 2008).One hundred replicates were performed for each comparison.
Finally, the niche similarity (background similarity) test estimates whether the observed overlap between the climatic niche of two taxa is different from the overlap between the observed climatic niche in the range of one taxon and a second climatic niche generated by random background points (i.e., whether the climatic niche of one taxon predicts the climatic niche of a second taxon better than expected by chance) (Peterson et al., 1999;Warren et al., 2008).In this case, niche points between pairs of taxa were simultaneously shifted using 100 replicates.All procedures were performed using the R package ecospat (Di Cola et al., 2017).
In general, we found low levels of observed heterozygosity for the three groups (H O < 0.06; Table 1).For populations within taxa, H O ranged from 0.04 to 0.22 (Table 1).Positive and significant values of F IS were found in all groups, with H. americana North being the one that presented the highest levels of heterozygote deficit.For populations within taxa, significant and positive F IS values for most of the populations were found.The RCh population for H. americana North was the exception with a negative F IS value (Table 1).

| Genomic differentiation between taxa
For dataset 1, the AMOVA showed a very high overall genetic structure among the three taxa (F ST = 0.89, p < .001)(Table 2).The highest pairwise F ST values were between H. calverti and H. americana South (F ST = 0.83, p < .0001)and between H. calverti and H. americana North (F ST = 0.78, p < .0001),while the F ST value between H. americana North and H. americana South was 0.57 (p < .0001).We also found that most loci presented high values of relative and absolute genetic differentiation values between pairs of taxa (Figure 2).STRUCTURE analysis supported the presence of two main genetic groups based on the ΔK method (K = 2) for dataset 1, separating H. calverti from a group merging H. americana North and H. americana South (Figure 1d).Nevertheless, for K = 3, it is evident that the three genetic clusters show very low admixture (Figure 1d) and H. americana North and South are clearly separated.This was also evident in the PCA analysis, in which the two first principal components explained more than 95% of the variation (Figure 1b).When we estimated the structure for data without H. calverti, we found that the most likely value of K was 2, separating H. americana North and H. americana South, with only a few admixed individuals in populations SN and EF (Figure S1).The first axis of the PCA analysis explained 30% of the variation and separated H. americana North and South.The second PC explained 15% of the variation and separated the SP population from the rest of H. americana North populations (Figure S1).
The DAPC analysis which included all individuals, retained 20 PCs (based on the lowest RMSE value from the cross-validation).
The plot showed three isolated clusters corresponding to the three taxa (Figure S2).Finally, the maximum likelihood phylogenetic tree showed different clades that corresponded to the taxa (Figure 1c) being H. calverti the sister species of H. americana North and South which are also monophyletic.Nevertheless, the position of SN and EF populations is not well supported either in the clade of H. americana North or in H. americana South.
Population pairwise F ST values showed that most population pairs have high and significant values of differentiation in H. calverti, with the highest differentiation observed between SG and Ct populations (F ST = 0.75, p < .01)(Figure S3).For H. americana North, the highest differentiation was between SP and EF populations (F ST = 0.78, p < .05).For H. americana South, the highest differentiation was between Ve and CC populations (F ST = 0.75, p < .05).
In contrast to the other taxa, in H. americana South most of the populations from the center of Mexico were not differentiated (Figure S3).
The relative and absolute genetic differentiation for each locus between populations within taxa also showed that some populations are highly differentiated but not all (Figure S4).For H. calverti, SG and Ct populations showed the greatest differentiation for F ST (0.59; s. d. = 0.37) and SG and Za for D XY (0.0016; s. d. = 0.0016).
On the opposite side, SM and Ve showed the lowest F ST value (0.14; s. d. = 0.17 The STRUCTURE analysis suggested that the most likely number of genetic clusters (K) for H. calverti was two.The first group included the north-east populations Ct and Za and the second group included the rest of the populations (Figure 3a).Nevertheless, the PCA analysis suggested that SG is a different group separated by the PC2 as is evident in the plot of K = 3 (Figure 3a).This is also supported by the maximum likelihood tree, which showed four different clades (Figure 3a).It is interesting to observe that population EF (Sinaloa; North-West) is more related to populations of the south (Chiapas; NR and CC) than to populations of the northeast (Coahuila; Ct, Za).
For H. americana North, STRUCTURE analysis suggested the presence of three clusters (K = 3) (Figure 3b).The first group included populations BHC, DC, LC and distributed in the centre-east of the USA.The second group included populations EF, SN and RCh, which are distributed in the north of Mexico.Finally, the last group included only the population SP, from the Baja California Peninsula (Figure 3b).The PCA showed a similar pattern which is also congruent with the maximum likelihood tree (Figure 3b).3c).The PCA showed the same pattern, with two highly differentiated groups, northeast and south, and an intermediate group of populations from the centre of Mexico (Figure 3c).The phylogenetic relationships between these populations suggested two main clades, the northeast and the centre-south groups (Figure 3c).

| Niche modelling
The most important variables according to the MAXENT results are presented in Table S4.The validation test values for the three groups America, there are also suitable areas in Guatemala, El Salvador, and Nicaragua (Figure 4c).
The area of potential overlap among the three Hetaerina taxa comprised, mainly, the coast region of Baja California Sur, Sonora and Sinaloa, as well as the south-centre of Mexico in the states of Jalisco, Michoacán and Guerrero (Figure 4d).The overlap between   The niche identity test indicated a significant niche divergence between H. calverti and H. americana North (p = .006).In the case of H. calverti and H. americana South, the climatic niche divergence was not greater than chance expectations (p = .39).Finally, for H.

H. calverti and H. americana
americana North and H. americana South, the niche identity test also indicated that the climatic niches of the two taxa are divergent (p = .003).
The niche similarity tests suggested that the climatic niches of all pairs of taxa are not more similar than random niches (p > .05for all comparisons).In general, the results indicate low climatic niche overlap among all Hetaerina groups, suggesting that each group occurs at specific conditions of temperature, precipitation and elevation.

| DISCUSS ION
4.1 | How many species are there in the Hetaerina americana complex?
The first remarkable result is the high degree of genetic differentiation between the three taxa for most loci.This indicates that gene flow is null among the three taxa even when there is co-occurrence at some localities (e.g., EF, Za, Ve, CC, NR; Figure 1).Furthermore, we found that H. calverti is the most divergent taxon and that, in agreement with previous findings (Vega-Sánchez et al., 2022), it shows complete genetic isolation from the other two.In contrast, the two groups of H. americana (North and South) showed evidence of a lower degree of divergence; however, the patterns of genetic differentiation allowed us to recognize them as two discrete taxa, except for two populations (EF and SN) that showed some mixed ancestry, which can be associated to incomplete or weaker isolation barriers between these two taxa, as are ecological barriers (Sánchez-Guillén et al., 2012).These patterns of genomic divergence indicate that the three taxa represent species in different stages of the speciation continuum (De Queiroz, 2007;Stankowski & Ravinet, 2021).
Morphological variation in caudal appendages is known to play an important role as a premating mechanical isolation barrier in odonates such as Enallagma, Ischnura and Hetaerina (McPeek et al., 2011;Sánchez-Guillén et al., 2012).The differentiation in the shape of the male caudal appendages is more conspicuous in H. calverti and it is considered the main morphological character that contributes to reproductive isolation in localities where it is sympatric with H.
americana South (Vega-Sánchez et al., 2022).In contrast, the isolation barriers between H. americana North and South seem to be different, since the divergence in the shape of caudal appendages is subtler (Vega-Sánchez et al., 2019).
Another premating barrier that has been reported in odonates with no courtship, is habitat isolation (Sánchez-Guillén et al., 2012).
In this sense, we found evidence of niche divergence between H.
americana North and South, and few areas of their distribution seem to be shared by the two taxa, suggesting habitat isolation.In other odonates, it has been described that temperature is the main environmental factor that shapes the distribution of these ectotherm animals (Dudaniec et al., 2018).In our case, for example, the minimum temperature of the coldest month (MnTCM) is a variable that ranges from −14 to 8°C, with an average of −1.stages are active and mating is constant during the whole year (Paulson, 2011a(Paulson, , 2011b)).These differences in the ability to deal with low temperatures, especially during immature stages, could be associated with local adaptation, but this needs to be tested.Also, the different reproductive phenology between H. americana North and South could be related to reproductive isolation (e.g., temporal isolation).
The adaptation to different habitats could limit the number of effective migrants and therefore, the gene flow between species.
Under this scenario of habitat isolation, we also could expect, due to the low divergence in the morphology of caudal appendages, incomplete or weak mechanical barriers and thus the existence of hybrid zones between H. americana North and H. americana South.
This could be especially probable in areas in the north of Mexico as in the EF and SN populations, where there is a transition to lower temperatures in the winter season.Thus, these areas could experience the occasional immigration of H. americana South individuals (in warmer months), resulting in some extent of hybridization.This is also supported by niche modelling and STRUCTURE analysis because it is in this area where a geographic overlap between the two taxa is predicted and where some individuals showed mixed ances-  SG and NR are also contrasting in environmental conditions related to precipitation and elevation and this could also be influencing the patterns of genetic variation.
The SP population of H. americana North also showed high genetic differentiation from all the other populations of this taxon, what could be associated with the drier conditions where this populations is located.Moreover, in this case, the Sonoran Desert could be also participating as a geographical barrier.

| Speciation by neutral or niche divergence processes?
In general, a multifactorial scenario (geographical barriers and environmental differences) could be related to the early stages of divergence in these organisms (e.g., between populations within species), but the presence of geographic barriers seems to be more important for a limited gene flow (e.g., the Sonoran Desert and Motagua-Polochic-Jocotán fault system).Based on this, we could suggest that the divergence between the three species probably started in an allopatric setting by vicariance processes.This agrees with findings in other damselflies in which "nonadaptive processes" drove the speciation (Wellenreuther & Sánchez-Guillén, 2016).Nevertheless, natural or sexual selection can act after secondary contact and maintain the genetic divergence if the cross between divergent individuals is maladaptive.For example, there is evidence in groups such as copepods and amphibians that suggests that divergent selection between thermal environments is often strong enough to maintain a bimodal genotype distribution, highlighting the relevance of thermal adaptation in ecological speciation processes (Keller & Seehausen, 2012).
Finally, understanding why in some species the divergence of morphological traits as the caudal appendages and body coloration has evolved and not in others is an outstanding question.An evolutionary hypothesis for the maintenance of reduced phenotypic divergence in Haeterina females is the catch-22 hypothesis (Drury et al., 2019), which attempts to explain why reproductive interference is common between species of this genus (Drury et al., 2015).
According to the catch-22 hypothesis, the females of different sympatric species are very similar in coloration patterns and thus the males are not able to distinguish among them, resulting in a high frequency of heterospecific sexual interactions.Since the cost of lost mating opportunities is higher than the costs of heterospecific interactions, there is no divergent selection in male mating recognition and consequently no divergent sexual selection in female phenotypic characters.The lack of phenotypic divergence in male coloration can also be explained by high levels of competitive interference; their competitor recognition system is based on wing coloration which is maintained because they share a common resource (females), nevertheless, this needs to be tested (Drury et al., 2015).
This hypothesis opens the door to the existence of more cryptic species in the genus, as is also shown by genetic data in other Hetaerina species such as H. titia and H. cruentata (Vega-Sánchez, 2013).
As conclusion, the H. americana species complex presents a very tangled evolutionary history.We have shown that this complex comprises three highly genetically differentiated taxa that still include at least one entity requiring taxonomic recognition as a distinct species.These biological entities are reproductively isolated by diverse mechanisms including morphological and ecological barriers.
Moreover, we found evidence that niche divergence and possible local adaptation related to temperature variation are associated to the genetic divergence between morphologically cryptic species, highlighting the importance of environmental variation in the speciation process in this order of insects, which has been suggested to have evolved by "nonadaptive" processes.Moreover, we found pronounced genetic differentiation between populations within species that could be also related to environmental variation and the presence of some geographical barriers.If such differentiation could suggest the presence of more cryptic species arises as a hypothesis.Nevertheless, better sampling and more detailed morphological analyses are necessary to evaluate this new hypothesis.

First
, we ran process_radtags to demultiplex, quality filter and trim the sequence data.Then, we used clone_filter on demultiplexed reads to identify PCR clones, using Illumina iTru5 barcodes.After filtering the data, five samples were removed due to low quality; the final number of individuals analysed was 90.After cleaning, we performed seven different SNP callings with different sample sets to enhance the number of SNPs obtained and to answer different questions with each dataset.To analyse genetic structure among the three taxa, dataset 1 included all individuals (n = 90) belonging to H. calverti, H. americana North and H. americana South.Secondly, to analyse genetic variation within each group, we carried out three SNP callings: dataset 2 included only individuals of H. calverti (n = 31), dataset 3 included individuals of H. americana North (n = 24) and dataset 4 included individuals of H. americana South (n = 35).Finally, to compare taxa by pairs, we performed three more SNP callings: dataset 5 included individuals of H. calverti and H. americana North (n = 55), dataset 6 included individuals of H. calverti and H. americana South (n = 66) and dataset 7 included individuals of H. americana North and South (n = 59).The parameters m (minimum depth coverage), n (distance allowed between catalogue loci)

F
Genetic structure patterns in the American Rubyspot species complex.(a) Sampling localities, with colours representing the three taxa studied.Circles with two colours indicate sympatric localities.Map background colours represent the variation in the Minimum Temperature of the Coldest Month (MnTCM) along the sampling area.See Table 1 for locality name abbreviations.(b) Principal component analysis (PCA) includes 90 individuals and 4722 SNPs, each circle represents an individual and includes the population label; colours represent the three taxa.(c) Maximum likelihood tree based on consensus sequences per population, bootstrap values are indicated by colours in the squares.(d) STRUCTURE analysis for K = 2 to K = 4.Each bar represents an individual and its assignment to different clusters, population labels are also shown.[Colour figure can be viewed at wileyonlinelibrary.com] populations within taxa.Moreover, we estimated the relative (F ST ) The average and standard deviation for H. calverti and H. americana North were F ST = 0.52 (s.d. = 0.37) and D XY = 0.0024 (s.d. = 0.0017); for H. calverti and H. americana South F ST = 0.58 (s.d. = 0.37) and D XY = 0.0027 (s.d. = 0.0017), and for H. americana North and H. americana South F ST = 0.36 (s.d. = 0.28) and D XY = 0.0017 (s.d. = 0.0012).
) and Ct and Za the lowest D XY value (0.0007; s. d. = 0.0009).For H. americana North the values ranged from F ST = 0.75 (s.d. = 0.34) between MC and SP to 0.19 (s.d. = 0.21) for EF and SN.The highest D XY value was 0.0019 (s.d. = 0.0019) for SP and EF populations and the lowest was 0.0005 (s.d. = 0.001) for BHC and DC populations.Finally, for H. americana South the highest value for F ST (0.65; s. d. = 0.35) was between Cm and Ve populations and the lowest was 0.07 (0.16) for Bo and Co populations.For D XY , populations CC and Ve showed the highest value (0.0016; s. d. = 0.0018) and Vi and Za populations showed the lowest value (0.0006; s. d. = 0.0009).

TA B L E 1
Genetic diversity values in three taxa of Hetaerina (H.calverti, H. americana North and H. americana South) and for populations within taxa.
indicated a good performance (H.calverti AUC = 0.95, H. americana North AUC = 0.94 and H. americana South AUC = 0.86).The potential distribution of H. calverti has a wide extension, from the southeast of the USA through Mexico to Guatemala and Nicaragua.In Mexico, the model predicts suitable conditions for the species mainly on the Pacific and Atlantic coasts in low-elevation areas (Figure 4a).The potential distribution of H. americana North is also wide, encompassing from the north of the United States to the northwest of Mexico, in the states of Sonora, Sinaloa, Chihuahua, Durango and Nayarit.Moreover, the distribution extends to the south in some regions of the Trans-Mexican Volcanic Belt (TMVB) (Figure 4b).Meanwhile, H. americana South has the most restricted potential distribution of the three taxa, comprising mainly the central part of Mexico, the west coast from Sonora to Chiapas states, and patchy areas in the Atlantic coast in Tamaulipas, Veracruz and Campeche.In Central North included the north-eastern region of Mexico and the Baja California Peninsula as well as the TA B L E 2 Analyses of molecular variance (AMOVA) evaluating overall differentiation among the three Hetaerina taxa (dataset 1) and populations within H. calverti (dataset 2), within H. americana North (dataset 3) and H. americana South (dataset 4).

Fixation indices F
IS = −0.03,F ST = 0.42**, F IT = 0.40** Note: Datasets refer to different SNP callings performed for each group of individuals.See text for details.**p < .0001.south-eastern region of the USA in Texas, Mississippi, Alabama and Georgia.On the other hand, the potential distribution of H. calverti and H. americana South showed the greatest area of overlap, mostly in the southwestern region of Mexico and Nicaragua.Finally, the overlap between H. americana North and H. americana South included Chihuahua, Sonora and Sinaloa states in the northwest of Mexico, as well as some regions of the TMVB in the states of Jalisco, Guanajuato and Nayarit.The niche overlap test values between H. calverti and H. americana North were low (I = 0.016, D = 0.003).Also, low values of overlap were found between H. calverti and H. americana South (I = 0.187, D = 0.054).Finally, the overlap values between H. americana North and South were the lowest (I = 0.005 and D = 0.0008).

F I G U R E 2
Genomic differentiation per locus between the three Hetaerina taxa.Frequency distributions of D XY and F ST values per locus are shown, as well as average values for each comparison.[Colour figure can be viewed at wileyonlinelibrary.com] individuals die during the winter season and only immature stages remain in latency until the next spring season when the temperature rises, while in southern populations adults and immature

F
Population structure within each of the three taxa.(a) H. calverti, (b) H. americana North and (c) H. americana South.For each taxon phylogenetic relationships among populations, sampling distribution map, plot of the two principal components and STRUCTURE analysis for K = 2 to K = 4 are shown.For each case, different populations are identified by colours.In the maps, blue shading indicates elevation.[Colour figure can be viewed at wileyonlinelibrary.com]Beyond the three species, we also found evidence of high genetic differentiation between some populations within species.Especially, populations Ct, Za and SG for H. calverti; SP population in H. americana North and Ve, Za and Vi in H. americana South showed the greatest D XY and F ST differentiation values as well as genetic distinctness patterns based on clustering analyses.This differentiation could be driven by neutral processes such as genetic drift as result of a limited dispersal and geographic isolation, or by adaptive processes.We found evidence of both processes in different populations of calverti.The SG and NR are populations that are close in distance (117 km in Euclidian distance), and despite the usually assumed high dispersal potential of these organisms, we found high genetic differentiation between them (F ST = 0.50).In contrast, the NR and EF populations are 2095 km apart and showed an F ST value of 0.2, suggesting that factors different to geographic distance are shaping genetic differentiation.Interestingly, the SG and NR populations are separated by an important physiographic feature, the Motagua-Polochic-Jocotán fault system (Gutiérrez-García & Vázquez-Domínguez, 2013), which is a barrier to gene flow in several animal species (Mendoza et al., 2019; Rodríguez-Gómez & Ornelas, 2014).

F
Potential distribution models for the three taxa.(a) Binary map for the potential distribution of H. calverti.(b) Potential distribution for H. americana North.(c) Potential distribution for H. americana South.(d) Overlap among the potential distribution of the three taxa.The red colour represents areas where the three taxa overlap their potential distribution.Orange-coloured areas represent the overlap between H. calverti and H. americana North.Yellow zones represent the overlap between H. calverti and H. americana South.Purple areas represent the overlap between H. americana North and South.Light green, light blue and light salmon areas represent the potential distribution for H. calverti, H. americana North and South, respectively.The sampling localities are shown.[Colour figure can be viewed at wileyonlinelibrary.com] and estimated the pairwise F ST values among taxa (i.e., H. calverti, H. americana North and H. americana South) and 52°C, in H. americanaNorth localities, and from 3 to 19°C, with an average of 10.3°C, in H. americana South localities, indicating that populations of H. amer- icana North are adapted to colder temperatures than populations of H. americana South.Actually, in northern populations, all adult