The Andes Mountains are a continuous chain of mountains along the western coast of South America with average altitudes of ~4,000 meters above sea level, whose formation has had a profound impact on the geography and climate of this continent (Hartley, 2003; Moreira-Muñoz, 2011). For instance, in Chile, these mountains define a unique, long, and narrow coastal environment, which exhibits a highly diverse climate due to the wide latitudinal variation (Moreira-Muñoz, 2011). The uplift of the Andes has had a profound effect on the climate of the Southern Hemisphere, creating an extensive barrier to atmospheric circulation and promoting biogeographic changes that have affected the evolutionary history of local fauna and flora (Pirie et al., 2006; Antonelli et al., 2009). The uplift of the Andes has exerted a significant influence on the dispersal and diversification of a variety of South American organisms, such as plants, mammals, insects and specifically in triatomines (Mauseth, 1996; Doan, 2003; Struwe et al., 2009; Monteiro et al., 2018; Pavan et al., 2021).

The genus Mepraia Mazza, Gajardo and Jörg, 1940 is a Triatominae Jeannel, 1919 endemic to arid and semi-arid regions (Frías-Lasserre et al., 2019). This previously considered monotypic genus (Lent et al., 1994), comprises three distinct species: Mepraia spinolai (Porter, 1933) Mepraia gajardoi (Frías-Lasserre et al., 1998), and Mepraia parapatrica (Frías-Lasserre, 2010). These species are of significant public health importance as they are involved in the sylvatic cycle of the etiologic agent (Chacón et al., 2016; Lambarri et al., 2018). M. gajardoi depends on seabirds as its main food source, as well as on lizards that inhabit coastal regions. M. parapatrica also feeds on these species but has been observed to include domesticated mammals in its diet when it invades human settlements (Frías-Lasserre, 2010). In contrast, M. spinolai prefers wild and domesticated mammals and exhibits a strong tendency to invade human habitats in Andean valleys, although it never feeds on seabirds (Frías-Lasserre, 2010).

Mepraia has unique attributes among the triatomines, such as diurnal feeding habits and wing polymorphism (Faúndez and Carvajal, 2012; Frías-Lasserre et al., 2019). In summary, all Mepraia females are micropterous (with reduced wings), which may be a new or apomorphic characteristic of Triatominae (Frías-Lasserre, 2010). On the other hand, all M. gajardoi males are brachypterous (with short wings) (Frías-Lasserre et al., 1998), and although M. parapatrica males are mostly brachypterous, some are macropterous (with long wings) (Frías-Lasserre, 2010). In contrast, M. spinolai has both macropterous and brachypterous males (Frías-Lasserre, 2010; Frías-Lasserre et al., 2019).

The studies on geographical distribution (Campos-Soto et al., 2020), molecular data (Campos-Soto et al., 2013; Campos-Soto et al., 2015; Campos-Soto et al., 2022), morphology (Campos-Soto et al., 2011; Frías-Lasserre et al., 2018; Frías-Lasserre et al., 2019), and cytogenetics (Frías-Lasserre et al., 1998; Frías-Lasserre et al., 2019) provide support for the three species of Mepraia. However, although phylogenies with few genes support the described taxa (Calleros et al., 2010; Campos-Soto et al., 2013), there is no consensus on the genus diversification process.

Two hypotheses have been proposed to explain the diversification of Mepraia. The first suggests that the genus has originated due to a vicariant process followed by dispersion and isolation, which was potentiated by the arid climate of the region (Frías-Lasserre, 2010). This hypothesis is based on morphological and cytogenetic data which suggest that after the emergence of the Andes, and their separation from their ancestor, a process of dispersion and speciation occurred in a North–South direction (Frías-Lasserre, 2010). According to this hypothesis, M. gajardoi is considered the most ancestral species, while M. parapatrica gave rise to M. spinolai after some cytogenetic alterations (Frías-Lasserre, 2010). This latter study suggests that the wing morphology supports the hypothesis that the macropterous wing structure of M. spinolai is derived from M. parapatrica and M. gajardoi, which have winged males with ancestral characteristics (Frías-Lasserre et al., 1998; Frías-Lasserre and Atria, 1998; Frías-Lasserre, 2010).

The second and most recent proposal suggests that the diversification of Mepraia was also driven by geological changes and climatic fluctuations over the past 5 million years (Campos-Soto et al., 2022). In this model, there were two geographically isolated clades (Campos-Soto et al., 2022). The northern group gave rise to M. parapatrica, M. gajardoi, and a new lineage, while the southern group formed M. spinolai. However, this hypothesis is based only on data from three mtDNA genes (Campos-Soto et al., 2022), and further studies are needed to confirm its validity.

Although there is still no definitive consensus on the processes that led to the diversification of Mepraia, recent research indicates that it likely occurred rapidly (Justi et al., 2016; Campos-Soto et al., 2020; Campos-Soto et al., 2022). Several lines of evidence, such as genetic and phenotypic differences between populations, geographic and reproductive barriers, and molecular dating, suggest that there was rapid species diversification (Calleros et al., 2010; Campos-Soto et al., 2022), though, currently, there are only a few studies that investigate levels of population isolation and gene flow between them. Resolving the phylogenetic relationships of groups that diverged recently can be challenging (Seehausen et al., 2014; Marques et al., 2016; Weng et al., 2020; Congrains et al., 2021) due to the retention of ancestral polymorphisms causing incomplete lineage sorting (ILS), introgression, and/or gene duplication (Mallet, 2005). This scenario can become even more complex if there has been ancient gene flow, which can produce conflicting phylogenetic signals throughout the genome (Zhou et al., 2016; Edelman et al., 2019). Multi-species coalescent methods have been useful in clarifying evolutionary relationships and mitigating the impact of deep coalescence (Degnan and Rosenberg, 2009; Young and Gillung, 2020).The combination of next-generation sequencing technologies and robust phylogenetic methods has effectively resolved ambiguous phylogenetic relationships for a variety of organisms (Rokas et al., 2003; Dunn et al., 2008; Regier et al., 2010; Irisarri et al., 2017) and may be useful in the study of Mepraia.

This study aimed to investigate phylogenomic relationships in the Mepraia genus for the first-time using transcriptome data from the heads and salivary glands of their three species. To do so, we also evaluated the impact of gene flow on the genus’ speciation by using different methodological approaches, which increases the reliability of the results. Concatenation-based approaches and multispecies coalescent methods were employed for this purpose, and the study also assessed the impact of ILS and introgression events on the inferred phylogenies.

This study contributed important genetic resources for the subfamily Triatominae, generating transcriptome of heads and salivary glands of M. gajardoi, M. parapatrica and M. spinolai. These species have epidemiological importance, as they are part of the wild cycle of T. cruzi in Chile (Campos-Soto et al., 2016). The 18 RNA libraries (six samples per species) generated transcriptome assemblies with similar quality metrics to other available Triatominae transcriptomes (Assumpção et al., 2011; Marchant et al., 2015; Carvalho et al., 2017; Brito et al., 2019). Transcriptome assemblies from the three species of Mepraia showed similar metrics, with N50 on average of 1,300 bp, which is in the range of N50 metrics of other available. Completeness and redundancy metrics are crucial for genomic studies, and in addition to N50, we used the BUSCO tool to estimate completeness using the Arthropoda database. Our transcriptomes contained between 71% and 90% complete orthologs, indicating a high percentage of complete copies for most predicted genes. The number of duplications and CDS annotations in downstream analysis correlated with the transcriptome quality.

Our transcriptomes provided a vast amount of data that allowed us to investigate the phylogenetic relationships of Mepraia. To identify strictly orthologous nuclear gene clusters, we employed the strategy described by Yang and Smith (2014) for non-model organisms, resulting in a total of 3,245 gene clusters. Our supermatrix has an overall occupancy of 72%, with about 28% of data missing from the matrix. Although the occupancy rate varied for each taxon, with some having a higher occupancy rate than others, there are an average of 3,200 loci per sample. The maximum likelihood (ML) tree inferred from the supermatrix produced a well-resolved topology with high support (bootstrap values >90), which is consistent with other studies (Justi et al., 2014; Campos-Soto et al., 2015; Campos-Soto et al., 2020), and corroborates the three lineages described in Mepraia, as well as infers M. gajardoi and M. parapatrica as sister groups and M. spinolai as more distant.

Based on mitochondrial genes, Bayesian analysis, and the Kimura two-parameter (K2P) model, we found low intra- and interspecific variability among the sequences, which was expected due to the close relationships among the species. The average genetic diversity rate for the 57 sequences was 0.20. It is worth noting that mitochondrial genes, such as COI, are widely used for species delimitation in triatomines (Hypša et al., 2002; Monteiro et al., 2003; Justi et al., 2014). However, some species, such as Triatoma maculata (Erichson, 1848), exhibit high levels of genetic distance, as evidenced by the CytB and ND4 genes (Gómez-Palacio et al., 2022). On the other hand, our findings are consistent with the low genetic distance observed in mtDNA between closely related species of the T. sordida subcomplex, as reported by Belintani et al. (2020).

The FBD method to estimate diversification time showed that the genus underwent recent divergence 2.3–1.2 Mya. Although the lack of morphological data on the fossil species P. hispaniolae used for calibration may have influenced our analysis, our findings are consistent with other available dating methods (Campos-Soto et al., 2020; Campos-Soto et al., 2022). Our dating also supports the hypothesis that Mepraia diverged from its Triatoma ancestors after the uprisal of the Andes mountains, which was estimated to have occurred around 5–2 Mya (Gregory-Wodzicki, 2000; Frías-Lasserre, 2010). The rise of the Andes and the resulting changes in climate have been identified as crucial factors that contributed to the diversification of the Triatomini and Rhodniini tribes in South America (Justi et al., 2016), which may have also influenced the diversification of Mepraia.

The phylogenetic analysis of Mepraia based on over 60,000 unlinked biallelic SNPs produced consistent phylogenies. Both support the relationships among species, with M. garjadoi and M. parapatrica inferred as sister taxa and M. spinolai as more distantly related. Although the phylogeny constructed with the quartet model has robust support, the Bayesian method used in the SNAPP software revealed a significant number of alternative topologies. The inferred alternative topologies may indicate that the evolutionary relationships between the species in question are uncertain or that there is a complex evolutionary history for these species that cannot be fully resolved with the available data set (Leaché and Oaks, 2017; McLean et al., 2022).

The coalescent-based species tree inferred using quartet inference in ASTRAL-III resulted in a topology consistent with the one produced by the concatenation method, as well as to that produced by mtDNA. However, despite the high bootstrap support values across the inferred tree, the normalized quartet scores were relatively low. The exact threshold for robust quartet support may vary depending on the study and the objectives of the analysis (Whidden and Burch, 2015). In general, quartet scores above 0.5 are considered to provide moderate support, while scores above 0.7 provide strong support for the topology. In our analysis, the inferred clades in the species tree have quartet values ranging from 43% to 56%, which prompted us to perform exploratory analyses to assess the potential conflict among gene trees across the Mepraia phylogeny.

A bipartition analysis of trees among different parts of the genome using PhyParts showed strong support for the general topology of species trees, while still indicating alternative topologies. This analysis is especially valuable as it offers a counterpoint to tree species, revealing that while we can identify a major topology, we can also clearly infer the existence of a significant number of incongruent topologies. Is worth noting that high posterior probability values obtained in ASTRAL-III can be misleading when evaluated alone and may obscure the underlying phylogenetic signal (Maddison and Knowles, 2006; Zhou et al., 2016).

Our results suggest that there are only a small number of concordant trees, indicating that several alternative topologies are being supported. This pattern is expected, as gene trees are rarely entirely concordant with species trees (Vargas et al., 2019; Meng et al., 2021; Xiao et al., 2022), especially in recently diverged taxa. Incongruence in the phylogenetic signal across distinct parts of the genome can be the result of various evolutionary phenomena, including incomplete lineage sorting, gene flow, hybridization, lateral gene transfer, recent divergence, convergent evolution, different evolutionary rates, or loss of information during evolution (Betancur et al., 2019). These phenomena shape the relationships among species, which, depending on the degree of conflict, can make phylogenetic inferences incredibly challenging. The conflict between the inferred species tree and the underlying gene trees indicates that some of the above-mentioned processes could have influenced the differentiation of Mepraia, especially incomplete lineage sorting and/or introgression.

Although phylogenies based on coalescent models at the species and gene level have produced similar relationship among Mepraia species, our data suggest a reticulated evolutionary history. As indicated above, several factors, especially ILS and introgression, can produce reticulate phylogenies, although distinguishing between the latter two can be challenging (Schumer et al., 2014; Wang et al., 2022). As ILS can occur even in the absence of introgression events, it is important to include it as null hypotheses when testing for introgression to properly investigate the causes of reticulate phylogenies (Mendes and Hahn, 2018; Hibbins and Hahn, 2022), which we did using different strategies.

Our study provides robust evidence of introgression between Mepraia species using four different methods. Firstly, an ABBA-BABA test considering only the most likely topology estimated a significant value of the D statistic (D = 0.158, Z-score 4.43), indicating introgression. Secondly, the HyDe test also revealed significant introgression between the three species, providing further support for introgression. Additionally, PhyloNet analyses showed strong evidence that M. parapatrica may result from introgression between M. gajardoi and M. spinolai. Finally, the JML test allowed identification of hybridization events in all comparisons (p-value <0.1).

Knowledge of current gene flow between Mepraia species is limited (Campos-Soto et al., 2013). Although hybridization has been observed in other Triatominae species (Panzera et al., 2021), there is still much to be discovered about the extent and effects of hybridization in Triatominae (Costa et al., 2009; Panzera et al., 2021; Pavan et al., 2021; Pinotti et al., 2021). To date, the only study that has investigated inter-specific crossing ability in the Mepraia genus, particularly between M. gajardoi and M. spinolai, found evidence of hybrid incompatibilities in their second backcross (Campos-Soto et al., 2016). Hybrid speciation involves the mixing of different lineages and gene flow through hybridization and backcrossing, resulting in the formation of a new species (Mallet, 2005; Mallet, 2007; Abbott et al., 2013). This process was previously thought to be rare in animals, but recent methodological and technical advances in the study of introgression (Guo et al., 2018; Lamichhaney et al., 2018) have shown its importance in evolution and prevalence in closely related species of several different organisms, such as plants, birds, fishes, and butterflies (Schwarz et al., 2005; Lamichhaney et al., 2018; Zhang et al., 2019).

The parapatric evolution model proposed for Mepraia suggests that cytogenetic changes were responsible for the emergence of M. parapatrica and later M. spinolai from the ancestral M. gajardoi (Frías-Lasserre, 2010). According to this hypothesis, character displacement among the three species in tension zones indicates absence of gene flow (Frías-Lasserre, 2010). However, the more recent speciation proposal by Campos-Soto et al. (2022) advocates a speciation model that involves two geographically distinct groups: the North, composed of M. garjadoi and M. parapatrica, and the South, composed of M. spinolai. The data presented in that study suggests an alternative model that may be compatible with ancient hybridization events or continued gene flow. In this sense, it is important to consider that the evolution of these species may not have occurred in a linear and bifurcating manner, and that the phylogenetic relationships between the species may not be described in a simplified way, but rather by diversification patterns that have been shaped by historical events of hybridization and geographic segmentation. This broader view of evolution may help to better understand the complexity of evolutionary processes and to reconstruct the evolutionary history of species more accurately.

Hybridization occurs when certain conditions are met, including the mixing of two parental lines, followed by reproductive isolation and a viable, genetically distinct hybrid line (Schumer et al., 2014). Based on our data and the current knowledge about Triatominae speciation (Pavan et al., 2021), we can speculate that incongruity patterns in gene trees cannot be explained only by incomplete lineage classification and stochasticity. Rather, to some extent, it may have been modulated by introgression, potentially due to ancestral or continuous gene flow. As discussed earlier, geoclimatic changes in Chile over the past 5 million years may have been crucial for the evolution of these Triatomine bugs (Campos-Soto et al., 2022). Similar diversification events have been observed in North American Drosophila (Campo et al., 2013), Anolis lizards on Caribbean islands (Mahler et al., 2010), and the Heliconius butterfly genus (Mavárez et al., 2006; Jiggins et al., 2008). In the Andean region of Chile, it is known that the Bombus dahlbomii bumblebee hybridized with other Bombus bee species in Chile (Arbetman et al., 2013).

This study successfully obtained and analyzed the transcriptomes of the head and salivary glands of three species of Mepraia, shedding new light on the relationships between M. gajardoi, M. parapatrica, and M. spinolai. Analyses of a large transcriptome dataset provided valuable information about the evolutionary history of the genus. We showed that while all methods produced concurrent phylogenetic relationships among species of the genus, a considerable level of incongruence is observed across gene trees. We demonstrated through several tests that such incongruence is not compatible only with ILS, rather, is an indication of introgression between lineages due to past hybridization events. Therefore, our data suggests that the diversification of Mepraia occurred recently and was significantly influenced by interspecific gene flow during its early stages of divergence. This study highlights the importance of considering introgression and hybridization in the evolution of organisms, especially in Chagas disease vector triatomines.

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/genbank/, PRJNA916468.

TB designed the study, wrote the entire text, discussed results, and performed all study analyses. CC assisted and reviewed all bioinformatics analyses, as well as reviewing the text, results, and discussion. HP and SC-E assisted in RNA extraction and reviewed the text. RB reviewed the text, assisted in the discussion, and provided servers for bioinformatics analyses. JO collected specimens, assisted in the discussion, and reviewed the text. DF-L collected specimens, assisted in the discussion, and reviewed the text. FF mapped and reviewed the test. JR designed the study, collected specimens, assisted in the discussion, reviewed, and obtained resources. All authors contributed to the article and approved the submitted version.