Phylogenomics for Chagas Disease Vectors of the Rhodnius Genus (Hemiptera, Triatominae): What We Learn From Mito-Nuclear Conflicts and Recommendations

We provide in this study a very large DNA dataset on Rhodnius species including 36 samples representing 16 valid species of the three Rhodnius groups, pictipes, prolixus and pallescens. Samples were sequenced at low-depth with whole-genome shotgun sequencing (Illumina technology). Using phylogenomics including 15 mitochondrial genes (13.3 kb), partial nuclear rDNA (5.2 kb) and 51 nuclear protein-coding genes (36.3 kb), we resolve sticking points in the Rhodnius phylogeny. At the species level, we confirmed the species-specific status of R. montenegrensis and R. marabaensis and we agree with the synonymy of R. taquarussuensis with R. neglectus. We also invite to revisit the species-specific status of R. milesi that is more likely R. nasutus. We proposed to define a robustus species complex that comprises the four close relative species: R. marabaensis, R. montenegrensis, R. prolixus and R. robustus. As Psammolestes tertius was included in the Rhodnius clade, we strongly recommend reclassifying this species as R. tertius. At the Rhodnius group level, molecular data consistently supports the clustering of the pictipes and pallescens groups, more related to each other than they are to the prolixus group. Moreover, comparing mitochondrial and nuclear tree topologies, our results demonstrated that various introgression events occurred in all the three Rhodnius groups, in laboratory strains but also in wild specimens. We demonstrated that introgressions occurred frequently in the prolixus group, involving the related species of the robustus complex but also the pairwise R. nasutus and R. neglectus. A genome wide analysis highlighted an introgression event in the pictipes group between R. stali and R. brethesi and suggested a complex gene flow between the three species of the pallescens group, R. colombiensis, R. pallescens and R. ecuadoriensis. The molecular data supports also a sylvatic distribution of R. prolixus in Brazil (Pará state) and the monophyly of R. robustus. As we detected extensive introgression events and selective pressure on mitochondrial genes, we strongly recommend performing separate mitochondrial and nuclear phylogenies and to take advantages of mito-nuclear conflicts in order to have a comprehensive evolutionary vision of this genus.

We provide in this study a very large DNA dataset on Rhodnius species including 36 samples representing 16 valid species of the three Rhodnius groups, pictipes, prolixus and pallescens. Samples were sequenced at low-depth with whole-genome shotgun sequencing (Illumina technology). Using phylogenomics including 15 mitochondrial genes (13.3 kb), partial nuclear rDNA (5.2 kb) and 51 nuclear protein-coding genes (36.3 kb), we resolve sticking points in the Rhodnius phylogeny. At the species level, we confirmed the species-specific status of R. montenegrensis and R. marabaensis and we agree with the synonymy of R. taquarussuensis with R. neglectus. We also invite to revisit the species-specific status of R. milesi that is more likely R. nasutus. We proposed to define a robustus species complex that comprises the four close relative species: R. marabaensis, R. montenegrensis, R. prolixus and R. robustus. As Psammolestes tertius was included in the Rhodnius clade, we strongly recommend reclassifying this species as R. tertius. At the Rhodnius group level, molecular data consistently supports the clustering of the pictipes and pallescens groups, more related to each other than they are to the prolixus group. Moreover, comparing mitochondrial and nuclear tree topologies, our results demonstrated that various introgression events occurred in all the three Rhodnius groups, in laboratory strains but also in wild specimens. We demonstrated that introgressions occurred frequently in the prolixus group, involving the related species of the robustus complex but also the pairwise R. nasutus and R. neglectus. A genome wide analysis highlighted an introgression event in the pictipes group between R. stali and R. brethesi and suggested a complex gene flow between the three species of the pallescens group, R. colombiensis, R. pallescens and R. ecuadoriensis. The molecular data supports also a sylvatic distribution

INTRODUCTION
The blood-sucking bugs (Hemiptera, Reduviidae, Triatominae) are vectors of the parasite Trypanosoma cruzi (Chagas, 1909) (Kinetoplastea, Trypanosomatida). Most species of live in natural habitats where they feed on a large variety of vertebrates. But some of them can be found in human dwellings responsible for the Chagas disease endemic to Latin America where about 6 million people are infected (PAHO, 2020). The Triatominae subfamily is distributed in five tribes, the most diverse among them being the Rhodniini and Triatomini, which represent more than 90% of the known species. The Rhodniini tribe is composed of two Rhodnius and Psammolestes genera, which show some morphological differentiations (Lent and Wygodzinsky, 1979) that has justified the erecting of two different genera. They also differ in their ecology, the Rhodnius species mainly associated with palm trees and the Psammolestes species with bird nests (Galvão and Justi, 2015). Whereas only three Psammolestes species were described, there are twenty-two Rhodnius species but one species, R. taquarussuensis (Rosa et al., 2017a), has been recently synonymized with R. neglectus (Nascimento et al., 2019) (see Table 1 with authors and year of description). The Rhodnius genus is divided into three major groups, pictipes, prolixus and pallescens (Hernández et al., 2020). The three species of the pallescens group (R. colombiensis, R. ecuadoriensis, and R. pallescens) are trans-Andean, distributed in the west of the Andes. The seven species of the pictipes group (R. amazonicus, R. brethesi, R. micki, R. paraensis, R. pictipes, R. stali, and R. zeledoni) are cis-Andean, distributed in the east of the Andes and in the Amazon region. The same distribution is observed for ten species of the prolixus groups (R. barretti, R. dalessandroi, R. domesticus, R. marabaensis, R. milesi, R. montenegrensis, R. nasutus, R. neglectus, R. prolixus, and R. robustus), the eleventh R. neivai having also some trans-Andean populations (Zhao et al., 2021).
Some species are difficult to differentiate from each other on morphological criteria thus generating taxonomic conflicts. Pioneer studies using allozymes had questioned the speciesspecific status of closely related species such as Venezuelan populations of R. robustus and R. prolixus (Harry et al., 1992;Harry, 1993). Concerning some newly described species, molecular data are sparse for R. barretti (Abad-Franch et al., 2013), not released in a database for R. marabaensis (Souza et al., 2016) or lacking for R. micki (Zhao et al., 2021). Recently, molecular studies on a newly described species R. montenegrensis (Rosa et al., 2012) suggested that it is not a true species but a part of the R. robustus variability (Abad-Franch et al., 2013;Monteiro et al., 2018). However, the two species R. montenegrensis and R. marabaensis were recognized by Castro et al. (2020) based on the geographical origin of their samples. Moreover, some authors consider R. milesi to be a R. neglectus variant from southeastern Amazonia (Abad-Franch et al., 2013). Otherwise, only a few specimens have been collected for R. paraensis first described as R. domesticus (Sherlock et al., 1977) and R. amazonicus close to R. pictipes but this latter species was revalidated by Bérenger and Pluot-Sigwalt (2002). Finally, the species-specific status for two species is questionable and difficult to deepen since R. zeledoni reported as morphologically close to R. paraensis was described from only one dead and dried male specimen (Jurberg et al., 2009), and R. dalessandroi is actually known only from its published description (Carcavallo and Barreto, 1976).
Some Rhodnius species are important vectors of T. cruzi, as R. prolixus, the main Chagas disease vector within the genus, which extends from Central America to the Andean countries and the Amazon basin. Control programs in Central America were thought to have achieved the elimination of this vector but it was reported again in Mexico in 2019 (Antonio-Campos et al., 2019). Three other species are domiciliated in some countries, namely R. ecuadoriensis in the northern zone of Peru and Ecuador, R. stali in Bolivia, and R. pallescens in Panama. In Brazil, R. neglectus and R. nasutus often colonize human environments (Galvão and Justi, 2015).
Phylogenetic studies performed on the Rhodniini tribe, mainly with mitochondrial genes, showed that Psammolestes was nested within the Rhodnius clade but discrepancies remained about its phylogenetic position within the prolixus group (reviewed in Hernández et al., 2020). Moreover, several molecular studies have also focused on the relationship between the three groups. Two tree topologies were obtained either pictipes and pallescens, or pictipes and prolixus groups as sister taxa, unrelated to the type of marker used, mitochondrial, nuclear, or both (reviewed in Hernández et al., 2020). Recently, the third configuration, prolixus and pallescens groups as sister taxa, was also reported by Kieran et al. (2021) using ultraconserved elements.
The principal sources of discrepancies for phylogenies may come from methodological pitfalls and/or reflect the biological reality. Methodological bias could therefore be induced by incomplete/different taxon samplings but also by the type and the number of markers used. However, discordance between mitochondrial and nuclear phylogenies are expected because the mitochondrial genome is haploid uniparentally inherited in most animals unlike the nuclear one and has faster rates of evolution. Incongruences between single gene phylogenies can be due to incomplete lineage sorting, deep coalescence, horizontal gene transfer, introgression, hybridization, hidden paralogy or lack of phylogenetic information (Toews and Brelsford, 2012;Campillo et al., 2019). For the Rhodniini tribe, only the complete mitogenome of R. pictipes (Zhao et al., 2019) and the complete nuclear genome of R. prolixus (Mesquita et al., 2015) are available. For all the phylogenies yet performed, the datasets were indeed limited in terms of the number of genes. So far, the most complete studies for the Rhodniini tribe were performed by Justi et al. (2014) using ten species and 1-4 genes depending on the species, Kieran et al. (2021) using from eight to sixteen species and ultraconserved elements and/or rDNA data, and Paula et al. (2021) using thirteen species and three genes.
In this study, we aimed to resolve major phylogenetic conflicts within the Rhodnius genus using phylogenomics from lowdepth whole-genome shotgun sequencing and to explore mitonuclear conflicts in order to have a better understanding of the evolution of this genus.
Data were obtained by Illumina technology for 36 samples of the Rhodniini tribe including 17 putative species identified using morphological characters. We performed both mitochondrial and nuclear phylogenies using a set 15 mitochondrial genes (13 protein-coding mitochondrial plus 2 rDNA genes) and two sets of nuclear markers (nu-rDNA genes and 51 protein-coding genes). The protein-coding datasets were tested for selection.
The outcome of this work led us to identified 16 valid Rhodnius species in our molecular dataset and to formulate recommendations for both taxonomic and phylogenetic issues.

Material and Genomic Data
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.

Mitogenomes and Nuclear rDNA Operon Assembly and Annotation
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.
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.

Genome Assembly and Protein-Coding Gene Annotation
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

Diversity Parameters
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.

Phylogenetic Analysis
All the individual mitochondrial genes and the three concatenated datasets namely the 15 mtG, the nu-rDNA, and the For all the samples used are given: the quantity in gigabase (Gb) of the data generating by the whole-genome shotgun sequencing by Illumina technology (Data), the corresponding genome depth (x), given that 0.7 Mb (R. prolixus genome reference) was used for the calculation, the length sequences in base pair (bp) for the mtG, the nu-rDNA and the nu-PCGs datasets and the corresponding sequence coverage in percentage. See Table 1 for the Rhodniini samples description. The Triatominae outgroups used were: T. brasiliensis (Tbra), T. dimidiata (Tdim), P. rufotuberculatus (Pruf), and B. colossus (Bcol).
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.

Detection of Selection Signatures
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.

Introgression Analysis
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 2 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 3 . 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.

Mitochondrial Phylogenetic Inferences
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 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 FIGURE 1 | Phylogenetic inferences for the concatenated mitochondrial (mtG) dataset for Rhodnius species with the Maximum Likelihood (ML) method. Statistical supports (1,000 replicates SH-aLRT) are indicated on the nodes. T. dimidiata and P. rufotuberculatus are used as outgroups. Introgressed species inferred from the nuclear phylogeny (see Figure 5) are indicated in purple, misidentified species in red, and species with incongruent position are highlighted in dark purple. 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 FIGURE 2 | Phylogenetic inferences for the concatenated mitochondrial (mtG) dataset for Rhodnius species using Bayesian Inference (BI) method. Statistical supports (specify Bayesian posterior probabilities in percent) are indicated on the branches. T. dimidiata and P. rufotuberculatus are used as outgroups. Introgressed species inferred from the nuclear phylogeny (see Figure 6) are indicated in purple, misidentified species in red, and species with incongruent position are highlighted in dark purple.

Nuclear rDNA Phylogenetic Inference
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.

Nuclear Protein-Coding Gene Phylogenetic Inferences
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 FIGURE 3 | Phylogenetic inferences for the concatenated nuclear rDNA (nu-rDNA) dataset for Rhodnius species with the Maximum Likelihood (ML) method. Statistical supports (1,000 replicates SH-aLRT) are indicated on the nodes. B. colossus was used as outgroup. Introgressed species inferred from the nuclear phylogeny (see Figure 5) are indicated in purple, and misidentified species in red. 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.

Genome-Wide Examination for pallescens and pictipes 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.

Phylogenetic Relationship Between the Three Rhodnius Groups
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(Paula et al., , 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.

Variation Rate Among Partitions
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.

Selective Pressure on Mitochondrial Genes
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.

Species Diversity in the Rhodnius Genus
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.

Discordance Between Mitochondrial and Nuclear Phylogenomic Inferences
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).

The pictipes Group
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, FIGURE 7 | Phylogenetic inferences and interpretation for the Cytb gene for Rhodnius species with the Maximum Likelihood (ML) method including Monteiro et al. (2003) samples. The statistical supports (1,000 replicates SH-aLRT) are indicated on the nodes. T. dimidiata was used as outgroup. For the dataset of this study, introgressed species inferred from the nuclear phylogeny (see Figure 5) are indicated in purple, misidentified species in red. For Monteiro et al. (2003) samples, based on the inferences made for our samples, we proposed putative introgressed species, highlighted in the tree and putative misidentifications for some others, indicated in red. For ROB I clade clustered with R. prolixus samples and composed of Venezuelan R. robustus samples, they putatively are R. robustus introgressed by mitochondrial DNA from R. prolixus. For ROB II clade, some samples clustered with the R. montenegrensis could be synonymized with R. montenegrensis while the other could be R. prolixus introgressed by mitochondrial DNA from R. montenegrensis. For the ROB III clade, specimens clustered with R. marabaensis and could be synonymize as R. marabaensis. For the ROB IV clade, some samples are clustered with R. robustus while others could be R. prolixus introgressed by mitochondrial DNA from R. robustus.
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.

The pallescens Group
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 underrepresentation 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.

The prolixus Group
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): (iii) specimens of the ROB III clade clustered with ROBB from Marabá and assimilated to R. marabaensis. So, we conclude that this clade could be synonymize as R. marabaensis has used by Castro et al. (2020). (iv) ROB IV clade is also split into two. One comprises samples including R. robustus from French Guyana and clustered with also our sample from this region. The same clustering was found by Barnabé et al. (2018) for an extensive survey of R. robustus from French Guyana. The other comprises samples clustered with our Venezuelan sample ProYRP for which we demonstrated the introgression of R. prolixus by mitochondrial DNA from R. robustus and the Hap3 also a R. prolixus introgressed by mitochondrial DNA from R. robustus and sharing the same haplotype than the roBR9 (Fitzpatrick et al., 2008). Then, we suggested that this part of samples from ROB IV clade are R. prolixus introgressed by mitochondrial DNA from R. robustus.
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.

CONCLUSION
In this study, we provide a very large dataset on Rhodnius species including 36 samples. With low-depth whole-genome shotgun sequencing, we resolve sticking points in the Rhodnius phylogeny using phylogenomic approaches based on 15 mtG (13.3 kb), nu-rDNA genes (5.2 kb), and 51 nu-PCGs (36.5 kb). Out the 17 putative species determined on phenotype, using molecular data, we identified 16 valid Rhodnius species (Table 1). We confirmed the species-specific status of R. montenegrensis and R. marabaensis and we agree with the synonymy of R. taquarussuensis with R. neglectus. We also invite to revisit the species-specific status of R. milesi that is more likely R. nasutus. We added R. marabaensis in our dataset, first identified as R. robustus. We propose the robustus complex species for including the four close related species R. marabaensis, R. montenegrensis, R. prolixus, and R. robustus. Moreover, we confirm the position of P. tertius into the Rhodnius genus and strongly recommended reclassifying P. tertius as R. tertius. For the newly described species R. micki and R. barretti, molecular data and especially nuclear ones are needed to assess their phylogenetic position. The same is true for R. paraensis rarely sampled since this species has been collected only from arboreal rodent nests in the Amazon region in Brazil (Sherlock et al., 1977) but also with light trap in French Guyana (Bérenger et al., 2009). For R. zeledoni and R. dalessandroi only known from their description, no molecular data can be obtained.
Both mitochondrial and nuclear data coherently support that the pictipes and pallescens groups are more related to each other than they are to the prolixus group. We can hypothesize that the pictipes and pallescens groups more likely diverged from an ancestral form when the Northern Andean uplift before the formation of the Pebas system around 23-10 MYA and that the latter event has induced the diversification of the prolixus group when this system started to be drained eastward into the Atlantic Ocean. This subsequently led to hydrological changes promoting allopatric speciation or favoring dispersals (Albert et al., 2018). However, robust time-calibrated phylogenies are now needed to infer the historical and spatial origin of the Rhodnius lineages.
By comparing the topologies obtained for mitochondrial and nuclear phylogenies but also by performing genome-wide analysis, our results demonstrated unexpected introgression events in all the different Rhodnius groups, in laboratory strains but also in wild specimens. In the pictipes group, we formulated the hypothesis of cytoplasmic introgression between R. brethesi and R. stali and in the pallescens group the occurrence of complex introgression. In the prolixus group introgression occurred between various species: (i) from R. robustus or R. montenegrensis mitochondrial DNA to R. prolixus,(ii) from R. robustus mitochondrial DNA to R. prolixus, (iii) from R. nasutus mitochondrial DNA to R. neglectus, (iv) the reverse introgression from R. neglectus mitochondrial DNA from to R. nasutus and (v) putatively from one species of the robustus complex to R. nasutus. The taking into account the introgression events makes it possible to explain the paraphyly observed for R. robustus in some mitochondrial phylogenies but leads to abolish it by considering nuclear genes and restores the monophyly for this species. The extensive introgression especially in the prolixus group does not question the species boundaries for the pairwise R. nasutus and R. neglectus nor for the species of the robustus complex, and even that complicates the identification of the species, it is challenging from an evolutionary point of view. Further studies are needed in order to date the introgressions found in field samples which can result from old events at the time of the speciation or from recent hybridization between sympatric species.
This study exemplifies that the evolutionary inferences from only mitochondrial markers notably by using only the cytb gene or incomplete dataset could be misleading. Because of the persuasive introgression which confuse not only the molecular issue but surely also the phenotypic one, we strongly recommend the use of both mitochondrial and nuclear markers prior any study on Rhodnius species and to perform separately mitochondrial and nuclear phylogenies in order to take advantages from mito-nuclear conflicts for having a comprehensive evolutionary vision.

DATA AVAILABILITY STATEMENT
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: NCBI (accession: PRJNA429761).

AUTHOR CONTRIBUTIONS
JF, MM, and HB: data analysis and writing. FM and EF-R: study design. CA: writing. J-MB: field collections and Rhodnius taxonomy. MH: study conception, DNA extractions, data analysis, and writing. All authors edited the manuscript and approved the submitted version.