Comparative Transcriptomic Analysis of Insecticide Resistant Aedes aegypti from Puerto Rico Reveals Insecticide-Specific Patterns of Gene Expression

1. Introduction

Chemical insecticides are widely used to control insects of public health significance. However, the intensive and widespread use of these insecticides has led to rapid and widespread evolution of insecticide resistance in major insect disease vectors, raising concerns about their sustained use in public health. Developing strategies to maintain efficacy of the limited number of chemical insecticides approved for public health use in the face of growing insecticide resistance is a major priority for the control of malaria. This effort could benefit from an in-depth understanding of the molecular basis of resistance to key insecticides used in public health vector control.

Mosquitoes are among the major disease vectors that have become increasingly resistant to the principal classes of chemical insecticides approved for vector control (Ranson and Lissenden 2016). Target site insensitivity and increased metabolic detoxification are the main mechanisms responsible for mosquito adaptation to insecticides, but other mechanisms such as behavioral avoidance of insecticide-treated surfaces and cuticular alterations have also been reported (Hemingway et al. 2004, Chareonviriyaphap et al. 2013, Liu 2015, Simma et al. 2019). Target site insensitivity occurs when a mutation of the protein targeted by the insecticide prevents the insecticide from binding to the molecule. These mutations are well characterized in mosquitoes and molecular tests for their detection are available for many major vectors (Weill et al. 2004, Bass et al. 2007, Mathias et al. 2011). In contrast, metabolic resistance is much less understood and is primarily conferred via gene amplifications, transcriptional enhancements, and coding mutations in detoxification enzymes, particularly cytochrome P450 monooxygenases, carboxylesterases and glutathione-S-transferases (Li et al. 2007). Because metabolic resistance can arise through multiple mechanisms and involves many gene families, understanding its molecular basis remains a complex public health challenge. Moreover, the specific detoxification genes involved in metabolic resistance can vary greatly between insecticides and across mosquito populations and species (David et al. 2013, Ishak et al. 2017, Riveron et al. 2017). These complexities hamper the development of universal molecular diagnostic markers for metabolic resistance and highlight the need for site-specific understanding of the molecular mechanisms responsible for insecticide resistance to facilitate development of evidence-based resistance management strategies.

Here, we use RNA-Seq to examine the gene expression profiles of malathion, alpha-cypermethrin, and lambda-cyhalothrin resistant populations of Aedes aegypti from two locations in Puerto Rico. This mosquito species lives near humans in urban settings and is the primary vector of several globally important arthropod-borne viruses (arboviruses) including dengue, chikungunya, Zika and yellow fever. Resistance of Ae. aegypti to multiple insecticide classes has been reported worldwide posing serious operational challenges towards its control (Braga et al. 2004, Vontas et al. 2012, Amelia-Yap et al. 2018, Schluep and Buckner 2021). Both target site insensitivity and metabolic resistance have been identified as important mechanisms of insecticide resistance in Ae. aegypti and a few previous studies have investigated the specific genes involved in metabolic detoxification of insecticides in this mosquito species (Strode et al. 2008, Lertkiatmongkol et al. 2010, Lima et al. 2011, David et al. 2014, Faucon et al. 2015, Aponte et al. 2019, Badolo et al. 2019, Lien et al. 2019, Saavedra-Rodriguez et al. 2021, Sene et al. 2021, Samal et al. 2022). This study seeks to disentangle the transcriptional basis of resistance to multiple insecticides in Ae. aegypti from two geographically distinct populations, with the aim of better understanding the specificity of metabolic resistance with respect to insecticide and geography.

2. Materials and Methods

2.1. Study area and mosquito collections

Three populations of Ae. aegypti were used in this study: the Rockefeller (ROCK) reference susceptible strain (obtained from the Malaria Research and Reference Reagent Resource Center, MR4) and two Ae. aegypti populations collected from the towns of Isabela (18.4704° N, 67.0242° W) and Manatí (18.4330° N, 66.4759° W) in Puerto Rico in spring of 2017. Eggs were collected from the adults and shipped to the U.S. Centers for Disease Control and Prevention (CDC) in Atlanta, USA. At CDC, the mosquitoes were reared at 27 ± 2°C and 70 ± 10% relative humidity, and a 14L:10D photoperiod. Adults were provided 10% sucrose ad libitum. All three strains were reared in parallel for each of the biological replicates for the RNA-Seq libraries.

2.2. Insecticide bioassays

Three-to-four-day old non-blood fed and unmated Ae. aegypti females from Isabela (ISA) and Manatí (MAN) were tested for resistance to the diagnostic doses of alpha-cypermethrin (10 μg/bottle), malathion (50 μg/bottle), and lambda-cyhalothrin (10 μg/bottle) using the CDC bottle bioassay (Brogdon and Chan 2010). Assays were carried out in four replicates, each containing approximately 25 individuals per bottle. For each test, a control bottle coated with only 1 ml of acetone was included. Mosquitoes alive after 30 minutes of insecticide exposure were considered resistant. Assays were repeated to obtain sufficient numbers of phenotyped mosquitoes for RNA extraction. For the susceptible ROCK strain, three-to-four-day-old adult non-blood fed and unmated female mosquitoes were killed by freezing and stored at -80˚C until RNA extraction was conducted. For ISA and MAN, mosquitoes that were alive (survivors) at the end of the 30-minute exposure period were separated and tallied, then frozen at -80 °C. In addition, individuals from ISA and MAN that were used in the control bottle and not exposed to insecticide were killed by freezing and stored at -80˚C for RNA extraction and sequencing.

2.3. RNA and library preparation methods

Three biological replicates with pools of 5 mosquitoes each were prepared for three groups of mosquitoes: the susceptible Ae. aegypti Rockefeller strain (ROCK), field-collected mosquitoes that were not exposed to insecticides during the bottle bioassay, and field-collected mosquitoes that survived exposure to insecticides at the diagnostic dose and diagnostic time during the bottle bioassay. RNA was extracted using the Applied Biosystems Arcturus PicoPure RNA isolation kit (Arcturus, Applied Biosystems, USA) and quantified using the Agilent RNA ScreenTape 4200 assay, according to the manufacturer’s protocols. Two micrograms of starting material were treated with Baseline-ZERO™ DNase (Lucigen, USA) and ribosomal RNA was removed using the Ribo-Zero™ Magnetic Core Kit and Ribo-Zero™ rRNA Removal kit (Illumina, USA), according to the manufacturer’s protocols. RNA-Seq libraries were prepared from each pool of extracted RNA using the ScriptSeq™ v2 RNA-Seq library preparation kit (Illumina, USA), using 12 cycles of PCR amplification, according to the manufacturer’s protocol. Libraries were purified using Agencourt AMPure XP beads (Beckman Coulter, USA) and assessed for quantity and size distribution using the Agilent DNA ScreenTape D5000 assay. Each library was barcoded, and equal amounts of each library were pooled and sequenced (2x125bp paired end) on an Illumina HiSeq 2500 sequencer, using v2 chemistry. Sequencing was done at the Biotechnology Core Facility at CDC, Atlanta, USA.

2.4. Measurement of gene expression by quantitative PCR

Five genes that were significantly differentially expressed in the RNA-Seq analysis were validated using quantitative PCR (qPCR). One microgram of RNA from 3 replicates of samples resistant to lambda-cyhalothrin from Manatí and Isabela, and samples resistant to alphacypermethrin from Isabela were used to synthesize cDNA using the High-Capacity cDNA reverse transcription kit (Applied Biosystems) with oligo-dT20 (NEB), according to the manufacturer’s instructions. RNA and cDNA synthesis was also prepared for the susceptible laboratory strain ROCK. There was insufficient material remaining from the unexposed samples from both locations hence this was not included in the qPCR validation. The primers used are listed in Table S1. Standard curves of Ct values for each gene were generated using a serial dilution of cDNA, allowing assessment of PCR efficiency. The qPCR amplification was carried out on a QuantStudio 6 Flex Real-Time PCR system (Applied Biosystems) using PowerUp SYBR Green Master Mix (Applied Biosystems). cDNA from each sample was used as a template in a three-step program as follows. Uracil-DNA glycosylase (UDG) activation at 50°C for 2 minutes, DNA polymerase activation at 95°C for 10 minutes, followed by 40 cycles of DNA denaturation for 15 seconds at 95°C, DNA annealing and extension for 1 minute at 60°C, and a last DNA extension step of 15 seconds at 95°C. The relative expression level and Fold Change (FC) of each target gene from resistant field samples relative to the susceptible lab strain were calculated using the 2−ΔΔCT method (Schmittgen et al 2008) incorporating the PCR efficiency. Housekeeping genes RPS3 (Smith et al 2018) and Actin (Mitra et al 2021) were used for normalization.

2.5. Data analysis

2.5.1. Quality control filtering and mapping

De-multiplexed Paired-End (PE) RNA-Seq reads for each sample were evaluated for quality using FastQC v0.11.5 (Andrews 2010). Subsequently, the PE reads were trimmed and filtered using fastp v0.21.0 (Chen et al. 2018) to remove adapter and low-quality reads according to the following criteria: minimum base quality score = 20, minimum length required = 25, polyG and poly tail trimming = True. The resulting trimmed and filtered read pairs (R1/R2) were aligned against the reference genome of Ae. aegypti (genome assembly version = AaegL5.1, GeneBank assembly identifier = GCA_002204515.1) (Giraldo-Calderón et al. 2015), using ‘subjunc’ v2.0.1, part of the subread aligner v2.0.1 (Liao et al. 2013), with default parameters. The resulting alignment was filtered to remove reads with low mapping quality (q < 10) and sorted successively using Samtools v1.10 (Li et al. 2009). Summary statistics for the QC filtering and sequencing alignments corresponding to malathion, alpha-cypermethrin, and lambda-cyhalothrin are reported in Table S2. Tags (a read pair or single, unpaired read) mapped to the sense orientation of the annotated Ae. aegypti gene set AaegL5.1 were quantified with ‘FeatureCounts’, as part of the subread-aligner package v2.0.1 (Liao et al. 2013), by using the following criteria: 1) count only read pairs that have both ends aligned; 2) count fragment instead of reads; 3) minimum number of overlaps required = 1; 4) feature type = exon; 5) attribute type = gene_id; and 6) strandness = sense. The FeatureCount analysis generated a tag count matrix table which was used as input to edgeR (Robinson et al. 2010) for differential gene expression analysis.

2.5.2. Differential Gene Expression Analysis

For each group, three pairwise comparisons of differential gene expression (DGE) analysis were conducted to identify genes associated with resistance to each insecticide: 1) resistant vs unexposed control (R-C), 2) resistant vs ROCK susceptible (R-S), and 3) unexposed control vs ROCK susceptible (C-S). The R-C comparison was used to detect transcription induced by insecticide exposure, while the R-S and C-S comparisons were used to distinguish genes that are constitutively associated with insecticide resistance as well as those induced by insecticide exposure.

The DGE analysis was conducted using the R package edgeR (Robinson et al. 2010). To remove the effect of noise and low-expressed genes, for each pairwise comparison, genes with a total tag count less than 50 across all libraries (any pairwise comparison) were filtered out before further analysis as previously suggested (Robinson et al. 2010, Chen et al. 2020). The function calcNormFactors was used to normalize tag counts among samples, by using the TMM (Trimmed Mean M-values) method. This approach normalizes the dataset for sequencing depth and library size, the two most important technical factors that influence DGE. The tag count was not normalized for gene length and GC content, as these values do not vary from sample to sample and would be expected to have little effect on differentially expressed genes (DEGs). The tag-wise and common dispersion of the read count distribution were estimated using the estimateDisp function from the edgeR package. The DEGs were selected after multiple testing using the decideTests function in the limma package (Ritchie et al., 2015). A critical value absolute log₂ fold-change |log2FC| = 1 (|FC| =2), and a (False Discovery Rate) FDR adjusted P value ≤ 0.01 were used to tag a gene as significantly differentially expressed.

2.5.3. Gene ontology annotation and functional enrichment analysis

The Ae. aegypti reference genome corresponding to genome assembly version ‘AaegL5.1’ and GeneBank assembly identifier ‘GCA_002204515.1’ (directly downloaded from Vectorbase, Release 51) was used. AaegL5.1 contains 14,718 protein coding genes, 4,704 ncRNA genes and 382 pseudogenes (corresponding to 19804 genes) (Matthews et al. 2018). However, functional annotation is provided for only 6,397 (36.68%), and gene ontology annotation is provided for only 10,910 (74.13%) genes in VectorBase. To help with the interpretation of our results, predicted genes of AaegL5.1 were functionally annotated using Blast2GO as follows. First, a local BLASTp (v2.9) search of the predicted protein coding sequences was conducted against the Arthropoda (taxid=6665) category of the nr protein NCBI database with maximum e-value cut-off 10^-3. Second, the protein sequences were searched against the InterPro database (Hunter et al. 2008), using InterProScan v5 (Jones et al. 2014). The Blastp and InterProScan outputs were simultaneously provided to Blast2GO command line as input, which map the RefSeq and InteProScan identifiers to the GO database as curated and updated in the Blast2GO database (July 2021). This analysis assigned putative descriptions to 13,302 (90.37%) protein coding genes and gene ontology annotation to 11,672 (79.30%), which could be considered a significant improvement. All annotation results (VectorBase and Blast2GO) are provided in Table S3. GO term enrichment analysis of up- and down-regulated genes was carried out using Goatools (Klopfenstein et al. 2018), based on the go-basic database (release 2021-02-01). The list of 11,672 blast2GO annotated genes of Ae. aegypti with their associated GO terms was used as the background reference set. The P-values used to evaluate significantly enriched GO terms were calculated based on Fisher’s exact test and corrected with the Benjamini-Hochberg multiple test correction method(Klopfenstein et al. 2018). Finally, we used an FDR adjusted P-value <0.05 to define statistically significantly overrepresented (enriched) and underrepresented (purified) GO terms associated with the list of DEGs.

3. Results

3.1. Bioassay results

The mosquitoes from both sites showed different patterns of resistance to the three insecticides tested. After 30 min of exposure to alpha-cypermethrin, lambda-cyhalothrin, or malathion the mortality rate of mosquitoes from Manatí was 46%, 85% and 56%, respectively, compared to 81%, 48%, and 25%, respectively, for mosquitoes from Isabela (Figure 1, Table S4).

3.2. RNA sequencing, quality control filtering and alignment rate

The Table S2 shows the statistics of the RNA-Seq sequencing results (before and after quality filtering), the alignment of the filtered reads to the reference genome and the read quantification. The Illumina sequencing generated approximatively 27.6-69.4 million of reads per sample. After adapter trimming and read quality filtering, an average of 41.52 ± 10.43 million reads per sample were retained for subsequent analysis. An average of 24.12 ± 9.43 million quality filtered PE reads per sample were mapped to the reference genome of Ae. aegypti, while an average of 7.2 ± 3.03 million tags (read pairs) per sample, representing 50% to 74% of the total alignment, were successfully assigned to the exonic regions of the gene set AaegL5.1 (Figure S1).

3.3. Differential gene expression analysis

3.3.1. Differential gene expression associated with malathion resistance

For ISA samples, 1132 genes (441 upregulated and 691 downregulated), 563 genes (256 upregulated and 307 downregulated), and 119 genes (33 up regulated and 86 down regulated) were significantly differentially expressed in R-S, C-S and R-C comparisons, respectively. For MAN samples, a total of 1,396 genes (739 upregulated and 657 downregulated), 653 genes (197 up regulated and 456 downregulated), and 215 genes (174 upregulated and 41) were significantly differentially expressed in R-S, C-S, and R-C comparisons, respectively (Table 1, Figure 2, Table S5).

The DEGs that overlapped 2 or more comparisons for each experiment (4 intersections in total) were extracted and inspected for important gene expression patterns (Figure S2). For the ISA samples, there were 69 DEGs (24 upregulated and 45 downregulated) that were shared between R-S and R-C groups (Figure S2A). The top 5 shared over-expressed genes between R-S and R-C groups included 2 odorant binding proteins (OBPs) (the duplicated OBP66), a kDa salivary protein, a zinc finger 512B-like, and an uncharacterized protein (Table S6). However, for MAN samples, there were 148 differentially expressed genes (126 upregulated and 22 downregulated) that overlapped between R-S and R-C groups (Figure S2B); the top 10 shared upregulated genes included 7 OBPs (including duplicated OBP66, and OBP67), serine protease snake-like, and two uncharacterized proteins (Table S6).

Focusing on R-S comparisons, we identified 753 genes (309 upregulated and 444 downregulated) that were differentially expressed in malathion survivors from both ISA and MAN relative to the susceptible ROCK strain (Figure 3A). The top 15 upregulated genes with retrievable annotations included 2 serine proteases, 3 uncharacterized proteins, 3 odorant binding proteins, a 39S ribosomal mitochondrial protein, a digestive organ expansion factor homolog, a kDa secreted -1, a nuclear valosin-containing -like, nucleolar dao-5, peritrophin-1-like, pre-mRNA-processing factor 39, and peritrophin-1-like protein coding genes (Table S7).

The detoxification genes annotated as cytochrome P450 monooxygenases were predominantly over-expressed, while the cuticular-related protein genes were predominantly downregulated in malathion resistant samples from both ISA and MAN (Figure 2A-B). 83% (20/24) and 82% (27/33) of differentially expressed CYP450s were upregulated in the R-S comparisons from ISA and MAN, respectively, while 78% (18/23) and 91% (39/33) of cuticular-related DEGs were significantly down regulated in R-S comparisons from ISA and MAN, respectively (Figure S3).

From the list of DEGs shared between the R-S groups from the two sites, we have identified a list of core detoxification genes that are associated with malathion resistance in Ae. aegypti in both populations. These core candidate detoxification genes include 22 CYP450s (19 upregulated: CYP6Z7, CYP6D5, CYP6BY1, CYP325Q1, CYP9J2, CYP6N9, CYP6Z6, CYP325V1, CYP6N12, CYP6BB2, CYP6M9, CYP9J23, CYP9J22, CYP9J9, CYP4G36 and four CYP9F2 orthologs; and 3 downregulated: CYP12F5, CYP325N2, and CYP6N15), 1 upregulated carboxylesterase and 3 glutathione-s-transferases (including 1 upregulated: GSTD6, and 2 downregulated: GSTS1 and GSTX1) (Figure 4, Table 2). Although ISA Ae. aegypti showed higher resistance to malathion than MAN (Figure 1), we did not detect any detoxification genes that were uniquely overexpressed in malathion-resistant ISA. However, several CYP450s (including CYP4D39, CYP28A5, CYP4J16, and CYP9J23) were overexpressed in malathion-resistant MAN, but not in malathion-resistant ISA. A detailed summary of all the DEGs shared between ISA and MAN is shown in Table S7.

3.3.2. Differential gene expression associated with alpha-cypermethrin resistance

For ISA, 1171 genes (642 upregulated and 529 downregulated), 478 genes (158 upregulated and 322 downregulated), and 433 genes (87 up regulated and 346 down regulated) were significantly differentially expressed in R-S, C-S and R-C comparisons, respectively. For MAN, a total of 1175 genes (605 upregulated and 570 downregulated), 811 genes (222 upregulated and 589 downregulated), and 602 genes (154 upregulated and 448 downregulated) were significantly differentially expressed in R-S, C-S, and R-C comparisons, respectively (Table 1, Table S5). The list of DEGs that overlapped 2 or more comparisons for each experiment (4 intersections) (Figure S2B–C), were manually extracted and are reported in Table S6. Unlike the malathion experiment, no DEGs were shared between the R-S and R-C groups.

Focusing on the R-S comparisons, a total of 801 differentially expressed genes (421 upregulated and 380 downregulated) were shared between the MAN and ISA alpha-cypermethrin survivors, relative to the susceptible strain (Figure 3B). The top 10 overlapping DEGs with retrievable annotation (FC=15.1-2101.6) included a serine protease, a peritrophin-1-like, a kDa secreted salivary protein, an uncharacterized protein LOC5569546 isoform, an uncharacterized protein LOC5569546 isoform X1, NPC intracellular cholesterol transporter, a nuclear valosin-containing -like, a V-type proton ATPase subunit d1, and 2 cytochrome P450s (CYP67 and CYP9F2) (Table S7).

Similar to the malathion analysis, the detoxification genes annotated as cytochrome P450 monooxygenases were predominantly over-expressed, while the cuticular protein genes were predominantly downregulated in both alpha-cypermethrin survivors (R-S) and unexposed groups (C-S) relative to the susceptible strain for both study sites (Figure 2C-D). In particular, 76% (22/29) and 74% (20/27) of differentially expressed CYP450s were upregulated in the R-S comparisons for ISA and MAN, respectively, while 97% (34/35) and 97% (33/34) of cuticular related DEGs were significantly down regulated in R-S comparisons for ISA and MAN, respectively (Figure 2, Figure S3). From the list of DEGs that overlapped the two R-S groups, a list of core detoxification genes associated with alpha-cypermethrin resistance was identified due to their consistent differential expression in both MAN and ISA survivors (Figure 4, Table 2). These candidate alpha-cypermethrin detoxification genes included 21 CYP450s, of which 15 were significantly upregulated (CYP6Z7, CYP9J2, CYP6N12, CYP6BB2, CYP9J22, CYP6N9, CYP6M9, CYP6D5, CYP325V1, CYP6Z6, CYP9J9, CYP4G36 and three CYP9F2 orthologs) and 6 were downregulated (including CYP12F5, CYP6Y3, CYP4AC1, CYP325N, CYP6N15, and CYP4H32) (Figure 4A). Additionally, we identified 4 GST genes (GSTS1, GSTD1, GSTD11 and GSTE4) that were downregulated in both ISA and MAN alpha-cypermethrin survivors as compared to the susceptible strain (Figure 4B).

3.3.3. Differential gene expression associated with lambda-cyhalothrin resistance

For ISA samples, 1261 genes (341 upregulated and 920 downregulated), 480 genes (158 upregulated and 322 downregulated), and 286 genes (30 upregulated and 256 downregulated) were significantly differentially expressed in R-S, C-S and R-C comparisons, respectively. For MAN samples, a total of 1413 genes (362 upregulated and 1051 downregulated), 811genes (222 up regulated and 589 downregulated), and 254 genes (67 upregulated and 197) were significantly differentially expressed in R-S, C-S, and R-C comparisons, respectively (Table 1, Table S5). The DEGs that overlapped 2 or more comparisons for each experiment were manually extracted and are shown in Table S6. Although 174 DEGs (8 upregulated and 166 downregulated) and 139 DEGs (12 upregulated and 127 downregulated) were identified in the intersection R-S/R-C, for ISA and MAN, respectively, no upregulated genes were shared between the resistant, susceptible and controls from each study site (Figure S2B-C). Among the over-expressed genes in the R-S/R-C intersection for MAN was a glutathione-s-transferase gene (GSTD6) with fold change expression of 10.56 and 2.73 in R-S and R-C, respectively. However, of the 8 upregulated genes that overlapped the R-S and C-S groups for ISA, none were previously reported to be associated with insecticide detoxification (Table S6).

Focusing on only the R-S groups, a total of 943 differently expressed genes (203 upregulated and 743 downregulated) were shared between MAN and ISA lambda-cyhalothrin survivors relative to the susceptible strain (Figure 3C). The top 10 shared upregulated genes (FC = 14.2-1901.7) included a serine protease, a peritrophin-1-like, a digestive organ expansion factor homolog, an uncharacterized protein, nucleolar dao-5, pre-mRNA-processing factor 39, kDa secreted salivary protein, nuclear valosin-containing -like, NPC intracellular cholesterol transporter 2 homolog a, and a CYP450 (CYP6Z7) (Table S7). Similar to the malathion and alpha-cypermethrin experiments, the detoxification genes annotated as cytochrome P450 monooxygenases were predominantly over-expressed, while the annotated cuticular protein genes were predominantly downregulated in both R-S and C-S groups (Figure 2). In particular, 69% (18/26) and 55% (12/22) of differentially expressed CYP450s were significantly overexpressed in the R-S comparisons from ISA and MAN, respectively, while 97% (34/35) and 97 % (33/34) of cuticular-related DEGs were significantly downregulated in R-S comparisons of ISA and MAN lambda-cyhalothrin survivors, respectively (Figure 2, Figure S3).

From the DEGs that overlapped the R-S comparisons of both study sites, several core insecticide detoxification genes were identified. These included 17 cytochrome P450s (10 upregulated: CYP6Z7, CYP9J2, CYP28A5, CYP6Z6, two CYP9F2 orthologs, CYP6M9, CYP6BB2, CYP307A1, CYP4C2-like; and 7 downregulated: CYP6P12, CYP6Y3, CYP4H29, CYP6S3, CYP6N15, CYP325N2, CYP4C1-like). Additionally, 5 GSTs (1 upregulated: GSTD6; and 4 downregulated: GSTX1, GSTE6, GSTD11 and GSTS1) were significantly differentially expressed in mosquitoes that survived lambda-cyhalothrin exposure relative to the susceptible strain for both ISA and MAN (Figure 4, Table 2).

3.3.4. Genes associated with resistance to multiple insecticides

A total of 234 genes (43 upregulated and 191 downregulated) were significantly differentially expressed in mosquitoes that survived either malathion, alpha-cypermethrin or lambda cyhalothrin exposure compared to the susceptible laboratory strain (Figure 5A, Table S8). Two of the top 10 upregulated DEGs that overlapped all the R-S groups included a putative serine protease (AAEL029072; FC=1197-6081) and a peritrophin-1-like (AAEL021372; FC=189-1078), reflecting their potential association with the multiple-insecticide resistance phenotypes (Figure 5B). Eight key genes coding for CYP450 enzymes including CYP6Z7, CYP28A5, CYP9J2, CYP6Z6, CYP6BB2, CYP6M9, and two CYP9F2 orthologs were upregulated across the six R-S comparisons. Two additional CYP450s (CYP325N2, CYP6N15) and one GST (GSTS1) were downregulated in all six resistant groups (Figure 5C). A total of 131, 325, and 240 genes were uniquely differentially expressed in malathion, alpha-cypermethrin, lambda cyhalothrin survivors, respectively, suggesting their association with a specific insecticide (Figure 5A). Two detoxification DEGs, a CYP450 (CYP6BY1) and a carboxylesterase (AAEL004022), were unique to malathion resistant samples (Figure 4A). Additionally, one downregulated detoxification gene (GSTD1) was only detected in alpha-cypermethrin-resistant samples. Unique downregulated detoxification genes associated with lambda-cyhalothrin resistance were CYP4H29, CYP307A1, CYP6P12, and GSTE6 (Table 2).

3.4. Gene Ontology enrichment analysis

Gene ontology enrichment analysis (GOEA) was performed on the list of DEGs associated with all the 6 R-S comparisons (Table S7). Not surprisingly, the differences seen in the transcriptomic profiles between the different comparisons was also reflected in the total number of enriched gene ontology (GO) terms associated with DE genes (Table 1, Table S9). GO analysis of the upregulated genes for the R-S comparisons, identified the enrichment of some relevant molecular functions that are strongly associated with multiple-insecticide resistance phenotypes, including “monooxygenase activity” (GO:0004497), “oxidoreductase activity” (GO:0016705), “heme binding” (GO:0020037), and “tetrapyrrole binding”(GO:0046906) (Figure 6). Additionally, the molecular functions “catalytic activity” (GO:0003824) and “iron ion binding” (GO:0005506) were enriched in 4 out of 6 upregulated gene lists analyzed, reflecting their association with multiple-insecticide resistance. Interestingly, the enrichment of some relevant GO terms was found to be specific to malathion resistant mosquitoes, including several GO terms associated with postsynaptic neurotransmitter activities (GO:1904315, GO:0099529, GO:0098960, GO:0030594), transmembrane transporter activities (GO:0022803, GO:002285, GO:0022857), ion channel activities (GO:0005216 , GO:0005261, GO:0015267, GO:0022848, GO:0099094, GO:0005231, GO:0022848) and protein receptor activities (GO:0004930, GO:0005102) (Figure 6). For the list of downregulated genes, several GO terms were overrepresented in 4 out of the 6 R-S comparisons. These included “catalytic activity” (GO:0003824), “carboxypeptidase activity” (GO:0004180), “metallocarboxypeptidase activity” (GO:0004181), “structural constituent of chitin-based cuticle” (GO:0005214), “chitin binding” (GO:0008061) and “exopeptidase activity” (GO:0008238) (Figure S4). GO terms associated with chitin binding activities were overrepresented in the list of downregulated genes of all the six R-S DEG comparisons, which was also clearly reflected in the functional annotation of the DEGs displayed in Figure 2, where genes associated with chitin and chitinase activities were mostly downregulated (Figure S4).

3.5. Validation of relative expression levels estimated by RNA-Seq with quantitative RT-PCR

Five transcripts differentially expressed in resistant samples (I-ACYP, I-LCT and M-LCT) including GSTD6, GSTE4, CYPBB2, CYP6M9 and CYP6N9, were used to validate the RNA-Seq results by quantitative RT-PCR. The threshold cycle (ct) value was not detected for 2 of the 15 qPCR assays, likely associated to the low expression level of these transcripts in the mosquitoes. This is evident in the low read count from the RNA-seq. The qRT-PCR results broadly supported the directionality of the changes in expression levels (75% of the essays with ct values), although for several genes, the magnitude of the expression difference was smaller than estimated by RNA-Seq (Figure S5).

4. Discussion

Insecticide resistance in Ae. aegypti continues to expand globally due the extensive use of insecticides for its control. In this study, we used whole transcriptomic approach to investigate the molecular basis of resistance in Ae. aegypti populations from two locations in Puerto Rico, Isabela and Manatí, that were resistant to an organophosphate (malathion) and two pyrethroids (lambda-cyhalothrin and alpha-cypermethrin).

Differential gene expression analysis of genes associated with resistance to all the three insecticides from the two locations showed similarities in overexpression of the detoxification genes belonging mainly to the cytochrome P450 gene family and downregulation in cuticular protein genes. The insect cuticle is comprised mainly of chitin and cuticular proteins, and modifications to the cuticle can lead to thickening of the cuticle hence slowing the penetration of insecticides (Balabanidou et al. 2018). Cuticular thickening has been associated with resistance to pyrethroids in Anopheles funestus (Wood et al. 2010) and with reduced penetration of deltamethrin in highly resistant strain of Anopheles gambiae (Balabanidou et al. 2016). The downregulation of cuticular proteins in this study suggests that resistance in these populations is predominantly driven by cytochrome P450-mediated detoxification of all three insecticides. A study on the transcriptomic profile of a resistant strain of Ae. aegypti from Vietnam similarly reported downregulation of cuticular proteins in resistant versus susceptible strains (Lien et al. 2019).

In all the resistant groups from all three insecticides, a serine protease (AAEL029072) was significantly upregulated with a high fold change (Figure 5B). Serine protease has been reported to be highly upregulated in An. gambiae resistant to DDT (Vontas et al. 2005) as well as in Culex pipiens pallens resistant to deltamethrin (Gong et al. 2005), where it was suggested to play a crucial role in the innate immunity in this species and other insects (El Moussawi et al. 2019, Wang et al. 2021, Shan et al. 2022). This finding suggests that the increase in protease activity may be important in modulating Ae. aegypti resistance to multiple insecticides, however this remains speculative as this has not yet been functionally proven.

Several CYP450s were also commonly differentially expressed in the three resistant groups. These included CYP6Z7, CYP28A5, CYP9J2, CYP6Z6, CYP6BB2, CYP6M9, and two CYP9F2 orthologs, suggesting that these CYPs may confer cross resistance to pyrethroids and organophosphates. The majority of these CYPs have been associated with resistance to insecticides in previous studies. CYP6BB2 has consistently been reported to be overexpressed in Ae. aegypti resistant to pyrethroids in Asia, Europe, and the Americas (Bariami et al. 2012, Kasai et al. 2014, Reid et al. 2014, Dusfour et al. 2015, Moyes et al. 2017). The overexpression of CYP6BB2 has also shown strong metabolic activity for permethrin in in vitro metabolic studies (Kasai et al. 2014), highlighting its importance as a candidate metabolic resistance marker for pyrethroids in Ae. aegypti populations. In the same study, CYP6Z7 and CYP6Z6 were also overexpressed in permethrin-resistant Ae. aegypti. In another study, CYP6Z7 and CYP6BB2 were among four of the only CYPs upregulated in three tested strains of Ae. aegypti resistant to deltamethrin while CYP9J2 also appeared in one of the resistant strains in the same study (Epelboin et al. 2021). Multiple studies have also reported the upregulation of CYP6Z6 in pyrethroid resistant Ae. aegypti populations around the world (Marcombe et al. 2009, Bariami et al. 2012, Saavedra-Rodriguez et al. 2012, Reid et al. 2014, Dusfour et al. 2015).

Despite some shared gene expression patterns across the three insecticides, some of the genes were only found to be differentially regulated in survivors of specific insecticides and not others. In the malathion resistant samples from both locations, CYP6BY1 was found to be upregulated but was not upregulated in the lambda-cyhalothrin or alpha-cypermethrin resistant samples. CYP6BY1 has previously been reported to be among only two overexpressed P450s in Colombian Ae. aegypti larvae resistant to temephos, which is also an organophosphate insecticide (Morgan et al. 2022). In samples resistant to alpha-cypermethrin, the glutathione S-transferase gene GSTD1 was found to be uniquely overexpressed in this group compared to survivors of the other two insecticides. GSTD1 has been implicated in resistance to DDT in both An. gambiae (Ranson et al. 2001) and D. melanogaster (Tang and Tu 1994). Two genes, CYP4H29 and GSTE6, were uniquely overexpressed in samples resistant to lambda-cyhalothrin. In a targeted RNA-seq study, GSTE6 was the only GST reported to be associated with deltamethrin resistance across eight Ae. aegypti strains from distinct geographical origins (Faucon et al. 2017). CYP4H29 has previously been associated with resistance to pyrethroids in Puerto Rican Ae. aegypti (Reid et al. 2014) The same study reported the overexpression of four CYP450s belonging to family four. Studies on the role of cytochrome P450 4 (CYP4) family in resistance are scarce, however, Reid et al showed increased survival of D. melanogaster expressing CYP4H29 when exposed to permethrin.

Another difference that stands out is the enrichment of GO terms belonging to the postsynaptic neurotransmitter receptor activity mainly in the alpha-cypermethrin resistant samples from both Isabela and Manatí. This supports the evidence that type II pyrethroids such as alpha-cypermethrin have been shown to mainly inhibit GABA_A receptor function, hence leading to higher excitation in the insect’s central nervous system (Eldefrawi and Eldefrawi 1989). However, the same GO terms did not appear enriched in the lambda-cyhalothrin resistant group which is also a type II pyrethroid. According to the WHO report 1998 (World Health 1998), alpha-cypermethrin was more toxic than other pyrethroids and was significantly more effective at killing anopheline mosquitoes than permethrin or lambda-cyhalothrin and this may explain the high activity in the enriched GO terms. In addition, the bottle bioassay results showed a higher frequency of resistance to alpha-cypermethrin in the Manatí population than in the Isabela population, potentially also explaining the higher enrichment score in this group.

6. Conclusions

This study investigated the gene expression profiles associated with pyrethroid and organophosphate resistance in Ae. aegypti populations from two locations in Puerto Rico. We identified several detoxification genes which exhibited similar differential over-expression patterns across all the insecticide resistant phenotypes from both populations, suggesting that these genes could be used as expression-based markers for multiple-insecticide resistance in Ae. aegypti. We also identified some detoxification genes that were uniquely overexpressed in mosquitoes that survived exposure to malathion, alpha-cypermethrin, or lambda-cyhalothrin, indicating that their overexpression is more closely associated with resistance to specific insecticides. Additionally, we report the association of significant overexpression of serine protease genes and downregulation of cuticular-related protein genes (chinase and chitin) with organophosphate and pyrethroid resistance in both Ae. aegypti populations. The genes identified in this study, particularly in the cytochrome P450 family, should be functionally validated to confirm their importance as molecular markers for resistance detection.