First Insights into the Comparative Transcriptomic Response of Field and Laboratory Aedes aegypti Strains to Partial-Mortality Concentration (<50%) Imidacloprid and Broflanilide Exposure

TL;DR

This study compares how wild and lab mosquitoes respond to two insecticides at the gene level, revealing differences in their molecular reactions that could help improve insecticide strategies.

Contribution

The study provides new insights into the transcriptomic responses of field and lab Aedes aegypti strains to sublethal insecticide exposure.

Findings

Both insecticides activated genes related to protein degradation, detoxification, and structural changes in mosquitoes.

The field strain showed more flexible gene responses, especially to imidacloprid, compared to the lab strain.

Broflanilide triggered similar early responses in both strains, while imidacloprid amplified existing differences.

Abstract

Aedes aegypti mosquitoes transmit viruses such as dengue, Zika and chikungunya and are mainly controlled with insecticides. However, many mosquito populations are becoming resistant, which reduces the effectiveness of these tools. In this study, we compared how two different populations of Aedes aegypti change their gene activity after being exposed to two insecticides: imidacloprid, a neonicotinoid, and broflanilide, a newer meta-diamide compound. One population was a field-derived strain from San Nicolás with a history of insecticide exposure, and the other was a fully susceptible laboratory strain called New Orleans. Using RNA sequencing, we measured which genes were turned on or off after a short, sublethal exposure to each insecticide. Both compounds caused strong changes in gene expression, mainly activating genes involved in protein degradation, membrane transport, detoxification and changes in the cuticle and cytoskeleton. We also found that the field population showed more flexible (plastic) gene responses than the laboratory strain, especially after exposure to imidacloprid. In contrast, early responses to broflanilide were more similar between populations. These results help us understand how different mosquito populations respond at the molecular level to new and repurposed insecticides and may support the design of better insecticide rotation strategies to delay resistance.

Insecticide resistance in Aedes aegypti (Linnaeus, 1762), the primary vector of several arboviruses, threatens vector control efficacy and motivates evaluation of current and candidate public health insecticides, such as imidacloprid and broflanilide, and their molecular impacts. Here, we used RNA sequencing (RNA-seq) to characterize the transcriptomic response to one-hour acute exposure to an operational partial-mortality concentration (<50%) of imidacloprid and broflanilide in two Ae. aegypti strains: a field-derived, pyrethroid-resistant population from San Nicolás and a susceptible laboratory strain (New Orleans). Adults were exposed for 1 h to partial-mortality concentration (<50%) doses of each insecticide or acetone control, and differential gene expression and Gene Ontology (GO) enrichment were assessed with DESeq2-based workflows. We detected pronounced baseline transcriptomic differences between strains and extensive activation of gene expression after insecticide exposure, with a strong bias toward up-regulation. A shared transcriptional core involving proteolysis, transmembrane transport, detoxification pathways, and structural remodeling of the cuticle and cytoskeleton was identified across contrasts. Despite these common elements, broflanilide elicited largely conserved early responses between strains, whereas imidacloprid amplified pre-existing divergence and produced marked population-specific transcriptional signatures. These findings suggest greater transcriptional changes in the field-derived strain, particularly in response to imidacloprid, and highlight the importance of integrating population-specific molecular information when designing insecticide rotation schemes and resistance management strategies targeting Ae. aegypti.

Figures

Full text

1. Introduction

Mosquitoes of the family Culicidae are among the most important vectors of human pathogens worldwide. Within this group, Aedes aegypti (Linnaeus, 1762) is the dominant urban vector of several arboviruses of public health concern in the Americas, including dengue, Zika, chikungunya and Mayaro [1,2]. Its close association with human dwellings, reliance on artificial containers, and strong dispersal capacity, as well as the ability of its eggs to withstand desiccation, all contribute to its successful establishment and persistence in a wide range of environments [3].

Current control programs for Ae. aegypti rely largely on insecticide-based interventions directed against larvae and adults. However, repeated use of a limited set of active ingredients has favored the spread of insecticide resistance in many populations, reducing the operational efficacy of existing pesticide formulations and complicating long-term vector management [3]. This scenario has stimulated interest in evaluating new compounds and alternative modes of action that could be incorporated into public health programs or used in rotation with existing insecticides. Messenger RNA (mRNA) represents the pool of protein-coding transcripts produced from the genome at a given time, and the complete set of these transcripts in a cell, tissue or organism under a specific developmental or physiological condition is referred to as the transcriptome [4,5]. Transcriptomic analyses, therefore, provide a powerful framework to quantify genome-wide changes in gene expression in response to insecticide exposure.

The continued spread of insecticide resistance in mosquito populations highlights the need to evaluate alternative compounds for vector control [6]. Imidacloprid is a neonicotinoid insecticide that has been proposed as an alternative in pyrethroid-resistant mosquito populations, with promising results reported for several culicid species, including Anopheles stephensi (Liston, 1901), Culex quinquefasciatus (Say, 1823) and Ae. aegypti [7]. Broflanilide is a recently developed meta-diamide insecticide that is currently registered mainly for agricultural use, but its novel mode of action and efficacy against other insect pests make it a promising candidate for inclusion in insecticide rotation schemes for vector control In insects, broflanilide acts on the γ-aminobutyric acid (GABA)-gated chloride channel at a distinct allosteric site, disrupting inhibitory neurotransmission and causing sustained neuronal hyperexcitation and death [8].

Understanding the molecular basis of insecticide resistance in Ae. aegypti is essential for designing effective resistance management strategies. Recent RNA-seq studies, such as that by Derilius et al., in 2023 [9], have characterized gene expression profiles associated with resistance to several widely used insecticides, including malathion, alpha-cypermethrin and lambda-cyhalothrin, in Ae. aegypti populations from Puerto Rico [9]. Similarly, Sun (2021) described how the transcriptome and proteome change in Ae. aegypti with respect to exposure to pyrethroid insecticides [10]. Furthermore, in the study conducted by Mack (2023), the transcriptome response of the mosquito to exposure to permethrin and different temperatures for 24 h was evaluated, generating data of great importance for understanding the reactions that affect the transcriptome in insect vectors [11]. However, the transcriptomic response of Ae. aegypti to neonicotinoids such as imidacloprid and to broflanilide remains largely unexplored, particularly when comparing field-derived pyrethroid-resistant and fully susceptible strains. Building on this approach, we used RNA-seq to investigate the transcriptomic response to partial-mortality concentration (<50%) exposure to two insecticides with distinct modes of action, imidacloprid and broflanilide, in two Ae. aegypti strains: a field-derived strain from San Nicolás (Nuevo León, Mexico) and a fully susceptible laboratory strain (New Orleans) free of known resistance mutations and previous insecticide exposure. By comparing the gene expression profiles of survivors across these strain–insecticide combinations, we aim to identify detoxification pathways and other stress-related genes that are differentially regulated in response to each compound. Generating this baseline information before the widespread use of these molecules in vector control could help anticipate potential resistance mechanisms and inform the design of more sustainable insecticide rotation strategies.

2. Materials and Methods

2.1. Mosquito Strains and Rearing

The field-derived strain from San Nicolás (SN) was established from Aedes aegypti eggs collected using ovitraps in the municipality of San Nicolás de los Garza, Nuevo León, Mexico (25.74167° N, −100.30222° W). After collection, egg papers were immersed in 5 L plastic containers filled with 4 L of dechlorinated water, and 0.3 g of brewer’s yeast was added to reduce dissolved oxygen and stimulate egg hatching. Newly hatched larvae were fed daily with ground commercial cat food suspended in water (5 g in 50 mL). Pupae were transferred to emergence chambers and adults were held in 30 × 30 × 30 cm entomological cages. The F_0_ adults were provided with 10% sucrose solution ad libitum, and females were offered heparinized bovine blood for egg production.

The susceptible laboratory strain New Orleans (NO), originally obtained from the Centers for Disease Control and Prevention (CDC), was reared under the same conditions for egg hatching, larval feeding and adult maintenance.

Both strains were maintained in the insectary at 24–26 °C, approximately 95% relative humidity and a 12:12 h light:dark photoperiod.

2.2. Experimental Design and Treatments

The experimental design for the gene expression analyses is summarized in Table 1. Adult mosquitoes from each strain were exposed to a partial-mortality concentration (<50%) of IMIDACLOPRID 18380 STD 10 mg (AccuStandard, distributed by CTR Scientific, Monterrey, Mexico) or BROFLANILIDE, 99% purity (Catalog No. E8093, Selleckchem, HO, USA), or to acetone only (solvent control). For each strain and treatment combination, 25 mosquitoes were used per technical replicate (separate exposure bottles from the same F_2_ cohort) per condition.

In Table 1, treatments for the San Nicolás strain are coded as SNB (SN exposed to broflanilide), SNI (SN exposed to imidacloprid) and SNC (SN exposed to acetone only). Treatments for the NO strain are coded as NOB (NO exposed to broflanilide), NOI (NO exposed to imidacloprid) and NOC (NO exposed to acetone only). Surviving mosquitoes from each condition were subsequently pooled for total RNA extraction and RNA-seq library preparation.

2.3. Insecticide Exposure Bioassays

Insecticide exposure bioassays were performed to obtain surviving mosquitoes after partial-mortality concentration (<50%) exposure to broflanilide or imidacloprid, with the aim of inducing sufficient toxicological stress to elicit transcriptomic changes, following the approach described by Trujillo et al. (2025) [12]. A total of 450 F_2_ adults females from each strain, 3–5 days post-emergence and fed only on 10% sucrose solution, were used in 250 mL glass bottles, with 25 mosquitoes per bottle and four replicate bottles per strain, plus a solvent control group.

Stock solutions of IMIDACLOPRID 18380 STD 10 mg (AccuStandard, CTR Scientific, Monterrey, Mexico) and BROFLANILIDE, 99% purity (Catalog No. E8093, Selleckchem, HO, USA), were prepared in molecular-grade acetone. Serial dilutions of each insecticide were tested in the NO strain to identify a partial-mortality concentration (<50%) that did not cause more than 50% mortality. Based on these preliminary assays, the internal surfaces of the bottles were coated with 1 mL of broflanilide solution (5 µg/mL) or imidacloprid solution (1 µg/mL). Control bottles were treated with acetone only. Mortality was recorded in the NO strain to confirm the partial-mortality concentration (<50%) effect. The same concentrations were then applied to the San Nicolás strain; LC_50_ values were not estimated for this population, as this was not the primary objective of the experiment. Using the same partial-mortality concentration (<50%) for both strains increased the number of surviving field mosquitoes, allowing the extraction of sufficient RNA.

Immediately after 1 h of exposure to broflanilide or imidacloprid, surviving mosquitoes were immobilized at −20 °C and then stored for total RNA extraction. For each experimental condition, groups of 12–15 surviving mosquitoes were pooled to generate each biological replicate for subsequent RNA-seq analyses.

2.4. RNA Extraction and Sequencing

For each strain–treatment combination, three biological replicates were generated by pooling approximately 12–15 surviving adult females per replicate, as described above. Total RNA was extracted from each pool using the TRIzol^®^ reagent protocol (Thermo Fisher Scientific, Waltham, MA, USA), following the manufacturer’s instructions. After extraction, RNA pellets were dried and stabilized in GenTegra^®^ RNAssure™ tubes (GenTegra LLC, Pleasanton, CA, USA) to prevent RNA degradation during storage and shipment for sequencing.

RNA quantity and integrity were assessed by the service provider, and RNA-seq libraries were prepared and sequenced by HaploX GeneTech (Hong Kong, China) using their transcriptome sequencing service (“library prep and PE150 sequencing, 6G package”) on an Illumina NovaSeq X Plus platform. Libraries were sequenced as 2 × 150 bp paired-end reads, targeting approximately 6 Gb of raw data per sample. After demultiplexing and initial quality control by the provider, FASTQ files were used for downstream mapping and quantification. Across the 18 libraries retained for analysis, between 9.2 and 33.5 million read pairs per library were assigned to genes, and all libraries met the planned depth and quality criteria.

2.5. Bioinformatic and Statistical Analyses

All bioinformatic and statistical analyses were performed in R [13] and relied primarily on the packages DESeq2 [14], tidyverse [15], ggplot2 [16], pheatmap [17], RColorBrewer [18], EnhancedVolcano [19], matrixStats [20], data.table [21], clusterProfiler [22] and patchwork [23].

2.5.1. Read Counting and Construction of the Gene Expression Matrix

Adapter- and quality-filtered RNA-seq reads were aligned against the Ae. aegypti reference genome (AaegL5; BioProject PRJNA318737). Gene-level read counts were obtained with featureCounts [24] using the corresponding GTF annotation, and only uniquely mapped reads were retained. The resulting featureCounts output files were imported into R and merged into a single count matrix, where rows represented genes (NCBI, GeneID) and columns represented biological replicates from each experimental condition.

For downstream analyses, condition codes for population–treatment combinations follow those defined in Table 1. Genes with very low expression were filtered out by requiring a minimum total count > 10 across all samples to reduce noise from sparsely expressed features.

2.5.2. Differential Expression Analysis

Differential expression analysis was carried out with DESeq2 [14]. The merged count matrix was used to construct a DESeqDataSet object with a design formula including a single factor, condition, representing the six population–treatment combinations. Size factors were estimated using the median-of-ratios method, and gene-wise dispersions were fitted following the default DESeq2 workflow. Wald tests were used to compare pairs of conditions, and Benjamini–Hochberg-adjusted p-values (padj) were reported.

We focused on the following biologically relevant contrasts: (i) insecticide versus control within each population (NOB vs. NOC, NOI vs. NOC, SNB vs. SNC, SNI vs. SNC) and (ii) population differences under specific treatments (SNC vs. NOC as the baseline difference between strains in the absence of insecticides, SNB vs. NOB under broflanilide, and SNI vs. NOI under imidacloprid). For each contrast, genes were considered differentially expressed (DEGs) when they showed padj < 0.05 and an absolute log_2_ fold change (|log_2_FC|) > 1. Lists of DEGs and full DESeq2 result tables were exported as comma-separated values (CSV) files for subsequent analyses and are included as Supplementary Material (Supplementary Material S1–S3).

To formally test whether insecticide responses differed between populations, we additionally fit a DESeq2 model including main effects of population and insecticide and their interaction (design = ~population + insecticide + population:insecticide). We then tested the interaction terms for broflanilide and imidacloprid (populationSanNicolas:insecticideBroflanilide and populationSanNicolas:insecticideImidacloprid). Genes with Benjamini–Hochberg-adjusted p-values (padj) < 0.05 for the interaction term were interpreted as showing statistically supported population-specific responses to the corresponding insecticide. Results were exported as Supplementary Tables (Supplementary Material S4 and S5).

2.5.3. Quality Control and Global Transcriptome Structure

To evaluate sample quality and global expression patterns, raw counts were normalized and transformed using the variance-stabilizing transformation (vst) implemented in DESeq2 [14]. Sample-to-sample Euclidean distances were computed from the vst matrix, and distance heatmaps were generated using pheatmap [17] with annotation tracks indicating population (NO vs. SN) and treatment (control, broflanilide, imidacloprid).

To highlight the most variable genes across all samples, row-wise variance was calculated using matrixStats [20], and the top 10 most variable genes were visualized as expression heatmaps (centered per gene) with pheatmap [17]. We deliberately restricted this visualization to 10 loci to keep the figure and interpretation tractable and to focus on the genes showing the strongest variance across samples, as an illustrative example rather than a comprehensive ranking. These quality control plots were used to verify clustering of biological replicates, detect potential outliers and characterize the relative contribution of population and insecticide treatment to the overall transcriptomic structure.

2.5.4. Visualization of Differential Expression

For each pairwise contrast, volcano plots were generated using the EnhancedVolcano package [19], displaying log_2_FC on the x-axis and −log_10_ (padj) on the y-axis. Genes passing the significance thresholds (padj < 0.05 and |log_2_FC| > 1) were highlighted, and the most strongly regulated loci were labeled by their NCBI GeneID. Multi-panel figures combining all contrasts were assembled using patchwork [23] to facilitate visual comparison of the magnitude and direction of transcriptional responses across populations and insecticides.

2.5.5. Gene Ontology Annotation and Enrichment Analysis

Functional enrichment analyses were based on Gene Ontology (GO) terms. A gene-to-GO mapping table for Ae. aegypti was constructed from the NCBI gene2go file, which was downloaded from the NCBI FTP site and imported into R with data. table [21]. Entries were filtered to retain only records with taxonomic identifier 7159 (Ae. aegypti), and separate TERM2GENE and TERM2NAME tables were generated for the three main GO ontologies: Biological Process (BP), Molecular Function (MF) and Cellular Component (CC).

For each DE contrast, we extracted the set of DEGs (padj < 0.05, |log_2_FC| > 1) with valid GO annotation. This combination of statistical significance threshold (padj < 0.05) and an effect-size threshold (|log_2_FC| > 1, i.e., at least a two-fold change in expression) was chosen to focus enrichment analyses on robust expression changes and to avoid GO terms being driven by minimal but statistically significant fold changes. We required at least 10 annotated genes per contrast to proceed.

The gene universe for enrichment corresponded to all genes tested by DESeq2 that also had GO annotation in the filtered gene2go table. Over-representation analyses for BP, MF and CC were performed using the enricher function from clusterProfiler [22], with the Benjamini–Hochberg procedure controlling the false discovery rate (q-value < 0.05). Enrichment results were exported as tab-separated files for each contrast and ontology.

Enriched GO terms were visualized as bar plots and comparative dot plots using ggplot2 [16] and RColorBrewer [18] for color palettes. Bar plots showed the top enriched terms ranked by adjusted p-value, while dot plots summarized enrichment across all contrasts, with dot color proportional to −log_10_(padj) and dot size reflecting the number of DEGs annotated to each term (gene count). These visualizations were used to identify core biological processes shared among treatments, as well as insecticide- and population-specific functional signatures.

3. Results

At the selected concentrations, mortality in the NO strain was approximately 52% for imidacloprid and 47% for broflanilide, whereas mortality in the SN strain was approximately 44% and 43%, respectively.

3.1. RNA-Seq Data Overview and Mapping

A total of 18 RNA-seq libraries were generated, corresponding to six strain–treatment combinations, with three biological replicates per condition. On average, each library yielded 21.3 million mapped read pairs (range: 9.2–33.5 million).

Across all samples, a total of 17,850 genes (NCBI GeneID) showed non-zero counts, of which 15,885 remained after filtering out lowly expressed genes (total count ≤ 10 across all samples). These filtered genes were used for all downstream normalization, differential expression, and functional enrichment analyses.

3.2. Global Sample Relationships and Quality

Pairwise Euclidean distances calculated from variance-stabilized counts were visualized as a sample-to-sample distance heatmap (Figure 1). Biological replicates of each condition (NOC, NOB, NOI, SNC, SNB, SNI) showed shorter distances to one another than to other samples, indicating good within-condition consistency and the absence of obvious outliers. When samples were ordered by population, libraries from New Orleans (NOC, NOB, NOI) formed a distinct block with relatively low within-block distances, and libraries from San Nicolás (SNC, SNB, SNI) formed a second block. Distances between these two blocks were consistently higher than those within each block.

This pattern indicates that population background exerts a major influence on global transcriptomic profiles, with additional variation associated with insecticide exposure superimposed onto this background. In other words, samples are most similar to those from the same population and treatment, and are more divergent from samples originating from the opposite population, providing a robust foundation for subsequent differential expression and functional enrichment analyses.

3.3. Differential Expression Within and Between Populations

Differential expression analysis revealed extensive transcriptional responses to both insecticides in the two Ae. aegypti populations (Figure 2). Because these analyses were conducted on surviving adult females after the 1 h exposure, the observed differences between exposed and control samples capture the net effect of both mortality-driven selection (removal of more susceptible individuals) and inducible transcriptional changes in the surviving subset. Within each population, contrasts of insecticide versus control (NOB vs. NOC, NOI vs. NOC, SNB vs. SNC, SNI vs. SNC) produced hundreds to thousands of DEGs at padj < 0.05 and |log_2_FC| > 1. In all four cases, volcano plots were dominated by genes with positive log_2_ fold changes, indicating a strong bias toward up-regulation after exposure to broflanilide or imidacloprid. Only a relatively small subset of loci showed significant down-regulation, suggesting that activation of transcriptional programs is the predominant response to insecticide challenge in both the New Orleans and San Nicolás populations.

Comparisons between populations highlighted substantial baseline and insecticide-specific divergence. Even in the absence of insecticide (SNC vs. NOC), we detected a large number of DEGs, again with a strong skew towards higher expression in one population, consistent with pronounced intrinsic transcriptomic differences between the San Nicolás and New Orleans strains. Under broflanilide exposure (SNB vs. NOB), however, almost no genes passed the significance thresholds, indicating that the global transcriptional response to this compound is largely conserved between populations. In contrast, imidacloprid exposure (SNI vs. NOI) yielded a substantial number of DEGs in both directions, revealing marked population-specific differences in the transcriptional response to this insecticide. Together, these patterns indicate that while both compounds strongly perturb gene expression, broflanilide elicits a broadly similar response in the two populations, whereas imidacloprid tends to amplify pre-existing transcriptional divergence between strains.

The DESeq2 interaction model identified a limited set of genes with statistically supported population-specific responses. Using FDR < 0.05 on the interaction terms, 88 genes showed a significant population × broflanilide interaction, whereas 20 genes showed a significant population × imidacloprid interaction (Supplementary Material S4 and S5). These results provide a formal statistical test for differential responses between populations, while leaving the main pairwise contrasts (insecticide vs. control within each population) as the primary framework for biological interpretation.

3.4. Expression Patterns of the Most Variable Genes Across Populations and Treatments

To illustrate exemplar transcriptional profiles, we visualized the 10 most variable genes across all samples using variance-stabilized, row-centered counts (Figure 3). For most of these loci, expression levels were consistently low in control samples and markedly higher in insecticide-exposed mosquitoes, particularly in the SNB and SNI conditions, highlighting strong induction after treatment. Several genes also showed pronounced population-specific baseline differences, with higher expression in New Orleans controls (NOC) than in San Nicolás controls (SNC), while remaining strongly up-regulated in both populations upon exposure. Biological replicates within each condition displayed highly similar expression patterns, reinforcing the robustness of the differential expression results and illustrating how population background and insecticide exposure jointly shape the transcriptional profiles of key responsive genes.

3.5. Functional Enrichment of DEGs

3.5.1. GO Biological Process Enrichment Across Contrasts

GO Biological Process enrichment revealed a shared core response to insecticide exposure together with compound- and population-specific features (Figure 4). “Proteolysis” and “transmembrane transport” were consistently enriched in all or most contrasts, indicating pervasive activation of protein turnover and membrane transporter remodeling in response to both broflanilide and imidacloprid. Imidacloprid treatments (NOI vs. NOC and SNI vs. SNC) additionally showed strong enrichment of “protein refolding” and related stress terms, consistent with engagement of proteostasis and unfolded-protein responses, whereas broflanilide exposure in San Nicolás (SNB vs. SNC) uniquely emphasized mitochondrial energy metabolism, including “proton motive force-driven ATP synthesis” and “tricarboxylic acid cycle.” Together, these results suggest that a common transcriptional core is overlaid by insecticide- and population-specific biological programs.

3.5.2. GO Cellular Component Enrichment Highlights Extracellular and Membrane-Associated Remodeling

At the Cellular Component level, enriched terms were dominated by extracellular and surface-associated structures (Figure 5). “Extracellular space” and “extracellular region” were significantly enriched in all contrasts, frequently with high gene counts, indicating that many DEGs encode secreted proteins or components of the extracellular milieu. Additional terms such as “extracellular matrix,” “external side of plasma membrane,” “chitin-based extracellular matrix,” “adherens junction” and “myosin complex” appeared in specific contrasts, suggesting remodeling of cuticle and extracellular matrix components, cell–cell adhesion structures and contractile machinery in response to insecticide exposure.

3.5.3. GO Molecular Function Enrichment Indicates Activation of Proteases, Detoxification Enzymes and Structural Proteins

GO Molecular Function enrichment revealed a coherent set of enzymatic and structural activities underlying the transcriptomic response (Figure 6). “Serine-type endopeptidase activity” and “serine-type endopeptidase inhibitor activity” were prominently enriched, particularly under imidacloprid exposure, highlighting activation of protease–antiprotease cascades. Detoxification-related functions, including “monooxygenase activity,” “heme binding,” “oxidoreductase activity” and “peroxidase activity,” were enriched in multiple contrasts, in line with cytochrome P450- and redox-based xenobiotic metabolism. Additional enriched functions such as “transmembrane transporter activity,” “symporter activity,” “chitin binding,” “structural constituent of chitin-based cuticle,” “actin binding” and “actin filament binding” indicate coordinated changes in membrane transport, cuticle and extracellular matrix composition and cytoskeletal organization.

3.5.4. GO Biological Process Enrichment in the New Orleans Population

In the New Orleans strain, both insecticides triggered a highly convergent transcriptional response at the Biological Process level (Figure 7). For broflanilide exposure (NOB vs. NOC), “proteolysis” was by far the most significantly enriched term, followed by “transmembrane transport” and several lipid-related categories, including “lipid metabolic process” and “fatty acid elongation” involving saturated and mono- and polyunsaturated fatty acids. Additional enriched processes such as “nervous system process,” “sodium ion transport,” “actin cytoskeleton organization” and “sphingolipid biosynthetic process” suggest that broflanilide impacts neuronal signaling, ion homeostasis and membrane/cytoskeletal remodeling.

Imidacloprid exposure (NOI vs. NOC) showed the same strong core of “proteolysis” and “transmembrane transport,” but was further characterized by enrichment of “protein refolding,” “cellular response to unfolded protein” and “tissue regeneration,” indicating prominent engagement of proteostasis and repair pathways. Terms related to ion and metabolite transport (“sodium ion transport,” “calcium ion transmembrane transport,” “lipid transport”) were also enriched, consistent with broad reorganization of membrane transport systems. Overall, these patterns indicate that in New Orleans, both insecticides elicit a shared response centered on protein turnover and membrane transport, with broflanilide more strongly associated with lipid and neuronal processes, and imidacloprid more strongly associated with proteostasis and tissue repair.

3.5.5. GO Biological Process Enrichment in the San Nicolás Population

In the San Nicolás strain, broflanilide and imidacloprid also induced distinct yet overlapping Biological Process profiles (Figure 8). For broflanilide exposure (SNB vs. SNC), the most significantly enriched terms were “proton motive force-driven ATP synthesis” and “tricarboxylic acid cycle,” followed by “mitotic cytokinesis,” “actin cytoskeleton organization,” “modulation of chemical synaptic transmission” and “mitochondrial electron transport, cytochrome c to oxygen.” These categories point to a strong mitochondrial energy response, increased cell division, and remodeling of the cytoskeleton and synaptic signaling, together with additional enrichment of “lipid metabolic process,” “chitin catabolic process” and “synaptic transmission, glutamatergic.”

Imidacloprid exposure in San Nicolás (SNI vs. SNC) showed a different emphasis. As in New Orleans, “proteolysis” and “transmembrane transport” were among the top enriched terms, but here they were accompanied by “tricarboxylic acid cycle,” “sodium ion transport,” “iron ion transport,” “carbohydrate metabolic process,” “carbohydrate transport,” “response to oxidative stress” and “monocarboxylic acid transport.” Together with enrichment of “actin cytoskeleton organization,” “chitin catabolic process” and “protein refolding,” these results suggest that imidacloprid in this population engages protein turnover, ion and metabolite transport, oxidative stress defenses and cuticle/extracellular matrix remodeling. Overall, broflanilide in San Nicolás is characterized by a pronounced mitochondrial and synaptic signature, whereas imidacloprid combines the proteolysis–transport core with oxidative stress and structural remodeling processes.

4. Discussion

In this study, we used RNA-seq to compare the transcriptomic response of a field-derived, pyrethroid-resistant strain of Ae. aegypti and a susceptible laboratory strain after partial-mortality concentration (<50%) exposure to two insecticides with distinct modes of action, imidacloprid and broflanilide. At the global level, our analyses revealed pronounced baseline transcriptomic differences between survivors of each strains; extensive transcriptional reprogramming after exposure to both insecticides, with a strong bias towards gene up-regulation; and contrasting patterns of population divergence between compounds, with broflanilide eliciting largely conserved responses and imidacloprid revealing marked population-specific transcriptional signatures. Functional enrichment analyses highlighted a shared core response involving proteolysis and transmembrane transport, overlaid with insecticide- and population-specific programs affecting energy metabolism, proteostasis, oxidative stress, cuticle remodeling and cytoskeletal organization.

The strong baseline differences we detected between the SN and NO strains, even in the absence of insecticide, are consistent with previous work documenting phenotypic divergence between laboratory and field Ae. aegypti populations. For example, López et al. [25] reported higher reproductive output and population growth in a laboratory-adapted colony compared with a wild population, indicating that colonization and long-term rearing can profoundly alter life-history traits. Vinauger and Chandrasegaran [26] further showed that genetic variability among Ae. aegypti populations is associated with differences in structural and behavioral traits, including host preference. In addition, differences in gut microbiota between field and laboratory mosquitoes have been linked to variation in phenotype and vector competence [27]. In the present work, our aim was exploratory and focused on describing global expression patterns and major functional pathways, rather than attributing specific DEGs primarily to microbiota or any other single factor. Nevertheless, it is plausible that part of the transcriptomic divergence we detect, particularly for genes related to immune function, digestion and epithelial or barrier tissue, reflects host responses to distinct microbial communities in the field-derived versus laboratory strain. Together, these studies support the idea that laboratory and field Ae. aegypti represent distinct biological entities at genetic, phenotypic and microbiological levels, which likely contributes to the extensive baseline transcriptomic divergence observed in our comparisons between NO and SN controls.

Our results also provide insight into the early transcriptional dynamics of broflanilide and imidacloprid from survivors. Despite strong insecticide-versus-control responses within each population’s survivors, almost no genes were differentially expressed between strains under broflanilide exposure (SNB vs. NOB), suggesting that the initial transcriptional response to this compound is largely conserved. This pattern is compatible with the known mode of action of broflanilide as a meta-diamide insecticide with delayed mortality, whose effects may unfold over several days [28]. Because our sampling was conducted 1 h after exposure, it is plausible that early-stage responses are dominated by a generic stress program that does not yet diverge strongly between populations; strain-specific trajectories might emerge at later time points. In contrast, imidacloprid exposure (SNI vs. NOI) revealed substantial transcriptional divergence between strains, indicating that this neonicotinoid not only triggers a robust response within each population but also amplifies pre-existing differences in gene regulation between them.

At the pathway level, GO enrichment results revealed a consistent core of Biological Process terms associated with “proteolysis” and “transmembrane transport” across most insecticide-versus-control contrasts, together with Molecular Function terms linked to protease activity, detoxification enzymes and transporters. These findings align with current models of insecticide resistance emphasizing transcriptional and post-transcriptional regulation of genes involved in xenobiotic metabolism and barrier functions [29]. The enrichment of proteolysis-related terms is compatible with increased protein turnover, processing of signaling molecules and remodeling of damaged proteins under toxic stress. Enrichment of transmembrane transport and symporter activity, together with functions such as monooxygenase activity, heme binding, oxidoreductase activity and peroxidase activity, points to coordinated regulation of detoxification pathways and membrane transporters that may modulate insecticide uptake, distribution or excretion [30,31]. Although fatty acid synthase and related lipid metabolic pathways have been implicated in development and stress responses in other insects [32], in our data, lipid-related processes were enriched only in specific contrasts, suggesting more context-dependent roles.

The comparison between NO and SN also revealed interesting differences in energy metabolism and stress responses. In the NO strain, both insecticides shared a strongly convergent Biological Process profile centered on proteolysis and transmembrane transport, with broflanilide more prominently associated with lipid metabolism and neuronal processes, and imidacloprid more prominently associated with protein refolding and tissue regeneration. In SN, broflanilide exposure was characterized by enrichment of mitochondrial energy metabolism (proton motive force-driven ATP synthesis, tricarboxylic acid cycle) together with actin cytoskeleton organization and synaptic transmission, whereas imidacloprid combined the proteolysis–transport core with additional processes related to oxidative stress, carbohydrate metabolism and ion transport. These patterns are consistent with a scenario in which both compounds engage a shared set of stress and detoxification pathways, while broflanilide elicits a stronger mitochondrial and synaptic signature and imidacloprid more strongly activates proteostasis, oxidative stress defenses and broader metabolic remodeling.

Inspection of the most variable genes across conditions provides further insight into the types of loci involved in these responses. Several highly variable transcripts corresponded to histidine-rich proteins or histidine-related functions. Histidine is not only an essential amino acid but also a precursor of histamine, a key neurotransmitter in insects [29,31,32]. The overexpression of histidine-rich or histidine-metabolism-related genes in imidacloprid-exposed mosquitoes may therefore reflect modulation of neuronal signaling cascades disrupted by this neonicotinoid’s action on nicotinic acetylcholine receptors. Although our data cannot establish a causal link, this pattern is compatible with previous proposals that amino acid metabolism, including histidine, can contribute to cellular responses to neurotoxic challenges [29,31].

Other highly expressed loci included genes related to actin and muscle function. For instance, actin genes expressed in the indirect flight muscle of Ae. aegypti are critical for myofibril formation and proper wing muscle development [31], and Myo-sex myosin heavy chain has been shown to be necessary for male flight [33,34]. Up-regulation of such muscle-related genes in insecticide-exposed groups, particularly in the field population, could reflect structural remodeling of flight musculature, increased locomotor activity or broader stress responses involving cytoskeletal reorganization. While it is tempting to interpret these changes as facilitating escape behavior from treated surfaces, such behavioral implications remain speculative without dedicated functional and behavioral assays.

We also detected elevated expression of cuticle-related genes, including cuticular protein 100A and chitin metabolism genes, particularly in SN under insecticide exposure. Cuticular proteins and chitin-based structures play important roles in cuticle formation, integrity and permeability, and have been implicated in penetration resistance to insecticides [35]. Their enrichment in our data is consistent with Cellular Component terms such as chitin-based extracellular matrix, extracellular matrix and external side of plasma membrane being over-represented in several contrasts [36]. These observations support the idea that structural remodeling of the cuticle and associated extracellular structures may be part of the response of Ae. aegypti to both broflanilide and imidacloprid, potentially contributing to reduced insecticide penetration or altered surface interactions.

Among enzymes involved in protein and peptide processing, we observed expression patterns consistent with leucine aminopeptidase activity, a zinc-dependent cytosolic enzyme that removes amino acids from the N-terminal end of proteins and peptides and has been implicated in protein metabolism, digestion, development and immune responses in insects [37,38]. In our data, some of these loci showed elevated expression not only in insecticide-exposed groups but also in at least one SN control replicate, suggesting a generally higher metabolic and/or immune tone in the field-derived strain. This reinforces the notion that wild populations may maintain a broader range of constitutive defenses compared with laboratory strains, which could influence both their intrinsic fitness and their responses to chemical stressors.

Taken together, our gene- and pathway-level analyses suggest that the field-derived SN population is more transcriptionally plastic than the laboratory NO strain, particularly in response to imidacloprid. These heightened changes may stem from prior exposure to heterogeneous environmental conditions and insecticide regimes, as well as differences in microbiota and life history [39,40]. From an applied perspective, the observation that broflanilide elicits largely conserved early responses in both populations, whereas imidacloprid accentuates strain-specific divergence, has potential implications for resistance management. Compounds that trigger strongly population-specific transcriptomic responses may be more likely to select for local adaptations and heterogeneous resistance trajectories, underscoring the importance of considering population background when designing insecticide rotation schemes.

Our study has several limitations that should be acknowledged. First, RNA-seq was performed only on mosquitoes that survived the 1 h exposure to a concentration causing partial mortality (<50%) in the susceptible strain. As a consequence, the observed differences between exposed survivors and unexposed controls may reflect a combination of selection (preferential survival of individuals with more tolerant baseline expression profiles) and acutely induce transcriptional changes in the surviving subset. With the current design, we cannot formally disentangle pre-existing tolerance-associated expression from inducible responses at the level of individual genes, and we therefore use the term “response” in a descriptive sense to refer to the transcriptomic profiles of surviving mosquitoes under each condition.

Second, we analyzed a single time point (1 h) after exposure. Many regulatory and physiological processes relevant to survival, recovery and resistance are likely to unfold over longer periods; thus, we may be capturing only the earliest phases of the response, particularly for broflanilide [41].

Third, although transcriptomic changes provide valuable hypotheses about candidate genes and pathways, they do not necessarily translate into functional resistance mechanisms. Functional validation through gene knockdown, overexpression or genome editing, as well as assays of insecticide uptake, metabolism and excretion, will be required to establish causal links between specific genes and resistance phenotypes.

Fourth, we did not explicitly integrate phenotypic, behavioral or microbiota data into our analyses; future work combining transcriptomics with life-history, behavioral, microbiological and vector competence assays would help clarify how the observed transcriptional changes impact mosquito biology and disease transmission.

Finally, GO enrichment analyses were based on fixed DEG thresholds (padj < 0.05 and |log_2_FC| > 1) without formal sensitivity analyses across alternative cutoffs. While this approach is widely used to focus on robust expression changes, different threshold choices could influence the statistical support of specific GO terms, and enrichment results should therefore be interpreted in this context.

Despite these limitations, our findings demonstrate that transcriptomic profiling of field and laboratory Ae. aegypti strains under controlled insecticide exposure can reveal both shared and population-specific components of the response to novel and repurposed compounds. By identifying core pathways such as proteolysis, transmembrane transport, detoxification, mitochondrial metabolism, cuticle remodeling and cytoskeletal organization, our study provides a foundation for mechanistic investigations of how neonicotinoids and meta-diamides affect mosquito physiology. Moreover, the marked divergence in imidacloprid responses between the SN and NO strains highlights the need to account for local population history and genetic background when implementing insecticide-based interventions. Integrating such population-specific molecular information into vector control programs may help anticipate and mitigate the evolution of resistance, contributing to more sustainable management of Ae. aegypti and the arboviral diseases it transmits.

5. Conclusions

In this study, we used RNA-seq to compare the early transcriptomic response of a field-derived, pyrethroid-resistant Ae. aegypti strain (San Nicolás) and a fully susceptible laboratory strain (New Orleans) after partial-mortality concentration (<50%) exposure to imidacloprid and broflanilide. We found pronounced baseline transcriptomic differences between strains and extensive activation of gene expression after insecticide exposure, with a strong bias towards up-regulation. A shared transcriptional core involving proteolysis, transmembrane transport, detoxification pathways and structural remodeling of cuticle and cytoskeleton was identified across contrasts.

Despite these common elements, the two insecticides differed in how they interacted with population background. Broflanilide elicited largely conserved early responses between strains, whereas imidacloprid amplified pre-existing divergence, producing marked population-specific transcriptional signatures. These findings suggest that field populations may display greater transcriptional changes than laboratory strains, particularly in response to neonicotinoids, and highlight the importance of incorporating population-specific molecular information when designing insecticide rotation schemes and resistance management strategies. Our results provide a transcriptomic baseline for future functional studies aiming to validate candidate genes and pathways involved in the response of Ae. aegypti to broflanilide and imidacloprid.