For this study, field specimens and laboratory strains were used with a special focus on the prolixus group for field specimens (Table 2). From our field collection, were selected specimens with doubts about their morphological identification (MILEP, NEGP), difficult to determine (INCP, NasG), or exhibited incoherent cytb sequences in a preliminary study (ROBB, ROBR, ROBQ). Some specimens were very old since they were also used by Harry (1993) and Harry (1994). We used laboratory strains from the Insetário de Triatominae da Faculdade de Ciências Farmacêuticas/Unesp/Araraquara, Brazil. Samples were either stored in absolute EtOH or air-dried.

DNA was extracted from legs and alary muscles from adult bugs, from 1 for field specimens to 6 for laboratory strains using the Qiagen DNEasy tissue kit. Genomic data was obtained from 36 samples of the Rhodniini tribe corresponding to 17 species identified using morphological characters. A whole-genome shotgun sequencing was performed using Illumina technology (100 bp paired-end, Imagif platform, Gif-sur-Yvette, France) generating from 3.7 to 21.1 Gb data per sample corresponding to a genome depth from 5 to 30x (Table 3).

Sequences have been deposited in the NCBI BioProject database via the accession number PRJNA429761.

In order to assemble the mitogenomes, the reads were subsampled (500 kb of reads for each sample) and assembled using the Trinity software with default parameters (Haas et al., 2013). This approach allows to assembly sequences that are highly abundant in the sequence pool like mitochondrial or ribosomal genes. In insect, mtDNA represents on average 0.42% of the total DNA in genome sequence project (Meng et al., 2019). As we used 500 kb of data/sample for the mitochondrial assembly which in average contains (500 × 0.42)/100 = 2,100 kb of mtDNA and given that the mitochondrial genome is about 15 kb, we obtain coverage of 140x for this genome. This assembly strategy also prevent to pick up some NUMS that are mitochondrial pseudogenes in nuclear chromosomes, because if these elements were present in the data, as they are nuclear, their coverage would be very low (from 5 to 30x according to the samples) compared to that of the (true) mtDNA and with a low probability to be retrieved from the data used to assemble the mtDNA. Contigs corresponding to the mitogenomes and the nuclear rDNA operon were identified by BLAST using as templates the Triatoma infestans mitogenome and the R. prolixus reference genome rDNA. In most cases, the mitogenome and the nuclear rDNA operon of each sample appeared as a single contig, sometimes as two or three contigs.

The annotation of the Rhodnius mitogenomes was achieved using the MITOS pipeline (Bernt et al., 2013). Each mitochondrial and nuclear rDNA datasets was independently aligned with MAFFT¹ (Katoh et al., 2019), visualized with Bioedit (version 7.2.6.1; Hall, 1999), and ambiguously aligned regions or gapped positions were manually corrected or removed. Two concatenated datasets were generated for the mitochondrial data. The first included the 15 mitochondrial genes (15 mtG), namely the 13 mtPCG (Nad2, COI, COII, Atp8, Atp6, COIII, Nad3, Nad5, Nad4, Nad4L, Nad6, Cytb and Nad1) and the 2 rDNA genes (16 rDNA and 12 rDNA). The second dataset including only the 13 mitochondrial protein-coding genes (mtPCG) was used for the selection signature analyses. This dataset was translated in order to test for the presence of pseudogenes (Supplementary Table 1).

For the nuclear rDNA genes (nu-rDNA), we retrieved from the data the complete 18S and 28S genes but the partial ITS1, 5.8S, and ITS2 genes. The alignment for this dataset is given in the Supplementary Table 1.

A second assembly was carried out for each of the 36 genomes with the SOAPdenovo2 software (Luo et al., 2012) with k-mers estimated using the Kmergenie program (Chikhi and Medvedev, 2014). Benchmarking Universal Single-Copy Orthologs (BUSCO) genes (Waterhouse et al., 2019) were searched in assemblies using the insecta_odb10 dataset comprising 1,367 genes universal to all insects. In a first step, we kept the 88 nuclear protein-coding genes (nu-PCGs) for which an entire or fragmented copy was found in each of the 36 Rhodniini genome assemblies and in the T. brasiliensis transcriptome (Marchant et al., 2015), as this species was used as an outgroup. Each set of genes was aligned with MAFFT (global homology option; Katoh et al., 2019) and then concatenated. The concatenated alignment was visualized using Bioedit (version 7.2.6.1; Hall, 1999). All regions without at least 29 aligned genomes were retrieved as well as gaps or misaligned sequences. This drastic final trimming resulted in 51 nu-PCGs with a sequence coverage of at least 80% except for R. brethesi (Table 3). In order to have R. prolixus samples validated without (or low) introgression, we also search the 51 nu-PCGs in both the R. prolixus used by Mesquita et al. (2015) for the reference genome (VectorBase: RproC3.4) and a Honduras strain reared in a French laboratory (Table 1). The corresponding BUSCO genes and the alignment are given in the Supplementary Table 1.

In order to compare the variation rate between the datasets, the nucleotide diversity (Pi) was estimated for each dataset using DnaSP 6 (Rozas et al., 2017). Pairwise distances were calculated between samples using MegaX (Kumar et al., 2018) using the Maximum Composite Likelihood model and all the substitutions. The similarity between pairwise samples was calculated in percentage using the ratio of the total number of conserved sites and the total number of both conserved and variable.

All the individual mitochondrial genes and the three concatenated datasets namely the 15 mtG, the nu-rDNA, and the nu-PCG alignments, were submitted to phylogenetic analyses. For each dataset, we used as outgroup the closest relatives identified in the sequence database. For the mitochondrial dataset, as some mitogenomes were sequenced for Triatominae, we used T. dimidiata (AF301594; Dotson and Beard, 2001) and Panstrongylus rufotuberculatus (NC_042682; Zhao et al., 2019). For the nu-rDNA dataset as only partial data is available for Triatominae, we used as outgroup the Reduviidae Brontostoma colossus (NC_024745; Kocher et al., 2014) and for the nu-PCGs, T. brasiliensis (Marchant et al., 2015). For the three concatened alignments, the sequence coverage is given for each sample in the Table 3. The optimal partitioning scheme and models of nucleotide substitution for each concatenated dataset were first analyzed using PartitionFinder v2.1.1 (Lanfear et al., 2012). For the concatenated mt-PCGs, alignments were also separated into single codon positions. Optimal partitioning and nucleotide substitution models were subsequently analyzed using the Bayesian Information Criteria (BIC) for all datasets. For the concatenated datasets two types of trees were generated. MrBayes (Ronquist et al., 2012) was used, to generate Bayesian Inference (BI) topologies using the partitioning scheme and models of evolution suggested by PartitionFinder and by performing two independent 10-million generation runs. Maximum Likelihood phylogenies (ML) were inferred using IQ-TREE (Minh et al., 2020) using the partitioning scheme and models of evolution suggested by PartitionFinder and the approximate likelihood ratio test with the non-parametric Shimodaira–Hasegawa correction (SH-aLRT) (Guindon et al., 2010) was used to estimate the support for each node with 1,000 replicates. For individual mitochondrial genes, Maximum Likelihood phylogenies (ML) were inferred using IQ-TREE (Minh et al., 2020) after searching the best substitution model with ModelFinder (Kalyaanamoorthy et al., 2017) and SH-aLRT branch test with 1,000 replicates. Phylogenetic trees were displayed and modified using iTol (Letunic and Bork, 2019). Gene partitions and best substitution models for the concatenated dataset are given in the Supplementary Table 1.

Two approaches were used to unravel potential selective pressures during the evolution of Rhodnius mtDNA genomes and more specifically to identify whether the specific groups, namely pallescens, pictipes and prolixus, evolved at different rates and were under unique signatures of selection.

First, we applied the McDonald–Kreitman test (McDonald and Kreitman, 1991) using DnaSP 6 (Rozas et al., 2017) to investigate the rates of adaptive evolution in mitochondrial and nuclear genes. We used the test to compare synonymous and non-synonymous variations within and between species groups. Under neutrality (estimated with the neutrality index NI and the Fisher exact test), the ratio of replacement to synonymous fixed substitutions between species groups should be the same as the ratio of replacement to synonymous polymorphisms within species groups.

Positive selection was also investigated using branch-site models to test for differential selection on the three species groups with CodeML (PAML package v4.9; Yang, 2007) using the graphical user interface EasyCodeML (Gao et al., 2019). The mt-PCG and the nu-PCG datasets were used for this test but also the four mitochondrial genes usually used for Rhodnius phylogenies, namely Cytb, COI, COII, Nad1.

To detect signatures of introgression in the pallescens and pictipes groups, we used the ABBA-BABA test (Green et al., 2010). The ABBA-BABA test considers that for phylogenetically-informative di-allelic sites in a four-taxon bifurcating tree, the most frequent configuration should be (A,(A,(B,B))), congruent with the species tree. Incomplete lineage sorting or homoplasy could lead to equal amounts of loci with (A,(B,(B,A))) and (B,(A,(B,A))) incongruent configurations. However, introgression could be detected if either one of the two configurations significantly exceeds the other. For the pallescens group, we used the three samples, R. pallescens R1J, R. colombiensis R10B and R. ecuadoriensis, ECUD with R. neivai NEII as the outgroup and for the pictipes outgroup, R. stali STAWZ, R. brethesi R9WX, R. pictipes R63K and R. amazonicus Ama1A as the outgroup. Raw Illumina reads of each species were mapped to the R. prolixus reference sequences² using minimap2 (v. 2.17-r941; Li, 2018) with default parameters. The mapped reads were filtered for a mapping quality of 20 with samtools (v. 1.9; Li et al., 2009). BAM files were cleaned by removing unmapped reads and sorted by coordinate using Picard v. 1.4.2³. The sorted BAM files were then combined in a mpileup file with samtools and finally converted to a synchronized file (one species per column) with custom scripts from PoPoolation2 (v. 1.201; Kofler et al., 2011). Using a customized perl script on the synchronized file, we counted the three configurations for all diallelic loci with a minimum coverage of 10 reads in the four species.

The mitochondrial dataset is composed of 36 samples from the Rhodniini tribe and the 15 mitochondrial genes. The concatenated alignment totalized 13.3 kb with a nucleotide diversity (Pi) of 0.12389. No pseudogene was retrieved from our dataset as attesting by the functional proteins obtained after the traduction (Supplementary Table 1). For the concatenated dataset mtG, highly congruent trees were obtained from both the Maximum-Likelihood (Figure 1) and the Bayesian Inference (Figure 2) methods and showed pallescens and pictipes groups as sister clades with high statistical support (ML: 98.8, BI:100). For the individual genes, various topologies were obtained (Supplementary Figures 1–15). For nine genes, the same topology as for the concatenated dataset was obtained, namely for Cytb (Supplementary Figure 8), Nad1 (Supplementary Figure 9), Nad5 (Supplementary Figure 14), Nad6 (Supplementary Figure 15) and also for Atp6 (Supplementary Figure 3), Atp8 (Supplementary Figure 4), COI (Supplementary Figure 6), COIII (Supplementary Figure 7), and Nad2 (Supplementary Figure 10), but for these five genes a paraphyly was observed for the prolixus group, R. neivaii, R. domesticus, and/or P. tertius, according to genes being outside the prolixus group. For four genes, the pictipes and the prolixus groups were sister clades, namely for 12S (Supplementary Figure 1), 16S (Supplementary Figure 2), Nad3 (Supplementary Figure 11), Nad4 (Supplementary Figure 12), and the Nad4L but for this gene the pictipes group is less supported (Supplementary Figure 13). The third topology, that is the pallescens and prolixus groups as sister clades, was displayed for the COII gene (Supplementary Figure 6) but the position of R. neivaii and R. domesticus were weak supported.

The analysis of rare insertion/deletion (InDel) patterns in the alignment indicated the presence of an InDel of 40 pb in the tRNA Cys that supported the genetic proximity between the pallescens and the pictipes groups and the split with the prolixus group (Supplementary Figure 16).

According to the concatenated mtG trees (Figures 1, 2), for the pictipes group, R. stali and R. brethesi were closer to each other than to R. pictipes. In the pallescens group, R. pallescens and R. colombiensis were closer to each other than to R. ecuadoriensis. The prolixus group comprised samples of R. domesticus, R. neivai, R. montenegrensis, R. taquarussuensis (syn R. neglectus), R. neglectus, R. nasutus, R. robustus, R. prolixus and also P. tertius. Indeed, the Psammolestes species was robustly included in the prolixus group independently of the considered mitochondrial dataset. However, in the two trees, ML (Figure 1) and BI (Figure 2), the genetic proximity between the three clades (R. montenegresis-R. robustus/R. nasutus-R. neglectus/R. prolixus) is not resolved. For R. amazonicus some discrepancies were observed according to the mitochondrial dataset. For the concatenated dataset, this taxon was related to the pictipes group at a basal position but for the individual datasets it was either related to the pallescens group or the pictipes group.

Three Brazilian specimens seemed to have an aberrant position in the trees. Specimens NEGP and MILEP, morphologically identified as R. neglectus and R. milesi respectively, were related to R. prolixus. The NASG specimen, identified as R. nasutus, is at a basal position in the ML tree (but with a weak bootstrap) and nested within the prolixus group but not related to the R. nasutus-R. neglectus clade in the BI tree. The non-identified INCP sample was related to R. nasutus.

For samples identified as R. robustus and R. prolixus, two clades were inferred, each of them comprising samples from the two species. One is composed of all the R. robustus samples except the field specimen RobQ from Brazil, the laboratory strain R. montenegrensis and one R. prolixus sample, YRP, a Venezuelan laboratory strain from Cojedes. The other grouped all other R. prolixus samples with the RobQ.

Concerning species with questionable species-specific status in addition to the topologies observed in the phylogeny, the pairwise distances (Supplementary Table 2) and the similarity index can be considered using the mitochondrial PCG concatenated dataset. R. neglectus (R5WU) and R. taquarussuensis (syn R. neglectus) (R3WT) appeared to be very closely related taxa (pairwise distance = 0.0012; similarity = 99.88%). Unexpectedly, in the phylogeny, R. milesi (MILE) was included in the same clade as the laboratory reference R5WU R. nasutus sample (MILE R7WU: 0.0069; 98.29%). More distant pairwise values were obtained between R. montegrenensis and R. robustus, the smallest being with the R. robustus samples from Peru (R8F RobR: 0.0149; 98.56%) and the highest with the Amazonian ones, the Brazilian ROBB (R8F ROBB: 0.0311; 97.13%) and the Guyanese Rob5s (R8F Rob5s: 0.0326; 97%). For comparison, the pairwise values between R. nasutus (RW7) and R. neglectus (R5WU) that are close but clearly distinct species, were smaller (0.0202; 98.09%) than those between R. robustus and R. montenegrensis.

The nuclear rDNA dataset included 36 taxa from the Rhodniini tribe for an alignment length of 5.2 kb and a Pi of 0.0089. ML and BI rDNA phylogenetic analysis confirmed the three major Rhodniini groups (pallescens, pictipes, and prolixus) but not supported by high bootstrap values (Figures 3, 4). However, an InDel of 19 pb in the 28S rDNA genes supported the genetic proximity between the pallescens and the pictipes groups (Supplementary Figure 16).

For the resolution within the Rhodnius groups, only the BI tree gave some reliable results (Figure 4). The pallescens group exhibited a good resolution level, with the two species R. pallescens and R. ecuadoriensis closer related than to R. colombiensis. For the prolixus group, only the R7WU R. nasutus is differentiated from the other R. neglectus samples including R. taquarenssensis (syn. R. neglectus) (R3WT). Two R. robustus specimens (RobR, Rob5s) were clustered with R. montenegrensis R8F. Note that Psammolestes was clustered with the R. prolixus species. The pictipes group was the least resolute.

The nu-PCG dataset included 36 taxa from the Rhodniini tribe. The trees were reconstructed using 51 PCGs for an alignment length of 36.3 kb with a Pi of 0.02648. The tree topologies between the reconstruction methods, ML (Figure 5) and BI (Figure 6) are fully congruent and robust.

As for the mitochondrial dataset, the three groups, pictipes, pallescens, and robustus, were robustly identified (MI and BI trees) with the same topology with pallescens and pictipes as sister groups. Moreover, the same topology was found for the pallescens group, R. pallescens and R. colombiensis being closer to each other than to R. ecuadoriensis. Discrepancies were found for the pictipes and the prolixus groups. Unlike the relationship revealed with mitochondrial data, R. stali and R. pictipes appeared with nuclear data closer to each other than to R. brethesi. Similar to what was observed with the concatenated mitochondrial dataset, R. amazonicus was related to the pictipes group at a basal position. For R. robustus and R. prolixus samples, two well supported related clades were observed. The Guyanese R. robustus (Rob5s), the Brazilian R. robustus (ROBB) and the R. montenegrensis (R8F) samples were grouped in the same clade. All the other R. robustus and R. prolixus samples were in another clade.

Concerning the pairwise distances (Supplementary Table 2) and the similarity index, R. neglectus (R5WU) and R. taquarussuensis (syn R. neglectus) (R3WT) were very closely related (pairwise distance = 0.00076; similarity = 99.82%), and R. milesi (MILE) was closest to R7WU R. nasutus (MILE R7WU: 0.00192; 99.81%). More distant pairwise values than for the two previous pairs of species were obtained between R. montegrenensis and R. robustus, e.g., RobB (R8F RobB: 0.00309; 99.62%) and Rob5s (R8F Rob5s: 0.00408; 99.59%) and were of the same order of magnitude as between the Guyanese R. robustus (Rob5s) and R. prolixus samples e.g., Pro100 (Rob5s Pro100: 0.00485; 99.52%). Similar distant pairwise values were obtained between R. nasutus and R. neglectus (R7WU R5WU: 0.0036; 99.64%) or between the R5WU R. neglectus and INCP (R5WU INCP: 0.004; 98.16%) or NASP (R5WU NASP: 0.00396; 99.65%). It should be of interest to note that NasG, placed outside the R. nasutus-R. neglectus clade for the concatenated mitochondrial dataset, was close related to R7WU R. neglectus (NASG R7WU: 0.00182; 99.82%).

The McDonald–Kreitman test detected significant positive selection in the concatenated mitochondrial dataset for the three group comparisons, prolixus-pallescens groups (NI = 0.696; p-value = 0.023), prolixus-pictipes groups (NI = 0.609; p-value = 0.003), and pallescens-pictipes groups (NI = 0.779; p-value = 0.04), revealing an excess of divergent non-synonymous mutations between the group tested. For the concatenated nuclear PCG dataset no selection was detected for any pairwise, prolixus-pallescens groups (NI = 0.754; p-value = 0.41), prolixus-pictipes groups (NI = 0.583, p-value = 0.199), and pallescens-pictipes groups (NI = 0.931; p-value = 0.792).

Selective pressure was also tested using PAML method (Table 4). For the concatenated mitochondrial dataset, selective pressure was detected in each Rhodnius group test against the others (LRT p-value = 2 × 10^–9 for pallescens vs. pictipes and prolixus groups; 0.0044 for pictipes vs. pallescens and prolixus groups; 9.6 × 10^–7 for prolixus vs. pictipes and pallescens groups). For the single mitochondrial genes tested, namely cytb, COI, COII, Nad1, all were under selective pressure for the prolixus group tested against the pictipes and pallescens groups (LRT p-value = 13 × 10^–6; 13 × 10⁶; 0.0009; 0.02 respectively) and two, COI (LRT p-value = 0.046) and COII, for the pallescens group tested against the pallescens and prolixus groups (LRT p-value = 0.04). For the concatenated nuclear dataset no selective pressure was detected for any of the groups (LRT p-value = 0.47; 0.10; 0.14 for pallescens, pictipes and prolixus groups, respectively, tested against the others groups).

We conducted the ABBA-BABA test of introgression in the pictipes group for which a discordance between the mitochondrial and nuclear PCG trees was noticed. R. stali and R. brethesi were indeed clustered together in the mitochondrial tree while R. stali was closer to R. pictipes in the nuclear tree.

The three species R. brethesi, R. stali and R. pictipes were tested for introgression and R. amazonicus as the outgroup. The ABBA-BABA results were in accordance with the nuclear phylogeny of the species and provided a strong support for a mitochondrial R. stali introgression in R. brethesi. Out of 6,437 bi-allelic SNPs, we counted 5,571 AABB, namely R. pictipes closer to R. stali, 409 ABBA, namely R. brethesi closer to R. pictipes and 457 BABA, namely R. brethesi closer to R. stali. No signal of further nuclear introgression was found between R. brethesi and R. stali since both ABBA and BABA counts were almost equal and far below the AABB (i.e., the nuclear tree topology) SNP counts.

We conducted the same analysis in the pallescens group, for which the mitochondrial and nuclear PCG trees are congruent showing R. colombiensis and R. pallescens clustered together, but in the rDNA phylogeny, R. pallescens and R. colombiensis are closer. We found a strong signal for pervasive and widespread phylogenetic incongruence between the three species, R. pallescens, R. colombiensis, and R. ecuadoriensis. Out of a total of 18,910 di-allelic SNPs, the ABBA-BABA test gave 7,666 sites with the configuration AABB, namely R. pallescens is closer to R. colombiensis, 4,084 with the ABBA configuration supporting a closer relationship between R. ecuadoriensis and R. colombiensis, and 7,160 with the configuration BABA supporting a closer relationship between R. ecuadoriensis and R. pallescens. Remarkably, the BABA counts even approached those of the AABB species tree counts, indicating that the exact position of R. pallescens could not be determined with certainty. In order to test if the ABBA or BABA configurations are restricted to specific genomic regions, the evaluated SNPs were mapped to the reference genome of R. prolixus. We did not observe any specific localization, the three configurations were widespread and overlapping on all scaffolds of the reference genome of R. prolixus.

In this study, our goal was to provide a large genomic dataset in order to elucidate key issues on phylogenetic relationships for yet the most taxonomically controversial Tritaominae tribe, the Rhodniini, using both mitochondrial and nuclear datasets analyzed separately. Previous Rhodnius phylogenetic studies used a limited number of genes, which could explain the divergent topologies found for the Rhodnius groups. Three topologies were recovered independently of the type of marker used. Considering studies with at least eight Rhodnius species and reliable phylogenetic reconstruction methods (Maximum Parsimony, ML, BI), the pictipes-pallescens clade was recovered using only mitochondrial markers, as cytb (Maia da Silva et al., 2007), as well as using a combination of mitochondrial and nuclear genes, as 16S + cytb + 28S (Monteiro et al., 2000). The pictipes-prolixus clade was obtained using also only mitochondrial markers: 16S (Paula et al., 2005, 2007) or 16S + cytb (Hypša et al., 2002), as well as mitochondrial and nuclear markers together, namely 18S + 28S + 16s + Cytb + COI + COII (Justi et al., 2014) or 18S + 28S + wg + 16S (Justi et al., 2016). It is to highlight that Paula et al. (2021) by combining the three genes, 16S, 18S, and wg, obtained contrasted results according to the type of analysis used, namely pictipes-prolixus clade with MP and pictipes-pallescens with ML. Using ultraconserved elements and/or rDNA data, Kieran et al. (2021) obtained the pallescens-prolixus clade. To limit the methodological bias due to a limited dataset, we used in this study a large genomic dataset of both mitochondrial (13.3 kb) and nuclear data (5.3 kb of rDNA; 51 PCGs, 36.3 kb). The large dataset was able to resolve the relationship between the three Rhodnius groups. Except for some single mitochondrial genes, both the concatenated mitochondrial and nuclear datasets showed pictipes and pallescens as sister groups, irrespective of the outgroup used.

The incongruence of phylogenetic signals between different genomic regions can result from the rate of variation among partitions. Our study indicates that the nuclear rDNA markers (including 18S, 5.8S, 28S, and ITS) led to a poorly resolved tree. Indeed, the nucleotide divergence of the rDNA marker dataset is low (max 1.99%) compared to the mitochondrial Rhodnius dataset (max 15.06%), and the nucleotide diversity (Pi) is 14 times slower than the mitochondrial dataset one. It results that slow-evolving rDNA genes are less suitable to infer phylogenetic relationships. Mitochondrial genes, which evolve more rapidly than the rDNA and the set of insect universal nuclear genes (BUSCO genes), seem to be candidates for more accurate phylogenetic reconstruction. In our dataset, based on Pi estimates, the mitochondrial genes evolve 4.7 times faster than the nuclear PCGs, which is in accordance with what was found in insects, the mitochondrial genes being estimated to evolve 2–9 times faster than the nuclear protein-coding genes (Lin and Danforth, 2004). Moreover, mitochondrial genes are easier to use for both technical and evolutionary reasons. Indeed, mitochondrial primers are widely available and as these mitochondrial genes are clonally inherited and non-recombining, they are not subjected to recombination and heterozygosity like the nuclear markers. Nevertheless, mitochondrial genes have some clear disadvantages (Lin and Danforth, 2004). Since all mitochondrial genes are linked to the same DNA unit, they cannot provide an independent estimate of phylogeny in the same way that unlinked single-copy, as the nuclear genes do (Harrison, 1989). Furthermore, the higher rate of substitution can be disadvantageous to resolve divergences between taxa of more than 5–10 million years, while they are useful to resolve relationships between closely related taxa that have diverged relatively recently.

In our study, we demonstrated that mitochondrial genes widely used for phylogenies, namely cytb and COI, were under selection for all the three groups, pallescens, pictipes and prolixus. By contrast, no selective pressure was evidenced for the nuclear PCGs. In the literature, there is increasing evidence that the genetic diversity of mtDNA may therefore be shaped not only by random genetic drift but also by natural selection given the functional importance of mitochondrial polymorphisms. Positive selection on the mitogenome has been shown for various species including vertebrates and insects (reviewed by Sun et al., 2018). This process may be related to energy production linked to migration but also aging, size, diet, neuronal function, and thermoregulation (see: Garvin et al., 2015; Sloan et al., 2017). For the Rhodnius species, selective pressure is likely linked to fitness-related traits that have to be investigated.

A delicate point of this work is to give molecular pertinent arguments about the species-specific status of some controversial or recently described species.

Concerning the status of the Psammolestes genus, some morphological differentiations with the Rhodnius genus has justified the erecting of two distinct genera as shorter head, stronger legs, wider femora and wider rostrum but also the lifestyle of the Psammolestes species in nests of birds while Rhodnius species live mainly in palm tree (Lent and Wygodzinsky, 1979). Our work clearly established using both large mitochondrial and nuclear datasets that P. tertius is included into the Rhodnius genus and the prolixus group. This position into the Rhodnius genus was also pointed out by previous studies (e.g., Lyman et al., 1999; Monteiro et al., 2000; Hypša et al., 2002; Justi et al., 2014; Paula et al., 2021). Thus, according to the cladistic principles and the rule of species name precedence, we strongly recommended reclassifying the P. tertius species as Rhodnius tertius. Complementary studies combining morphological and molecular approaches will be useful to know to what extent the morphological characters put forward in Psammolestes are part of some adaptive processes acting on Rhodnius species from differentiated environments. It is interesting to note that among the species of the Rhodnius genus, R. paraensis which presents a particular ecology because found in the nests of arboreal rodents, is a species of small size as Psammolestes with also short head, stout legs, short and stout antennae that led Lent and Wygodzinsky (1979) to suggest that these characters “make this one of the most plesiomorphic species of Rhodnius, if the progressive elongation of the head and appendages found in most other species is considered as derived.” Phylogenetic studies including R. paraensis and Psammolestes species are required to determine the polarity (ancestral vs. derived) of the character states. Since the three Psammolestes species are morphologically and ecologically close (Lent and Wygodzinsky, 1979), share cytogenetic characteristics indicating a chromosomal homogeneity (Oliveira et al., 2018), and are phylogenetically closely related (rDNA, Kieran et al., 2021), we recommended to re-examine also the taxonomic status of the two other Psammolestes species, P. arthuri and P. coreodes.

The species R. amazonicus was initially described by Almeida et al. (1973). It was then invalidated by Lent and Wygodzinsky (1979) because of the description made on only one female collected in the municipality of Manaus (Amazonas, Brazil) with several R. pictipes specimens and was then seen as a poorly preserved specimen of R. pictipes. The species was then revalidated by Bérenger and Pluot-Sigwalt (2002) from both male and female specimens collected in French Guyana and more recently recorded by Rosa et al. (2017b) in Pará, Brazil. We provided the first molecular data for this species and confirmed by both our robust mitochondrial and nuclear phylogenies its species-specific status as bona species. Moreover, this taxon belonged to the pictipes group at a basal position.

In contrast, the large molecular data set of our study definitively closes the question about the species-specific status of R. taquarussuensis (syn R. neglectus), described by Rosa et al. (2017a) from a female, which invaded a rural dwelling in the city of Taquarussu, Mato Grosso do Sul, Brazil. Because of its very close genetic similarity for both mitochondrial and nuclear genes, we demonstrated that this sample belongs undoubtedly to the R. neglectus species. This result corroborated the revision performed by Nascimento et al. (2019) using two mitochondrial and four nuclear genes leading the authors to state that the morphological and constitutive heterochromatin pattern differences between R. taquarussuensis and R. neglectus are likely intraspecific polymorphism of R. neglectus and consequently to synonymize R. taquarussuensis with R. neglectus.

In our dataset, we included R. milesi, a less studied species. This species was described from Pará, Brazil and morphologically similar to R. dalessandroi (Valente et al., 2001). Yet, for this latter species, only the published description and illustration are available leading Lent and Wygodzinsky (1979) to not include it in their taxonomical work. Doubts about the species-specific status of R. milesi had been issued by Monteiro et al. (2018) in their review. The authors concluded to the proximity of R. milesi with R. neglectus based on molecular cytb and ITS sequences but the accession number of the R. milesi sequences are not provided. In our dataset, for both the nu-PCG and mitochondrial phylogenies, MILE is congruently related to the R7WU R. nasutus sample with low nuclear genetic pairwise values (<0.002) but is clustered with R. neglectus in the rDNA phylogeny. A possible explanation is that the MILE specimen used in this study is a R. nasutus with a probable introgression with some nuclear DNA from R. neglectus. This molecular finding is in accordance with the morphological observation made by Carcavallo et al. (in Valente et al., 2001) observing that R. milesi male genitalia have a second phallosome process, only found in R. nasutus. With respect to our study, we conclude that R. milesi should be considered as synonym to R. nasutus and we recommend reconsidering the species-specific status of this species with a larger sampling.

A crucial point concerns two species closely related to R. robustus, namely the two newly described R. montenegrensis (Rosa et al., 2012) and R. marabaensis (Souza et al., 2016). Rhodnius montenegrensis was described from a strain established from eight specimens collected in a palm tree in the municipality of Monte Negro, Rondonia, Brazil and morphologically identical to specimens previously caught in the same locality (Rosa et al., 2012). The authors assessed their description on morphological and cytb sequence divergences with R. robustus. But, previously, using a large set of R. prolixus and R. robustus samples, Monteiro et al. (2003) based on cytb phylogeny have described four clades for R. robustus and proposed to see R. robustus as a complex of species and concluded for a paraphyly for R. robustus. Thus, the validity of the newly described species was pointed out by Abad-Franch et al. (2013) by declaring that “R. montenegrensis likely represents one of the R. robustus lineages of Monteiro et al. (2003).” Using cytb and ITS2 genes from transcriptomic data, Brito et al. (2019) evidenced that R. montenegrensis and R. robustus clade II are likely the same species but this does “not invalidate the former as a separate species.” Concerning R. marabaensis, this species was described from Marabá, Pará, Brazil (Souza et al., 2016). Unfortunately, no DNA sequence is available in public database for the described species, the authors only stated that the cytb sequence shown 99% of identity with that of R. robustus. In this study, based on the nuclear genes, we evidenced similar genetic distance values between the three pairwise involving Rob5s R. robustus from field French Guyana, R. montenegrensis (R8F), and ROBB. Similar distant pairwise values were also obtained between the two species R. nasutus and R. neglectus. Indeed, if we recognized R. nasustus and R. neglectus as valid species, the three entities must be too, especially since they displayed congruent position in all trees, mitochondrial and nuclear ones. Thus, we state that the sample Rob5s is representative of R. robustus, R8F of R. montenegrensis and the specimen ROBB has also a species-specific status. It is worth noting that ROBB was caught in a palm tree near the Benfica Field Station (5°16′S, 49°50′E) at 50 km from Marabá, Pará, Brazil. Since R. marabaensis was described from Marabá (Souza et al., 2016), even if molecular data for this are not publicly available, a parsimonious hypothesis could be formulated that ROBB is putatively R. marabaensis. The same assumption was made by Castro et al. (2020) assimilated the R. robustus clade III to R. marabaensis based on the geographical origin of the samples. For the four close related species R. marabaensis, R. montenegrensis, R. robustus, and R. prolixus, we proposed to use the terminology of “robustus species complex.”

In order to corroborate the inference about the R. robustus clades defined using cytb sequences (Monteiro et al., 2003), despite we demonstrated that this gene was under selective pressure and that Rhodnius phylogeny is not well resolved using this marker, we nonetheless performed an ML phylogeny including the R. robustus clade sequences and also cytb sequences from R. barretti (Figure 7). In the cytb phylogeny: (i) R. montenegrensis R8F was robustly included in the ROB II clade, (ii) ROBB in the ROB III clade assimilated to R. marabaensis and (iii) R. robustus Rob5s in the ROB IV clade.

Note that the new species R. barretti described from both Colombia and Ecuador and morphologically close to R. robustus (Abad-Franch et al., 2013), showed on the cytb phylogeny in Monteiro et al. (2018) a basal position to R. neglectus, R. nasutus, and the robustus species complex. In the cytb phylogeny (Figure 7), the NASG has a position similar to R. barretti. If only mitochondrial genes were analyzed, due to the high divergence observed, the proposal of a specific status for this specimen could arise. But the nuclear data are unequivocal, NasG is relative to R. nasutus but having a mitochondrial introgression likely by DNA from a species of the robustus complex (see below). If the position of R. barretti in the cytb is irrefutable, more molecular data and especially nuclear ones are needed to discard putative introgression and definitively assess its phylogenetic position.

If the mitochondrial genes used for phylogenetics suffer from disadvantages, discordance between mitochondrial and nuclear phylogenies could highlight evolutionary processes and be a potential source of additional information about the evolutionary history of a recent clade. Such mito-nuclear discordance is not rare and was reported in numerous taxa. Incomplete lineage sorting (ILS) could explain this discordance by the conservation of ancestral polymorphisms in multiple lineages after species splits (Degnan and Rosenberg, 2009). However, polymorphisms in mtDNA are expected to be lost relatively rapidly through genetic drift because of the reduced effective population size in mitochondrial genomes. This leads ILS to unlikely be the predominant cause of mito-nuclear discordance, notably when a large number of nuclear loci are congruent in phylogenies but in conflict with mitochondrial ones (Good et al., 2015). Another more likely cause of discordance could be introgression, highlighted between otherwise well-defined species by differentially affecting nuclear vs. mitochondrial DNA. In many cases, mitochondrial introgression occurs with few or no detectable movement of nuclear genes (Sloan et al., 2017). Nonetheless, evidence for mitochondrial introgression with accompanying nuclear genes usually selected against in the recipient species was described (Rubinoff and Holland, 2005; Matute et al., 2020). A least plausible explanation is hybrid speciation, though more common in nature than previously thought (Mavárez and Linares, 2008; Schumer et al., 2014; Payseur and Rieseberg, 2016). In this study, for the mitochondrial genes, misleading phylogenetic reconstructions could be due to sequence convergence as we found the occurrence of selective pressure on these genes, but also to introgression, with mitochondrial introgression being the most likely process rather than the nuclear one (Good et al., 2015).

In our study, we aim to test signatures of introgression using the ABBA-BABA test (Green et al., 2010) for the pictipes group for which a discordance between the mitochondrial and nuclear PCG trees was noticed. However, the mapping on R. prolixus, the only reference genome available for the Rhodnius species, could be a limitation for this study. Withal, the genome-wide analysis supports in the pictipes group the hypothesis of a mitochondrial introgression between R. brethesi and R. stali as the source of incongruence. Based on the nuclear phylogenies, we can suggest an introgression of R. stali with mitochondrial DNA from R. brethesi. This assumption is supported by the present distribution of the two species, R. stali being further south (Bolivia, Brazil: Mato Grosso and Acre) than R. brethesi (Colombia, Venezuela, Brazil: Pará and Amazonas) and also with the morphological similarity of the two species. However, an unsampled (or undiscovered) species of the pictipes group as the source of the introgressed mitogenome cannot be excluded. The recent discovery of a new species R. micki (Zhao et al., 2021) from Bolivia could be a putative candidate as the source of mitochondrial introgression since this species is morphologically similar to R. pictipes and R. stali but R. brethesi was not included in the morphometric study and molecular data from R. micki is still lacking. A selective sweep of mitochondrial DNA linked to an expending endosymbiont could also be considered. In this way, a study of Wolbachia genomic data (Filée et al., unpublished) showed that the same Wolbachia strain is shared by species of the pictipes group (R. brethesi, R. amazonicus, R. stali, and R. pictipes) but that it differs from that of the prolixus group (R. robustus, R. prolixus, R. neglectus, and R. nasutus). Indeed, this endosymbiont does not seem to be involved in this mechanism.

For the pallescens group, we tested signatures of introgression because if the mitochondrial and nu-PCG trees are congruent showing R. colombiensis and R. pallescens clustered together, in the rDNA phylogeny, R. pallescens and R. colombiensis are closer. Such divergent topologies were already obtained in the literature, indicating either a close relationship between R. pallescens and R. colombiensis (Hypša et al., 2002) or between R. ecuadoriensis and R. colombiensis (Justi et al., 2016). However, the first topology is coherent with the vector distribution since the Chocó rainforests along the Colombian Pacific coast is assumed to separate R. ecuadoriensis from R. pallescens–R. colombiensis (Abad-Franch et al., 2009). Although the tests were carried out using R. prolixus for the mapping, three hypotheses could explain the numerous SNPs supporting either each grouping (ABBA- or BABA distribution): Incomplete Lineage Sorting (ILS), hybrid speciation and introgression. ILS could be the result of a fast radiation of the three species that led to the sharing of ancestral polymorphisms. However, this hypothesis cannot explain alone the under-representation of the ABBA configuration, i.e., less allele sharing between R. ecuadoriensis and R. colombiensis than between R. pallescens and R. colombiensis. A least conventional yet plausible explanation is hybrid speciation. The ambiguous positioning of R. pallescens as closer to either ecuadoriensis or colombiensis could be the result of crosses between R. ecuadoriensis males and R. colombiensis females leading to the third species R. pallescens, which will reflect the mitochondrial tree topology. But hybrid speciation remains to be demonstrated in Rhodnius. Introgression, which results from the occurrence of gene flow after speciation, could explain the ABBA under-representation if R. ecuadoriensis and R. colombiensis were genetically isolated but each still hybridizing with R. pallescens. Indeed, Díaz et al. (2014) demonstrated in the laboratory the crossbreeding between R. pallescens and R. colombiensis although the laboratory hybrids were infertile. We did not observe ABBA or BABA configurations restricted to specific genomic regions, the three configurations were widespread and overlapping on all scaffolds of the reference genome of R. prolixus used for the read mapping. Nonetheless, we cannot rule out the possibility that multiple rearrangements have taken place since the divergence of the three species and R. prolixus. At least, nuclear R. ecuadoriensis introgression in R. pallescens could explain the ABBA under-representation since the nuclear rDNA showed a closest genetic proximity between these two species. Therefore, genomic studies with large population sampling in the R. pallescens species group wherein different demographic models could be tested are strongly needed in the future.

In this study, the robust nu-PCG phylogeny splits the prolixus group into two major clades, the neglectus-nasutus clade and the robustus species complex clade (R. marabaensis, R. montenegrensis, R. robustus, and R. prolixus). In these two clades, we observed several discrepancies between mitochondrial and nuclear trees.

The wild Piauí-Brazil NasG specimen showed an incongruent position between the mitochondrial and nuclear PCG trees. This specimen showed a basal position for the two clades in the mitochondrial tree but it is clustered with the R. neglectus samples in the nu-rDNA tree and with R7WU R. nasutus in the nu-PCGs tree. Considering both its morphology is similar to R. nasutus and the substantial data for the nu-PCG dataset (36.3 kb), the mitochondrial introgression is the most likely process rather than the nuclear one. The nature of the mitochondrial introgression remains speculative. Unexpectedly, an introgression with a member of the robustus complex species cannot be ruled out. Moreover, some nuclear introgression for this R. nasutus sample could be hypothesized probably by rDNA from R. neglectus.

Mitochondrial introgression was also observed for two other field specimens from Piauí. The INCP sample exhibited a dark morph and an atypical phenotype which cannot be assigned to R. nasutus or R. neglectus. The molecular study evidenced that this sample is related to R. neglectus in the nuclear tree but to R. nasutus in the mitochondrial one, thus INCP could be considered as a R. neglectus introgressed by mitochondrial DNA from R. nasutus. The reverse introgression is exhibited by the NASP sample, that is a R. nasutus introgressed by mitochondrial DNA from R. neglectus. Introgressions between R. nasutus and R. neglectus was not yet reported, including from field specimens but are not surprising due to the close genetic proximity between the two species as revealed by genetic distance values but also to the geographic distribution of the species. Results obtained by Abad-Franch et al. (2009) via geometric morphometrics showed that R. neglectus and R. nasutus are sympatric in the Brazilian Cerrado-Caatinga transitional area. There is also some evidence of populations of R. neglectus occurring also in the true Caatinga biome where R. nasutus is distributed (Dias et al., 2011), which could favor hybridization. Moreover, chromatic variation in R. nasutus was already observed and adaptive phenotypic plasticity was suggested related to microhabitat features (Abad-Franch et al., 2009). It will be of interest using molecular markers to test the various R. nasutus chromatic forms for introgression with R. neglectus, that could explain the darker morph observed.

For the robustus species complex, three field specimens NEGP, RobQ and MILEP, all originated from Pará Brazil, and morphologically identified as R. neglectus, R. robustus and R. milesi respectively, a misidentification is likely since all the trees (mtG, rDNA, nu-PCGs) are congruent with R. prolixus species. This exemplifies the difficulty to identify the species of the prolixus group and to admit the occurrence of sylvatic R. prolixus. If the occurrence of R. prolixus in Venezuela is now well documented (i.e., Fitzpatrick et al., 2008), in Brazil the controversy persists but it is actually accepted that R. prolixus occurs basically in the states of Amazonas and Pará (Pinho et al., 1998). Indeed, the discovery of R. neglectus in this region (Lent, 1954) but also the widely distribution of R. robustus and the description of R. milesi from Pará cast doubt on past identifications of R. prolixus and guided the future ones. By instance, it was suggested that R. neglectus, and not R. prolixus, was the species invading houses in central Brazil (Gurgel-Gonçalves et al., 2008). In our study, we then evidence by molecular data the presence of sylvatic R. prolixus in the Brazilian Pará state.

For two strain samples, R4H and RobR, a mito-nuclear conflict was evidenced. In the nuclear tree, these two samples were clustered with R. prolixus samples but in the mitochondrial tree, they were closely related to R. montenegrensis and likely have experienced mitochondrial introgression with this species. In addition, for the RobR sample, the rDNA tree suggested an introgression by R. robustus. The domiciliary Trujillo-Venezuelan R. prolixus specimen (ProYRP) was also introgressed by R. robustus but for mitochondrial DNA. Similar results were obtained for Venezuelan populations by Fitzpatrick et al. (2008) using both the mitochondrial cytb gene and the D2 nuclear gene, a population of domiciliated R. prolixus identified as R. robustus was introgressed with mitochondrial DNA from R. robustus (Hap3, Figure 7). Introgression events of R. prolixus with R. robustus in Venezuela could explain the results obtained by Harry (1994) with morphometric data in which no clear-cut differentiation was evidenced between Venezuelan R. prolixus and R. robustus samples from Trujillo.

Considering the cytb phylogeny including Monteiro et al. (2003) samples, the apparent paraphyly of R. robustus pointed out by the authors is resolved by the consideration of introgression events (Figure 7):

In this study, we demonstrated that introgression could be a biological event as field samples were introgressed. Misplacement in phylogeny for some Rhodnius lab strains could also result from a past biological event and not only be the fact of mixed or misidentified samples as pointed out by Brito et al. (2019). Indeed, at this time, it is impossible to state if lab strains with such pattern were accidentally introgressed or were founded with specimens naturally introgressed.