Transcriptomic modulation in response to an intoxication with deltamethrin in Triatoma infestans, vector of Chagas disease

Background Triatoma infestans is the main vector of Chagas disease in the Southern Cone. The resistance to pyrethroid insecticides developed by populations of this species impairs the effectiveness of vector control campaigns in wide regions of Argentina. The study of the global transcriptomic response to pyrethroid insecticides is important to deepen the knowledge about detoxification in triatomines. Methodology and findings We used RNA-Seq to explore the early transcriptomic response of T. infestans after intoxication with deltamethrin. We were able to assemble a complete transcriptome of this vector and found evidence of differentially expressed genes belonging to diverse families such as chemosensory and odorant-binding proteins, ABC transporters and heat-shock proteins. Moreover, genes related to transcription and translation, energetic metabolism and cuticle rearrangements were also modulated. Finally, we characterized the repertoire of previously uncharacterized detoxification-related gene families in T. infestans and Rhodnius prolixus. Conclusions and significance Our work contributes to the understanding of the detoxification response in vectors of Chagas disease. Given the absence of genomic information from T. infestans, the analysis presented here constitutes a resource for molecular and physiological studies in this species. The results increase the knowledge on detoxification processes in vectors of Chagas disease, and provide relevant information to explore new potential insecticide resistance mechanisms in these insects. Author summary Chagas disease affects millions of people worldwide. In the Southern Cone, the development of pyrethroid resistant populations from T. infestans is related to vector persistence and affects the efficiency of vector control campaigns. Several studies have explored the causes of insecticide resistance in T. infestans populations. However, the global transcriptomic response after insecticide treatment has not been analyzed in this species so far. In this study, we obtained transcriptomic information which allowed us to characterize important gene families despite the absence of genomic information. Furthermore, we performed a quantitative analysis of gene expression after deltamethrin intoxication. The results provided here increase the knowledge on detoxification processes in vectors of Chagas disease, which is essential for the design of new vector control strategies.


Introduction
Chagas is a neglected tropical disease that can provoke disability and death.
Six million inhabitants of Latin America are infected, whereas a fifth of the population in this region remains at risk. In the last decades, human migrations spread Chagas all over the world (https://www.paho.org/es/temas/enfermedad-chagas). The causative agent of the disease is the protozoan Trypanosoma cruzi, primarily transmitted to humans by feces of triatomine bugs deposited during blood-feeding.
Vectorial transmission of T. cruzi occurs in wide zones of the Southern Cone, where Triatoma infestans is the primary vector [1]. In spite of its medical relevance, the genomic sequence of T. infestans has not been published to date. The transcriptomic data available was obtained from particular body structures [2][3][4], or from a normalized library [5].
Given the absence of vaccines and efficient treatments during the chronic stage of Chagas disease, the control of triatomine populations is crucial for reducing vectorial transmission. Pyrethroid insecticides have been used for the control of triatomines for more than 20 years, given their efficiency and favorable toxicological properties [6]. Other types of compounds were tested with dissimilar success.
However, no commercial formulations with low toxicological risk to replace pyrethroids are available to date. In this context, the existence of pyrethroid resistant populations of T. infestans in regions from Argentina and Bolivia represents a public health problem [6]. The higher resistance ratio levels described in this species are associated with mutations in the insecticide target site, the voltage-gated sodium channel [7,8]. Other resistance-associated mechanisms are reduced penetration of insecticide due to cuticular alterations [9] and enhanced detoxification metabolism [10][11][12][13][14][15]. A focus located in Güemes Department (Chaco Province, Argentina) in the Gran Chaco ecoregion is of particular relevance given its extension and high levels of insecticide resistance detected in some triatomine populations. In this region, vectorial transmission of Chagas could not be stopped, and a mosaic of T. infestans populations presenting either susceptible, low or high pyrethroid resistant profiles coexist, probably due to a non-homogenous use of pyrethroids in the area [16].
The quantitative study of the modulation of gene transcription in response to an insecticide could give cues for understanding both detoxificant response and insecticide resistance/susceptibility. In this sense, high-throughput RNA-Seq represents an approach to identify expression changes that could be involved in detoxification mechanisms, e.g. comparing transcription profiles of susceptible and resistant populations (as for example in [17][18][19]). Recently, the analysis of transcriptomic databases allowed the identification of new relevant components of the detoxifying response to insecticides in mosquitoes [20][21][22]. Studies addressing the transcriptomic response after a treatment with an insecticide in hemipterans are scarce (see for example [23,24]) and none of them have been performed in a vector of Chagas disease to date.
Studies on insect detoxification have been mainly focused on three protein superfamilies: cytochromes P450 (CYPs), carboxyl-cholinesterases (CCEs) and glutathione transferases (GSTs). These families have been previously characterized in triatomine species [13,25]. However, other protein superfamilies involved in stressful and/or toxic stimuli responses were not comprehensively analyzed in vectors of Chagas disease. These superfamilies include chemosensory proteins (CSPs), which have been studied in triatomine insects [26 -28] but not in the context of detoxification, Heat-Shock Proteins (HSPs) and ATP-binding cassette (ABC) transporters. Research on the role of these proteins in the response to xenobiotics can provide interesting information on susceptibility and resistance to insecticides in triatomines.
Chemosensory proteins are small soluble proteins present only in arthropods [29]. Until very recently, the assigned role of insect CSPs was reduced to olfaction through the transport of hydrophobic odorant molecules, but recent evidence points to other functions (reviewed in [30]). Among them, a possible role in insecticide binding was suggested for these proteins [21,31], and growing evidence points to a CSP role in the response to xenobiotics: an overexpression of CSPs has been observed in several insect orders after a toxic stimulus [22, [32][33][34] and in a pyrethroid-resistant mosquito population [21]. Moreover, a role of these proteins in insecticide susceptibility [35,36] and resistance [21] was recently proposed. All this recent evidence strongly indicates that, in addition to their involvement in the transport of odorant molecules, a role in insecticide susceptibility and/or resistance should be considered for CSPs.
Heat-shock proteins, also known as molecular chaperones, are involved in essential processes in the cell such as protein homeostasis. Their expression can be constitutive or induced by a wide variety of stressors, and are classified according to their sequence homology, molecular weight and function (reviewed in [37]). Among them, HSP70s can be classified as canonical and atypical, having the latter activity as co-chaperones [38]. A modulation in the expression of HSP70s in response to different toxic stimuli has been observed in insects [22, [39][40][41] and an upregulation of HSP70 expression was observed in a pyrethroid resistant mosquito population [42].
A role of these proteins in resistance to starvation was proposed in T. infestans [43].
A similar observation was made in R. prolixus, where a role of these proteins in survival under starvation but also after feeding was proposed, and the expression of HSP70 transcripts was induced after thermal stress and feeding [44].
ABC transporters are membrane proteins which mediate the movement of substrates using ATP. Two types of ABC transporters have been described: 1) Full transporters (FT) are functional and possess two nucleotide binding domains (NBDs) and two transmembrane domains (TMDs); 2) Half transporters (HT) that contain one NBD and one TMD domain, and require homo or heterodimerization to be functional (reviewed in [45]). This superfamily is further classified into eight subfamilies (A-H) based on the homology of their NBDs [46]. A new subfamily (ABCJ) was recently proposed in Ae. aegypti [47]. Despite their nomenclature, some subfamilies such as 8 ABCE and ABCF lack transporting roles as they only encode for NBD domains. ABC superfamily, which is much less studied in arthropods than in other organisms such as bacteria and vertebrates, has been associated to a wide variety of physiological functions including the excretion of toxic compounds (reviewed in [48]). In insects, these transporters have been linked to insecticide transport and/or resistance, specially those belonging to subfamilies B, C and G [49].
In the present study, we used RNA-Seq to analyze the transcriptomic response of T. infestans after topical application with the pyrethroid insecticide deltamethrin. This allowed us to find, through the evidence of their differential expression, possible candidates to be involved in the detoxification response of this insect, related to a wide variety of processes such as metabolism, cuticular rearrangements, transport and transcription/translation-related processes. We also performed a comprehensive characterization of some gene superfamilies that could be involved in this response. The obtained results extend the genetic knowledge about T. infestans and provide information of the initial transcriptomic response to a pyrethroid. To our knowledge, this is the first RNA-Seq experiment to study the detoxificant response in a Chagas disease vector, and the first quantitative transcriptomic study in T. infestans.  [51] was used in the paired-end mode to remove low quality bases from 5' and 3' ends with the parameters TRAILING: 5 and LEADING: 5.

Sublethal intoxication assay and sample preparation
Besides, the SLIDING-WINDOW parameter was set as 4:18 and the adapters were removed using the parameters ILLUMINACLIP:TruSeq3-PE-2.fa:2:30:10. Only reads longer than 50 bp were maintained.
Transcriptome de novo assembly, prediction of coding regions and quality assessment The resulting trimmed and quality-filtered paired-end reads of the eight samples were used to construct a de novo transcriptome assembly using the software Trinity v2.10.0 [52,53] with default parameters. The basic statistics of the resulting assembled transcriptome were assessed with the TrinityStats.pl script.
De novo transcriptome assemblies contain thousands of assembled transcripts, putative isoforms and chimeric transcripts that can affect subsequent analyses, like transcript identification and differential expression analysis. For this reason, and in order to reduce the redundancy and complexity of our database, we decided to extract the predicted coding sequences (CDS) and to cluster them according to their sequence identity. Coding regions of the assembled transcripts were predicted using TransDecoder software v5.5.0 (http://transdecoder.github.io).
Open reading frames (ORFs) with a product of at least 100 amino acids were predicted using the TransDecoder.LongOrfs script. BLASTp (v2.5.0+) and HMMscan (v3.2.1) searches were performed on the predicted ORF database using the complete UniProtKB/Swiss-Prot and Pfam-A databases as queries, with an e-value cutoff of 1e -50 . The results of these searches were used in the TransDecoder.Predict script to obtain the predicted CDS. To avoid redundant CDS, the resulting sequences were then analyzed with the software CD-HIT v4.8.1 in cd-hit-est mode [54] to cluster them considering a sequence identity threshold of 0.95. This non-redundant CDS database was used to perform the following analyzes. Finally, to assess the completeness of this database, BUSCO v4.1.4 [55] was used in protein mode against the hemiptera_odb10 lineage data set.

Analysis of voltage-gated sodium channel transcripts
To identify the transcripts encoding the fragment of the voltage-gated sodium channel (the target site of pyrethroids) where pyrethroid-resistance associated mutations have been described [7,8], a blastp search was performed in the translated non-redundant CDS dataset. The protein fragments of this sequence previously identified in T. infestans [8] were used as query, with an e-value cut-off set at 1e-5. To assess the presence of mutations in the fragment of interest on the identified sequences, the trimmed paired-end reads from each sample were mapped to the non-redundant CDS dataset using STAR software v2.7.1a [56] (with outFilterMultimapNmax set at 50, and outFilterScoreMinOverLread and outFilterMatchNminOverLread set at 0.3). The resulting bam files (sorted by coordinate) were merged and indexed using Samtools software v1.10, and visualized along with the reference using the software IGV v2.9.2. Finally, a nucleotide alignment was constructed using Clustal Omega [57] to compare the identified fragment along with previously reported sequences.

Transcript quantification and differential expression analysis
The non-redundant CDS database was used as a reference to quantify transcript abundance in every sample with Salmon software v1.4.0 [58], using the align_and_estimate_abundance.pl Trinity script and the trimmed paired-end reads 13 from each sample. The matrix of counts per transcript in each sample was reported with the abundance_estimates_to_matrix.pl Trinity script. The obtained matrix was filtered using the HSTFilter v.1.28 [59] in RStudio to eliminate transcripts with low expression and/or high variation. Following, the filtered matrix was used to perform differential expression analysis using DESeq2 v.1.28 [60] with the lfcShrink function. This is based on the "Approximate Posterior Estimation for Generalized Linear Model'' method that utilizes an adaptive Bayesian shrinkage estimator to obtain more accurate log2 fold-change values [61,62]. Differentially expressed (DE) genes were identified based on a False Discovery Rate (FDR) < 0.05 and submitted to BLASTx searches against: 1) UniProtKB/Swiss-Prot, 2) the predicted protein sequences from D. melanogaster (release FB2021_02, downloaded from FlyBase (flybase.org)); 3) the predicted proteins form R. prolixus (release 51, downloaded from VectorBase [26,63]). Additionally, HMMscan searches on protein sequences were performed using the Pfam-A database. The output of these searches was integrated using Trinotate v3.2.1 (https://trinotate.github.io), with an e-value 1e-5 to extract the positive hits. Additionally, BLAST searches against the non-redundant (nr) database from NCBI were made, identifying four sequences in the DE set with hits on noninsect sequences that were discarded from the analysis.
A heatmap of the DE set was created with pheatmap v.1.0 in RStudio, with count data transformed using the rlogTransformation function included in DESeq2 [60]. A dendrogram was plotted with hierarchical clustering among genes based on Euclidean distances and ward.D2 method for clustering.

Characterization of gene families and sequence analysis
HMMscan searches (with an e-value threshold for reporting models set at 1e-5) were performed on both the translated non-redundant CDS database from T.
infestans and in the annotated proteins from R. prolixus, available at VectorBase (version 51; [26,63]). Additional BLAST searches were conducted with an e-value cutoff set at 0.01, using tBLASTn to analyze the non-redundant CDS database from T. infestans and BLASTp to analyze the R. prolixus predicted proteins. The following PFAM seed alignments (used as queries for BLAST) and PFAM profile HMMs (used for hmmscan searches) were used for each family: PF00005 for ABC transporters,

PF00012 for HSP70 proteins and PF03392 for CSPs (searched only in T. infestans).
Iterative searches were conducted using resulting sequences as queries until no new sequences were identified. Finally, BLAST searches against the non-redundant database from NCBI were made and transcripts with hits on non-insect sequences were discarded from the analysis.
Partial transcriptomic sequences were manually analyzed and corrected if evident assembly and/or sequencing errors were detected, using sequences from other insects for comparisons. These manually curated sequences were only used for the sequence and phylogenetic analyses. To avoid redundancy, CD-HIT v4.8.1 was used to analyze partial sequences and only those with less than 90% of amino acid identity with other sequences in the family were kept. The CSP sequences previously reported for T. infestans [5] were manually curated as described and clustered with the set found in this study using CD-HIT (in cd-hit-est mode, using a sequence identity threshold of 0.95) to obtain a non-redundant set for further analysis. In the case of ABC transporters, sequences from A to H subfamilies were analyzed.
The final sequence set was analyzed with the following softwares: TOPCONS2 [64] for the prediction of transmembrane helices in the transporter proteins; Interproscan v85.0 [65] to analyze protein domains, and SignalP v5.0 [66] to assess the presence of signal peptide. In the case of CSPs, Clustal Omega [57] was used to align Drosophila melanogaster CSP sequences along with the predicted proteins to assess the presence of conserved residues and PSIPRED v4.0 [67]

Phylogenetic analysis
The protein sequences identified in T. infestans and R. prolixus were aligned with MAFFT [68] along with sequences previously reported in D. melanogaster for ABC transporters [46] and HSP70 (FlyBase gene groups FBgg0000497 and FBgg0000498). In the case of CSP family, sequences previously reported for R.
prolixus [26,28] and T. brasiliensis [27] were aligned along with T. infestans predicted proteins. For this, G-INS-i strategy with the following settings was used: unaligned level=0.1; offset value= 0.12; maxiterate= 1000 and the option leave gappy regions. The resulting alignments were trimmed using trimAl v1.2 [69] with default parameters except for the gap threshold, which was fixed at 0.3. Following trimming, IQ-tree v1.6.12 [70] was used to build phylogenetic trees based on the maximum-likelihood approach. The branch support was estimated using both the approximate Likelihood Ratio Test based on the Shimodaira-Hasegawa (aLRT-SH, [71]) and ultrafast bootstrap approximation (UFBoot) [72]. The best-fit amino acid substitution models, selected by ModelFinder [73] within IQ-tree and chosen according to the Bayesian Information Criterion, were LG+I+G4 for CSPs, LG+G4 for HSPs and LG+F+R7 for ABCs. The phylogenetic trees were visualized and edited with iTol online tool [74] and rooted at midpoint.

RNA-Sequencing, de novo assembly and prediction of coding regions
The sequencing output of the eight samples generated more than 185 million paired-end reads of 150 bp length of raw data ( Table 1 in S1 file). The resulting de novo assembled transcriptome contained a total of 431,511 transcripts and 292,177,510 assembled bases, with a GC percent of 34.18.
A total of 63,890 predicted CDS were identified in the de novo assembled transcriptome and 33,321 of them were kept after removing redundancy. The BUSCO searches revealed more than 96% of complete BUSCOs (S1 Fig), showing the completeness of our database.

Analysis of target-site mutations associated to pyrethroid resistance
Although T. infestans population from CC used in this work was classified as low-pyrethroid resistant, it is surrounded by localities where both susceptible and resistant populations were detected [16]. It is known that high pyrethroid resistance in T. infestans is associated with the fixation of target site (kdr) mutations in a population [7,8]. Hence, we hypothesized that the low resistance levels reported in CC population could be explained by the presence of a kdr mutation with middle or low frequency, as a consequence of insect migrations from neighbouring localities.
To address this hypothesis, we identified all the reads that mapped on transcripts encoding a fragment of the domain II of the voltage-gated sodium channel (S3 Fig), which is the target site of pyrethroids. Two point mutations associated with pyrethroid resistance (named L925I and L1014F) are located in this region of the gene [7,8]. All the detected reads have the wild-type (susceptible) sequence, indicating that the low resistance levels observed in CC are not explained by target site mutations. A similar result was found in Fronza et al. (2020) [15], where both substitutions were absent in CC insects analyzed. Interestingly, a silent substitution from thymidine to cytosine was detected in 63% of the mapped reads (base pair 183 in S3 Fig). This substitution was also found in high proportion in a resistant population from a neighbouring locality, which carries the non-silent kdr mutation L925I, and in low proportion in the susceptible population [8]. In addition, a microsatellite study performed in the zone proposed some degree of gene flow between CC and a near high-resistant locality [75]. This reinforces the hypothesis that the ancestral T.
infestans populations from the region had a common genetic background. The emergence of populations with different levels and mechanisms of resistance could have been promoted by discontinuities in the insecticide spraying and enhanced by the effect of environmental factors [76].

Differential expression analysis
Given the absence of target-site mutations, other detoxifying responses could be responsible for the low level of resistance observed in the population under study.
The quantitative analysis of the transcriptomic response to a low dose of deltamethrin could help to identify genes and gene families potentially involved in detoxification processes. Among all eight sequenced samples, more than 70 million reads mapped back to the coding regions of the assembled transcriptome ( Table 2 in S1 file). The differential expression analysis showed that 77 genes significantly changed their expression levels after deltamethrin treatment (FDR < 0.05); 61 of them were upregulated and 16 were downregulated (Fig 1, Table 3 in S1 file and S2  Table 4 in S1 file), indicating that they are putative proteins and not artifacts of the bioinformatic predictions.
Several cellular processes were affected after the insecticide application: among others, DE genes related to transcription and translation processes, rearrangement of the cuticle and transmembrane transport were identified ( Table 4 in S1 file).
Transcription and translation processes: More than 20% of the DE genes encode proteins related with transcription and translation processes and their regulation ( Table 4 in S1 file). From them, three upregulated genes possess a basic leucine zipper (bZIP_1) domain ( Table 4 in S1 file). Proteins belonging to bZIP family have been related to detoxification and insecticide resistance in other insects [77][78][79][80][81][82].
One of the transcripts that contains this domain is TRINITY_DN1498_c3_g1_i1.p1, which seems to be an homologue of the Cyclic-AMP response element binding protein B (CREB) (Fig 1 and Table 4 in S1 file).
Interestingly, the over-expression of this transcription factor was also observed in an imidacloprid resistant strain of Bemisia tabaci, and it was proposed that its phosphorylation has a key role in the resistance to this insecticide through the upregulation of a cytochrome P450 [80]. Another upregulated gene (TRINITY_DN1194_c5_g1_i1.p1), which encodes an helix-loop-helix DNA-binding domain (PF00010) (Fig 1 and Table 4 in S1 file), seems to be related to the D.  Table 4 in S1 file). An upregulation of chitinase-encoding genes was previously observed in Sogatella furcifera after treatment with different types of insecticides, including deltamethrin [23]. The transcriptional changes seen in our study indicate that the rearrangement of the cuticle could be relevant in the response of T. infestans to the pyrethroid. This is in agreement with a role of the integument in detoxification and resistance to insecticides in T. infestans [9,14]. 20 Transmembrane transport: The abundance of six transcripts related to transmembrane transport was affected by the insecticide administration (Table 4 in S1 file). One of the upregulated transcripts belongs to ABC transporters family, which will be further analyzed here (TRINITY_DN3318_c1_g1_i9.p1) (Fig 1 and Table 4 in S1 file). Two of them (TRINITY_DN1232_c0_g1_i3.p1 and TRINITY_DN2432_c0_g1_i2.p1), whose expression decreases after treatment, are solute carriers that belong to the Major Facilitator Superfamily (MFS) (Fig 1 and Table 4 in S1 file). The underexpression of some of these transporters was also observed in Ae. aegypti larvae after intoxication with an essential oil [22]. These results, along with previous evidence [84,85], support the hypothesis that MFS could participate in the response to xenobiotics.
Energetic metabolism: A significant difference in the abundance of some transcripts potentially related to lipid and carbohydrate metabolism was found in this study (Fig 1 and Table 4 in S1 file). Consistently, transcriptional changes related to lipid metabolism have been reported in Ae. aegypti larvae in response to toxics [22,86] and alterations in carbohydrate metabolism after pyrethroid treatment has been observed [87]. The regulation of energetic metabolism after intoxication in triatomines deserves to be studied.
In addition to these processes, a significantly higher abundance of a transcript (TRINITY_DN22143_c0_g1_i8.p1) that encodes for an odorant binding protein (OBP) was found after treatment with deltamethrin (Fig 1 and Table 4 in S1 file). In agreement with our results, the overexpression of some OBPs was observed after treatment with dissimilar toxic compounds in insects from different orders (see for example [22,33,86]). Moreover, the role of OBPs in the defense against xenobiotics was demonstrated in Tribolium castaneum [35,88]. Three members of T. infestans CSP family were also found DE (see details below). Our results reinforce the idea of the involvement of some OBPs and CSPs in the response to xenobiotics. Finally, the overexpression of one HSP belonging to the HSP70 family was observed (Fig 1 and Tables 3 and 4 in S1 file); a complete analysis of these proteins is presented below.

Characterization of selected detoxification-related gene families
The analysis of transcription modulation after intoxication carried out in this study did not reveal a differential expression of CYP, GST or CCE genes, although the expression of some members of these families was reported to be induced after insecticide treatment (see for example [23, [89][90][91]). Quantitative PCR showed an induction of CYP genes in both fat body [12] and integument [14] in T. infestans after intoxication with deltamethrin. The fact that a differential expression of these genes is not observed in this study may be due to differences in the populations, doses and tissues analyzed. Previous transcriptomic studies on Ae. aegypti larvae did not detect changes in the expression of GST, CYP or CCE exposed for 48 h to sublethal doses of the pyrethroid permethrin [86]. However, a time dependent modulation of transcripts belonging to these gene families was observed for Anopheles stephensi larvae at different times post-exposition to permethrin [92] and for a resistant population of Anopheles coluzzii after sublethal exposure to deltamethrin [93].
Moreover, genes belonging to the CYP superfamily were found overexpressed in pyrethroid resistant populations, including some of T. infestans [12][13][14]. The results

Chemosensory proteins
The analysis of the predicted CDS revealed 18 complete CSPs in T.
infestans, which encode proteins between 113 and 136 amino acid residues in length ( Table 5 in S1 file). The characteristic features of these proteins (such as signal peptide, four conserved cysteine residues and six α-helices, reviewed in [94]) were identified in most of these sequences, confirming the identity and completeness of the assembled CSP sequences. In addition, four partial sequences were detected; one of them encoded a very short protein (54 amino acids) that was excluded from phylogenetic analysis (  [95]. The results obtained here along with previous evidence indicate a lower variation in CSP numbers in triatomines.
The phylogenetic analysis revealed that CSP family is conserved in the subfamily Triatominae: with the exception of Triin-CSP20 and Triin-CSP21, all T.
infestans CSPs have an orthologue in R. prolixus and/or T. brasiliensis. This suggests an evolutive pressure to conserve both the structure of the family and the sequence of CSP across the kissing-bug evolution (Fig 2). Whereas  Triatoma infestans CSP4 and CSP21 were found upregulated, whereas CSP16 was downregulated in response to deltamethrin (Fig 1), pointing to a differential modulation of expression within the CSP family. Even though previous studies showed that CSPs tend to be upregulated during the detoxificant response (see for example [21,22]), a member of the CSP family was also downregulated in response to imidacloprid in Ae. aegypti [86]. This could be an adaptive response to Chemosensory proteins were clustered in 5 groups based on their expression patterns (Fig 3). The group containing the upregulated sequences CSP4 and CSP21 (along with CSP5 and CSP2) show the highest expression in both treated and control groups. Interestingly, CSP2, CSP4 and CSP21 are phylogenetically related (Fig 2). The high conservation at the sequence level and a similar expression pattern could indicate that genes encoding these CSPs are organized in genomic clusters in T. infestans, as previously observed for CSPs in insect genomes [29]. The group integrated by CSP7, CSP12, CSP16 and CSP17 is the second most expressed,

Fig 2. Phylogenetic tree of the CSP family from T. infestans, T. brasiliensis
and R. prolixus. The maximum-likelihood tree was constructed using IQ-Tree and protein sequences obtained from the transcriptomes of T. infestans (Triin-, red, available in Table 5 in S1 file) and T. brasiliensis (Tribr-, black, obtained from [27], and from the R. prolixus genome (Rhopr-, blue, NTE: N-terminus missing in gap; FX: gene model repaired based on de novo transcriptome assemblies; obtained from [28]). Sequences from T. infestans were named according to their R. prolixus orthologues, with the exception of Triin-CSP20 and Triin-CSP21 which have no orthologues in R. prolixus. Support values are shown for branches with both aLRT-SH and UFBoot values > 80. The resulting tree was rooted at midpoint.
followed by the group integrated by CSP1, CSP8, CSP11 and CSP15 and the group integrated by CSP10, CSP13, CSP19 and CSP20. The group containing CSP6, CSP9 and CSP14 shows the lowest expression levels within the family in both analyzed conditions (Fig 3).

Heat-shock proteins 70
We found 10 complete ORFs identified as HSP70 on the T. infestans transcriptome, and one of them, named Triin_HSP70-9, was significantly upregulated (Tables 3 and 4 in S1 file). We also performed an analysis of the HSP70 family in R. prolixus, where they have not been comprehensively characterized so far despite the availability of its genome. In this species, 13 HSP70 candidates were detected; six of them with complete CDS, and seven partial sequences (including one chimeric protein, RPRC003073-PA), probably due to errors in the automatic gene prediction process ( Table 6 in S1 file). The number of HSP70 found in T.
infestans and R. prolixus is similar to the numbers found in mosquitoes, which ranged between 7 and 11 [42].
Most T. infestans HSP70 proteins possess the HSP70 domain, the nucleotide binding domain, the peptide-binding domain and the C-terminal domain, with the exception of Triin_HSP70-1 and Triin_HSP70-2 ( Table 6 in S1 file). The presence of HSP70 domain was confirmed in all R. prolixus sequences, with the exception of the Expression data of T. infestans CSP3, CSP18 and CSP22 is not available because these sequences were obtained from [5] and they were not detected in the transcriptome generated in this work.
chimeric protein and the sequence RPRC002179-PA, which encodes a protein with less than 100 amino acids. The sequences RPRC011930-PA and Triin-HSP70-2 represent a complete ORF but lack the peptide-binding and the C-terminal domains ( Table 6 in S1 file).
The phylogenetic analysis included T. infestans and R. prolixus sequences, with the exception of the R. prolixus chimeric sequence and those shorter than 250 amino acids. It shows 6 clades that represent 1:1 orthologies between R. prolixus and T. infestans sequences (Fig 4).  Table 6 in S1 file), and Dmel-Hsc70Cb (HSP110) is grouped with Triin_HSP70-3 from T. infestans (Fig   4 and Table 6 in S1 file). 27 Expression analysis revealed that Triin_HSP70-5 and Triin_HSP70-7 show the highest expression within the family, followed by Triin_HSP70-6 ( Fig 5).

ABC transporters
The ABC transporters superfamily has not been analyzed in detail in triatomine insects, with the exception of a recent analysis that reported the amount of ABC genes in R. prolixus [96]. In this work, we found 65 and 58 candidate sequences for ABC transporters in T. infestans and R. prolixus, respectively ( Table 7 in S1 file), belonging to A-H subfamilies. The recently proposed ABCJ subfamily [47] was not analyzed here. A total of 21 candidates (10 from T. infestans and 11 from R. The maximum-likelihood tree was constructed using IQ-Tree and protein sequences obtained from T. infestans transcriptome (Triin-, red, available in Table 6 in S1 file) and from R. prolixus (RPRC-, blue, VectorBase ID is shown) and D. melanogaster (Drome-, green, FlyBase gene name and isoform are shown) genomes. Rhodnius prolixus chimeric sequence and those shorter than 250 amino acids were excluded from the analysis ( Table 6 in S1 file). Triatoma infestans sequences were numbered according to their position in the tree. Support values are shown for branches with both aLRT-SH and UFBoot values > 80. The resulting tree was rooted at midpoint. prolixus) were excluded from further analysis given that the PF00005 domain, characteristic of ABC transporters, was absent. According to our sequence and phylogenetic analyses, the ABC repertoire of T. infestans was composed of 55 sequences (45 complete and 10 partial) while R. prolixus presented 47 ABC transporters, being 10 complete and 37 partial/misannotated sequences ( Table 7 in S1 file). A recent work identified a similar number of ABC genes in R. prolixus (49) [96]; this subtle difference between studies may be due to the different search strategies employed.
Subfamily G is the most abundant in both triatomines analyzed, with almost one third of the total ABC sequences (16 in T. infestans and 14 in R. prolixus, Fig 6 and Table 7 in S1 file). Subfamilies C and H are also abundant, with the latter being more abundant and subfamily C being less abundant in T. infestans (8 ABCCs and 15 ABCHs) in comparison to R. prolixus (9 members both in ABCCs and ABCHs).
Interestingly, a recent analysis of the ABC complement of >150 arthropod species found that the H subfamily shows expansions in hemipterans [96]. This is consistent with the number of Triin-ABCH transcripts found in this work, where these sequences represent almost 30% of the entire family, although genomic information on the species is needed to rule out variants or products of alternative splicing.
Regarding the other families, ABCB accounts for more than 10% of the sequences within the family in both triatomine species (7 ABCB in T. infestans and 5 in R. prolixus, being HTs more abundant than FTs) and the subfamilies A, D and E have the same number of transcripts in both triatomine species (3, 2 and 1 members respectively, Fig 6 and Table 7 in S1 file).
The phylogenetic analysis shows a T. infestans orthologue for most R.  Table 7 in S1 file) and from R. prolixus (RPRC-, blue, VectorBase ID is shown) and D. melanogaster (Drome-, green, FlyBase gene name and isoform are shown) genomes. Only those sequences complete enough for phylogenetic characterization were included in the tree (remaining sequences were given a preliminary classification shown in brackets in Table 7 in S1 file). Triatoma infestans sequences were numbered according to their position in the tree. Support values are shown for branches with both aLRT-SH and UFBoot values > 80. The resulting tree was rooted at midpoint. levels of expression of its members: the most contrasting case is that of H subfamily, which has members that show very low expression in both conditions, while others are highly expressed (Fig 7). These differences could be due to the expression of members restricted to particular tissues.

Concluding remarks
Here we present a de novo assembled transcriptome and the first quantitative RNA-Seq study performed on T. infestans, which is the main vector of Chagas disease in the Southern Cone. Four hours after the exposure to a sublethal dose of deltamethrin, an activation of the transcription machinery was observed, leading to significant changes in expression of genes belonging to ABC transporters, HSP70, CSP and OBP families. Transcripts related to transcription and translation processes, energetic metabolism and cuticle rearrangements were also modulated after intoxication. In the absence of genomic information, the transcriptomic data provided here constitute a key resource for molecular and physiological studies in T.
infestans. This information will be necessary in further studies directed to particular genes, such as RNAi-mediated gene silencing, expression modulation in particular tissues and/or at different time points, as well as comparisons of gene expression in susceptible and resistant populations. The study of detoxification and insecticide resistance in T. infestans will be useful for the rational design of an integrated vector control strategy.  Supporting information S1 File. Table 1: Sequencing results and number of reads after processing raw data. Table 2: Number of reads mapped to the non-redundant CDS database. Table 3: Differentially expressed genes between acetone (control) and deltamethrin-treated samples. Table 4: Sequence and putative annotation of transcripts with a significant change in their abundance after treatment with deltamethrin. Table 5: Triatoma infestans CSP sequence analysis. Table 6: Triatoma infestans and R. prolixus HSP70 sequence analysis.