
==== Front
Genome Biol Evol
Genome Biol Evol
gbe
Genome Biology and Evolution
1759-6653
Oxford University Press UK

39031593
10.1093/gbe/evae158
evae158
Article
AcademicSubjects/SCI01130
AcademicSubjects/SCI01140
Transcriptomic Response to Pyrethroid Treatment in Closely Related Bed Bug Strains Varying in Resistance
https://orcid.org/0000-0002-7371-9177
Haberkorn Chloé CNRS, VetAgro Sup, UMR 5558, Laboratoire de Biométrie et Biologie Évolutive, Universite Lyon 1, Villeurbanne, France
IZInovation, 13 Rue des Émeraudes, Lyon 69006, France

https://orcid.org/0000-0003-2027-8504
Belgaïdi Zaïnab CNRS, VetAgro Sup, UMR 5558, Laboratoire de Biométrie et Biologie Évolutive, Universite Lyon 1, Villeurbanne, France

Lasseur Romain IZInovation, 13 Rue des Émeraudes, Lyon 69006, France

https://orcid.org/0000-0003-0909-2936
Vavre Fabrice CNRS, VetAgro Sup, UMR 5558, Laboratoire de Biométrie et Biologie Évolutive, Universite Lyon 1, Villeurbanne, France

https://orcid.org/0000-0002-2100-1542
Varaldi Julien CNRS, VetAgro Sup, UMR 5558, Laboratoire de Biométrie et Biologie Évolutive, Universite Lyon 1, Villeurbanne, France

Betancourt Andrea Associate Editor
Corresponding author: E-mail: chloehbk@gmail.com.
Conflict of interest The authors declare no competing interests.

8 2024
20 7 2024
20 7 2024
16 8 evae15812 7 2024
06 8 2024
© The Author(s) 2024. Published by Oxford University Press on behalf of Society for Molecular Biology and Evolution.
2024
https://creativecommons.org/licenses/by-nc/4.0/ This is an Open Access article distributed under the terms of the Creative Commons Attribution-NonCommercial License (https://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact reprints@oup.com for reprints and translation rights for reprints. All other permissions can be obtained through our RightsLink service via the Permissions link on the article page on our site—for further information please contact journals.permissions@oup.com.

Abstract

The common bed bug, Cimex lectularius, is one of the main human parasites. The world-wide resurgence of this pest is mainly due to globalization, and the spread of insecticide resistance. A few studies have compared the transcriptomes of susceptible and resistant strains; however, these studies usually relied on strains originating from distant locations, possibly explaining their extended candidate gene lists. Here, we compared the transcriptomes of 2 strains originating from the same location and showing low overall genetic differentiation (FST=0.018) but varying in their susceptibility to pyrethroids, before and after insecticide exposure. In sharp contrast with previous studies, only 24 genes showing constitutive differential expression between the strains were identified. Interestingly, most of the genes with increased expression in the resistant strain encoded cuticular proteins. However, those changes were not associated with significant difference in cuticular thickness, suggesting that they might be involved in qualitative changes in the cuticle. In contrast, insecticide exposure induced the expression of a multitude of genes, mostly involved in detoxification. Finally, our set of transcriptome candidate loci showed little overlap with a set of loci strongly genetically differentiated in a previous study using the same strains. Several hypothesis explaining this discrepancy are discussed.

Cimex lectularius
insecticide resistance
RNA-seq
transcriptome
CIFRE 2019/0800 Scientific Breakthrough Project Micro-be-have Universite de Lyon, within the program “Investissements d’Avenir” ANR-11-IDEX-0007 ANR-16-IDEX-0005
==== Body
pmcSignificance

Insecticide resistance in bed bugs has been studied by comparing distantly related strains, leading to the detection of a very large number of overexpressed genes, which is difficult to exploit for improving insect control. Here, 2 closely related strains were compared, showing that differences in expression between insecticide resistant and susceptible strains were mainly due to cuticular genes. A greater number of genes were detected as responding to insecticide treatment, delineating a common plastic response in both strains. This analysis, therefore, places cuticular resistance as a major phenomenon by which bed bugs protect themselves against insecticides.

Introduction

Cimex lectularius, also known as the common bed bug, is an obligate blood-feeding parasite that mostly feed on humans. Physical and psychological disorders induced in their human preys can range from allergic reactions (O’Donel Alexander 1984), to psychosis and paranoia (Goddard and Deshazo 2009). The current resurgence of this species, starting late 90s (Potter 2011), is likely to be due to globalization, together with a widespread second-hand market (Doggett et al. 2004) and growing insecticide resistance (Davies et al. 2012). Pyrethroids are the main insecticides used to fight this pest and are, therefore, suspected to be the main pressure selecting resistance (Romero et al. 2007).

Pyrethroid resistance in C. lectularius is thought to be due to several mechanisms, common in insects. Cuticular resistance can first impair insecticide penetration in the body. Studies based on RT-qPCR assays showed that cuticular protein genes, such as chitin synthase (CHS) or cuticle protein, had higher transcript levels in resistant bed bug strains (Mamidala et al. 2012; Koganemaru et al. 2013). Following penetration of the cuticular barrier, detoxification metabolism may be mobilized to degrade or excrete the insecticide. Indeed, an increased activity of several enzymes such as cytochromes P450s (Romero et al. 2009) and esterases (Lilly et al. 2016a) has been observed in pyrethroid resistant strains. In addition, administration of piperonyl butoxide (PBO), a primary inhibitor of some cytochrome P450 monooxygenases, was associated with a significant decrease in resistance, further suggesting their involvement in resistance (Romero et al. 2009). Increased expressions of glutathione-S-transferases (GSTs) were also detected in resistant juvenile bed bugs (Mamidala et al. 2011). Finally, mutations can affect genes involved in the nervous system functioning, and more specifically in the voltage-gated sodium channels (VGSC), targeted by pyrethroids (Dang et al. 2014, 2015; Akhoundi et al. 2015; Balvín and Booth 2018). A single nonsynonymous mutation can alter VGSC conformation and hinder pyrethroids binding, thus conferring knock-down resistance (kdr mutations). The kdr mutant L925I (leucine 925 to isoleucine) has been identified in most bed bug populations (88% of the 117 strains tested in United States in Zhu et al. 2010, 100% of individuals in France in Durand et al. 2012), although some bed bugs had an additional V419L mutation (valine 419 to leucine, 40.9% in Zhu et al. 2010) or, more rarely, V419L alone (2.7% in Zhu et al. 2010). To provide insights into the genomic sequences underlying insecticide resistance, we recently contrasted allele frequencies between resistant and susceptible strains of C. lectularius through a DNA pool-sequencing approach (Haberkorn et al. 2023). A 6 Mb superlocus showing high genetic differentiation between resistant and susceptible strains was identified. This region was enriched for SNPs showing strong genetic differentiation between strains, and for the presence of structural variants. This genomic region also contains the major QTL involved in pyrethroid resistance, as identified by Fountain et al. (2016). Additionally, several much shorter peaks of genetic differentiation were identified throughout the genome. These data thus suggest that the 6 Mb superlocus and possibly other smaller regions are involved in insecticide resistance, although the mechanisms and exact loci responsible for resistance remain to be confirmed.

In this work, using the same strains as in our genomic study, we compared the transcriptome of susceptible and resistant bed bugs exposed or not to insecticides. In order to identify genes possibly involved in insecticide resistance, both constitutive and plastic responses were explored. Transcriptomic approaches are particularly relevant for identifying genes involved in detoxification, since their efficiency is often related to their level of expression (Li et al. 2007). A few studies, all published before the first bed bug genome was released, have addressed this question (Adelman et al. 2011; Bai et al. 2011; Zhu et al. 2013). However, since they all used 454 technology, quantification of each transcript was not possible (because of the relatively low throughput of this technology). Instead, a small set of transcripts/genes was further investigated using qRT-PCR. Since then, to our knowledge, only one resistance study using RNA-seq was conducted, leading to the identification of a huge quantity of transcripts differentially expressed (DE) between resistant/susceptible strains (15,540 out of 51,492 expressed sequence tags (ESTs), Mamidala et al. 2012). Although these studies provided significant insights into mechanisms possibly involved in resistance in bed bugs, they all suffered from several caveats. First, all the comparative studies mentioned so far were conducted on strains having very different genetic backgrounds. Indeed, most of them used the Harlan strain (sampled in Fort Dix, NJ in 1973) as the susceptible reference strain, and compared it to strains sampled from distant locations (800 km away for Columbus, OH in Mamidala et al. 2012 and Bai et al. 2011, and 440 km for Richmond, VA in Adelman et al. 2011) which probably differ in many traits unrelated to resistance phenotype. Consequently, it is unclear whether the genes or transcripts identified in these studies are DE as a result of selection by insecticides or other selective factors, or as a result of nonadaptive genetic differentiation (due to genetic drift). This may explain why so many transcripts were detected as DE in Mamidala et al. (2012). Additionally, these studies focused on constitutive differences between strains, since no exposure to insecticide was performed. This may limit our power to detect important genes involved in resistance, since resistance genes are often induced upon insecticide exposure (Poupardin et al. 2008; Guedes et al. 2017).

In the present study, we analyzed the whole protein-coding transcriptome of 2 bed bugs strains differing in resistance to pyrethroids. The 2 strains, that were either exposed or not to pyrethroids, were genetically very similar, since the overall index of differentiation (FST) was only 0.018 (Haberkorn et al. 2023). Expression levels were then compared in order to identify (i) genes showing constitutive differences between the 2 strains, (ii) genes whose expression is altered after pyrethroid exposure, and (iii) genes whose expression is differentially altered after insecticide exposure depending on the strain (interaction term). Whereas category (iii) can highlight both constitutive difference between strains or plastic response, category (ii) shows a shared plastic response across strains. We then crossed these results with our previous whole-genome analysis on the very same strains (Haberkorn et al. 2023), in order to test whether the set of candidates obtained through the genomic and transcriptomic datasets significantly overlapped. As we pinpointed several cuticular genes overexpressed in the resistant strain with numerous nonsynonymous mutations in close proximity, we explored the possibility of a cuticular thickening in this strain, as it has been observed in other resistant strains (Lilly et al. 2016b).

Results

Assessment of Resistance Phenotype

London Lab (LL) and London Field (LF) strains were first exposed to deltamethrin (1 ng per insect), in order to assess their pyrethroid resistance status. This dose was chosen based on a previous analysis showing significant differences in mortality between the strains, i.e. a resistance ratio of 17X, computed on Lethal Doses for 95% of individuals (see Table 1 in Haberkorn et al. 2023). As expected, the mean mortality was significantly higher for LL compared to LF (59.4% versus 31.3%, n=24 for each strain, fisher test with P-value=0.044), although the difference was expected to be greater.

Analysis of Differential Gene Expression

RNA-seq was performed to compare the transcriptomes of the 2 strains both in the absence of insecticide (controls) and after insecticide exposure (n=48). The read counts table obtained after RNA-seq was first analyzed by principal component analysis (PCA; Fig. 1). The sum of these 2 axes explained 45% of the variance (53% by adding PC3). No clear pattern was observed, neither between strains, nor between treated and untreated samples; this suggests that the transcriptomes were overall relatively similar between strains and between treated and untreated samples.

Fig. 1. Projection of LL and LF populations on the top 2 principal components using PCA. Read counts obtained from DESeq were processed with vst, which computes a variance stabilizing transformation. Individuals are separated within each population between survivors of insecticide treatment and untreated (n=48). Each point represents the top 500 most variable genes of a sample (default parameter for plotPCA function).

Variance in read counts was partitioned using a model including a strain effect (LL and LF), a treatment effect (treated and untreated) and their interaction. Genes were considered as DE when the fold change exceeded 1.5 (in absolute value) and when Padj<0.05.

Constitutive Differences between Resistant and Susceptible Strains

The “strain” effect was first analyzed in order to identify genes showing constitutive differences in expression between the 2 strains. Out of 1,2187 genes, 20 genes were detected as up-regulated and 4 as down-regulated in LF compared to LL (Fig. 2).

Fig. 2. Volcano plot using enhanced volcano on strain effect (ngenes=12,187, Padj<0.05, LFC>or<|0.58|). Genes coding for cuticular proteins are labeled with a star.

Among the genes overexpressed in LF, 6 were putatively involved in insecticide resistance which represents a highly significant overrepresentation for these genes (χ2 with Yates correction for continuity, χ2=33.21, P-value=8.29e-09). Strikingly, all 6 genes were coding for cuticular proteins (Fig. 3, Table S1, Supplementary Material online), which also represents an overrepresentation of this specific category (χ2 with Yates correction for continuity, χ2=146.01, P-value<2.2e-16). Among overexpressed genes, the one showing simultaneously the highest log2 fold-change (LFC) and lowest Padj encoded a putative leucine-rich repeat-containing protein DDB_G0290503 (106669337). To date, this gene is not known to be involved in insecticide resistance.

Fig. 3. Cuticular proteins detected as significantly overexpressed in LF. Expression is given in log2 of normalized counts.

Insecticide-Altered Genes

The “treatment” effect was then analyzed, i.e. the difference between untreated and treated survivors. Genes were similarly filtered on LFC and adjusted P-values, leading to, respectively, 375 genes significantly up-regulated and 388 down-regulated after insecticide exposure (Fig. 4).

Fig. 4. Volcano plot using enhanced volcano on treatment effect (ngenes=12,187, Padj<0.05, LFC>or<|0.58|). Genes coding for acetylcholinesterase-like proteins are labeled with a star.

Among the genes overexpressed in treated survivors, 20 were putatively involved in insecticide resistance (Fig. 5, Table S2, Supplementary Material online), which does not constitutes a significant enrichment (χ2 with Yates correction for continuity, χ2=2.89, P-value=0.09). However, 2 out of the 6 resistance categories among overexpressed genes were enriched, namely ABC transporters (χ2=7.70, P-value=0.006), and other detox (χ2=5.76, P-value=0.02).

Fig. 5. Resistance genes by categories detected as significantly overexpressed in treated survivors. Expression is given in log2 of normalized counts.

Most of the induced genes identified in a resistance category were involved in the detoxification metabolism: 6 ABC transporters (4 G4, 1 G1, and 1 G23), 6 various other detox (including transcription factor cap’n’collar and MAF, and 3 sulfotransferases), 2 GST, and 2 P450 (9f2 and 6k1). Three acetylcholinesterase-like and 1 cuticular protein were also detected. All acetylcholinesterase-like genes detected (labeled with a star on Fig. 4) were in the top 15 of genes having the highest LFC. The up-regulated gene showing the highest LFC encoded a peroxidasin-like protein (106671105), whereas the up-regulated gene showing the lowest Padj encoded an elongation of very long chain fatty acids protein 7-like (106664652).

Of the down-regulated genes, 15 were detected in an insecticide resistance category (Table S3, Supplementary Material online), which is not more than what is expected under the null hypothesis of no association between resistance and down-regulation (χ2 with Yates correction for continuity, χ2=0.04, P-value=0.83). Similarly, no significant enrichment in any resistance gene category was detected (χ2 with Yates correction for continuity and FDR correction for multiple comparisons, all Padj>0.05). Among the 15 genes identified, 7 were coding for cuticular proteins, 2 for ABC transporters (both on LG 2), 2 for insecticide target and nervous system (1 sodium channel protein Nach and 1 acetylcholinesterase-like), 1 from other detox category (sulfotransferase), 2 P450, and 1 Red/Ox gene (superoxide dismutase [Cu-Zn]-like). The down-regulated gene showing both the lowest LFC and lowest Padj was an uncharacterized protein (106673864).

Interaction between Treatment and Strain

Only one gene showed a significant interaction between treatment and strain: a leucine-rich repeat protein soc-2. This gene had opposite expression patterns in the 2 strains upon insecticide exposure (106661213, Fig. 6). Leucine rich repeat (LRR) motif are found in many immune receptors of animals, from Homo to Drosophila (Dushay and Eldon 1998).

Fig. 6. Gene coding for a LRR associated with a significant interaction term between strain and treatment effects. Expression is given in log2 of normalized counts.

Combining Transcriptomic Results with the Genomic Data

In a previous population genomic study conducted on the same 2 strains (Haberkorn et al. 2023), we identified a set of SNPs and structural variants (inversions, duplications) showing high between-strain genetic differentiation. These loci may, therefore, delineate genomic regions involved in the resistance phenotype. This genomic analysis revealed a clustering of candidate SNPs in some genomic regions, and particularly in 3 adjacent scaffolds constituting a 6 Mb region. This “superlocus” was furthermore significantly enriched in putative resistance genes, and in particular in GSTs and P450s, suggesting a clustering of those detoxification genes involved in insecticide resistance. One of the aims of the present study was to cross the population genomic analysis with the present transcriptomic analysis in order to identify some overlapping gene sets. Genes identified by both analyses are expected to be strong candidates for being involved in the phenotypic differences observed between the strains. In the following sections, we will, therefore, test whether the genes identified here are distributed evenly across the genome or whether they are clustered in some genomic regions, in particular in the highly differentiated genomic regions previously identified.

Distribution of DE Genes

We first tested whether the DE genes for treatment or strain effects were distributed evenly along the genome, at the linkage group (LG) or scaffold levels. Our hypothesis was that genes constitutively overexpressed in resistant individuals could be found in the superlocus previously identified. However, no pattern of enrichment was observed for DE genes of both conditions, at any scale (binomial test with Bonferroni corrections, all P-values>0.23). Hence, contrary to the SNPs that were clustered in some genomic locations, DE genes for treatment or strain effects were scattered along the genome. Nonetheless, 6 significantly DE genes were identified in the 6 Mb superlocus, all of them being induced upon insecticide exposure. One of them was putatively involved in resistance, namely a GST on scaffold NW_019392763.1 (106664023, Table S2, Supplementary Material online).

Genetic Differentiation of DE Genes

FST mean values were computed by gene, and compared between DE genes and the non-DE pool. We hypothesized that FST values should be higher for DE genes, since both the genomic and the transcriptomic approaches are expected to be enriched for genes involved in insecticide resistance. Although no link was observed between FST value and differential expression for the strain effect (t-test, P-value=0.926), there was a marginally albeit nonsignificant effect for the treatment effect (t-test, P-value=0.059): FST values were slightly higher for DE genes upon insecticide exposure compared to non-DE genes (FSTDE=0.01356>FSTnon-DE=0.00826). This would suggest that DE genes upon insecticide exposure are located in genomic regions that are slightly more differentiated between susceptible and resistant individuals than the rest of the genome. However, the same test conducted within each individual LG did not reveal any pattern.

Exploration of SNPs Associated with Transcriptomic Differences between Strains

Out of the 24 genes showing constitutive transcriptomic differences between LL and LF, 12 genes carried at least 1 outlier SNP, including 4 resistance genes—all being overexpressed cuticular genes (Table S4, Supplementary Material online). Two of them carried outlier SNPs only in intergenic regions (106662986 and 106664047), whereas the others also had outlier SNPs in promoter (106664492) or intronic region (106666062). Finally, the only gene showing interaction between strain and treatment, i.e. the leucine-rich repeat protein soc-2 (106661213), was associated with 5 outlier SNPs (Table S4, Supplementary Material online). Two of those SNPs were localized in intergenic regions, one in a promoter and 2 in introns.

Exploration of Structural Variants Associated with Transcriptomic Differences between Strains

We first tested whether the genes identified in the present transcriptomic analysis showed a higher LFC value between strains when they contained a structural variant. The rationale was that structural variants may underlie transcriptomic changes, either through gene duplications or through rearrangements in promoter sequences. Overall, this was not the case (t-test, P-value=0.901). However, among the 20 overexpressed genes in LF compared to LL, 5 were detected inside SVs (Table S5, Supplementary Material online), including 2 cuticle protein A3A-like within a single putative inverted duplication located on LG5 (106664492 and 106664495). The 3 other genes were not part of any resistance category. Two genes out of the 4 under-expressed in LF compared to LL were part of SVs, but none of them had annotation suggesting their involvement in insecticide resistance. Finally, no structural variants were detected as overlapping with the only gene detected with an interaction between treatment and strain effects, namely LRR soc-2.

Exploration of Nonsynonymous Outlier SNPs Associated with Transcriptomic Differences after Insecticide Exposure

Genes that are similarly induced upon insecticide treatment in both susceptible and resistant strains may still underlie variance in resistance if their coding sequences differ. We thus checked whether the DE genes for treatment (375 up-regulated, 388 down-regulated in LF) carried nonsynonymous SNPs identified as outliers in our previous population genomic data. This analysis yielded 11 genes, including 7 up-regulated and 4 down-regulated (Table S6, Supplementary Material online). However, none of them belonged to a putative resistance category.

No Significant Difference in Cuticular Thickness between Resistant and Susceptible Strains

As our transcriptomic analysis revealed a clear enrichment for genes encoding cuticle proteins among the constitutively DE genes, we sought to compare cuticle thickness between the resistant and susceptible strains. Mean cuticular thickness were thus measured on bed bug legs extracted from both strains but no significant difference was observed (repeated measures ANOVA, P-value=0.572, Fig. 7).

Fig. 7. Analysis of Cimex lectularius cuticle. a) Example of a transverse section of a bed bug middle leg tibia in electron microscopy with the circle pattern and subsequent 12 cuticle thickness measurements (darker segments, letters ‘a’ to ‘l’). b) Cuticle thickness (in μm) of LL (n=10) and LF (n=10) strains. Mean values are represented with dots.

Discussion

The resurgence of C. lectularius likely stems from insecticide resistance evolution (Davies et al. 2012). Previous transcriptomic analyses focused on resistant versus susceptible strains with different genetic backgrounds. Consequently, it is unclear to what extent the genes or transcripts identified so far are truly relevant to insecticide resistance or other stresses, or whether they have any adaptive significance. Here, we provide a comprehensive transcriptomic analysis by comparing 2 bed bug strains with very similar genetic backgrounds (average genome-wide FST=0.018, Haberkorn et al. 2023) but differing in their resistance phenotypes. In addition to baseline constitutive gene expression, we analyzed gene expression before and after deltamethrin treatment in order to identify induced responses. We interpret our results in the light of a previous genome-wide analysis performed on the same strains, which identified genomic islands of exceptional differentiation (Haberkorn et al. 2023).

A Small Set of Constitutively DE Genes

We first compared the 2 strains in order to detect constitutive differences between strains that would be related to resistance. Only 24 genes (out of 12,187) were detected as differentially expressed, despite the insecticide assay confirmed the difference in resistance between the strains. However, this difference was lower than expected, which might be explained by individual variations in resistance level of the random insects selected.

Among the DE genes, 20 genes were overexpressed in the resistant strain, with a single category of putative resistance genes (cuticular genes). These results are in sharp contrast with previous comparative whole transcriptome studies performed in bed bugs, where several thousands of DE transcripts were identified (Mamidala et al. 2012), for multiple categories of putative resistance genes—esterases, GSTs and P450s (Adelman et al. 2011; Bai et al. 2011) Although part of the difference can be explained by the fact that all these previous studies did not have a reference genome at the time (and thus had to rely on de novo assembly of the transcripts, which can inflate the number of inferred genes), another important explanation lies in the heterogeneity of the strains used. Indeed, in all 3 cases, the strains were collected in distant locations (>440 km), meaning that they may differ in their allelic composition just because of other selective pressures experienced in the different environment they encountered (apart from insecticide), or because of genetic drift due to independent evolution. On the contrary, from our previous genetic analysis (Haberkorn et al. 2023), it appeared that LL and LF strains were the most closely related of the 4 strains we analyzed with an overall low genetic differentiation (1.2% relative to the maximum possible value). Genetic differences (that should underlie the constitutive expression differences observed) are thus more likely to be the consequence of selection for the major difference they have, i.e. the difference in pyrethroid resistance. This result thus challenges the view that transcriptomic changes involving various detoxification mechanisms (esterases, GSTs, P450) have evolved in pyrethroid-resistant bed bugs populations. Our rather stringent filters helped us to detect a small set of genes that, in our opinion, would deserve to be functionally tested, to advance our understanding of previously unknown insecticide resistance mechanisms. Our result suggests that among the genes that experience quantitative transcriptomic changes, cuticular genes seem to play a predominant role. The cuticular barrier is the first encountered by a molecule in contact with a bed bug, and its construction is developed from the earliest larval stage, making any modification constitutive.

A Key Role for Cuticular Genes?

Six cuticular genes were detected as constitutively differentially expressed, including 3 A3A larval cuticle proteins, 1 protein A1A and 2 uncharacterized cuticle proteins. Cuticular resistance has been observed in many pyrethroid-resistant insects, such as Helicoverpa armigera (Ahmad et al. 2006), Anopheles gambiae (Yahouédo et al. 2017), Culex pipiens (Pan et al. 2009), or bed bugs (Lilly et al. 2016b; Balabanidou et al. 2018). In particular, the larval cuticle protein A3A has also been detected among up-regulated transcripts in resistant insects (Gao et al. 2018; Zhang et al. 2022), including bed bugs (see Supplementary Table 6 in Mamidala et al. 2012).

Interestingly, 2 A3A proteins detected as up-regulated in LF strain were overlapping with a putative inverted duplication (Haberkorn et al. 2023). Whether the small difference in the estimated frequencies for this inverted duplication (12.3% in LF versus 5.7% in LL) underlies the observed difference in A3A transcription is unclear but surely deserves additional investigations.

Cuticle-based resistance may arise thanks to an increase of its thickness, and/or thanks to modifications of the composition of the cuticle (Balabanidou et al. 2018). A previous morphological study showed a positive correlation between cuticle thickness and resistance status of bed bug strains (Lilly et al. 2016b). ‘Resistant’ bed bug cuticles were significantly thicker than ‘intolerant’ ones. In our present analysis, no significant difference was found when comparing leg cuticle thickness between resistant and susceptible strains. Our data thus favor the scenario of a change in cuticle composition rather than an increased thickness.

Genetic Variants Associated with Differences in Expression between Strains

Overall, the set of genes showing constitutive difference in expression was not particularly enriched for genetic variants identified by our previous population genomic analysis. However, a few outlier SNPs were detected as located in flanking regions, sometimes annotated as promoters, or introns. These SNPs were part of the most differentiated loci detected in the whole-genome analysis we conducted previously (Haberkorn et al. 2023). Both flanking and intronic mutations can lead to transcriptomic changes, since they may serve as enhancers that ultimately recruit the RNA polymerase II on the promoter (Cannavò et al. 2016), thus being involved in the constitutive expression difference observed. In Nilaparvata lugens, the up-regulation of the P450 genes CYP6AY1 and CYP6ER1 were associated with resistance to imidacloprid together with several mutations in their respective promoters (Liang et al. 2018). Thus, we may speculate that one or several of the mutations we highlight here are involved in the transcriptomic difference observed.

The unique gene showing an interaction between the strain and treatment effect encoded a leucine-rich repeat soc-2 protein. Although soc-2 is not directly known to be involved in insecticide resistance, it seems to play a role in nicotinic acetylcholine receptor (nAChR) sensitivity, the target site of neonicotinoids (Gottschalk et al. 2005). High levels of neonicotinoid resistance have been observed in bed bug strains and attributed partially to enzymatic activities, without exploring the involvement of nAChR (Romero and Anderson 2016). Whilst we did not test the resistance to neonicotinoids in the present study, it is reasonable to think that the recently collected LF strain (2008) has more chances of having been exposed to neonicotinoids in the past rather than the ancient LL strain (<1980) that was collected before the deployment of neonicotinoids in the 1990s. Five SNPs showing high differentiation between LL and LF strains were detected in this gene, including 1 lying within a genomic region annotated as promoter. This polymorphism may thus be the genetic basis for the transcriptomic difference observed. More studies should be done to test the putative involvement of soc-2 in insecticide resistance phenotype.

Genes Induced upon Insecticide Exposure

Next, we aimed at identifying genes whose expression is modified upon insecticide exposure. In total, 375 genes were jointly induced in both strains. Among the 20 resistance genes up-regulated upon insecticide exposure, both ABC transporter and other detox categories of resistance genes were detected as significantly enriched (Table S2, Supplementary Material online). Almost all of the resistance genes overexpressed upon insecticide exposure are detoxification enzymes, showing a strong generalist plastic response in both strains. Detoxification genes have been shown to play a role in the defence of susceptible populations of insects against insecticides, with changes of expression upon insecticide exposure (Mastrantonio et al. 2017). Transcriptional plasticity of ABC transporters in particular has previously been identified to enable adaptation of efflux capacity in Tribolium castaneum to get rid of insecticides (Rösner and Merzendorfer 2020). This result also supports the hypothesis that resistance genes have been acquired over millions of years of arms race to counter toxins produced by plants in a distant stinging-sucking past, still helping to excrete insecticides nowadays (Després et al. 2007).

The protein-coding genes which are similarly transcriptionally induced after insecticide exposure in both strains are not expected to explain a difference in resistance phenotype between the strains, unless they encode different proteins. Among the multitude of genes that are similarly induced in both resistant and susceptible strains, we identified a few showing this pattern (Table S6, Supplementary Material online). Although none of these genes belonged to a resistance category, we considered that 1 deserved attention. A gene encoding a serpin B3 protein was associated with 3 nonsynonymous outlier SNPs, with the putative derived allele being fixed in the resistant strain. Additionally, this gene was associated in total with 28 nonsynonymous mutations (including nonoutlier SNPs) in our previous genomic analysis (Haberkorn et al. 2023), part of the top 1% of genes having the highest ratio between nonsynonymous and synonymous mutations (πn/πs=4.67), and highly differentiated between resistant and susceptible strains (FST=0.15, top 5%). Overall, this strongly suggests that this gene has experienced positive selection in the recent past, leading to an increased frequency of a putative derived allele. Although to our knowledge this gene has never been associated with pyrethroid resistance, the fact that it is also induced upon deltamethrin exposure renders its possible involvement in insecticide resistance plausible. The function of this particular gene is unclear, but serpin’s transcription has been shown to correlate with insecticide resistance in C. pipiens mosquitoes (Li et al. 2016).

No Evidence for Differentially Expressed Genes to Be Clustered in a Superlocus

DE genes, both between strains and between untreated and treated conditions, were dispersed across putative chromosomes (LG) and scaffolds of the bed bug assembly. More generally, there was no difference in average genetic differentiation between DE or non-DE genes. This is in sharp contrast with the high level of clustering of candidate SNPs observed in the bed bug genome, where a highly differentiated 6 Mb genomic region was identified (Haberkorn et al. 2023).

We consider 3 alternative explanations for this discrepancy: (i) most of the genomic regions identified previously are not associated with insecticide resistance, (ii) transcription factors located inside the superlocus might drive expression differences elsewhere in the genome, (iii) the resistance difference between the strains is mediated by quantitative or qualitative changes in messenger RNAs that were not measured in this study, or (iv) the resistance difference between the strains is not mediated by changes in mRNAs but rather by other mechanisms.

In support of hypothesis (i), although the 2 strains are highly comparable (FST=0.018), we cannot exclude the possibility that part of the genetic differentiation observed is unrelated to insecticide resistance, i.e. related to other environmental factors, or even unrelated to adaptation, but rather due to genetic drift.

However, as mentioned in hypothesis (ii), distant regulatory elements can also drive gene expression (West and Fraser 2005). We can, therefore, hypothesize that modifications in the sequences of transcription factors located in the superlocus detected could affect the regulation of genes detected in the present study, elsewhere in the genome.

In support of hypothesis (iii), we acknowledge that insecticide resistance may be triggered by between-strains differences in mRNA concentrations that are not quantified in this work. First, our transcriptomic analysis is focusing on adult insects, and any transcriptomic change relevant for insecticide resistance that would occur before adulthood would obviously not be detected. Indeed, previous studies have shown that the transcriptomic induction of some detoxification genes such as ABC transporter or P450 may confer insecticide resistance in the larval stages (Stevens et al. 2000; Wang et al. 2023). Exploring the transcriptomes of different stages of resistant bed bugs could be very informative on the timeline of resistance. Additionally, genes targeted by insecticide resistance such as VGSC may present nonsynonymous changes underlying resistance, which are not accompanied by transcriptomic changes. They may be detected in the genomic differentiation analysis (this is the case for the VGSC locus in Haberkorn et al. 2023) but surely not by the present transcriptomic analysis.

Finally, in support of (iv), a growing recent literature has brought to light examples of insecticide resistance mediated by noncoding RNAs that were not quantified in the present work. A recent study demonstrated that a long noncoding RNA (lncRNA) modulated the expression of a GST detoxification enzyme by competing with micro RNA (miRNA) binding, thus mediating cyflumetofen resistance (Feng et al. 2020). In Aphis gossypii, RT-qPCR and RNAi studies assessed the role of lncRNAs in acetyl-CoA carboxylase regulation, the target site of spirotetramat insecticide (Peng et al. 2021). miRNAs are another type of ncRNAs that can affect insecticide resistance by themselves, as shown in Spodoptera frugiperda, where injection of mi-RNAs-190-5p antagomir enhanced the expression of CYP6K2, improving insecticide tolerance (Yang et al. 2022). Hence, future studies should quantify the contribution of these mechanisms to the overall resistance pattern observed in C. lectularius.

To conclude, we identified a small set of genes (20 overexpressed in resistant strain, and 1 in interaction between resistance status and insecticide exposure) representing very strong candidates for insecticide resistance that should be functionally tested in the future.

Materials and Methods

Insects

The 2 strains used in this study were provided by Cimex Store Ltd (Chepstow, United Kingdom). The susceptible strain, LL (collected in London, Great Britain), was collected before the massive use of insecticide and raised in laboratories for more than 40 years. On the contrary, LF (collected in 2008 in London, in Great Britain) was moderately resistant to pyrethroids.

Insects were kept isolated before imaginal molt in 24-wells Petri dishes containing accordion-folded blotting papers, serving as harborage. Bed bugs were maintained at 25 ∘C, 40% relative humidity (RH), and a photoperiod of 12:12h. Since males perform the so-called traumatic insemination to copulate, which can be costly for their partners (Stutt and Siva-Jothy 2001), only 7-days-old unfed virgin adult females were used for both insecticide treatment and control.

Topical Assay

Insects were treated with a discriminating dose of 1 ng/μL of deltamethrin (1.98 μM), as assumed to reflect the LD50 of CimexStore resistant strains (Haberkorn et al. 2023). Insecticidal assays were carried out with deltamethrin (98% purity, Cluzeau, Sainte-Foy-La-Grande, France), a pyrethroid, as it remains the most used insecticide family recently for bed bug control. For each strain, 4 replicates of 8 insects were performed both for control and 1 ng/μL insecticide exposure (n=64). Bed bugs were previously immobilized by placing them in a Petri dish on ice for 5 min. Topical applications were then made onto the ventral surface of the thorax, between the coxae, with a 50 μL glass syringe attached to a repeating dispenser (Hamilton Co., Reno, NV). Treated insects were exposed to 1 μL of insecticide powder diluted in acetone, whereas control insects received 1 μL of acetone only.

Mortality was assessed after 24 h by flipping each insect on dorsal side with a featherweight forceps, to see if it was able to reverse on the ventral side (alive) or if its movements were not coordinated enough to do so (moribund, considered as dead). A fisher test was used to assess the significance of mortality difference between pairs of strains.

RNA Extraction and Sequencing

For each replicate, 3 individuals out of the survivors among the 8 individuals per replicate were randomly sampled for each condition (4 replicates for LL insecticide-exposed, LL control, LF insecticide-exposed, and LF control, n=48), immediately flash-frozen in liquid nitrogen and stored at −80∘C. Bed bugs were individually extracted using the RNAeasy Mini Kit (Qiagen, Hilden, Germany) with Turbo DNAse treatment on column (Thermofisher, Waltham MASS, USA). RNA concentration was measured using Nanodrop, and Qubit with RNA HS Kit (Agilent, Santa Clara, CA, USA).

Both reverse-transcription and sequencing were performed by Macrogen Europe (Amsterdam, Netherlands). Individual libraries (n=48) were prepared using TruSeq stranded mRNA kit (Illumina, San Diego, CA, USA), and sequencing was performed on an Illumina Novaseq6000 machine producing 40 million reads 2*100PE. The raw data have been submitted to the sequence read archive (SRA) database of NCBI under BioProject PRJNA832557.

Transcripts Quantification

The whole pipeline with the detail of parameters used is available on GitHub (https://github.com/chaberko-lbbe/clec-rnaseq). Adapters have already been removed by Macrogen, and no trimming was performed considering the very high quality of the reads (Q30>92%), assessed using Fastqc (Andrews 2015). Reads were then mapped on the bed bug genome using the STAR v2.7.3a workflow (Dobin et al. 2013). We used the Clec_2.1 and associated GTF annotation file (GCF_000648675.2). Parameters used for the mapping were as follow: “sjdbOverhang” of 100 (read size—1), a computed optimal “genomeSAindexNbases” of 13, “genomeChrBinNbits” of 18 and “quantMode” on “GeneCounts.” This last option allowed STAR to count the number of reads per gene while mapping, using the GTF annotation provided above (file “ReadsPerGene.out.tab”). For all 48 samples sequenced (24 samples per strain), the proportion of uniquely mapped reads was high (on average 84,09%, minimum 75,56%). Multiple mapping reads were discarded when over 10 occurrences (default parameter).

Mapping Scaffolds on Linkage Groups

Fountain et al. (2016) constructed a genetic map for C. lectularius, using RAD-seq markers on an F2 recombined population. However, as the genome assembly of C. lectularius has since been updated (v1.0 to 2.1), we used Fountain’s RAD-seq data and R scripts to generate a genetic map based on the latest genome assembly, as described in Haberkorn et al. (2023). This allowed us to obtain 14 putative autosomes, called LGs.

Differential Expression Analysis on London Strains

Reads counts by gene (except tRNA) were analyzed using the R package DESeq2 v1.34.0 (Love et al. 2014) with “lab” and “untreated” as the reference levels. A model with interaction between conditions treatment and strain was built. Only genes with at least 10 reads (summed over all samples) were kept, leading to a total of 12,187 genes out of 13,208 (total number of genes without tRNA).

Genes were considered as significantly DE when the adjusted P-value was below 5% (Padj<0.05) and when the fold change in expression was above 1.5, i.e. (LFC)>0.58 (up) or <−0.58 (down), as in Nardini et al. (2012).

Functional Analysis

In order to test whether some relevant biological functions were enriched in the set of candidate genes, we defined 10 categories of genes potentially involved in resistance based on the reference genome annotation, up to a total of 431 genes (as in Haberkorn et al. 2023). Six of these categories corresponded to metabolic resistance : “Binding/Sequestration” (6 genes) , “GST” (14), “CCE” (39), “UDPGT” (39), “ABC transporter/MRP” (55), and “P450” (56). Other categories were separated as: “Other detox” (64), “Cuticle” (113), “Insecticide target and nervous system” (31), and “Redox homeostasis” (14 genes) complete list with annotations is available in Table S7, Supplementary Material online, and all categories are further detailed in Haberkorn et al. 2023). Analysis of enrichment of candidate genes in each resistance category were performed by comparing their proportions in genomic regions of interest with the whole genome, using the function prop.test in R/stats default package (χ2 with Yates correction for continuity, with FDR correction for multiple comparisons).

Associating Differentially Expressed Genes with Outliers SNPs and Structural Variants

Genomic data were obtained from our previous analysis, using the whole-genome pool-seq data deposited on the SRA database of NCBI under BioProject PRJNA826750. Briefly, reads were mapped on the C. lectularius reference genome (Clec_2.1 assembly, Harlan strain), SNPs were then called, and differentiation indexes (FST) were computed for each SNP. The whole pipeline is available on GitHub (https://github.com/chaberko-lbbe/clec-poolseq). Outlier SNPs were detected based on the intersection of 3 criteria: (i) the differentiation between the 2 strains originating from London should be high (top 5% FST), (ii) alleles should be in a derived state, if one considers that the alleles of the reference genome (the insecticide susceptible Harlan strain) are ancestral (at least for loci involved in insecticide resistance), (iii) (derived) alleles should be in higher frequency in LF than in LL. With these criteria applied using only the 2 London strains, 100,941 outlier SNPs were detected (versus 576 in Haberkorn et al. 2023 that included another set of strains). Among these SNPs, 59,537 were distributed in 4,987 genes (out of the 12,187).

Using a combination of read-pair orientation, insert size information and read depth, the genomic analysis also revealed 650 inverted duplications, 509 tandem duplications, and 1,338 inversions in LF compared to the reference genome, with higher frequencies in LF compared to LL. These putative structural variants were overlapping with 4,118 genes (Haberkorn et al. 2023).

Analyses of enrichment in DE genes across LG or scaffolds were performed using R/stats default package (binomial test with Bonferroni correction for multiple testing). We also checked whether DE genes had higher values of genetic differentiation (computed as the mean value of FST per gene, i.e. the mean value of all SNPs’ FST detected in a gene), and the fold change value in constitutive DE genes whenever containing a structural variant or not, by performing t-tests with alternative “greater.”

Preparation of Cuticular Samples for Scanning Electron Microscope

Cuticular resistance, mentioned above, could be due to an increased thickness of the bed bug cuticle (Balabanidou et al. 2018). We measured this trait on 10 control specimens (exposed with acetone only) of each strain picked up from the RNA-seq experiment. These samples were preserved dry at −80∘C. The left middle leg of each specimen was separated from the body at the apical region of the femur in fresh PBS buffer. The next steps were performed by the “Centre Technologique des Microstructures” (CTMU) at University Lyon 1. Briefly, fixation was first carried out for 2 h with 2% glutaraldehyde in 0.1 M sodium cacodylate buffer. Samples were then rinsed 3 times for 15 min in 0.2 M sodium cacodylate buffer. Postfixation was performed in 1% osmium tetroxide in 0.1 M sodium cacodylate buffer. After a rapid rinsing with distilled water, gradual dehydration were made onto increasing baths from 30% to 100% alcohol (30, 50, 70, 80, 95, 100) for 10 min each, followed by 2 15 min baths in propylene oxide. Progressive impregnation with EPON resin were then performed, followed by polymerization at 56 ∘C (i.e. resin hardening in an oven). Blocks were then surfaced on the middle part of each tibia at room temperature with a diamond knife, and copper-plated with a 10 nm deposit.

Electron Microscope Imaging and Image Analysis

Observations were made by CTMU laboratory with a Zeiss Merlin VP Compact SEM at 3 kV. Each sample was individually observed and a picture was captured with a working distance of 5 mm, scan speed of 9, aperture size of 60 μm, and magnification of 729X.

We then analyzed raw images using Fiji in ImageJ2 (version 2.9.0/1.53t) (Rueden et al. 2017). The scale given in microns on the raw image was converted into a number of pixels to measure the cuticle thickness. A circle pattern with 12 radii was affixed to each image, enabling measurements to be taken from the point where the circle intersects the outside of the cuticle to the inside of the cuticle perpendicularly, whilst avoiding obvious debris, damage, or setae as made by Lilly et al. (2016b). Hence, 12 measurements were made for each cuticular cross-section sample. The experimenter was blind to the resistance status of the sample to avoid any bias. Cuticle thickness was compared between LL and LF strains in a repeated measures ANOVA framework.

Supplementary Material

evae158_Supplementary_Data

Acknowledgments

Bioinformatics analyses were performed using the computing facilities of the CC LBBE/PRABI. We thank the Symbiotron platform from the FR3728 BioEEnViS, especially Angelo Jacquet, and the Equipex+ InfectioTron [ANR-21-ESRE-0023] for facilities and equipment for rearing and experimentation on bed bugs. Electron microscopy studies have been done at the “Centre Technologique des Microstructures”—Claude Bernard University of Lyon, with special thanks to Lucie Geay. We also thank Dr Rike Stelkens (Stockholm University) for her gracious proofreading and valuable insights.

Supplementary Material

Supplementary material is available at Genome Biology and Evolution online.

Funding

This work was supported by a CIFRE grant (2019/0800), funded by Association Nationale de la Recherche et de la Technologie, from the collaboration between the Laboratoire de Biométrie et Biologie Évolutive and IZInovation SAS; and the Scientific Breakthrough Project Micro-be-have (Microbial impact on insect behavior) of Universite de Lyon, within the program “Investissements d’Avenir” (ANR-11-IDEX-0007 and ANR-16-IDEX-0005).

Data Availability

The data that support the findings of this study are openly available on SRA under BioProject PRJNA832557.
==== Refs
References

Adelman ZN , KilcullenKA, KoganemaruR, AndersonMAE, AndersonTD, MillerDM, HansenIA. Deep sequencing of pyrethroid-resistant bed bugs reveals multiple mechanisms of resistance within a single population. PLoS One. 2011:6 (10 ):1–9. 10.1371/journal.pone.0026228.
Ahmad M , DenholmI, BromilowRH. Delayed cuticular penetration and enhanced metabolism of deltamethrin in pyrethroid-resistant strains of Helicoverpa armigera from China and Pakistan. Pest Manag Sci. 2006:62 (9 ):805–810. 10.1002/ps.1225.16649192
Akhoundi M , KengneP, CannetA, BrenguesC, BerengerJ-M, IzriA, MartyP, SimardF, FontenilleD, DelaunayP. Spatial genetic structure and restricted gene flow in bed bugs (Cimex lectularius) populations in France. Infect Genet Evol. 2015:34 :236–243. 10.1016/j.meegid.2015.06.028.26140960
Andrews S . 2015,FastQC: A Quality Control Tool for High Throughput Sequence Data [Online], URL http://www.bioinformatics.babraham.ac.uk/projects/fastqc/.
Bai X , MamidalaP, RajarapuSP, JonesSC, MittapalliO. Transcriptomics of the bed bug (Cimex lectularius). PLoS One. 2011:6 (1 ):1–10. 10.1371/journal.pone.0016336.
Balabanidou V , GrigorakiL, VontasJ. Insect cuticle: a critical determinant of insecticide resistance. Curr Opin Insect Sci. 2018:27 :68–74. 10.1016/j.cois.2018.03.001.30025637
Balvín O , BoothW. Distribution and frequency of pyrethroid resistance-associated mutations in host lineages of the bed bug (Hemiptera: Cimicidae) across Europe. J Med Entomol. 2018:55 (4 ):923–928. 10.1093/jme/tjy023.29562293
Cannavò E , KhoueiryP, GarfieldDA, GeeleherP, ZichnerT, GustafsonEH, CiglarL, KorbelJO, FurlongEEM. Shadow enhancers are pervasive features of developmental regulatory networks. Curr Biol. 2016:26 (1 ):38–51. 10.1016/j.cub.2015.11.034.26687625
Dang K . Detection of knockdown resistance (KDR) in Cimex lectularius and Cimex hemipterus (Hemiptera: Cimicidae). In: Müller G, Pospischil R, Robinson WH, editors. The 8th International Conference on Urban Pests in Zurich, Switzerland. 71 (February 2018); 2014. p. 914–922.
Dang K , ToiCS, LillyDG, BuW, DoggettSL. Detection of knockdown resistance mutations in the common bed bug, Cimex lectularius (Hemiptera: Cimicidae), in Australia. Pest Manag Sci. 2015:71 (7 ):914–922. 10.1002/ps.3861.25046700
Davies TGE , FieldLM, WilliamsonMS. The re-emergence of the bed bug as a nuisance pest: implications of resistance to the pyrethroid insecticides. Med Vet Entomol. 2012:26 (3 ):241–254. 10.1111/j.1365-2915.2011.01006.x.22235873
Després L , DavidJP, GalletC. The evolutionary ecology of insect resistance to plant chemicals. Trends Ecol Evol. 2007:22 (6 ):298–307. 10.1016/j.tree.2007.02.010.17324485
Dobin A , DavisCA, SchlesingerF, DrenkowJ, ZaleskiC, JhaS, BatutP, ChaissonM, GingerasTR. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013:29 (1 ):15–21. 10.1093/bioinformatics/bts635.23104886
Doggett SL , GearyMJ, RussellRC. The resurgence of bed bugs in Australia: with notes on their ecology and control. Environ Health. 2004:4 (2 ):30–38.
Durand R , CannetA, BerdjaneZ, BruelC, HaouchineD, DelaunayP, IzriA. Infestation by pyrethroids resistant bed bugs in the suburb of Paris, France. Parasite. 2012:19 (4 ):381–387. 10.1051/parasite/2012194381.23193523
Dushay MS , EldonED. Drosophila immune responses as models for human immunity. Am J Hum Genet. 1998:62 (1 ):10–14. 10.1086/301694.9443887
Feng K , LiuJ, WeiP, OuS, WenX, ShenG, XuZ, XuQ, HeL. lincRNA_Tc13743.2-miR-133-5p-TcGSTm02 regulation pathway mediates cyflumetofen resistance in Tetranychus cinnabarinus. Insect Biochem Mol Biol. 2020:123 :103413. 10.1016/j.ibmb.2020.103413.32534987
Fountain T , RavinetM, NaylorR, ReinhardtK, ButlinRK. A linkage map and QTL analysis for pyrethroid resistance in the bed bug Cimex lectularius. G3. 2016:6 (12 ):4059–4066. 10.1534/g3.116.033092.27733453
Gao Y , KimK, KwonDH, JeongIH, ClarkJM, LeeSH. Transcriptome-based identification and characterization of genes commonly responding to five different insecticides in the diamondback moth, Plutella xylostella. Pestic Biochem Physiol. 2018:144 :1–9. 10.1016/j.pestbp.2017.11.007.29463402
Goddard J , DeshazoR. Bed bugs (Cimex lectularius) and clinical consequences of their bites. J Am Med Assoc. 2009:301 (13 ):1358–1366. 10.1001/jama.2009.405.
Gottschalk A , AlmedomRB, SchedletzkyT, AndersonSD, YatesJR, SchaferWR. Identification and characterization of novel nicotinic receptor-associated proteins in Caenorhabditis elegans. EMBO J. 2005:24 (14 ):2566–2578. 10.1038/sj.emboj.7600741.15990870
Guedes RNC , WalseSS, ThroneJE. Sublethal exposure, insecticide resistance, and community stress. Curr Opin Insect Sci. 2017:21 (31 ):47–53. 10.1016/j.cois.2017.04.010.28822488
Haberkorn C , DavidJ-P, HenriHlne, DelpuechJ-M, LasseurR, VavreF, VaraldiJ. A major 6 Mb superlocus is involved in pyrethroid resistance in the common bed bug Cimex lectularius. Evol Appl. 2023:16 (5 ):1012–1028. 10.1111/eva.13550.37216030
Koganemaru R , MillerDM, AdelmanZN. Robust cuticular penetration resistance in the common bed bug (Cimex lectularius L.) correlates with increased steady-state transcript levels of CPR-type cuticle protein genes. Pestic Biochem Physiol. 2013:106 (3 ):190–197. 10.1016/j.pestbp.2013.01.001.
Li C-x , GuoX-x, ZhangY-m, DongY-d, XingD, YanT, WangG, ZhangH-d, ZhaoT-y. Identification of genes involved in pyrethroid-, propoxur-, and dichlorvos- insecticides resistance in the mosquitoes, Culex pipiens complex (Diptera: Culicidae). Acta Trop. 2016:157 :84–95. 10.1016/j.actatropica.2016.01.019.26802491
Li X , SchulerMA, BerenbaumMR. Molecular mechanisms of metabolic resistance to synthetic and natural xenobiotics. Annu Rev Entomol. 2007:52 (1 ):231–253. 10.1146/annurev.ento.51.110104.151104.16925478
Liang Z-K , PangR, DongY, SunZ-X, LingY, ZhangW-Q. Identification of SNPs involved in regulating a novel alternative transcript of P450 CYP6ER1 in the brown planthopper. Insect Sci. 2018:25 (5 ):726–738. 10.1111/1744-7917.12472.28459131
Lilly DG , DangK, WebbCE, DoggettSL. Evidence for metabolic pyrethroid resistance in the common bed bug (Hemiptera: Cimicidae). J Econ Entomol. 2016a:109 (3 ):1364–1368. 10.1093/jee/tow041.27018436
Lilly DG , LathamSL, WebbCE, DoggettSL. Cuticle thickening in a pyrethroid-resistant strain of the common bed bug, Cimex lectularius L. (Hemiptera: Cimicidae). PLoS One. 2016b:11 (4 ):6–16. 10.1371/journal.pone.0153302.
Love MI , HuberW, AndersS. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014:15 (12 ):550. 10.1186/s13059-014-0550-8.25516281
Mamidala P , JonesSC, MittapalliO. Metabolic resistance in bed bugs. Insects. 2011:2 (1 ):36–48. 10.3390/insects2010036.26467498
Mamidala P , WijeratneAJ, WijeratneS, KornackerK, SudhamallaB, Rivera-VegaLJ, HoelmerA, MeuliaT, JonesSC, MittapalliO. RNA-Seq and molecular docking reveal multi-level pesticide resistance in the bed bug. BMC Genomics. 2012:13 (1 ):6. 10.1186/1471-2164-13-6.22226239
Mastrantonio V , FerrariM, EpisS, NegriA, ScuccimarraG, MontagnaM, FaviaG, PorrettaD, UrbanelliS, BandiC. Gene expression modulation of ABC transporter genes in response to permethrin in adults of the mosquito malaria vector Anopheles stephensi. Acta Trop. 2017:171 :37–43. 10.1016/j.actatropica.2017.03.012.28302529
Nardini L , ChristianRN, CoetzerN, RansonH, CoetzeeM, KoekemoerLL. Detoxification enzymes associated with insecticide resistance in laboratory strains of Anopheles arabiensis of different geographic origin. Parasit Vectors. 2012:5 (1 ):113. 10.1186/1756-3305-5-113.22676389
O’Donel Alexander J . Infestation by Hemiptera. In: Arthropods and human skin. London: Springer London; 1984. p. 57–74.
Pan C , ZhouY, MoJ. The clone of laccase gene and its potential function in cuticular penetration resistance of Culex pipiens pallens to fenvalerate. Pestic Biochem Physiol. 2009:93 (3 ):105–111. 10.1016/j.pestbp.2008.12.003.
Peng T , PanY, TianF, XuH, YangF, ChenX, GaoX, LiJ, WangH, ShangQ. Identification and the potential roles of long non-coding RNAs in regulating acetyl-CoA carboxylase ACC transcription in spirotetramat-resistant Aphis gossypii. Pestic Biochem Physiol. 2021:179 :104972. 10.1016/j.pestbp.2021.104972.34802522
Potter MF . The history of bed bug management- with lessons from the past. Am Entomol. 2011:57 (1 ):14–25. 10.1093/ae/57.1.14.
Poupardin R , ReynaudS, StrodeC, RansonH, VontasJ, DavidJ-P. Cross-induction of detoxification genes by environmental xenobiotics and insecticides in the mosquito Aedes aegypti: impact on larval tolerance to chemical insecticides. Insect Biochem Mol Biol. 2008:38 (5 ):540–551. 10.1016/j.ibmb.2008.01.004.18405832
Romero A , AndersonTD. High levels of resistance in the common bed bug, Cimex lectularius (Hemiptera: Cimicidae), to neonicotinoid insecticides. J Med Entomol. 2016:53 (3 ):727–731. 10.1093/jme/tjv253.26823499
Romero A , PotterMF, HaynesKF. Evaluation of piperonyl butoxide as a deltamethrin synergist for pyrethroid-resistant bed bugs. J Econ Entomol. 2009:102 (6 ):2310–2315. 10.1603/029.102.0637.20069862
Romero A , PotterMF, PotterDA, HaynesKF. Insecticide resistance in the bed bug: a factor in the Pest’s Sudden Resurgence? J Med Entomol. 2007:44 (2 ):175–178. 10.1603/0022-2585(2007)44[175:iritbb]2.0.co;2.17427684
Rösner J , MerzendorferH. Transcriptional plasticity of different ABC transporter genes from Tribolium castaneum contributes to diflubenzuron resistance. Insect Biochem Mol Biol. 2020:116 :103282. 10.1016/j.ibmb.2019.103282.31740345
Rueden CT , SchindelinJ, HinerMC, DeZoniaBE, WalterAE, ArenaET, EliceiriKW. ImageJ2: ImageJ for the next generation of scientific image data. BMC Bioinformatics. 2017:18 (1 ). 10.1186/s12859-017-1934-z.
Stevens JL , SnyderMJ, KoenerJF, FeyereisenR. Inducible P450s of the CYP9 family from larval Manduca sexta midgut. Insect Biochem Mol Biol. 2000:30 (7 ):559–568. 10.1016/s0965-1748(00)00024-2.10844248
Stutt AD , Siva-JothyMT. Traumatic insemination and sexual conflict in the bed bug Cimex lectularius. Proc Natl Acad Sci USA. 2001:98 (10 ):5683–5687. 10.1073/pnas.101440698.11331783
Wang L , TianS-H, ZhaoW, WangJ-J, WeiD-D. Overexpression of ABCB transporter genes confer multiple insecticide tolerances in Bactrocera dorsalis (Hendel) (Diptera: Tephritidae). Pestic Biochem Physiol. 2023:197 :105690. 10.1016/j.pestbp.2023.105690.38072545
West AG , FraserP. Remote control of gene transcription. Hum Mol Genet. 2005:14 (suppl_1 ):R101–R111. 10.1093/hmg/ddi104.15809261
Yahouédo GA , ChandreF, RossignolM, GinibreC, BalabanidouV, MendezNGA, PigeonO, VontasJ, CornelieS. Contributions of cuticle permeability and enzyme detoxification to pyrethroid resistance in the major malaria vector Anopheles gambiae. Sci Rep. 2017:7 (1 ):11091. 10.1038/s41598-017-11357-z.28894186
Yang Y , ZhangY, WangA, DuanA, XueC, WangK, ZhaoM, ZhangJ. Four MicroRNAs, miR-13b-3p, miR-278-5p, miR-10483-5p, and miR-10485-5p, mediate insecticide tolerance in Spodoptera frugiperda. Front Genet. 2022:12 (January ):1–11. 10.3389/fgene.2021.820778.
Zhang C , GuoX, LiT, ChengP, GongM. New insights into cypermethrin insecticide resistance mechanisms of Culex pipiens pallens by proteome analysis. Pest Manag Sci. 2022:78 (11 ):4579–4588. 10.1002/ps.v78.11.35837767
Zhu F , GujarH, GordonJR, HaynesKF, PotterMF, PalliSR. Bed bugs evolved unique adaptive strategy to resist pyrethroid insecticides. Sci Rep. 2013:3 :1–8. 10.1038/srep01456.
Zhu F , WiggintonJ, RomeroA, MooreA, FergusonK, PalliR, PotterMF, HaynesKF, PalliSR. Widespread distribution of knockdown resistance mutations in the bed bug, Cimex lectularius (Hemiptera: Cimicidae), populations in the United States. Arch Insect Biochem Physiol. 2010:73 (4 ):245–257. 10.1002/arch.20355.20301216
