Historical Biogeography and the Evolution of Hematophagy in Rhodniini (Heteroptera: Reduviidae: Triatominae)

The Rhodniini tribe is one of the five tribes in the subfamily Triatominae and is notorious for its domestic blood-sucking pests and vectors of Trypanosoma cruzi across Latin America. The human and economic costs of the Chagas disease in the American tropics are considerable, and these insects are of unquestionable importance to humans. We used mitochondrial rDNA (16S), nuclear ribosomal RNA (28S) and wingless (Wg) sequences to perform phylogenetic analysis to derive trees based on parsimony and maximum likelihood. Nucleotide sequences were used in molecular-clock analyses to estimate time divergence between species of Rhodniini. The potential distribution of each species was modeled and compared with Kappa statistic. Multivariate niches with bioclimatic variables were used to describe differences between the species using discriminant analysis. The results of this study indicate that the Rhodniini originated 17.91 Mya ago. Rhodnius domesticus is the oldest species having its origin at 9.13 Mya. Rhodniini are closely related to Salyavatinae that are specialist termite predators and diverged from this subfamily 30.43 Mya. Most species are clearly allopatric and have distinct bioclimatic niches. The colonization of bromeliads, palms trees and bird nests represent important events for the speciation of these taxa. The hematophagous habit can be described as a scenario where Rhodniini’s ancestor could be pre-adapted for the invasion of bromeliads, palm trees, and bird nests where they would find significant water availability and thermal damping. These environments are widely used by vertebrate inquilines that would be the source of food for the species of Rhodniini. Lastly, our results show an alternative position of Psammolestes in the phylogenetic tree.


INTRODUCTION
The Rhodniini tribe is one of the five tribes in the subfamily Triatominae and is notorious for its domestic blood-sucking pests and vectors of Trypanosoma cruzi (Chagas, 1909) (Kinetoplastida, Trypanosomatidae) across Latin America. The tribe has a wide geographical distribution ranging from Central America to the Southern Cone (Hernández et al., 2020). The human and economic cost of the Chagas disease in the American tropics is considerable, totaling the annual burden of $627.46 million in health-care costs (Lee et al., 2013). There are about six million chronic disease carriers in the Americas and these insects are of unquestionable importance to humans (PAHO/WHO, 2020). Because of this, Rhodniini has been extensively studied due to its epidemiological importance (Schofield and Galvão, 2009). The tribe is considerably diverse, with 24 described species (21 of the genus Rhodnius Stål and three of Psammolestes Bergroth) (Zhao et al., 2021). The members of the genus Rhodnius have been classified into three groups-pallescens (trans-Andean), pictipes (cis-Andean), and prolixus-based on geographical distribution and morphology (Schofield and Galvão, 2009;Hernández et al., 2020). Although the taxonomy of Rhodnius has been widely studied, there are controversies regarding the number of species, the classification of these species into groups, and the phylogenetic relationships and monophyletic status of these groups that deserve further studies to elucidate Rhodnius taxonomic status. This is the case of R. zeledoni Jurberg et al., 2009, which seems to be very similar to R. domesticus Neiva and Pinto, 1923(Galvão, 2014. The unique specimen was found very damaged in Sergipe state, Brazil, a region included in the distribution range of R. domesticus, therefore, examination of further material is essential to confirm whether or not R. zeledoni is a valid species. Following Monteiro et al. (2018), R. marabaensis dos Santos Souza et al., 2016 has been suggested to be a genetic variant very close to R. robustus Larrousse, 1927 (specific status recently considered valid by Castro et al., 2020) and R. taquarussuensis (Rosa et al., 2017), was synonymized with R. neglectus Lent, 1954by Nascimento et al. (2019 based on interspecific crosses and molecular markers. Issues in the phylogenetic relationships in Rhodnius have arisen from the conflicting results of morphometric analyses and molecular markers (Schofield and Galvão, 2009;Monteiro et al., 2018). Monteiro et al. (2003) and Pavan et al. (2013) proposed the paraphyly of R. robustus with respect to R. prolixus Stål, 1859, despite not including other populations of the species of the prolixus complex-R. nasutus Stål, 1859 and R. neglectus. Although rarely studied and recognized, the cross-contamination between R. prolixus and R. robustus colonies is fairly frequent. Mesquita et al. (2015) and Monteiro et al. (2018) proposed that the colony material should always be genotyped and compared with the established taxonomic standards, particularly when the research involves taxa with known cryptic variation. Misidentification of cryptic species must also be considered.
Triatominae species are almost entirely hematophagous, although there are cases of some species that feed on other invertebrates (Sandoval et al., 2004(Sandoval et al., , 2010. Through evolutionary processes, Triatominae species have specialized morphological, physiological, and behavioral adaptations for feeding on blood. Among them are mouthparts, specifically saliva composition and enzymes, along with digestive symbionts that are the most noticeable (Otálora-Luna et al., 2015). The ecological, morphological and molecular changes that are involved in the transition from a predatory lifestyle to hematophagous Triatominae are still unclear due to the lack of robust and densely sampled phylogenies, preventing the development of significant hypothesis tests on these issues (Monteiro et al., 2018). The transition from predatory habits in Reduviidae to hematophagous habits in Triatominae could have occurred by colonization of vertebrate nests. Furthermore, the evolution of blood feeding in the tribes Triatomini and Rhodniini + Cavernicolini may have occurred independently (Hwang and Weirauch, 2012).
de Paula et al. (2005) presented the first molecular evidence for the polyphyletic origin of the subfamily Triatominae and showed a separation between Rhodniini and Triatomini tribes from predatory groups. In that study, it is suggested that the Reduviinae subfamily is a sister group of Triatomini, and the subfamilies Salyavatinae and Harpactorinae are sister groups of Rhodniini. According to Hwang and Weirauch (2012), Triatominae are paraphyletic with respect to the Opisthacidius Berg, 1879 (Reduviinae), and this contrasts prior hypotheses that found Triatominae to be monophyletic or polyphyletic. In addition, the same authors, as well as Patterson and Gaunt (2010) and Justi et al. (2016) proposed Zelurus Hahn, 1826 and Opisthacidius as sister groups to Triatominae in their phylogenetic analysis. Despite that, several issues remain underlying on the hypotheses of the origin of Rhodniini. The most common involves the period in which one group has arisen and the geologic and biogeographic events that have conditioned it. Understanding the mechanisms and processes underlying patterns of species distributions remains a central question in historical biogeography (Lomolino et al., 2016).
Molecular dating studies involving Triatominae species included Gaunt and Miles (2002); Hwang andWeirauch (2012), andJusti et al. (2016). Until the work of Justi et al. (2016), the Rhodniini tribe had been analyzed in biogeography studies by de Paula et al. (2007) and Abad-Franch et al. (2009). Other biogeographical hypotheses for Rhodniini's evolution have also been tested (Hypša et al., 2002;Hwang and Weirauch, 2012). Justi et al. (2016) presented the first comprehensive phylogenetic study aiming to understand the geological events that led to Triatominae diversification using fossils as calibration points for time estimates and correlating the geological changes in the Neotropics to Triatominae evolution. They showed that the major geological changes in the Neotropics since the Eocene have played a role in the diversification of Triatominae, including the effects of the Andean uplift in the Amazonian area, comprising the formation of the Pebas system (23-10 Mya); the formation of the Acre system (10-7 Mya); the formation of the connection between the Amazon and the Atlantic Forest; the GAARlandia (GAAR = Greater Antilles + Aves Ridge) land bridge; the period of biodiversity exchange resulting from the closure of the Isthmus of Panama (10-2.7 Mya); the Florida and Gulf Coast inundations during the Miocene; and dispersal across the Bering Land Bridge during the Eocene.
The goal of the study reported here was to expand on our analysis (de Paula et al., 2007) by providing an update on the biogeographic hypotheses to explain the modern geographic distribution of Rhodniini species. We also performed molecular dating analysis for each species and presented an evolutionary hypothesis for the evolution of hematophagy. A fundamental question that we propose in this work is to evaluate the relationship between the Rhodniini tribe and the subfamilies Reduviinae and Salyavatinae. Salyavatinae are specialist termite predators (Cobben, 1978;van Doesburg and Forero, 2012) while Reduviinae are generalist predators (Hwang and Weirauch, 2012). Consequently, there are two possible scenarios to explain the haematophagic habit in Rhodniini that need to be elucidated. Lastly, we used bioclimatic data to show potential areas of distribution of each species.

Sequence Alignment, Phylogenetic Analysis, and Molecular Dating
The study makes use of mitochondrial rDNA sequences (16S), nuclear ribosomal RNA (28S) and wingless (Wg) sequences currently available in GenBank. We chose these markers because they are available for the outgroup and for most of the Rhodniini species studied.
Nucleotide sequences have been the most common form of data used in molecular clock analyses. Several models using molecular information have been developed and implemented to infer divergence times, although they do not assume a strict molecular clock, and yet are commonly applied to empirical data. The divergence time estimate analysis was conducted using BEAST 2.6.2 (Bouckaert et al., 2019), which goal is to estimate the joint posterior probability of the branch rates and times given a set of sequences and calibration information. The gamma model 1 https://www.ebi.ac.uk/Tools/msa/muscle/ was set to general time reversible (GTR), and base frequencies and lognormal relaxed clock (Uncorrelated) were estimated. We adjusted the scale parameter to reduce the prior density (α = 2 and β = 0.25). The ucldMean and ucldStdev were set to exponential. Age estimates from fossil organisms are the most common form of divergence time calibration information. Fossils are part of the same diversification process (i.e., birth-death model) that generated the extant species and these data were used as age constraints on their putative ancestral nodes (Ho and Phillips, 2009;Drummond, 2010, 2012;Heath, 2012). Recent analyses using relaxed clock models calibrated with fossils converge estimates on the origin of Triatominae around 35 Mya to just over 40 Mya (Hwang and Weirauch, 2012;Ibarra-Cerdeña et al., 2014;Justi et al., 2016). We calibrated the origin of Salyavatinae to 42.0 Mya (Hwang and Weirauch, 2012) and lognormal prior distribution (M = 1.0 and S = 1.25). By default, BEAST sets the number of generations to 10,000,000 and length of chain to 1,000,000, thus these parameters were applied. The first 10% of the samples in each tree file were discarded using TREEANNOTATOR 2.6.2 (Bouckaert et al., 2019). The remaining trees were used to produce the maximum clade credibility tree that was visualized in FIGTREE 1.3.1 (Rambaut, 2018).

Ecological Niche Model
The potential distribution of each species was modeled using the MAXENT 3.4.1 (Phillips et al., 2018). The distribution data for Rhodniini was obtained from Ceccarelli et al. (2018) (see Supplementary Table 2). Georeferenced data were checked for inconsistencies using DIVA-GIS 7.5 (Hijmans et al., 2012), and duplicated data were removed. WorldClim bioclimatic variables version 2 (Fick and Hijmans, 2017) were used to perform the analyses. The rasters had a 2.5 min resolution, which represents the value of the 19 bioclimatic variables in the study area (see Supplementary Table 3). We used 10 replicate analyses based on cross validation. The distributional data for each species were separated into two sets: one for model calibration (70% of points) and one for model evaluation (30% of points) (Elith et al., 2010;Merow et al., 2013;de Paula and Barreto, 2020). The cloglog format was used to output the results in ASCII format, and the average output grids were imported into DIVA-GIS 7.5 (Hijmans et al., 2012) for assessments and analyses. To define the potential distribution for Rhodniini species, we used a 10% training presence threshold (the probability value at which 90% of the presence points fall within the potential area) to generate binary rasters with potential distribution areas for each individual species (Phillips et al., 2006Phillips, 2017). We used the Area Under Curve (AUC) to evaluate models' performances. AUC values range from 0 (a model worse than random) to 1 (perfect discrimination) with a cut-off of 0.5, meaning the model does not have a better fit than random prediction models.
We compared potential distribution maps of related species with Kappa statistic tool in the MAP COMPARISON KIT (MCK) 3.2.1 (Visser and de Nijs, 2006; Research Institute for Knowledge Systems BV, 2010). The Kappa comparison method is based on a straightforward cell-by-cell map comparison, which considers whether each pair of cells on the two maps are equal or not. This results in a comparison map displaying the spatial distribution of agreement. This comparison method does not require any parameters. The kappa statistic ranges from −1 to +1, where +1 indicates perfect agreement and values of zero or less indicate a performance no better than random (Cohen, 1960;van Vliet et al., 2011).

Species Discrimination
Linear discriminant function analysis (i.e., discriminant analysis) performs a multivariate test of differences between groups. This classification technique introduced by Fisher (1936) is based on the idea of discriminant function analysis to provide a set of weightings that allow the groups to be distinguished. It requires that the individuals be divided into groups and it is used when there are observations from pre-determined groups with two or more response variables recorded for each observation. It then generates a linear combination of variables that maximizes the probability of correctly assigning observations to their predetermined groups. Linear discriminant function analysis was used to test the efficacy of the bioclimatic variables to classify the species of Rhodniini. Multivariate niches with all 19 bioclimatic variables were used to describe differences between the species (Supplementary Table 4). The values of bioclimatic variables were extracted from the WorldClim 2 layers (Fick and Hijmans, 2017). Twenty locations for each species were randomly selected, except for R. brethesi Matta, 1919 (12), R. domesticus (15), and R. stali Lent et al., 1993 (10), which selection comprised all locations. The discriminant analysis was performed using the stepwise method, Willks' Lambda (λ), and computing F values for pairwise distances. Wilks' λ is the ratio of within-groups sum of squares to the total sum of squares and it is designed to indicate whether a particular variable contributes significantly to explaining additional variance in the dependent variable. There are an F and a p-value associated with Wilks' λ that indicates the level of significance (Quinn and Keough, 2002;McLachlan, 2004). The analysis was performed in SPSS SOFTWARE 20 (IBM, 2011). The criteria for entry sets were F values of 3.84 to enter a new variable and 2.71 to remove an already entered variable.
These F values represent significance levels of approximately 0.05 and 0.10, respectively. We defined F values from groups size option for prior probabilities for classification-the groups are weighted by the proportion of the number of cases in each group.

Sequence Alignment, Phylogenetic Analysis, and Molecular Dating
The parsimony analysis combining 16S + 28S + Wg sequences retained five trees with CI = 0.838, −ln L ranging from 6101.613 to 6110.204, and AIC from 12245.226 to 12262.408 ( Table 1).
The tree in Figure 1A shows the topology with the best scores FIGURE 1 | The topology of the retained trees in the parsimony analysis combining 16S + 28S + wingless sequences retained five trees with CI = 0.838, -ln L ranging from 6101.613 to 6110.204, and AIC from 12245.226 to 12262.408 (see Table 1). (A) the topology with the best scores for the test for significance of likelihood-score differences-Kishino-Hasegawa test (KH). (B) bootstrap support values employing a heuristic search with 1,000 bootstrap replicates using 100 random stepwise addition (tree-bisection-reconnection, TBR). We defined values for bootstrap support > 90% as strongly supported, 90-70% as well-supported/moderate support, and <70% as weakly supported. (C) the topology showing the alternative position of the genus Psammolestes. It was used PAUP 4.0 in the phylogenetic analysis to derive trees based on parsimony.
Frontiers in Ecology and Evolution | www.frontiersin.org for the test for significance of likelihood-score differences-Kishino-Hasegawa test (KH) ( Table 1). The bootstrap values varied between 62 and 100, although the species R. domesticus + Psammolestes had no bootstrap support ( Figure 1B). The stronger bootstrap value was seen for the pictipes group (100%). The pallescens group had weakly supported bootstrap values (65%) and the prolixus group was well-supported (87%). One of the five trees had a phylogenetic relationship where Psammolestes species-P. tertius and P. coreodes-were not retained in the clade R. domesticus + prolixus group (Figure 1C, −ln L = 6104.380, AIC = 12250.760). This topology was not showed in any publication prior to our study. The result of the maximum likelihood analysis (ML) is shown in Figure 2 and showed better supports than the trees retained in the parsimony analysis (−ln L = 5895.714, AIC = 11833.428). The main difference in the topology of this tree is the split of the tribe into two evolutionary lineages-pallescens group + pictipes group, and R. domesticus + Psammolestes + prolixus group, contrasting the trees retained in the analysis of parsimony that showed three clades-pallescens group, pictipes group, and R. domesticus + Psammolestes + prolixus group. The results of our phylogenetic analyses did not show the species Z. alcides and O. chinai as closely related to Rhodniini, as previously found (Patterson and  Gaunt, 2010;Hwang and Weirauch, 2012;Justi et al., 2016). The results support the hypothesis proposed in de Paula et al. (2005), where Salyavatinae was the sister clade to Rhodniini. The BEAST analysis (Figure 3) showed that O. chinai and Z. alcides are more recent than the origin of the Rhodniini and originated 12.91 Mya (95% credibility interval: 0.02; 0.67). This more recent dating for the species of the subfamily Reduviinae in our results reinforces the hypothesis that Salyavatinae is closely related to Rhodniini. Four main clades appear in this figurepallescens group, pictipes group, R. domesticus + Psammolestes, and prolixus group. The posterior probability values are strong for most of the species of Rhodniini, ranging from 1 to 0.84.

Ecological Niche Model
All ecological niche models had Training AUC values above 0.8 and Test AUC values above 0.8, which confirms the high sensitivity of the analysis (Supplementary Table 5).
The distribution of the species Rhodnius ecuadoriensis, R. colombiensis, and R. pallescens are shown in Figures 4A-C, and these species form the pallescens group (trans-Andean). All these species have combinations of Kappa values below 0 (performance no better than random). The Fraction correct was above 0.95. These results indicate speciation by vicariance, where the species are clearly allopatric. The pictipes group includes Rhodnius brethesi (Figure 4D), R. stali (Figure 4E), and R. pictipes ( Figure 4F) (cis-Andean). Rhodnius brethesi and R. stali are clearly allopatric (Kappa = −0.083, Supplementary Table 6). Rhodnius pictipes distribution overlapped with that of R. brethesi (Kappa = 0.218) to the north of its distributional area. Rhodnius pictipes and R. stali are closely related but had little overlap in their areas of potential distribution (Kappa = −0.119). Based on this, we can assume that R. stali and R. brethesi originated by vicariance, and R. pictipes would have also originated by vicariance followed by post-speciation dispersal.

Species Discrimination
Multivariate niches with all 19 bioclimatic variables were used to describe differences between the species. Discriminant dimensions for all Rhodniini species (Figure 7) were significantly different in the results of pairwise F-ratios for testing group differences (Supplementary Table 7). In total, 75.4% of original grouped cases were correctly classified, and 66.5% of crossvalidated grouped cases were correctly classified. Wilks' λ showed p < 0.05 for all bioclimatic variables, indicating that they are indeed significant predictors for Rhodniini species (Tests of Equality of Group Means, Supplementary Table 8). Variables in species discrimination analysis after 12 steps (Supplementary Table 9) included Temperature seasonality, Precipitation seasonality, Isothermality, Mean temperature of driest quarter, Annual mean diurnal range, Annual precipitation, Precipitation of driest quarter, Precipitation of warmest quarter, Precipitation of wettest quarter, Precipitation of wettest month, Annual temperature range, and Mean temperature of wettest quarter. These results indicate that Rhodniini species clearly show distinct multivariate niches.

Historical Biogeography
The software TREEMAP was used by de Paula et al. (2007) to deduce taxon-area associations to propose a biogeographical hypothesis of the Rhodniini in the Neotropics. In this study, tanglegrams were used to show the relationship between the Neotropical areas of endemism and Rhodniini species. TREEMAP generated 12 optimal solutions with the lowest total cost. These optimal reconstructions required six vicariance events, 20 duplications (sympatry), at least three dispersals, and at least one extinction event. In our study, though, results indicate a parsimonious scenario comprising vicariance events. The origin of the Rhodniini tribe occurred 17.91 Mya (Figure 3), and this result is more recent than the study of Hwang and Weirauch (2012), which defined the origin of Rhodniini as 22.18 Mya. This divergence estimate is much younger than the 107 Mya age proposed by Patterson and Gaunt (2010). Although molecular and biogeographic divergences support the diversification of the species of the Rhodniini tribe (as discussed below), karyotype and chromosomal homogeneity demonstrate that species of this tribe have not suffered structural changes in chromosomes during speciation events Ravazi et al., 2021).
The formation of the Pebas system (23-10 Mya, Hoorn et al., 2010) split the clades that originate the pallescens (trans-Andean) and pictipes groups (cis-Andean) (17.91 Mya). The rapid uplift of the Eastern Cordillera between 5 and 2 Mya (Hoorn et al., 2010) may be related to the origin of the species R. colombiensis and R. pallescens (1.09 Mya). The tectonic uplifts of the Northern Andes, called Quechua 2 (9-8 Mya), may be related to the origin of R. ecuadoriensis (7.93 Mya). Hydrological connections between the Araguaia-Amazon basins in Early and Middle Miocene (20.4-9.0 Mya, Albert et al., 2018) would be a possible additional event that split the cis-Andean + trans-Andean clades of the R. domesticus + Psammolestes + prolixus group (15.2 Mya). This hydrological connection incorporated the Pantanal system in the Upper Miocene-Lower Miocene (10.0-4.5 Mya, Albert et al., 2018), resulting in the separation of Pantanal from the Atlantic Forest. This event may have contributed to the formation of Psammolestes species when P. coreodes and P. tertius were then originated (4.98 Mya).
New species establish themselves in new habitats under new climate regimes and can persist in a changing environment. Rhodnius nasutus occurs in the Borborema Province (NE Brazil), which consists of a complex mosaic of rocks composed mainly by Proterozoic metasedimentary belts that amalgamate Archeanto-Paleoproterozoic gneiss-migmatite complexes (Oliveira and Medeiros, 2018). The distribution of R. nasutus is bounded to the west by the Parnaíba River, which is trapped between the Amazonian craton and the Borborema orogenic belt (Daly et al., 2014). Arid biomes mainly expanded after 15 Mya ago by pronounced cooling events during the Miocene (10 Mya). Nearly all the dated shifts into the arid biome occurred after this time (Crisp et al., 2009). Rhodnius nasutus appeared 6.92 Mya (Figure 3) and occurs in Caatinga Steppe, an arid biome. We here propose that R. nasutus originated during cooling events, possibly when it colonized carnaúba palm trees [Copernicia prunifera (Miller)].
Rhodnius prolixus originated 3.68 Mya and its distribution area may be related to the hydrological connection between the Llanos-Amacuro basins that occurred between 10 and 4.5 Mya (Albert et al., 2018). Rhodnius robustus and R. neglectus were dated at 2.26 Mya and would probably be related to the hydrological connection between the Amazonas-Araguaia basins that occurred between 4.5 and 0 Mya (Albert et al., 2018). Lastly, the transcontinental Amazon Stage 2 (Albert et al., 2018) occurred since middle Pliocene to Recent (4.5-0 Mya), and may be related to the origin of R. stali (1.01 Mya). Schofield (1988) proposed that different species and genera of triatomines have different ecological preferences and postulated that Triatominae may not represent a clade, but different hematophagous lineages that could have evolved independently of predatory ancestors that shared microhabitat. According to Hwang and Weirauch (2012), the evolution of bloodfeeding may have occurred once or twice independently among predatory assassin bugs. All prey specialists evolved from generalist ancestors, with multiple evolutionary origins of termite and ant specializations. The hematophagous habit seen in Rhodniini is associated with major morphological, physiological and behavioral adaptations that have accumulated throughout the evolutionary history of the various lineages of bloodsucking arthropods (Lazzari et al., 2013). Hwang and Weirauch (2012) showed that the generalist predatory feeding strategy is ancestral for Reduviidae, and that all prey specialists evolved from generalist ancestors. The transition from predatory to hematophagous habits could have occurred by the colonization of vertebrate nests. Ambrose (1999) suggested that the reduviids could be broadly divided into two groups based on whether or not they possessed tibial pads. Those without tibial pads live in tropical forest ecosystems and are known as timid predators that do not use their forelegs to capture prey, instead impaling prey items with their long rostra. Rain forest reduviids may have developed tibial pads and other features that made them more efficient predators when they migrated to deciduous scrub forest and other semiarid habitats. Although the forelegs are more simplified in Triatominae, some species exhibit similar extendible adhesive organs (spongy pads) in both the fore and middle tibia of the adults, which allow them climbing smooth surfaces (Weirauch, 2007), a useful adaptation of hemipterans that fly and live in trees, as do species of the Rhodnius genus. Ambrose (1999) also considered the members of Salyavatinae ancestral to the subfamilies Triatominae and Ectrichodiinae. This hypothesis is supported by de Paula et al. (2005).

Evolution of Hematophagy
The pantropical Salyavatinae, subfamily used as outgroup in this study, are suspected specialist termite predators (Cobben, 1978;van Doesburg and Forero, 2012), and show morphological and camouflaging behavior that may represent adaptations to FIGURE 7 | Linear discriminant functional analysis for Rhodniini. All species were significantly different in the results of pairwise F-ratios for testing group differences (see Supplementary Table 7). 75.4% of original grouped cases were correctly classified. 66.5% of cross-validated grouped cases were correctly classified.
termite feeding (Weirauch, 2006). Considering Salyavatinae and Rhodniini closely related (de Paula et al., 2005(de Paula et al., , 2007, we propose a scenario where the Rhodniini's ancestor could be pre-adapted to the invasion of bromeliads and palm trees, where they would find significant water availability and thermal damping. Rhodniini are closely related to Salyavatinae that are specialist termite predators and diverged from this subfamily 30.43 Mya (Figure 3). Dias et al. (2011) suggested that microclimate plays an important role in establishing a stable relationship between Rhodniini species and palm trees, further influencing the adaptation of the bugs to different palm genera. Rhodniini expresses a clear thermo and hygropreference, such as amenable microclimates in terms of temperature and relative humidity (Lazzari et al., 2013). Furthermore, Rhodniini species need to find adequate shelter for protection from daylight, which can make them vulnerable to potential predators. Besides that, food source detection in Rhodniini is achieved via the detection of air currents carrying distinctive odors, water vapors and heat (Lazzari et al., 2013). Natural habitats (burrows, nests, hollow tree-trunks) provide an abundance of shelters, although many of them may not offer suitable conditions. Shelters can be subject to high temperatures, humidity and intense illumination or may not be inhabited by colonies of conspecifics. We believe that the transition from the predatory lifestyle to hematophagy may have been a gradual process that involved specialist termite predators who sought shelter in vertebrate nests and then opportunistically fed on arthropods that inhabit nests. In laboratory conditions, R. prolixus nymphs showed signs of cannibalism and cleptohematophagy (Piñero et al., 1988).
The species of the tribe Rhodniini included in this study have clearly distinct bioclimatic niches (Figure 7 and Supplementary Table 7). Consequently, the colonization of bromeliads, palm trees and nests represent important events for the speciation of these taxa. Rhodnius domesticus can be found in bromeliads, rodents and marsupial refuges, hollow trees, and under barks (Lent and Wygodzinsky, 1979;Monteiro et al., 2018). The results of our study indicate that R. domesticus is the oldest species of Rhodniini, having its origin at 9.13 Mya (Figure 3). In this scenario, the invasion of palm trees would have occurred only once, with subsequent colonization of bromeliads and bird nests in the R. domesticus + Psammolestes lineage. These clades had their origins at 15.2 Mya (pictipes), 7.93 Mya (pallescens), 11.01 Mya (prolixus + R. domesticus + Psammolestes), respectively (see Figure 3).
Psammolestes species are adapted to explore bird nest microhabitats, and seem to use palms opportunistically (Abad-Franch et al., 2015). Psammolestes species are wild and found in bird nests, mostly from birds in Furnariidae family (Lent and Wygodzinsky, 1979). For instance, the nest of Furnarius rufus (Gmelin) (Furnariidae) is a domed mud structure, with a clear separation between the breeding chamber and the outside. In general, each couple builds one nest per year, but each nest is used for one clutch or two consecutive clutches in the same breeding season. Some nests can remain in the field for 2 or 3 years, but some have a longer permanence (up to more than 8 years). Thus, the nests are widely used by other vertebrate inquilines too (Turienzo and di Iorio, 2010), which would be the source of food for the species of Psammolestes, in addition to birds. The species of Psammolestes, P. coreodes and P. tertius, have an origin dated 4.98 Mya (see Figure 3).
The results of this study indicate that the Rhodniini tribe originated 17.91 Mya ago, with R. domesticus being the oldest species, originated 9.13 Mya. The biogeographic history of the studied species comprised vicariance events. Rhodniini species clearly have distinct bioclimatic niches, and the colonization of bromeliads, palm trees and bird nests also represent important events for the speciation of these taxa. Their hematophagous habit can be described as a scenario where the ancestor of Rhodniini could have been pre-adapted to the invasion of bromeliads, palm trees, and bird nests, where it would have had available water and thermal damping. These habitats are used by vertebrates that the ancestor would have fed on. These environments are widely used by vertebrate inquilines that would be the food source for Rhodniini.

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 in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
ASDP: conceptualization, methodology, formal analysis, and writing. CB: methodology, formal analysis, and writingŮreview and editing. MT: methodology, formal analysis, and writing. LD and CG: conceptualization and writing. All authors contributed to the article and approved the submitted version.

ACKNOWLEDGMENTS
Special thanks to Dr. Danon Clemes Cardoso (Laboratório de Genética Evolutiva e de Populações-UFOP) for critical reading and suggestions in the first draft of the manuscript. We also thank Matthew Meehan (Western University) for final comments, and JADR and AR for their constructive criticism and time. CB and MT thank UFOP for the support of their undergraduate thesis research.