Molecular evidence of hybridization in sympatric populations of the Enantia jethys complex (Lepidoptera: Pieridae)

Hybridization events are frequently demonstrated in natural butterfly populations. One interesting butterfly complex species is the Enantia jethys complex that has been studied for over a century; many debates exist regarding the species composition of this complex. Currently, three species that live sympatrically in the Gulf slope of Mexico (Enantia jethys, E. mazai, and E. albania) are recognized in this complex (based on morphological and molecular studies). Where these species live in sympatry, some cases of interspecific mating have been observed, suggesting hybridization events. Considering this, we employed a multilocus approach (analyses of mitochondrial and nuclear sequences: COI, RpS5, and Wg; and nuclear dominant markers: inter-simple sequence repeat (ISSRs) to study hybridization in sympatric populations from Veracruz, Mexico. Genetic diversity parameters were determined for all molecular markers, and species identification was assessed by different methods such as analyses of molecular variance (AMOVA), clustering, principal coordinate analysis (PCoA), gene flow, and PhiPT parameters. ISSR molecular markers were used for a more profound study of hybridization process. Although species of the Enantia jethys complex have a low dispersal capacity, we observed high genetic diversity, probably reflecting a high density of individuals locally. ISSR markers provided evidence of a contemporary hybridization process, detecting a high number of hybrids (from 17% to 53%) with significant differences in genetic diversity. Furthermore, a directional pattern of hybridization was observed from E. albania to other species. Phylogenetic study through DNA sequencing confirmed the existence of three clades corresponding to the three species previously recognized by morphological and molecular studies. This study underlines the importance of assessing hybridization in evolutionary studies, by tracing the lineage separation process that leads to the origin of new species. Our research demonstrates that hybridization processes have a high occurrence in natural populations.

Introduction indistinguishable from E. albania (cryptic species), but forming a sister clade of E. mazai, suggesting the possibility of recent speciation or even hybridization. The existence of hybridization is supported by field expeditions that have permitted the observation of two separate events of interspecific mating: one between male E. jethys and female E. mazai, and another between male E. albania and female E. jethys (personal observation, Castañeda-Sortibrán and Jasso-Martínez in this paper), suggesting that hybridization may be a relatively frequent process in the Enantia jethys complex.
In this study, we employed the molecular markers commonly used in phylogenetic studies (mitochondrial and nuclear sequences), as well as inter-simple sequence repeat (ISSR) markers to examine contemporary genetic exchange among the butterflies of the Enantia jethys complex. Identification of fast molecular markers, such as ISSR, are the principal tools to study hybridization processes, and are recommended instead of using microsatellites [18]. Hybridization and introgression events have been successfully studied through ISSRs in many taxonomic groups [19][20][21][22], including butterflies [23]. The ISSR-PCR method is a molecular technique used to screen a large part of the genome without prior knowledge of the sequences. This method provides highly reproducible results and generates abundant polymorphisms in many systems [20]. The ISSRs are dominant molecular markers, where the absence of a band is interpreted as the loss of a locus/allele through either the deletion of the SSR site or a chromosomal rearrangement [24]. The use of ISSR to study species of Lepidoptera is relatively recent, yet abundant (not an exhaustive list: [25][26][27][28][29][30][31]).
In summary, considering the hypothesis of hybridization events in the Enantia jethys complex, we used a set of molecular markers (one mitochondrial and two nuclear sequences, and ISSR) to: i) confirm the hypothesis of the hybridization process in the butterfly complex, and, in the event our hypothesis is correct, to evaluate the extent of this process; ii) identify the directionality of the introgression; iii) determine if genetic diversity presents variation between admixed and non-admixed individuals (i.e., hybrids and non-hybrids) for each morphospecies; and iv) evaluate, using different molecular markers, the monophyly state of the butterfly complex. Finally, we discuss different hypotheses to explain the evolutionary history of this complex.

Ethics statement
Species involved in this study are not endangered or protected in the study area, and no specific permission is required for scientific research. No specific permission was required for our fieldwork in any location visited.  Table 1). All butterflies were collected in mountain cloud forests during 2015 (April-September), where the mean monthly temperature over this period varied from 23˚C to 28.5˚C according to locality. Butterflies were collected during consecutive days from 9:00 am to 4:00 pm using entomological nets. All individuals from the Enantia jethys complex were kept separately in glassine bags, labeled and geo-referenced (mobile navigation application: Locus Map). Other species of butterflies captured in the nets were released. All samples were brought back to the "Laboratorio de Genética y Evolución" (Genetics and Evolution Lab) of the "Universidad Nacional Autónoma de México" (UNAM) for their morphological identification; subsequently, abdomens from all individuals were cut, placed in absolute ethanol, and finally conserved at 4˚C until DNA extraction. Voucher specimens were deposited at the "Museo de Zoología de la Facultad de Ciencias" (Science Faculty Museum of Zoology) UNAM (MZFC-UNAM) in Mexico City.

Butterfly samples
Species identification (Fig 2) was based on color wing phenotypes using the descriptions from Llorente-Bousquets and colleagues [15]. Considering the difficult and sensitive identification of females (from E. jethys and E. mazai), the following steps were taken to ensure correct identification. First, we identified all males, which have wing pigmentation patterns that differ among the three species. Subsequently, females of E. albania were easily identified (presence of a diagnostic mark on the posterior wings), and separated from the remaining females. Finally, females of E. jethys and E. mazai were separated by careful observation with a stereoscopic microscope (Stemi DV4, Zeiss) and photographing each specimen. Once morphospecies (males and females) were identified, all individuals were correctly labelled and deposited (physically and in photograph) in the MZFC-UNAM collection.
The PCR mix consisted of 14.15 μl (RpS5) and 14.4 μl (Wg) ddH 2 O, 5 μl 5X Green Buffer, 0.5 μl dNTPs mix, 0.15 μl (RpS5) and 0.1 μl (Wg) Taq Polymerase 5 U/μl (Promega), 2.2 μl (RpS5) and 2 μl (Wg) of MgCl 2 , 2 μl of template DNA, and 0.5 μl of each primer (10 μM). Amplifications were carried out under specific conditions: initial denaturation at 95˚C for 6 min (RpS5) and 3 min (Wg), 40 cycles of denaturation at 95˚C for 30 s, primer annealing temperature at 51˚C for 30 s (RpS5) and 60 s (Wg), an extension at 72˚C for 90 s, and a final extension at 72˚C for 10 min (RpS5) and 6 min (Wg). The protocol for visualizing the PCR products was identical to mitochondrial DNA. The PCR products were sent for purification and sequencing to MACROGEN INC. (Seoul, Korea). GenBank accession numbers are presented in S1 Table. ISSR nuclear gene. We tested a total of 24 different ISSR markers. The selection of markers for this study was based on the identification of a particular band profile to facilitate the identification of each morphospecies, in addition to presenting good resolution and a high number of bands. Two ISSR markers were retained for our study: (AG) 8 Y and (GA) 8 C ( Table 2). Percentage of guanine and cytosine content (%GC), melting temperature (Tm), annealing temperature (Ta), total number of bands per primer over all localities (N bands), size range of the DNA fragments for each primer (size). The designation Y (C or T) was used for degenerated sites.
PCR amplifications were performed in 15 μl reaction volume containing~20 ng of template DNA, 1.5 μl 5X Green Buffer (Promega), 200 μM dNTP (dNTP mix; Promega), 3 mM MgCl 2 (Promega), 1 μM of primer (Integrated DNA Technologies), and 1.25 U GoTaq Flexi DNA Polymerase (Promega); finally, the volume was adjusted with ultrapure water. Amplifications were conducted in a T100 Thermal Cycler (Bio-Rad™): initial denaturation step at 94˚C for 4 min, 39 cycles of denaturation at 94˚C for 45 s, annealing temperature (Ta) 54˚C or 56˚C depending on the ISSR primer (Table 2), extension temperature at 72˚C for 2 min, and a final extension at 72˚C for 10 min. DNA banding patterns were visualized by electrophoresis, performed with 3 μl of amplified products on a 2% agarose gel using 1X TAE buffer and poststaining with GelRed™ (Biotium), at 110 V for 2 h. A 100 bp DNA Ladder (Promega) was used to estimate amplification product lengths. Fragment (band) patterns were visualized and digitized using an imaging system (PhotoDoc-it, UVP1).

DNA sequence analyses
Chromatograms of the forward and reverse sequences were edited and assembled using GEN-EIOUS R9 [41], and sequences were then aligned with GENEIOUS R9 and MAFFT online [42].
Population genetic parameters were determined for each of the three morphospecies (E. albania, E. jethys, and E. mazai) for mitochondrial and nuclear genes. Program DNASP 5.10.1 [43] was used to determine the number of haplotypes (nHap), haplotype diversity (h, also called genetic diversity), and nucleotide diversity (π). To determine the genetic structure among species, an analysis of molecular variance (AMOVA; 10,100 permutations) was performed using the ARLEQUIN program 3.5 [44] under the Tamura and Nei (TrN) model, the best-fit model for nucleotide substitution, as selected by the corrected Akaike Information Criterion (AIC C ) in JMODELTEST 2.1.9 [45][46]. Haplotype network analyses were performed for nuclear and mitochondrial sequences with Network software [47] though the probabilities method with the algorithm Median Joining.
For each matrix of DNA sequences, Bayesian inference was conducted in MRBAYES 3.2 [48] over the CIPRES web platform [49]. We used the reversible jump technique for selecting the best substitution model [50] and considered a gamma distribution with four categories and a percentage of invariant sites; the MCMCMC chains were run for 15 million generations, and the first 25% of the trees (burn-in) were discarded. A consensus tree was calculated after burn-in, and the posterior probabilities summarized in the MRBAYES consensus tree were indicated on the nodes.
We performed a molecular clock analysis among species using a relaxed clock approximation with BEAST 1.8.2 [51]. The analysis was performed over the web platform CIPRES [49] using mean and standard deviation of the COI marker per Papadopoulou and colleagues [52]. We ran the analysis for 100 million generations. The ESS values were revised on TRACER software [51]. Trees sampled after burnin were used to construct the Maximum Clade Credibility Tree in the TREEANOTATOR software included in the BEAST package.

ISSR analysis
The ISSRs were treated as dominant markers; amplified fragments were scored as 1 (presence of a band represents a dominant allele) or 0 (absence represents a recessive allele). With this principle, a binary data matrix was generated for all individuals and the two ISSR markers used. Only bands that could be scored consistently among localities, and individuals that presented genetic information for all primers, were used for analysis. This can explain the low discrepancy between the numbers of collected individuals, and the number of individuals used in our analysis.
All individuals were classified in one of the three morphospecies based on their genetic ISSR profile; after which, the number of individuals of each morphospecies was compared by wing color phenotype and ISSR profile identification. For each morphospecies, we determined the parameters using GENALEX 6.5 [53][54] and POPGENE 1.31 software [55], considering all localities together and individually: total number of bands (TB), number of rare bands (RB), number of private bands (PB), percentage of polymorphic loci (%P), and Nei's gene diversity (h). Furthermore, the percentage of polymorphic loci (%P) and Nei's gene diversity (h) were determined for "no-admixture" individuals and "admixture" individuals, considering all localities together and individually, using POPGENE 1.31 software. A one-way analysis of variance (ANOVA) was used to evaluate the differences in genetic diversity (h) among the three morphospecies, considering all localities together and individually, among no-admixture and admixture individuals for each morphospecies, and for each morphospecies among localities, using STATISTICA7.0. To visualize the relationship among individuals from the three morphospecies, a principal coordinate analysis (PC o A) was performed in GENALEX 6.5.
The level of genetic differentiation among morphospecies was evaluated through the level of gene flow (Nm) and PhiPT parameters (Φ PT ), an analogue to standardized F ST for binary data, determined for pairs of morphospecies with 9999 permutations via AMOVA analysis. Analyses were performed using POPGENE 1.31 and GENALEX 6.5 programs.
To evaluate the level of genetic homogeneity for each morphospecies, we performed a Bayesian analysis implemented in STRUCTURE v2.3.3 [56][57][58]. This method was designed to identify K (unknown) groups (or clusters), characterized by a set of allele frequencies for each locus, followed by assigning each individual, with a probability (q i ), to one group or more, if group genotypes indicate that they are admixed. In this study, we did not look for the optimal number of populations (K), but we worked with a K = 3, which corresponds to the a priori three morphospecies, identified by wing phenotypes. Equally, considering the a priori morphospecies, we ran STRUCTURE software with the no admixture as the ancestry model, and the allele frequencies as the correlated models. The program was run 10 times for K = 3 to verify the homogeneity of results; for each run, the Markov chain Monte Carlo (MCMC) algorithm was run with a burn-in period of 100,000 steps, followed by 100,000 steps. Each genotype individual was assessed and allocated to one of the three clusters based on values of membership probability (q i ). Furthermore, each individual was assigned to only one cluster if q i > 0.90 (genetically "pure" individuals considered as no-admixture), otherwise they were associated with two or more clusters (admixture individuals) [19].
Finally, we performed a mean distance analysis (minimum evolution) using a heuristic search for an optimal tree, carried out with tree bisection and reconnection branch swapping. Distance analysis was performed using PAUP version 4.0b10 [59] and the tree was displayed using TREEVIEW 1.5 [60]. Negative branch lengths were allowed, but set to zero for tree-score calculation. The steepest descent options were not in effect. Starting tree(s) were obtained via neighbor joining, and bootstrap values were calculated under the same criteria.

Results
A total of 342 butterflies in the Enantia jethys complex were collected and classified according to their morphological characteristics, resulting in 76 E. albania (44 females and 32 males), 123 E. jethys (53 females and 70 males), and 140 E. mazai (45 females and 95 males) ( Table 1). A total of 100 randomly selected individuals were sequenced, out of which 83 provided highquality DNA sequences for the COI, 75 for RpS5, and 77 for Wg genes. Furthermore, 339 samples amplified correctly with ISSR markers and were used for analysis.

DNA sequences
A variable number of haplotypes (nHap) was observed following molecular markers, demonstrating the lowest number for COI (from 2 to 4) and highest for RpS5 (14 and 26) ( Table 3). The COI marker showed a very low variability of nHap among the three species, in contrast to the two nuclear genes. The RpS5 gene and Wg gene showed the highest nHap for E. jethys and E. albania, respectively. Haplotype networks (S1 Fig) of COI and RpS5 showed a clear separation among the three morphospecies, whereas for the Wg haplotype network, E. jethys and E. mazai share its haplotypes. Generally, haplotype diversity is high and relatively homogeneous among molecular markers and species, ranging from 0.5 to 0.9 (Table 3). Only E. jethys presented a very low value of haplotype diversity with respect to the COI gene. Nucleotide diversity did not show a systematic pattern that followed molecular marker or species. All values ranged from 0.002 to 0.005, with a tendency for the lowest values with the COI gene.
Each molecular marker, based on sequences (Wg, RpS5, and COI), provided different trees (Fig 3), leading to different phylogenetic hypotheses. In the Wg tree (Fig 3A), two clades were observed. In the first clade, all individuals were identified as E. mazai (blue color) and E. jethys (green color) morphospecies; they were clustered with a very high posterior probability (PP) from Bayesian inference (BI). The second clade was composed only of individuals that were identified as E. albania morphospecies (red color). Both molecular markers (RpS5 and COI) define three different clades corresponding to the three morphospecies, as previously reported [15][16][17]; however, the relationship among clades for each tree was different. The RpS5 tree ( Fig  3B) suggests that E. albania is the sister clade of E. jethys, whereas E. mazai constitutes a clade located at the base of the tree. Furthermore, the COI tree (Fig 3C) shows that E. jethys and E. mazai are sister groups, whereas E. albania is located at the base of the tree.
The concatenated analysis of the three molecular markers (S2 Fig) was consistent with our COI results (Fig 3C). There was good separation among the three species, where E. albania is a sister clade of the two other species. The molecular clock analysis (S3 Fig) identified the recent origin of the Enantia jethys complex species, starting around 1 Mya.
The AMOVA results varied greatly between mitochondrial and nuclear genes (S2 Table). The mitochondrial marker showed a high and significant level of differentiation among the three morphospecies, whereas both nuclear markers showed a lower but significant level of differentiation.

ISSR
Of the 342 Enantia butterflies collected, 339 samples amplified correctly and were used for the ISSR analysis. All ISSR data are available in S3 Table. We obtained a high number of individuals for each morphospecies (76 to 140; Table 1) and a good number at each locality (Table 1), except for Puebla, where we obtained only 15 individuals belonging to two morphospecies (used to characterize species; Table 4); however, we identified some band patterns (generally two bands nearly always together; Fig 4) that allowed us to classify individuals as one of the morphospecies. The observation of ISSR band patterns enabled us to identify many individuals with mixed patterns (admixture individuals) (Fig 4). The use of the ISSR band pattern allowed the correct identification of the majority of individuals (93.3%). Out of a total of 76 E. albania collected, 90.8% (n = 69) were confirmed by molecular ISSR technique. From the 123 E. jethys individuals, 97.6% (n = 120) were confirmed, and from the 140 E. mazai individuals, 88,6% (n = 123) were confirmed. Generally, polymorphism is very high, demonstrating values above 90% for each morphospecies (Table 4). If we consider all localities together, the morphospecies E. jethys presents a significantly higher genetic diversity than the other morphospecies, whereas E. mazai demonstrates a significantly lower genetic diversity (ANOVA test: F 195,2 = 13.8, P < 0.001) ( Table 4). For each locality, values of polymorphism and genetic diversity of each morphospecies remained similar to those observed for all localities together, even with a decrease in the sampling size. Genetic diversity was not significantly different among localities for the three morphospecies (E. albania: F 261,2 = 0.19, P = 0.82; E. jethys: F 261,2 = 0.68, P = 0.51; E. mazai: F 260,3 = 1.49, P = 0.22).
Graphical visualization of the relationship between the three morphospecies, through PC o A analysis, revealed good separation among morphospecies; however, there was an overlapping zone regarding individuals (Fig 5). After identification through the STRUCTURE analysis, hybrid Both genetic differentiation parameters (Nm and Φ PT ) showed the same tendency: a lower level of differentiation between E. albania and the two other morphospecies (E. jethys and E. mazai), than between E. jethys and E. mazai (Table 5).
Genetic structure results obtained with the software STRUCTURE identified individuals that represent a discrepancy between morphospecies and genetic profile (Fig 6). Based on these results, all individuals were classified as "non-hybrids" if 100% of its membership coefficient (q i ) corresponded to the same morphospecies. They were classified as "admixture" when individuals presented a genetic profile belonging to other morphospecies, or when individuals had two or more genetic components. Enantia albania had the lowest level of hybrid individuals (17%), followed by E. jethys with 24%, whereas over half of E. mazai individuals (53%) were classified as hybrids. Each morphospecies presented individuals with a genetic profile belonging completely to another morphospecies, and individuals that presented a genetic profile with two, or sometimes three, different genetic components (Fig 6). Notably, a unique profile that had not been observed was the combination between the E. jethys and E. mazai profiles (blue and green together; Fig 6). Polymorphism and genetic diversity of hybrids were lower than for non-hybrids of E. albania and E. jethys; however, significant differences were only demonstrated for E. jethys. Enantia mazai hybrids showed a significantly higher genetic diversity than non-hybrids (Table 4). At the locality, E. albania did not present significant differences in genetic diversity between hybrids and non-hybrids. Enantia jethys showed a systematic, significant decrease in genetic diversity of hybrids, while E. mazai presented a systematic increase in genetic diversity of hybrid individuals (Table 4). If we consider non-hybrid and hybrid individuals from the three morphospecies at each locality (Fig 7), the three localities from Veracruz present a very similar proportion of hybrids: Colonia Á lvaro Obregón 27%; Camino a la Cascada Texolo 33%; and Finca Mariposa: 35%. The Puebla locality presented a very high proportion of hybrids (87%).  From the total number of hybrid individuals (n = 117), two kinds of individuals could be distinguished: 1) individuals with only one genetic component (but different from its morphospecies), which represent 37% (n = 43); and 2) individuals with two, or three genetic components, representing 63% (n = 74). Different kinds of genetic profile combinations were found (Table 6). One individual presented a mix of the genetic profile from the three morphospecies, but interestingly no individuals combining E. jethys and E. mazai were found. Notably, almost all hybrid individuals presented an E. albania genetic component (96%, n = 112), whereas the genetic components of E. jethys and E. mazai were found in low and similar proportions (38% and 30%, respectively).
The results of the distance analysis cluster (Fig 8) showed a clear separation of the three morphospecies; however, all bootstrap values were very low. The general topology obtained by cluster analysis corresponded to the classification of individuals obtained by STRUCTURE analysis (Fig 6).

Genetic diversity in the Enantia jethys complex
The levels of polymorphism found through ISSRs in the species of the Enantia jethys complex were equal or higher than those reported by other authors for butterflies [28,[30][31]61]. Heterogeneous environments and gene exchange among populations (high connectivity) may lead to high levels of genetic variability (i.e., polymorphism) [62][63], for example P = 100% [26]. However, species with lower dispersal rates may attain lower values of polymorphism such as P 85% [29,31]. Species of the Enantia jethys complex have limited flight capacity [15], which implies low dispersal rates among populations, and therefore probably cannot explain our high level of genetic diversity. Even if no formal studies on the abundance and population density of the Enantia jethys complex exist, we collected around 1,500 Enantia individuals in the field. This suggests a high population density, which could explain the high level of polymorphism.
High genetic diversity in butterflies has been associated with a wide geographic territory and the possibility of populations being isolated from others throughout evolutionary time. For example, Papilio machaon (Lepidoptera: Papilionidae) shows high haplotype (h) and nucleotide diversity (π) (h = 0.856, π = 0.0084) [64], as does Aglais urticae (Lepidoptera: Nymphalidae) (h = 0.9647, π = 0.00983), using mtDNA [65]. The intermediate levels of genetic diversity observed for species of the Enantia jethys complex could reflect their relatively small distribution ranges.

Hybridization in the Enantia jethys complex
In contrast to DNA sequences, ISSR molecular markers provided evidence on the hybridization process in the Enantia jethys complex. Although some studies demonstrated hybridization with  [66][67], these hybridization processes have an ancient origin. The choice of molecular markers in hybridization studies is fundamental, considering that genetic isolation is a characteristic of some regions of the genome and not of the entire genome [68]. Some markers may hybridize and/or introgress further/faster than others [68]. This can explain the differences in results between ISSR and DNA sequences. ISSR allows us to observe multilocus variation at many independent loci (random primer amplifying nuclear noncoding region), while DNA sequences show variation at a single locus (with results potentially more affected by natural selection and stochastic effects, and generally showing lower variability [69]) [70]. Consequently, ISSR molecular markers are more useful in identifying hybrids [71] and studying biological processes (e.g., diversification, dispersion, and hybridization) for species of relatively recent origin [72][73][74], like the hybridization processes observed in the Enantia jethys complex, which could be considered contemporary. Llorente-Bousquets [15] suggested that the Enantia jethys complex is a very recent group. Our molecular clock analysis supports this hypothesis, suggesting that the species of the Enantia jethys complex form a clade with a recent divergence starting around 1 Mya. Our molecular results suggest a directional pattern for hybridization from E. albania to E. jethys and E. mazai. Some field observations (e.g., mating between species) combined with molecular information support this hypothesis. First, in the field, interspecific mating events were observed between E. albania and E. jethys, and between E. mazai and E. jethys, but no genetic profiles corresponding to a combination of E. jethys and E. mazai were identified among the 117 hybrid individuals detected in this study. Second, very few hybrids were observed in the E. albania group; and third, all hybrid individuals from the two other morphospecies contained a genetic component of E. albania. These results could suggest some reproductive barriers among species of the complex and differences in the viability of hybrids, with the existence of a prezygotic barrier (occuring after mating or gametic contact, thereby reducing the probability of fertilization [75]) between E. mazai and E. jethys, while mating involving E. albania does not appear to have a strong prezygotic barrier and hybrids are probably fertile. Prezygotic isolation is favored by natural selection as a process against unfit hybrids, as shown for Agrodiaetus butterflies [76]. Many studies on a variety of taxa ( [75] for review) have shown that prezygotic barriers are important barriers to prevent gene flow between closely related species. Other types of reproductive barriers are reported in other butterflies, such as Aricia and Polyommatus, two non-sister cryptic butterfly species for which barriers seem to be precopulatory [77], or between two sympatric sister species (Leptidea sinapis and L. reali) that present premating reproductive isolation [78].
Many cases of interspecific hybridization have been reported for Lepidoptera (non-exhaustive list: [66][67]79]). It is not surprising to observe hybridization in the Enantia jethys complex, considering that closely related species are more likely to hybridize [80], and that hybridization appears to be facilitated by wider sympatric distributions of parental species [81]. Furthermore, species of the Enantia jethys complex are subject to the same environmental pressures as they share the same oviposition and larval host plants [15], which may also facilitate opportunities for interspecific mating.
A study on Lepidoptera showed that hybrid individuals may present intermediate morphology [82]. process in natural populations may act in different ways, such as increasing genetic variation and new gene combinations, which favor novel adaptations [12,24], or reducing hybrid fitness [13] because of a postzygotic barrier [75]. Our study demonstrates an interesting case with different effects of hybridization on genetic diversity. First, E. albania hybrids showed a lower, although not significant, genetic diversity value (h) than non-hybrids, probably a reflection of sample size. Second, hybrid individuals of E. jethys presented a low but significant decrease in genetic diversity value (h), which could reflect a reduction in hybrid fitness. Finally, hybrid individuals of E. mazai showed an important and significant increase in genetic diversity value (h). This is not surprising because non-hybrid E. mazai individuals presented the lowest value of genetic diversity. Therefore, when they cross with individuals with higher genetic diversity, hybrids presented higher genetic diversity values (h), which could result in novel adaptations. Nevertheless, additional ecological, behavioral, chromosomal, and genomic studies will be necessary to test our hypothesis.

Phylogenetic relationships in the Enantia jethys complex
A previous study of the Enantia jethys complex using morphological data supported the existence of three distinct species: E. albania, E. jethys, and E. mazai, subsequently subdivided in subspecies E. mazai mazai and E. m. diazi [15]. The first study that investigated this species complex through DNA sequences (COI barcode) [17] validated the existence of the species E. albania, E. jethys, and E. mazai, without any evidence for subspecies or deep intraspecific differentiation. However, these authors discovered a sister clade of E. mazai composed of specimens that were morphologically indistinguishable from E. albania, suggesting the existence of a potential cryptic species originating in a specific geographic area (Ahuaxentitla, Puebla, Mexico). Using a multilocus approach (mitochondrial and nuclear markers), we confirmed the existence of three clades in the Enantia jethys complex in agreement with the results reported by Jasso-Matínez and colleagues [17]. The different topology recovered through the Wg marker (E. mazai and E. jethys in the same clade and shared haplotypes, Panel A in S1 Fig) could be explained by incomplete lineage sorting, as observed in other Lepidoptera studies [83][84]. Incomplete lineage sorting frequently results in different patterns between nuclear and mitochondrial markers [85], which is common in species of recent divergence [84]. To understand this phenomenon, a multiple independent loci analysis was used to generate a robust phylogenetic hypothesis (S2 Fig), as recommended by Edwards and Bensch [86]. Overall, our results confirmed the existence of three well-defined groups, corresponding to the three different morphospecies previously described [15][16][17].