Ae. aegypti Field Collection and Colonization

Ae. aegypti were collected from the Colombian municipalities of Bello and Cúcuta (Fig. 1). These municipalities were previously shown to be distinct in the burden of Ae. aegypti borne disease, socioeconomic status, and climate (Morgan et al. 2021). Mosquito collections took place in 2016 (Bello) and 2017 (Cúcuta) with the assistance of personnel from biology and control of infectious diseases research groups (University of Antioquia) and vector-borne disease program staff within each municipality. Immature Ae. aegypti were collected from deployed oviposition traps (ovitraps) and reared to adults under standard conditions at Universidad de Antioquia, Colombia. Standard rearing conditions were 28 ± 1°C, 80 ± 5% relative humidity and a 12 hr light: 12 hr dark photoperiod. Reared adults were offered a bloodmeal and the eggs were collected for the establishment of colonies. Upon establishment of colonies, eggs were collected and sent to Edge Hill University, UK for insecticide resistance profiling.

Larval Bioassays

Larvae for use in insecticide bioassays were reared under standard conditions within Edge Hill University Vector Research Group insectaries. Standard conditions were 27°C and 70% RH with an 11-hr day/night cycle with 60-min dawn/dusk simulation periods, using a lighting system of 4× Osram Dulux 26W 840 lights. Eggs were submerged in a hatching broth of 350 ml distilled H₂O (dH₂O), 0.125 g nutrient broth (Sigma–Aldrich, Darmstadt, Germany), and 0.025 g brewer's yeast (Holland & Barrett, Ormskirk, UK) for 48 hr (Zheng et al. 2015). Larvae were fed ground fish food (AQUARIAN advanced nutrition) and raised until third to fourth instar. Larval bioassays were conducted following WHO standard test procedures (WHO 2016). Preliminary testing was conducted to identify the activity range of temephos to larvae from each of the study municipalities and a susceptible laboratory strain (New Orleans). The activity ranges, yielding mortality of 10–95%, for each Ae. aegypti population as identified by preliminary testing are displayed in Table 1. At least four replicates, each consisting of 20 third to fourth instar larvae, were conducted for each temephos concentration. Fresh insecticide solutions were made for each replicate using temephos (93.7%; Pestanal, Sigma–Aldrich Darmstadt, Germany) and acetone (□≥99.9%; Sigma–Aldrich, Darmstadt, Germany). Bioassays were conducted inside the insectaries under standard conditions and mortality was recorded after a 24-hr exposure period. Following WHO guidelines, moribund larvae were counted as dead. Controls were exposed to the acetone solvent only.

Fig. 1.

The location of the two study sites. Mosquito collections took place in these two locations, Bello and Cúcuta, within Colombia. Departments are the largest units of local government. (A) Department of Antioquia governs Bello. (B) Department of Norte de Santander has as its capital Cúcuta, a city to the East of this department on the border with Venezuela. Map base layers were obtained from  https://data.humdata.org/dataset/colombia-administrative-boundaries-levels-0-3 covered by a Creative Commons Attribution 4.0 International (CC BY) License ( https://creativecommons.org/licenses/by/4.0/legalcode). Map base layers were modified by the addition of colors.

RNA Extraction and cDNA Synthesis

Aedes aegypti Sample Groups

Following the larval bioassays, field samples were categorized as resistant (Cúcuta) or susceptible (Bello). In the latter category were also samples from the lab adapted reference strain New Orleans. The Cúcuta temephos resistant samples were further divided into two groups: one exposed to temephos (for 24 hr) immediately before sampling for RNA extraction, and a control group of no temephos exposure (unexposed) (Fig. 2). For each group (field susceptible (FS), lab susceptible (LS), field resistant exposed (FRE) and field resistant unexposed (FRU) RNA extractions were carried out from four different larvae batches considered here as biological replicates.

Table 1.

Temephos activity range yielding between 10-95% mortality to each mosquito population. Results of preliminary testing to identify the activity range of temephos to each mosquito population. Activity ranges displayed in parts per million (ppm)

Temephos Exposure Assays for Resistant Larvae

Larvae in the resistant exposed group were exposed to temephos at the LC₅₀ of 0.06 ppm (Fig. 2). Larvae in all unexposed groups were exposed to the equivalent volume of acetone (the solvent used in temephos solutions). After the 24 hr bioassay larvae were taken for RNA extraction. For each experimental group there were four independent replicates, conducted using eggs from different batches and rearing and extraction conducted on different days. Egg submersion, feeding, bioassays, and RNA extraction on all replicates were all conducted at the same times of day.

Fig. 2.

Block diagram of our experimental approach. This study concerned the larval stages of the mosquito Ae. aegypti. Field samples denoted as Resistant originated from the Cúcuta, Colombia population. The Susceptible samples had dual origin: Field samples from Bello, Colombia (Field Susceptible) and the New Orleans reference lab strain (Lab Susceptible). The total RNA sequenced and mapped (Data analysis) originated from four different experiments (biological replicates) from each population. The gene expression levels of the Resistant samples compared against the Lab Susceptible and Field susceptible had at least 503 differentially expressed genes (DEG). The Resistant samples of larvae transiently exposed to temephos had 13 DEGs in comparison to the unexposed larvae. The functional annotation for the DEG sets was carried with several different repositories: VectorBase, Gene Ontology (GO), and KEGG Enrichment Analysis.

Data Analysis

RNA-Seq Data Quality Control, Mapping, and Differential Gene Expression

Analyses of RNA-Seq data were conducted using the Linux command line (Ubuntu 18.04) and R statistical software (version 4.0.3) (R Core Team 2020). Sequence read quality was assessed using FastQC (version 0.11.3). Reads with quality scores less than 30 and lengths less than 50 bp were trimmed using cutadapt (version 2.10) with Python (version 3.8.3). Read quality was then reassessed using FastQC to ensure only high-quality reads remained. Cleaned reads were mapped to the Ae. aegypti LVP_AGWG reference genome (version AaegL5, GenBank: NIGP00000000.1) using Rsubread (version 2.2.6) (Liao et al. 2019). The resultant BAM files were sorted and indexed using samtools (version 1.11). Alignment quality metrics from Rsubread were visualized using R's plotting function. Gene count tables were generated using Rsubreads. Read counts were normalized using the trimmed mean of M values method (Robinson and Oshlack 2010) in edgeR (Version 3.30.3) which accounts for library size and expression bias in RNA-Seq datasets (Robinson et al. 2010). TMM normalization was conducted between groups in each comparison rather than globally across all samples to ensure within-group variability could be detected. Differential gene expression was then calculated using a Quasi-likelihood negative binomial generalized log linear model (edgeR). Quasi-likelihood error family was selected due to its ability to account for uncertainty in dispersion. Counts per million (CPM) were calculated in edgeR and reads per kilobase million (RPKM) were calculated in edgeR using transcript lengths obtained from Enseml Metazoa (LVP AGWG (aalvpagwg_eg_gene) dataset) using biomaRt (version 2.44.4). Transcripts with fold-change >2 and FDR <0.05 were selected for gene ontology and KEGG pathway enrichment analyses, these thresholds were used to select genes with significant and larger differential expression.

Table 2.

LC50 and LC95 of Bello, Cúcuta, and New Orleans Ae. aegypti larvae to temephos.Temephos bioassays showing LC50, LC95, and their 95% confidence limits, calculated using probit analysis. Resistance ratios (RR) were calculated as a ratio of the lethal concentration (LC50 and LC95) of each population compared to the lab susceptible (New Orleans) strain. Bello and New Orleans are both susceptible to temephos whilst larvae from Cúcuta are resistant. SE = standard error

Temephos Susceptibility of Ae. aegypti Field Isolates and the Reference Strain

The lethal concentrations (LC₅₀) of temephos were 0.019 ppm (95% Cl 0.016–0.029) and 0.060 ppm (Cl 95% 0.052–0.070) for Bello and Cúcuta, respectively. This corresponded to resistance ratios of 2.6 and 8.0 when compared to the New Orleans susceptible laboratory strain (LC₅₀ = 0.008; Cl 95% 0.006–0.012). The LC₉₅ for Bello was 0.055 (Cl 95% 0.034–0.292), 2.8-fold higher than New Orleans (LC₉₅ = 0.020; Cl 95% 0.012–0.21). Cúcuta had an LC₉₅ of 0.182, 9-fold higher than New Orleans (Table 2). Following the WHO guidelines (WHO 2016) larvae from Bello were considered susceptible to temephos and denoted as the field susceptible (FS) population whist the larvae from Cúcuta were moderately resistant to this larvicide and denoted as the field resistant (FR) population (Table 2). New Orleans is referred to as the lab susceptible (LS) population.

RNA-Seq Mapping Summary

Three different populations of Ae. aegypti larvae were profiled using RNA-Seq: temephos field resistant (FR), field susceptible (FS), and the reference lab susceptible strain New Orleans (LS). The FR population was also split into two further groups as either exposed or not (unexposed) to temephos to determine further expression changes associated with insecticide exposure in already resistant populations (Fig. 2). Mosquito larvae from each of these four populations were grown in four different batches, one to four weeks apart, and were considered here four biological replicas. This generated a total of 16 sequenced samples.

Extracted total RNA was used to generate the Illumina RNASeq libraries that produced 65.7–120 million reads per sample with quality scores >30 and lengths >50 bp and 84.4% and 85.6% of those reads in each sample were successfully aligned to the reference genome (Table 3 and Materials and Methods). Using the most current gene model available for Ae. aegypti, LVP_AGWG reference genome which contains 19,381 open reading frames (version AaegL5, GenBank: NIGP00000000.1) the number of genes with successfully aligned reads ranged from 15,704 and 16,607 genes across all samples, corresponding to 81–86% of the total genes in the reference genome.

Table 3.

RNA-Seq sequencing data summary. The total number of obtained reads after quality control and the percentage of reads that mapped to the Aedes aegypti reference genome. Quality control removed reads with quality <30 and lengths <50 bp

Overview of Differential Gene Expression

Both field samples, susceptible and resistant (FS, FR) were equally and evidently distant, by RPKM number, to the samples of the lab strain New Orleans (LS) (Fig. 3). Crucially the temephos susceptible samples did not cluster together nor did the two field samples (Fig. 3). Field and lab strain samples were similarly distant in terms of differentially expressed (DE) transcripts: FS versus LS = 5324 DE transcripts (Fig. 4D and E), FR versus LS = 5579 (Fig. 4B and E). However, when comparing field samples (resistant and susceptible) the number of DE transcripts visibly lowered by four-fold to 1,454 (Fig. 4A and E). Therefore, the common practice in gene expression studies of comparing mosquito field samples to lab strains (e.g., Marcombe et al. 2009, Grisales et al. 2013, Dusfour et al. 2015, Ishak et al. 2017), would have generated an approximately 4-fold overestimation in the number of DE transcripts detected in the field resistant samples.

We also sought to investigate the gene expression changes in insecticide resistant larvae under transient exposure to insecticide. There were only 19 transcripts significantly differentially expressed in larvae within the resistant population which were exposed to temephos when compared to unexposed larvae from the same population (Fig. 4E). All 19 of those transcripts were overexpressed in the exposed group with no significant down regulated gene expression detected (Fig. 4C).

We addressed the issue of potential misrepresentation of gene expression metrics by triangulating both, the differentially expressed gene (DEG) sets and the RPKM counts between the field resistant (FR) samples against the field susceptible (FS) as well as the lab susceptible (LS) samples. The DEG set obtained contained transcripts which were found to be significantly differentially expressed with a fold change of >2 and a false-discovery rate (adjusted P value) of <0.05 in both comparisons. Under this threshold, a total of 623 (down from 1,454 transcripts in only the field-to-field comparison) transcripts covering 503 genes, were found differentially expressed in the field resistant population when compared to both the field and lab susceptible populations (Fig. 4E). This set of 503 genes comprised 301 overexpressed genes and 202 under expressed genes ( Supp Table 1 (tjab179_suppl_supplementary_table_s1.xlsx) [online only]). Of the 301 significantly overexpressed genes 239 were found in the category of protein coding genes: 88 annotated and 151 hypothetical genes. In the significantly under expressed gene set 166 were protein coding with 75 annotated and 91 hypothetical genes ( Supp Table 1 (tjab179_suppl_supplementary_table_s1.xlsx) [online only]). The significant differentially expressed genes also included 55 overexpressed genes encoding for long noncoding RNA (lncRNA) and 30 under expressed lncRNA genes in the temephos resistant larvae ( Supp Table 1 (tjab179_suppl_supplementary_table_s1.xlsx) [online only]).

Fig. 3.

Distribution of the transcriptomic profiles by sample groups.The multidimensionality of the RPKM values calculated for each mapped transcript per sample was reduced by principal components analysis (PCA). The field samples were seen linearly distant from each other across one component whist the reference lab strain NO cluster (Lab Susceptible) separated from both field samples. The orthogonal dispersion of these samples allowed for the triangulation of the data as described in the main text.

Fig. 4.

Differential gene expression in the samples from field susceptible, lab susceptible, and field resistant populations. (A–D) Differential expression of all transcripts including those over expressed, under expressed, and with no significant differential expression in the field resistant unexposed population compared to the field susceptible population (A) and susceptible lab population (B), the resistant temephos exposed population when compared to the resistant unexposed population (C) and in the field susceptible population when compared to the lab susceptible population (D). (E) The number of transcripts significantly differentially expressed (FC > 2, FDR < 0.05) between each of the experimental groups. The number of differentially expressed transcripts shown here (E) includes both significantly over and under expressed transcripts. The comparison groups and sample notation as detailed in Fig. 2.

Fig. 5.

GO terms and KEGG pathways enriched in the resistant population when compared to both field and lab susceptible populations. GO terms and KEGG pathways found to be significantly enriched (P < 0.05) following Benjamini Hochberg correction in the significantly over (A) and under (B) expressed transcripts. Enrichment percentage was calculated as the number of differentially expressed transcripts in each category/pathway divided by the total number of transcripts in the same category/pathway. Number in bars indicates the number of differentially expressed transcripts in each category.

Gene Ontology and KEGG Pathway Enrichment

All 623 transcripts with fold-change >2 and FDR <0.05 were selected for gene ontology (GO) and KEGG pathway enrichment analyses. GO categories and KEGG pathways with corrected p values (FDR) <0.05 were considered significantly enriched. Gene ontology and KEGG enrichment analyses conducted on the 301 significantly over expressed genes identified eight significantly enriched GO categories; one involved with biological processes, oxidation–reduction processes (GO:0055114), and seven associated with molecular functions (GO:0045735, GO:0016705, GO:0004100, GO:0004022, GO:0005506, GO:0016491, GO:0047938) (Fig. 5). KEGG enrichment analysis identified two significantly enriched KEGG pathways in the over expressed genes; insect hormone biosynthesis (path:00981) and ubiquinone and other terpenoid-quinone biosynthesis (path:00130) (Fig. 5).

GO and KEGG enrichment analysis were also conducted on the 202 significantly under expressed genes identifying 12 significantly enriched GO terms and one KEGG pathway (path:00981) (Fig. 5). The enriched GO terms include three terms involved in biological processes (GO:0090150, GO:0046416, GO:0006470), one involved with cellular components (GO:0005615), and seven associated with molecular functions (GO:0004866, GO:0004721, GO:0003884, GO:0032977, GO:0047938, GO:0016641, GO:0071949) (Fig. 5).

The Overexpressed Transcriptome of Field Temephos Resistant Aedes aegypti Larvae

The transcriptomic overview provided by the GO and KEGG enrichment models was interrogated by quantifying the represented genes with CPM and RPKM metrics. The expression profiles of the 88 annotated overexpressed protein coding genes in the resistant population were visualized using heatmaps of gene expression as log₂ values of counts per million (CPM) (Fig. 6) as well as bar plots of reads per kilobase million (RPKM) values of the represented genes (Fig. 7). The former allowed the visualization of the data's granularity by comparing gene expression across all 16 samples individually rather than just across groups (Fig. 6). Variation in expression between samples from the same experimental group can be seen across all genes (Fig. 6) highlighting the importance of biological replication in gene expression experimentation. The heatmap also showed the importance of using field susceptible populations in addition to lab susceptible comparator populations. Differences in gene expression (e.g., CYP6B1 and PGRPLA) can be overrepresented when comparing expression profiles between field and lab populations (here FR and LS) rather than between field to field (e.g., FR and FS). Differences in gene expression were also visualized using RPKM bar plots which enable ranking of genes specifically overexpressed in resistant samples, the majority (151 of the 239 protein coding genes) of those did not have a functional annotation in the current repository for VectorBase (Fig. 7,  Supp Table 1 (tjab179_suppl_supplementary_table_s1.xlsx) [online only]).

The over expressed annotated protein coding genes included detoxification enzymes; two cytochrome P450s (CYP12F6 – AAEL002005 and CYP6BY1 – AAEL017539), a carboxy/cholinesterase (CCEAE4C – AAEL003187), a glutathione S-transferase (AAEL006818), two glucosyl/glucuronosyl transferases (AAEL002688 and AAEL003076) and an aldehyde oxidase (AAEL014493). The cuticular biosynthesis enzyme chitin synthase (AAEL002718) and the digestive enzymes, putative trypsin genes, AAEL007102, AAEL014579, and AAEL003308, were also present.

Other over expressed genes included the hydrocarbon biosynthesis pathway enzyme acetyl-CoA dehydrogenase (AAEL014452) (Jones et al. 2013), glutamate decarboxylase (AAEL010951) which catalyzes biosynthesis of GABA through glutamate decarboxylation (Richardson et al. 2010), sarcosine dehydrogenase (AAEL014936), a mitochondrial glycine synthesizing enzyme (Huang et al. 2015), leucine aminopeptidase (AAEL006975), a proteolytic enzyme that hydrolyses amino acids with roles in toxin biosynthesis (Matsui et al. 2006), a manganese-iron (Mn-Fe) superoxide dismutase (MNSOD1 – AAEL004823), a mitochondrial antioxidant associated with increased life span in insects (Clancy et al. 2001, Kang et al. 2008) and two mannose-binding C-Type Lectins (CTLs) AAEL011612 and AAEL000533, ubiquitous proteins in multicellular organisms that provide the pattern recognition required for the initial phase of an immune response (Ourth et al. 2005, Phillips and Clark 2017).

The Under Expressed Transcriptome of Field Resistant Aedes aegypti Larvae

The transcript profiles of the 75 annotated protein coding genes significantly under expressed genes in the resistant population were also visualized in heatmaps of gene expression (log₂ values of CPM) (Fig. 8) as well as bar plots of RPKMs (square root values) of the represented genes (Fig. 9). The former showing variation between the 16 samples and the latter displaying variation between genes and between resistance status.

Genes encoding detoxification enzymes were also presented in this set of under expressed genes. Those included the cytochrome P450 CYP314A1 (AAEL010946) and a glucosyl/glucuronosyl transferase (AAEL003098). A cytochrome oxidase biogenesis protein (oxa1 mitochondrial – AAEL009183), essential for full expression of cytochrome c oxidase was also under expressed. Other genes significantly under expressed in the resistant population include a putative pupal cuticle protein (AAEL011444) and transferrin (Tf1 – AAEL015458) a regulator of iron metabolism with roles in mosquito innate immunity (Yoshiga et al. 1997). The mdg4-binding protein ortholog gene in Ae. aegypti (AAEL010576: Modifier of mdg4 [Mod(mdg4)]), responsible for chromosome remodeling was also represented in this group of under expressed genes.

The expression of several ion and solute membrane transporters was also down regulated. These included the sodium-coupled cation–chloride cotransporter AAEL009886 (aeCC3), the sodium/chloride dependent amino acid transporter AAEL000298, the sodium/solute symporter AAEL001198, and the sugar transporter AAEL010348. In the group of under expressed genes were also 30 lncRNA genes in the temephos resistant larvae ( Supp Table 1 (tjab179_suppl_supplementary_table_s1.xlsx) [online only]).

Gene Expression Profile of Temephos Exposed Larvae from the Resistant Population

Gene expression in the field resistant population following the controlled exposure to temephos was compared with gene expression of samples from the same population without insecticide exposure. The exposed samples showed 19 significantly (FC > 2, P < 0.05) overexpressed transcripts (Fig. 10A &  Supp Table 2 (tjab179_suppl_supplementary_table_s2.xlsx) [online only]) in comparison to the nonexposed samples. These 19 transcripts were mapped to 13 genes (Fig. 10). The products of the over expressed genes include a sodium/chloride dependent amino acid transporter (AAEL003619), an alkyl dihydroxyacetone phosphate synthase (AAEL007793), cathepsin-1 (AAEL011167), trypsin-1 (AAEL016975), and a serine protease stubble (AAEL020367). The remaining eight overexpressed genes had uncharacterized products in Ae. aegypti ( Supp Table 2 (tjab179_suppl_supplementary_table_s2.xlsx) [online only]).

Metabolic Detoxification Genes

Metabolic detoxification of insecticides is one of the most reported insecticide resistance mechanisms in mosquitoes. Abundant overexpression of detoxification genes, most commonly cytochrome P450 monooxygenases (P450), glutathione S-transferases (GST), and carboxylesterases (CE) has frequently been associated with insecticide resistance in mosquitoes (Vontas et al. 2020). Here we reported the overexpression of only two P450s (CYP12F6 – AAEL002005 and CYP6BY1 – AAEL017539), one GST (AAEL006818) and one CE (CCEAE4C [AAEL003187]), in the resistant population compared to both field and lab susceptible populations. CYP12F6 has previously been shown to be overexpressed in a permethrin resistant population of adult Ae. aegypti from Mexico albeit compared with a lab susceptible population only (David et al. 2014). GO terms associated with insecticide detoxification (oxidoreductase activity (GO:0016491) and oxidation–reduction process [GO:0055114]) were also found to be enriched in the temephos resistant larvae. Thus, by cross examining the data in field-to-field and field-to-lab population comparison, we observed genes representing these three forms of insecticide deactivation in much reduced number compared to what is commonly reported (Marcombe et al. 2009, Grisales et al. 2013, Dusfour et al. 2015, Ishak et al. 2017). To illustrate the above, if the resistant population had been compared with the lab susceptible population only a total of 49 cytochrome P450s, six GTSs, and 11 CEs would have been reported as differentially expressed ( Supp Table 3 (tjab179_suppl_supplementary_table_s3.xlsx) [online only]). This suggests that large overexpression of detoxification genes may be partly related to differences between field and lab mosquitoes rather than associated with the insecticide resistant phenotype. Large overexpression of detoxification in mosquitoes may also only be observed in mosquitoes when they have high levels of resistance rather than the moderate resistance reported here (Dusfour et al. 2015, Ishak et al. 2015).

Chitin Biosynthesis

The thickness and composition of the cuticle have been identified as a critical determinant of insecticide resistance due to its role in reducing insecticide penetration (Balabanidou et al. 2018). Over expression of genes associated with formation and maintenance of the cuticle have been reported in insecticide resistant populations of medically relevant species including An. gambiae (Awolola et al. 2009, Balabanidou et al. 2016, Yahouédo et al. 2017), An. funestus (Giles, Diptera: Culicidae) (Gregory et al. 2011) and Culex pipiens pallens (Linnaeus, Diptera: Culicidae) (Pan et al. 2009, Fang et al. 2015). The cuticle has also been associated with resistance in Ae. aegypti including in larvae (David et al. 2010, Riaz et al. 2013). The over expression of the chitin biosynthesis enzyme AAEL002718 and the enrichment of chitin synthase activity (GO:0004100) in temephos resistant Ae. aegypti larvae reported in this study further highlights the potential role of the cuticle in the development of insecticide resistance in Ae. aegypti larvae, even when resistance is moderate. Chitin, a biopolymer of N-acetylglucosamine, is a major constituent of the mosquito cuticle (exoskeleton [epidermal cuticle], tracheal cuticle, and eggshell) providing it with both strength and rigidity and is also found in midgut peritrophic matrices (Merzendorfer and Zimoch 2003). The use of chitin synthesis inhibitors (CSI), a type of insect growth regulators (IGRs) which interfere with the synthesis and deposition of chitin on the exoskeleton (Tunaz and Uygun 2004), has been highlighted as a potential approach to control Ae. aegypti with some promising findings in laboratory studies (Martins et al. 2008, Fontoura et al. 2012).

Pattern Recognition and Innate Immunity

Temephos resistant Ae. aegypti larvae were shown to express high levels of the mannose-binding C-type lectins (CTLs) AAEL011612 and AAEL000533 which are predominantly produced in the salivary glands of adult female Ae. aegypti (Ribeiro et al. 2016, Adelman and Myles 2018). Lectins are ubiquitous proteins in multicellular organisms that provide the pattern recognition required for the initial phase of an immune response (Ourth et al. 2005, Phillips and Clark 2017). C-type lectins are a group of calcium-dependant carbohydrate binding proteins (Zelensky and Gready 2005). In mosquitoes CTLs are primarily involved in facilitating viral infection (e.g., dengue, Rift Valley fever and Japanese encephalitis viruses (Liu et al. 2017, Licciardi et al. 2020)) through the enhanced viral entry, acting as bridges between flaviviruses and host cell receptors (Cheng et al. 2010, Liu et al. 2017). However, these proficient pattern recognition proteins seem to have evolved to mediate multiple multicellular processes beyond mosquito immune response including lifespan and reproductive capability (Li et al. 2020) as well as maintenance of gut microbiome homeostasis (Pang et al. 2016). Transferrin (Tf1 –AAEL015458), found to be under expressed in temephos resistant larvae in the current study, also has roles in the innate immune response to arbovirus infection (Wessling-Resnick 2018, Licciardi et al. 2020) and has previously been reported to be down regulated in CHIKV and DENV infected mosquitoes which may favor viral replication (Tchankouo-Nguetcheu et al. 2010). Transferrin expression has also been related to insecticide resistance in Culex pipiens with increased expression reported in mosquitoes with target-site resistance to pyrethroids and organophosphates, the biggest difference in transferrin expression was observed in adults (Tan et al. 2012, Vézilier et al. 2013).

Cell Membrane Transport

The expression of several ions coupled solute membrane transporters was down regulated in temephos resistant larvae: the sodium-coupled cation–chloride cotransporter AAEL009886 (aeCC3), the sodium/chloride dependent amino acid transporter AAEL000298, the sodium/solute symporter AAEL001198, and the sugar transporter AAEL010348. aeCCC3 is a larvae specific membrane transporter abundant in the anal papillae responsible for the absorption of external ions (Piermarini et al. 2017) which belongs to a family of cation-coupled chloride cotransporters (CCCs) which contribute to ion homeostasis by undertaking electroneutral transport of Na⁺, K+, and Cl^– (Pullikuth et al. 2003). A similar role is expected from the ion-coupled transporters AAEL000298 and AAEL001198 in the homeostasis of ion content, particularly in the midgut and Malpighian tubes where they are most abundant (Li et al. 2017).

The aquatic life of the Ae. aegypti larval stages demands an ion exchange homeostasis that differs from that of the adult mosquitoes. Due to their freshwater habitat Ae. aegypti larvae must excrete water gained by osmosis, reabsorb salt before excreting urine, and absorb salt from their surroundings (Ramasamy et al. 2021). Whilst the opposite is true in adults where water retention is needed due to constant loss through evaporation. A key process in this is Na+-dependent co-transport which is typically down the large inward (extracellular to intracellular) Na+ gradient generated by the Na+/K+-ATPase (Payne 2012). We speculate that the ion homeostasis changes caused by the reduced expression of the three CCCs transporters AAEL009886, AAEL000298, and AAEL001198 could reduce the exposure of larvae to temephos by reducing net uptake of the molecule, protecting the organs where they are commonly expressed (e.g., midgut and Malpighian tubes). Transcriptome studies of insecticide resistant mosquito populations tend to overlook the potential role of down regulated genes in favor of overexpressed genes, but this finding demonstrates the importance of investigating reduced expression when studying potential mechanisms of insecticide resistance. The potential role of CCC transporters in reducing insecticide uptake and therefore facilitating resistance warrants further investigation.

Chromosomal Remodeling

The mdg4-binding protein ortholog gene (AAEL010576: modifier of mdg4 [Mod(mdg4)]), responsible for chromosome remodeling was also significantly under expressed in the resistant Ae. aegypti larvae. Originally described as a protein binding the transposon mdg4 (Bayev et al. 1984), Mod(mdg4) gene encodes for a family of proteins due to at least 30 different alternative splicing variants in Diptera and Lepidoptera (Krauss and Dorn 2004, Gabler et al. 2005, Tikhonov et al. 2018). Mod(mdg4) variants bind a variety of insulators (DNA domains involved in nuclear organization and higher order chromatin structures) (Ghosh et al. 2001, Melnikova et al. 2004, Golovnin et al. 2007) and have been involved in regulating numerous traits of the insect embryonic progression such as synapsis structure (Gorczyca et al. 1999), chromosome Y-linked testis development (Branco et al. 2013), and midgut maturation (Cai et al. 2012). Changes in expression of Mod(mdg4) have been reported in Drosophila Kc cells treated with deltamethrin (Liu et al. 2020). The downstream targets of mod (mdg4) have also been shown to affect the down regulation of proteins associated with protection from various stress conditions, therefore acting as modulators of cytotoxic damage in Drosophila (Lin et al. 2021). The identification of under expression of the mdg4-binding protein in temephos resistant larvae suggests a further role for this protein in mediating insecticide resistance with a potential role in regulating the stress response.

Long Noncoding RNA

There were 55 over expressed and 30 under expressed genes encoding for long noncoding RNA (lncRNA) in the temephos resistant larvae ( Supp Table 1 (tjab179_suppl_supplementary_table_s1.xlsx) [online only]). Noncoding RNAs (ncRNA) are abundant cellular effectors of great prolific functionality (Eddy 2001) and long ncRNA are defined as transcripts, more than 200 nucleotides long, that are produced by RNA polymerase II and are not translated into proteins (Bonasio and Shiekhattar 2014). In Aedes lncRNAs, mainly involved in regulating gene expression, are multifunction with roles including sex differentiation (Xu et al. 2019), embryo-genesis (Azlan et al. 2019), and suppression of viral replication in DENV infected mosquitoes (Etebari et al. 2016). Long ncRNAs have also been associated with insect's response to xenobiotics, with reports of differential lncRNA expression in resistant populations of Plutella xylostella (Etebari et al. 2015). The findings of 85 differentially expressed lncRNAs reported here in resistant populations of Ae. aegypti supports the potential roles that lncRNAs could have in the development of insecticide resistance. Whilst gene expression studies have focused primarily on differential expression in protein coding genes, the development of next generation techniques has now provided an opportunity to also study noncoding RNA. Whilst work has been conducted into identifying lncRNAs in medically relevant mosquito species including An. gambaie (Jenkins et al. 2015) and Ae. aegypti (Azlan et al. 2019) there have been no studies that have aimed to investigate the role of lncRNAs in insecticide resistance in mosquitoes. Previous RNA-Seq studies on insecticide resistant populations of Culex pipiens pallens have also identified differential expression of lncRNAs (Shi et al. 2020), however, an in-depth discussion of their role in insecticide resistance has been neglected.

Differential Gene Expression in Resistant Larvae Following Temephos Exposure

In the study, we also tracked gene expression in insecticide resistant larvae following direct response to temephos exposure. Thirteen genes were found to have a significantly increased expression following a controlled exposure to temephos. Among those 13 genes were two serine proteases: trypsin –1 (AAEL016975) and serine protease stubble (AAEL020367), a cysteine protease: cathepsin-1 (AAEL011167), a sodium/chloride dependent amino acid transporter (AAEL003619), and an alkyl dihydroxyacetone phosphate synthase (AAEL007793). Serine proteases are a group of enzymes with a variety of known functions including digestion, metamorphosis, oogenesis, blood coagulation, and viral immune response (Terra and Ferreira 1994, Mesquita-Rodrigues et al. 2011, Licciardi et al. 2020). Cathepsin-1 (AAEL011167), a cysteine proteinase, is also a multifunctional digestive enzyme (Terra and Ferreira 1994, Liu et al. 2006). Upregulation of serine proteases has been previously reported in insecticide resistant mosquito populations (Bonizzoni et al. 2012, Reid et al. 2012, David et al. 2014, Zou et al. 2016), including temephos resistant Ae. aegypti from Cúcuta (Grisales et al. 2013). Serine proteases have also been shown to degrade insecticides through hydrolysis within the insect digestive tract, however, so far evidence of this is limited to pyrethroids such as deltamethrin (Yang et al. 2008a, b; Xiong et al. 2014; Wilkins 2017). Overexpressed proteases in this current study support the findings of previous studies that proteases may have a role in the metabolism of other insecticide classes besides pyrethroids. The changes responsible for resistance are often associated with modification of physiological processes that can lead to decreased performance and fitness disadvantage.

Deleterious effects of insecticide resistance can affect a wide range of life-history traits (e.g., longevity, biting behavior, and vector competence) (Alout et al. 2014, 2017). Although the cost of resistance genes is believed to gradually decrease due to subsequent modifier mutations (Raymond et al. 2001). With the relatively limited diversity of insecticide targets (Swale 2019), the gene expression patterns that resistant mosquitoes further undergo when exposed to the insecticide could be a source for novel assets for vector control. The study of such targets for insecticide development is a strategy that, to our knowledge, has not yet been explored.