Chemosensory Gene Expression for Two Closely Relative Species Rhodnius robustus and R. prolixus (Hemiptera, Reduviidade, Triatominae) Vectors of Chagas Disease

Two closely related species, Rhodnius prolixus and Rhodnius robustus, are the vectors of Trypanosoma cruzi, which is the causative agent of Chagas disease, but clearly exhibit clear-cut differences in their ecological behavior. R. prolixus is considered as a domiciliated species, whereas R. robustus only sporadically visits human houses in Amazonia. We performed a chemosensory gene expression study via RNA-sequencing (RNA-seq) for the two species and also included a laboratory introgressed R. robustus strain. We built an assembled transcriptome for each sample and for both sexes and compiled all in a reference transcriptome for a differential gene expression study. Because the genes specifically expressed in one condition and not expressed in another may also reflect differences in the adaptation of organisms, a comparative study of the presence/absence of transcripts was also performed for the chemosensory transcripts, namely chemosensory proteins (CSPs), odorant-binding proteins (OBPs), odorant receptors (ORs), gustatory receptors (GRs), and ionotropic receptors (IRs), as well as takeout (TO) transcripts because TO proteins have been proposed to be associated with chemosensory perception in both olfactory and taste systems. In this study, 12 novel TO transcripts from the R. prolixus genome were annotated. Among the 199 transcripts, out of interest, annotated in this study, 93% were conserved between R. prolixus and the sylvatic R. robustus. Moreover, 10 transcripts out of interest were specifically expressed in one sex and absent in another. Three chemosensory transcripts were found to be expressed only in the reared R. prolixus (CSP19, OBP9, and OR89) and only one in sylvatic R. robustus (OR22). A large set of transcripts were found to be differentially expressed (DE) between males and females (1,630), with a majority of them (83%) overexpressed in males. Between environmental conditions, 8,596 transcripts were DE, with most (67%) overexpressed in the sylvatic R. robustus samples, including 17 chemosensory transcripts (4 CSPs, 1 OBP, 5 ORs, 1 GR, 4 IR, and 2 TO), but 4 genes (OBP19, OR13, OR40, and OR79) were overexpressed in the reared samples.

Two closely related species, Rhodnius prolixus and Rhodnius robustus, are the vectors of Trypanosoma cruzi, which is the causative agent of Chagas disease, but clearly exhibit clear-cut differences in their ecological behavior. R. prolixus is considered as a domiciliated species, whereas R. robustus only sporadically visits human houses in Amazonia. We performed a chemosensory gene expression study via RNA-sequencing (RNA-seq) for the two species and also included a laboratory introgressed R. robustus strain. We built an assembled transcriptome for each sample and for both sexes and compiled all in a reference transcriptome for a differential gene expression study. Because the genes specifically expressed in one condition and not expressed in another may also reflect differences in the adaptation of organisms, a comparative study of the presence/absence of transcripts was also performed for the chemosensory transcripts, namely chemosensory proteins (CSPs), odorant-binding proteins (OBPs), odorant receptors (ORs), gustatory receptors (GRs), and ionotropic receptors (IRs), as well as takeout (TO) transcripts because TO proteins have been proposed to be associated with chemosensory perception in both olfactory and taste systems. In this study, 12 novel TO transcripts from the R. prolixus genome were annotated. Among the 199 transcripts, out of interest, annotated in this study, 93% were conserved between R. prolixus and the sylvatic R. robustus. Moreover, 10 transcripts out of interest were specifically expressed in one sex and absent in another. Three chemosensory transcripts were found to be expressed only in the reared R. prolixus (CSP19, OBP9, and OR89) and only one in sylvatic R. robustus (OR22). A large set of transcripts were found to be differentially expressed (DE) between males and females (1,630), with a majority of them (83%) overexpressed in males. Between environmental conditions,

INTRODUCTION
In Latin America, bloodsucking bugs (Hemiptera, Reduviidae, and Triatominae) are the vectors of the parasite Trypanosoma cruzi (Kinetoplastida and Trypanosomatidae), which is the cause of the Chagas disease affecting about 6 million people (PAHO/WHO, 2020). At present, 151 extant and 3 fossil species of Triatominae are described (Zhao et al., 2021), including 136 species restricted to Latin America. Most of the species living in natural ecotopes feed on a large variety of vertebrates, but some of them found in anthroposystems can feed on humans and transmit them the parasite via the infected feces.
Among the Chagas disease vectors, Rhodnius prolixus is one of the most famous vectors because of its epidemiological importance and its broad distribution in both Central and Northeastern Latin America. Chemical control campaigns have considerably reduced the prevalence of the disease due to this species vector. Therefore, in 2011, all previously endemic countries of Central America have been certified as free of Chagas disease transmission occurring due to R. prolixus (Hashimoto and Schofield, 2012), but it has been reported again in Mexico in 2019 (Antonio-Campos et al., 2019). R. prolixus is considered as fully domiciliated, but some researchers described the sylvatic populations of R. prolixus in Venezuela (Sanchez-Martin et al., 2006;Feliciangeli et al., 2007;Fitzpatrick et al., 2008) and Colombia (Guhl et al., 2009;Rendón et al., 2015).
Rhodnius robustus is closely related to R. prolixus and is widely distributed in the Amazon region. This species only sporadically visit human dwellings, which is being attracted by artificial light sources, as described in Peru (Cabrera et al., 2013), Guyana (Bérenger et al., 2009), Venezuela (Feliciangeli et al., 2002), and Brazil (Aguilar et al., 2007). Its sporadic visits of domiciles have favored the transmission of T. cruzi to humans (Castro et al., 2010) especially because this species exhibits a high prevalence of natural T. cruzi (Miles et al., 1983). Moreover, the mode of transmission involving an oral route through the contamination of fruit storages or palm juice (açaí juice) by the infected sylvatic insects, including R. robustus, is progressively increasing in the Amazon region (Valente et al., 1999;Coura et al., 2002;Monteiro et al., 2010).
According to Lent and Wygodzinsky (1979), R. prolixus and R. robustus-morphologically very closer to each other-are distinguished only by the leg color of the oldest larvae and some differences in the basal plates of the male genitalia. In addition, the two species are differentiated by wing morphometrics (Matias et al., 2001;Márquez and Saldamando-Benjumea, 2013) and external female genitalia (da Rosa et al., 2014). The lack of diagnostic allozyme markers for Venezuelan populations of R. prolixus and R. robustus (Harry et al., 1992;Harry, 1993a) and the absence of differences for the male median process of pygophore made Harry (1993b) to suggest that Venezuelan populations of R. prolixus and R. robustus were not distinct taxonomic units. Using cytochrome b sequences, the paraphyly of R. robustus was highlighted as the clade composed of R. robustus from the Orinoco region (Venezuela) was more closely related to R. prolixus clade in comparison to other R. robustus clades from the Amazon region. Then, R. robustus was described as a complex of five operational taxonomic units (Monteiro et al., 2003;Abad-Franch and Monteiro, 2007). In laboratory strains, the evidence of the introgression between R. prolixus and R. robustus has been reported as 11 of the 15 presuming R. prolixus colonies that are tested for the whole R. prolixus genome sequencing project, exhibited past introgression events with R. robustus (Mesquita et al., 2015).
Irrespective of their way of life (domiciliary or sylvatic), the survival, reproduction, and social communication of individuals rely on the chemosensory system. The detection and integration of environmental chemical signals, including smell and taste, allow organisms to detect food, hosts, and predators. This system plays a key role in ecological adaptation to new or changing hosts or habitats (Benton, 2015). The recognition of chemical molecules as signals depends on the activation of specific proteins, including odorant-binding proteins (OBPs), chemosensory proteins (CSPs), odorant receptors (ORs), and gustatory receptors (GRs). A functional OR comprises a heterodimer formed by a specific OR and an ubiquitous coreceptor, named Orco. The soluble and secreted OBP and CSP proteins are proposed to bind and transport chemical molecules to membrane receptors (Pelosi et al., 2014). ORs and GRs recognize specific ligands and send a signal to the brain, translating the chemical signal into an electrical signal. This chain reaction leads to the insect response (Jacquin-Joly and Merlin, 2004;Sánchez-Gracia et al., 2009;Bohbot and Pitts, 2015). In R. prolixus, the silencing of Orco using the RNA interference (RNAi) technique affects the physiology and behavior of insects by decreasing their ability to engorge on the blood, the ecdysis, and oviposition rates and by increasing the mortality rate (Franco et al., 2016).
Changes in olfactory specificity and/or sensitivity could give insects an opportunity to recognize new environmental signals. This adaptation may be driven by fast evolution via gene duplication/loss events of chemosensory genes (e.g., Guo and Kim, 2007;Vieira et al., 2007;Xiao et al., 2013;Engsontia et al., 2014;Suzuki et al., 2018) but also by polymorphism (e.g., Leary  Brand et al., 2015;Wada-Katsumata et al., 2018;Ferreira et al., 2021) or a variation in gene expression (e.g., McBride et al., 2014;Purandare and Brisson, 2020). In a previous study on Triatoma brasiliensis (Marchant et al., 2016a), another Chagas disease vector, we showed that chemosensory genes, including OBPs and CSPs, were differentially expressed (DE) according to the populations. We also evidenced variations in the expression of takeout (TO) genes. TO proteins have been proposed to be associated with chemosensory perception in both olfactory and taste systems (Fujikawa et al., 2006) and have been involved in the circadian cycle linking the temporal cycle and feeding status information So et al., 2000;Meunier et al., 2007). This study focuses on the chemosensory transcriptome of the two closely relative species R. prolixus and R. robustus. Wild R. robustus was collected in the field in French Guyana to ensure the integrity of its taxonomic status as this country seems to be outside the geographic distributional area of R. prolixus (Bérenger et al., 2009). Moreover, R. robustus Amazon lineage seems to be a monophyletic clade (Abad-Franch and Monteiro, 2007). Because there were a few sylvatic R. prolixus populations with uncertain status, only the rearing strains of R. prolixus obtained from domiciliary colonies were used. We also used the strain of R. robustus and tested it for its potential introgression with R. prolixus. To perform expression analysis, we built a reference transcriptome containing transcripts from all samples. The reads from each sample were then mapped against this reference to quantify gene expression. A comparative study was then carried out based on the differential expression study and on the presence/absence of chemosensory transcripts [OBP, CSP, OR, GR, and an ionotropic receptor (IR)] and also TO.

Sampling, Sequencing, and RNA Extraction
Samples are given in Table 1. The first letter of each sample name indicates the species (P for prolixus and R for robustus); the second refers to sex (M for male and F for female), and the last two letters indicate the insect origin: EF for the strains from Brazilian rearing, ES for the strains from Swiss rearing, and G for the sylvatic samples from French Guyana. For R. prolixus, we used two reared strains, one from the Insetário de Triatominae da Faculdade de Ciências Farmacêuticas/Unesp/Araraquara, Brazil, including males and females (PMEF and PFEF) and another with males only (PMES) from the Swiss Biology Institute, University of Neuchâtel, tested by Marchant et al. (2016b) for the transcriptome assembly strategy. In this study, the reared samples of R. robustus (RFEF and RMEF) from the same Araraquara Brazilian laboratory in the same conditions like R. prolixus were included. The strain samples were fed on ducks once a month. Once insects reached the adult phase, they were bloodfed on ducks two times in the 1st week. The strain samples (males and females) were sacrificed 7 days after their last bloodmeal and their head stored in RNA later R (Thermo Fisher Scientific, Waltham, MA, USA) before RNA extraction. Sylvatic R. robustus from French Guyana (RFG1, RFG2, and RMG) was collected using light traps and stored in RNA later R before RNA extraction.
The taxonomic status of all samples was checked using cytb sequences (data not shown) as the cytb of the reared R. robustus samples was not differentiated from that of the reared R. prolixus, we considered that the reared R. robustus was introgressed by R. prolixus. The sylvatic R. robustus samples (RMG1 and RFG3) were clustered in the robustus clade IV according to the definition of Monteiro et al. (2003). For each species, we first extracted RNA from entire heads to obtain a large panel of genes in the transcriptomes (PMES and RFG1 samples). The eyes were removed from the specimens prior to RNA extraction to discard the pigments that could interfere with the use of chemical reagents during the estimations of extraction and the RNA amount. Secondly, we extracted RNA from the isolated rostrums and antennae (PMEF, PFEF, RFEF, RMEF, RFG2, and RMG samples) to enrich the samples in chemosensory transcripts as some are known to be expressed at very low levels, especially chemosensory receptors. For these tissues, we grouped for each sex 2-4 individuals for the sylvatic R. robustus and to have some variabilities for the inbreeding strain samples, 18-21 individuals.
The extraction was carried out by TRIzol R . The samples were then sequenced using the Illumina paired-end technology by the IBENS genomics platform (Institute of Biology of the Ecole Normale Supérieure, Paris, France) according to the protocol described in Marchant et al. (2016b). Raw data are available at the European Nucleotide Archive (https://www.ebi.ac.uk/ena), Accession: SAMEA3900098 to SAMEA3900105.

Transcriptome Assemblies
Transcriptomes from each sample were first assembled using the procedure developed by Marchant et al. (2016b) for the laboratory R. prolixus strains, which is clearly the most effective one compared to de novo strategies alone (i.e., Trinity software). The quality of individual assemblies was evaluated using several criteria, including contig lengths, N50, the total number of bases assembled, and contig number. To estimate the completeness of the transcriptome in terms of the expected gene content, we used the Benchmarking Universal Single-Copy Ortholog (BUSCO V3) assessment tool (Waterhouse et al., 2018).
A common reference transcriptome that was needed to compare the expression between samples was built using the Cuffmerge script in the Cufflinks package and contigs from all samples irrespective of the species/sex/condition. This transcriptome (RefT-PR-PA) was used for the presence-absence analysis. For differential expression analysis, a non-redundant transcriptome (RefT-PR-DE) was derived from the RefT-PR-PA transcriptome, selecting only the most expressed isoform per gene based on the indices of the fragments per kilobase per million mapped reads (FPKM).

Annotation of Transcripts of Interest
Transcripts of interest, namely chemosensory (OBPs, CSPs, ORs, and GRs) and TO transcripts, were annotated. To identify chemosensory transcripts, we searched the contigs assembled from the reads that are mapped onto the annotated genes in the R. prolixus genome (accessible at www.vectorbase.org). We completed the OBP annotation using OBP28 identified in Marchant et al. (2016b). For TO genes that were not yet annotated in the R. prolixus genome at the beginning of this study, transcriptome annotation was performed by running tblastn with an e-value threshold of 10 −6 and using sets of queries of insect proteins from the family of interest retrieved from the GenBank database. The selected contigs were then blasted in a third validation step to the non-redundant protein database (Blastx with an e-value threshold of 10 −6 ). Only contigs for which the best blast hit was a protein from the searched family were kept.

Gene Phylogenies
For the CSP family, open reading frames in the identified transcripts were inspected for all the samples. We used R. prolixus (Mesquita et al., 2015) and T. brasiliensis (Marchant et al., 2016a) as reference sequences. Amino acid sequence alignments were performed using MAFFT (Katoh et al., 2019), and the peptide signals were deleted before tree reconstruction. In the final tree, only the isoforms detected in at least three sequences from the data set were retained to discard misassembled sequences.
Maximum likelihood trees were built using IQ-TREE (Nguyen et al., 2015), with the search of the best-fit model according to the Bayesian information criteria (BIC) in ModelFinder (Kalyaanamoorthy et al., 2017). Node support was assessed using SH-aLRT parallelized over the bootstrap samples to maximize load balance and efficiency. Phylogenetic trees were displayed using iTOL (Letunic and Bork, 2006).

Presence/Absence Analysis of Transcripts of Interest
We used FPKM information provided by Cuffdiff from Cufflink. Positive FPKM indicated that reads are mapped to the gene revealing gene expression while null FPKM reflected a lack of read mapping on the given gene revealing the absence of expression. We used two positive FPKM thresholds, 20 and 100, that defined three intervals (0<FPKM<20, 20<FPKM<100, and FPKM>100) to evaluate gene expression.

Differential Expression Analysis
We performed a chemosensory gene expression study via RNAseq (Trapnell et al., 2010). Read pairs from each sample were mapped against the RefT-PR-DE transcriptome using BWA with default options (Li and Durbin, 2009). Reads with multiple hits were excluded to avoid a counting bias. Count data were therefore produced using SAMtool facilities . PMES and RFG1 samples were not used in this analysis because their RNAs were from whole heads and thus not comparable to the samples obtained from antennae and rostrum only. Furthermore, transcripts showing an expression limited to one environment were discarded from the data set before the differential expression study. Contigs with a low coverage were filtered using HTSfilter (Rau et al., 2013) using the DESeq normalization method and a threshold of five mapped reads on average. The sex status (male vs. female irrespective of their environment) and the environmental conditions (sylvatic vs. reared irrespective of their sex) were investigated as two independent factors in the model. The analyses were carried out by using DESeq2, turning off the independent filtering (Love et al., 2014). The Benjamini-Hochberg correction (Benjamini and Hochberg, 1995) was applied, and contigs were considered as DE when the p-adjusted value was <0.05. The contigs identified as DE were annotated by Blast on the non-redundant protein database (NR, version May 2015).

Transcriptome Assembly
Transcriptome assemblies showed N50, which comprises between 907 and 2,415 pb according to the recovery of samples

Annotation of Transcripts of Interest
As discussed in the previous publication of R. prolixus genome, Mesquita et al. (2015) annotated 19 CSPs, 27 OBPs, 106 Ors, 33 Irs, and 28 Grs, but some gene models were incomplete and/or annotated as pseudogenes. In the RefT-PR-PA transcriptome, all the 19 CSP transcripts were retrieved, but only 26 OBPs (lacking OBP2), 76 ORs, 27 IRs, and 26 GRs. The OBP28 annotated in Marchant et al. (2016b) was also retrieved in this transcriptome. Overall, 39 of the 213 genes annotated as chemosensory in the R. prolixus genome were not identified in our transcriptome. Their expression level may be lower than our detection threshold. They may also be expressed in other chemosensory organs such as legs or at different developmental stages or physiological conditions. However, their genomic annotation may also be erroneous as it remains speculative and no functional data are available for chemosensory genes in Triatominae. In addition to chemosensory genes, we annotated 23 TO transcripts (Supplementary Table 2). As 15 TOs were annotated from the antennal transcriptome of R. prolixus by Latorre-Estivalis et al. (2020), we used their nomenclature. From our data sets, we retrieved 12 TOs annotated by Latorre-Estivalis et al. (2020). TO1, TO4, and TO6 were absent from our data sets, but we described 12 new TOs. In total, both of the studies identified 27 TO transcripts.
In total, 199 transcripts were annotated and studied. All these transcripts were retrieved from the chemosensory samples (antennas + rostrum) but only 184 from the head samples, some ORs being not detected in this tissue. These results validated our strategy targeting chemosensory organs to enrich the transcriptome in chemosensory genes, especially chemosensory receptors, known to be expressed at very low levels.

CSP Phylogenies
Because all 19 CSPs annotated from the R. prolixus genome were found to be expressed in the transcriptomic samples, a phylogeny was performed for this protein family (Figure 1). The best model test for the reconstruction according to BIC was LG+G4. The CSP sequences were strongly conserved between the samples irrespective of the species/origin as no sequence variation in amino acids was observed for the most expressed CSPs. However, we found isoforms for the three CSPs (CSP4, CSP15, and CSP16). The CSP phylogeny confirms the cluster of the three CSPs, CSP2-4, found in the same R. prolixus genomic contig. Other clades were identified as CSP10, CSP11, CSP18, CSP19, and CSP12-15.
The phylogenetic CSP tree was reconstructed using 13 Hemiptera species, namely 3 Triatominae (T. brasiliensis, R. prolixus, and R. robustus), 3 Aphidae (M. persicae, A. gossypii, and A. pisum), 3 Miridae (A. suturalis, A. lucorum, and A. lineolatus), 3 Delphacidae (N. lugens, L. striatellus, and S. furcifera), and 1 Pentatomidae (N. viridula). The hematophagous Phthiraptera P. humanus was also included (Figure 2). It is noteworthy that four generalist CSP clades comprised at least 10 species attesting for the conserved orthologous sequences probably derived from an ancestral copy and suggesting that these CSPs may have a common function. By contrast, four other clades comprised three to five CSPs sequences for all species belonging to only one family. These clades represent CSP family extension, two for the Triatominae, and one for both the Miridae and the Delphacidae. No CSP extension for Aphididae was found with the used data. Moreover, one clade clustered only some CSPs from the hematophagous species, suggesting that these CSPs had probably a common function related to hematophagy. FIGURE 1 | Phylogenetic tree of chemosensory proteins (CSPs). Maximum likelihood tree in amino acids from the transcriptomic Rhodnius robustus and Rhodnius prolixus data sets (see Table 1 for sample labeling). Sylvatic R. robustus is indicated in green. The predicted R. prolixus protein genes (Mesquita et al., 2015) are indicated in purple and Triatoma brasiliensis transcripts (Marchant et al., 2016a) are indicated in red. A bootstrap > 70% is indicated on the node.

Presence/Absence of Transcripts of Interest
We used the count data (FPKM) to estimate the presence/absence of transcripts in our data set (Supplementary Table 3). Among the 199 transcripts of interest, 16 were not expressed in the head samples (3 CSP, 1 OBP, 9 OR, 1 IR, 1GR, and 1 TO) compared to the antennas and rostrum samples. It was noted that 24 ORs were not expressed in the R. prolixus head samples (PMES) neither 18 in the sylvatic R. robustus 1 (RFG1), and 15 were only expressed in the former and 9 in the latter.  Considering the antennas and rostrum samples (Figure 3), 10 transcripts of interest were exclusively expressed according to sex and 19 to the environment (reared vs. sylvatic samples). Four transcripts seemed to be expressed specifically in females (OBP9, OR72, OR89, and IR75b) and six in males (OBP7, OBP18, OR20, OR21, GR24, and TO23). A comparison between environments (reared vs. sylvatic samples) highlights 1 transcript (OR22) specifically expressed in sylvatic samples and 14 specifically detected in the reared samples, two of them (OBP19 and GR16) having a FPKM higher than 20.
Considering the reared and sylvatic R. robustus samples, both of the samples lacked the three transcripts (CSP19, OBP9, and OR89) that are found to be only expressed in R. prolixus samples but shared the presence of three specific transcripts, OR60, OR75, and IR75b, the latter had a high FPKM (>100). It was noted that GR13 transcript was absent in the reared R. robustus while the GR16 in sylvatic R. robustus.
Expression data are also very useful to confirm gene identification and therefore annotation. Mesquita et al. (2015) reported that the genomic sequences of Gr2, Or33, Or52, and Or106 diverged importantly from their counterparts from the same families, with evidence for an additional intron in the first exon. All four genes were expressed in our transcriptome and showed high RPKM values in our data (ranging from 22 for GR2 to 804 RPKM for OR106) validating their corresponding gene models (Supplementary Table 4). Figure 4 shows the ACP made from count data normalized with the rlog function in DESeq2. The reared samples were closer with respect to each other irrespective of their sex and species, but the sylvatic R. robustus samples were much more dispatched between sexes. This high divergence of sylvatic samples could reflect a sampling bias as these samples were composed of only 2-4 individuals, while the reared samples consisted of a pool of about 20 highly inbred individuals per sample. Similarly, sylvatic samples are highly variable due to uncontrolled environmental conditions compared to the reared samples obtained from controlled rearing.

DE Genes
The performance of sample clustering according to the expression data of the 30 most DE contigs in the samples confirmed the observations provided by the ACP. In particular, both sylvatic R. robustus males and females were welldiscriminated from other samples. Within the reared clade, sex was more discriminating than species condition (Figure 5).
When the sex condition was considered, numerous contigs (1,630) were found as DE between males and females, with a majority of them (83%) overexpressed in males (Table 3). When reared vs. sylvatic was considered irrespective of the sex, many more contigs (8,596) were found as DE with most (67%) overexpressed in sylvatic samples (Table 3).

DISCUSSION
Because the insect chemosensory system plays a key role in the survival, reproduction, foraging, social communication of individuals, and ecological adaptation, we focused on FIGURE 4 | ACP based on count data normalized by the function of rlog DESseq2. see Table 1 for a sample description. chemosensory transcriptomes of males and females of R. prolixus and R. robustus samples.

Divergence of Expression: Males vs. Females
In this study, a high number of transcripts (1,630) common to both R. prolixus and sylvatic R. robustus were identified as DE in the chemosensory transcriptome between the two sexes, with a majority of contigs (83%) overexpressed in males that probably reflects sexual dimorphism or is associated with specific sexual behaviors. Male chemosensory genes could be implicated in the perception of odors or pheronomes emitted by females. In R. prolixus, it was shown that odors from metasternal glands were more frequently emitted by females in greater quantities compared to males and may be involved in sexual communication (Pontes et al., 2008). Moreover, R. prolixus females release a pheromone that promotes flight initiation by males (Zacharias et al., 2010). However, most of the DE genes are poorly annotated and their precise function is unknown. Concerning the chemosensory transcripts, if none were DE, 10 were specifically expressed in one sex, OBP9, OR72, OR89, IR75b in females, OBP7, OBP18, OR20, OR21, GR24, and TO23 in males. One (IR75b) was expressed only in R. robustus females (reared and wild), two (OBP9 and OR89) only in R. prolixus females, and three (OBP7, OR21, and TO23) only in the reared male samples (R. prolixus and R. robustus). Similar expression profiles for antennal transcriptomes for males and females were notified by Latorre-Estivalis et al. (2017), whereas a significant differential expression for diverse chemosensory genes was observed between larvae and adults for R. prolixus. Using quantitative PCR analysis on R. prolixus, a significant expression in male antennae was observed for OBP26 and OBP27  and also for OR1, OR3, OR74, and OR80 . Concerning TO transcripts, we evidenced that TO23 was exclusively expressed in reared males. For R. prolixus, TO2 was highlighted by Latorre-Estivalis et al. (2020) as significantly downregulated in male antennae are compared to those of larvae. It is, however, difficult to compare previous R. prolixus transcriptomic studies with our study because the data and the methods differed, namely antennae and rostrum in this study vs. antennas only in the other cited studies, and RPKM vs. quantitative PCR (qPCR).
In other insects, several studies have also shown significant differences in the expression of chemosensory genes between sexes, as for Camponotus floridanus and Harpegnathos saltator ant species (Zhou et al., 2012), the mosquito Anopheles gambiae (Pitts et al., 2011), or the moth Sesamia nonagrioides (Acín et al., 2009). For another Chagas disease vector, T. brasiliensis for which chemosensory genes were annotated (Marchant et al., 2015), in our previous study, we found three CSP (CSP14, CSP15, and CSP16) and one OBP (OBP4) overexpressed in males (Marchant et al., 2016a). Interestingly, the CSP16 clustered with CSPs of other Hemiptera and could suggest a common function related to sex, while the other two CSPs (CSP14 and CSP15) were in a CSP expansion specific to Triatominae. In Drosophila, Gr32a and Gr68a genes are expressed only in males and are necessary for effective courtship behavior (Bray and Amrein, 2003;Miyamoto and Amrein, 2008). OBP99b gene of Drosophila (also called tsx) was expressed in the adipose tissue of males but not in females while OBP99a gene was specifically expressed in females (Fujii, 2002). These two sex-specific OBPs genes are involved in the perception of pheromones and the modulation of reproductive behavior in Drosophila. In Bombyx mori, the OR-1 was exclusively expressed in male antennae and was involved in the perception of the sex pheromone (Sakurai et al., 2004). In Drosophila, several TO proteins have also been identified as expressed exclusively in male heads and may be related to sex (Dauwalder et al., 2002).
Functional experiments are required to investigate if the transcripts evidenced in our study as exclusively or DE in one sex, could be involved in courtship behavior, mating control, and/or other sex-specific behaviors such as searching for oviposition sites in females. Functional studies of chemosensory genes in R. prolixus are still sparse. Oliveira et al. (2018) by silencing R. prolixus OBP27 using RNAi and examining the sexual behavior of the phenotype suggested that this OBP is likely involved in the reception of sex pheromones. Another study by Franco et al. (2018) focused on four ORs using the Xenopus oocyte expression system, but none ORs responded to sex pheromones.

R. robustus and Environmental Factors
While R. robustus is considered sylvatic but sporadically visits houses, R. prolixus is widely domiciliated. As the species-specific status of R. robustus is debatable for populations of the Orinoco FIGURE 5 | Sample-to-sample distances. Heatmap performed using DESeq2 shows the Euclidean distances between sylvatic (red) and reared samples (blue) as calculated from the regularized log transformation expression data of the 30 most differentially expressed (DE) contigs between environment-species. see sample information (Table 1).
region (see Harry et al., 1992;Abad-Franch and Monteiro, 2007), we used for this study R. robustus from the Amazon region, i.e., French Guyana clearly differentiated from R. prolixus (Abad-Franch and Monteiro, 2007). For the laboratory strains used in this study, the introgression between R. prolixus and R. robustus has been evidenced at both mitochondrial and nuclear levels, a phenomenon also reported in other studies (Mesquita et al., 2015).
Among the chemosensory genes studied, 93.33% (182/199) were conserved between R. prolixus and the sylvatic R. robustus that support the close genetic relatedness between the two species. In the comparison between the sylvatic R. robustus and the reared R. prolixus, two factors overlap: the ecology of the vectors and their taxonomic status. However, the comparison "+" indicates contigs overexpressed in the first condition and "-" indicates contigs under expressed in the first condition.
with the reared R. robustus samples could be of interest. While introgressed with R. prolixus, the reared R. robustus samples did not express all the R. prolixus chemosensory genes like the wild R. robustus. Contrariwise, some genes were only In addition to these uncontrolled environmental factors, the sylvatic/domiciliary status is likely to influence gene expression through host adaptation. Previous studies showed differences in the expression of chemosensory genes between insect species or populations adapted to different hosts. This is the case for the two moth species, one adapted to cotton (Helicoverpa armigera) and another to tobacco (H. assulta), which differ in the expression of OBP (Li et al., 2013). The differential expression analysis between the Kenyan population of S. nonagrioides sampled in sylvatic Poaceae plants (Typha domingensis) and the French population sampled in corn crops also revealed differences in the expression of an OBP transcript and one OR (Glaser et al., 2015). In a previous study on the Chagas disease vector, T. brasiliensis (Marchant et al., 2016a), we highlighted a significant overexpression of many chemosensory genes in the sylvatic samples of T. brasiliensis compared to the domiciliary ones. The phylogenetic reconstruction identified that OBP and CSP are the orthologs in T. brasiliensis and R. prolixus. Interestingly, one OBP (OBP23) and two CSPs (CSP1 and CSP11) were overexpressed in sylvatic samples from both T. brasiliensis and R. robustus. In the phylogenetic tree, the Triatominae CSP1 were in a specific clade, but the closest CSPs were from Delphacidae, including NlugCSP3. Binding characteristics studies of this N. lugens CSP and behavioral assays revealed that plant metabolites and flavoring agents (as non-adecane and 2-tridecanone) have high binding affinities in fluorescence competition-binding assays and displayed strong attractiveness to N. lugens (Waris et al., 2020). Concerning Triatominae CSP11, it is clustered with CSPs of the other Hemiptera in a generalist clade, including SfurCSP5. This S. furcifera CSP presented an antennae-enriched expression pattern, suggesting its potential olfactory function . Using ligandbinding experiments in vitro, Chen et al. (2018) suggested that SfurCSP5 plays a key role in host odorant perception as SfurCSP5 bound several rice plant volatiles. Such a function in odorant perception could be therefore suggested for the Triatominae CSP1 and CSP11.
Apart from chemosensory genes, two contigs annotated as TO seemed to be over-expressed in the sylvatic samples compared to the reared/domiciliary ones. Similar results were obtained for T. brasiliensis (Marchant et al., 2016a). TO proteins are involved in the circadian cycle, food metabolism, and foraging behavior So et al., 2000), and in the regulation of locomotor activity during foraging activity (Meunier et al., 2007). In the case of Chagas disease vectors, the overexpression of TO transcripts in the sylvatic R. robustus samples compared to the reared ones may be related to dispersal behavior differences vs. aggregation. The implication of this gene family in behavioral changes in migratory Locusta migratoria was also noticed (Guo et al., 2011). These locusts adapt their behavior in response to population density changes, moving from a gregarious to a solitary migration phase and a TO gene LmigTo1 was overexpressed in solitary individuals. By contrast, in Danaus plexippus (monarch butterfly of North America) an under-expression of TO protein in migrants compared to sedentary summer was shown (Zhu et al., 2008). This change in the expression of TO could be related to the sunshine duration variation used by migrants to help their navigation to their wintering areas.
Differences in expression between samples from different environments could be attributed to the adaptation of bloodsucking bugs to different hosts and thus utilizing different olfactory cues. In mosquitos, a comparison between anthropophilic Aedes aegypti and sylvatic Aedes aegypti formosus showed chemosensory changes, leading to a preference for human scent (McBride et al., 2014), and these changes were correlated with increased expression of the AaegOr4 that binds a molecule presented at high levels in human scent.

CONCLUSIONS
Rhodnius prolixus and R. robustus, two closely related species vectors of Chagas disease, exhibited a large common panel of chemosensory genes. This study showed species-specific expression of some chemosensory genes and the overexpression of some other genes according to sex or environment (sylvatic/rearing). These differences in expression may induce differences in chemosensory sensitivity by a regulation at the perireceptor level, via the differential expression of OBPs and CSPs, supposed to transport odorants through the sensillum lymph to the chemosensory receptors, and at the receptor level, via the regulation of the OR and GR expression. Further electrophysiological studies comparing chemosensory sensitivity between sex and the insect origin may support this hypothesis.
Further studies are also needed at the genus scale, including different domiciliary and sylvatic Rhodnius species to link the differences in the chemosensory expression of chemosensory genes to environmental conditions that could reflect the adaptation of bloodsucking bugs to different environments and/or hosts, and thus utilizing different olfactory cues. Moreover, studies are required to functionally characterize the genes specifically or DE in the domiciliated species to identify the key cues involved in environmental preferences. Such types of research will open up promising perspectives in vector regulation based on the manipulation in chemosensory behavior.

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 Materials.

AUTHOR CONTRIBUTIONS
AM and FM: data analysis and writing. EJ-J: study design and sample extractions. CA: study design. DB and J-MB: field collections. JR: strain collections. MH: study conception, data analysis, and writing. All authors edited the manuscript and approved the submitted version.

FUNDING
This study was funded by the French Agence Nationale de la Recherche (ADAPTANTHROP project, ANR-097-PEXT-009) and supported by the LabEx BASC (University Paris Saclay, France) and the Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP, process number 2016/08176-9, 2017/50329-0). AM was funded by the Idex Paris Saclay, France. Financial support had a role neither in study design, collection, analysis, and interpretation of data nor in manuscript writing.