Principle Component Analysis Highlights Major Divergence in 24-Hour Post-Permethrin Treatment Expression Profile:
To assess transcriptomic changes across the first 24 hours after exposure to permethrin, 35 libraries were created from samples collected before exposure, and 6, 10, and 24 hours after exposure (summarized in Fig. 1a). These time points were chosen based on a microarray study in Anopheles gambiae that showed significant changes in response at these intervals10. All 35 libraries produced high quality 3’ Tag-Seq11 single end reads, with an average library size of 4,092,249 reads (min:3,493,811—unexposed replicate 2, max: 4,640,468—24 hrs post permethrin exposure replicate 3). Across samples, 96.7% of reads mapped to the reference genome at least once on average. Gene annotations were extracted from Vectorbase and genes with an “unspecified product” annotation were searched with BlastP against a custom database of the Culicidae and Drosophila databases to search for orthologs of unannotated genes. Of the 19,804 genes annotated in the Aedes aegypti genome, 11,155 (56.3%) were expressed in this dataset.
To assess how samples grouped based on overall expression profile, a principal component analysis was performed to ensure grouping was consistent with treatment conditions. PC1 explained 35.4% of the variation, while PC2 explained 16.4% of the variation (Fig. 1b). Because only 51.8% of the variation in the dataset was explained by PC1 and PC2, further clustering analyses were performed.
To further define the differences and similarities between groups, a kmeans cluster analysis was performed and the gap statistic was calculated to determine the optimal number of clusters within the data set12. Through this analysis, it was determined that the data could be categorized into 4 clusters as follows: baseline and 24 hours post-acetone treatment, 6- and 10-hours post-acetone treatment, 6- and 10-hours post permethrin treatment, and 24 hours post permethrin treatment (Fig. 1B). Based on the clustering, the acetone exposed insects appear to return to baseline 24 hours after exposure, indicating that the use of this reagent as a control in insecticide resistance testing experiments is adequate.
Because replicate samples from the 24 hours post-permethrin treatment demonstrate the greatest variance relative to the other treatment groups, we further investigated the genes driving this separation using a biplot. Biplot analysis (Supp. Figure 1), reveals that the gene most strongly influencing this treatment group is catalase (AAEL013407), which codes for a protein that mediates oxidative stress by catalyzing the breakdown of hydrogen peroxide into two water molecules. Catalase has undergone thorough investigation in mosquitoes and is associated with multiple physiological processes incluing insecticide resistance, fecundity, pathogen infection, diapause, and maintenance of homeostasis after blood feeding13–16. All of these processes result in high levels of reactive oxygen species (ROS) production, driving increased catalase expression to neutralize these molecules and mediate oxidative stress. Some studies in Anopheles species have found that insecticide resistant individuals live longer than their susceptible counterparts, likely due to an increase in catalase, though this phenotype seems to be strain specific17–20. Further study on the role of catalase in the strain specific longevity trends observed in resistant populations would improve our understanding of life-history tradeoffs due to resistance. For this reason, exploring chemical agents that exploit ROS production or knock down catalase activity could be viable for mosquito control.
Redox Homeostasis Processes Exhibit Significant Expression Changes Over Time After Insecticide Exposure:
We next characterized the time course expression response to pyrethroid exposure to holistically assess differential expression. This analysis provides information on the relationship between time points and expression values, and assess which genes are experiencing significant changes in expression across time. To accomplish this, we used a cubic regression spline curve with 3 degrees of freedom to assess expression trends across the time course. The analysis revealed significant responses over time after insecticide exposure with 383 genes upregulated and 200 downregulated.
We performed a gene ontology (GO) analysis on the upregulated and downregulated genes to identify enriched functions within these groups. Of the 383 upregulated genes, we found enrichment of 5 primary GO categories. Among these are acetyl-CoA biosynthesis, glucose metabolism, lipid metabolism, redox homeostasis processes, and ATP production (Fig. 2).
Redox homeostasis processes are the primary response component to many stressors, including insecticides. Pyrethroid exposure causes oxidative stress, resulting from an imbalance of reactive oxygen species (ROS) and enzymes capable of neutralizing them21–24. Free radicals are formed by many cellular processes including the electron transport chain, fatty acid beta oxidation, and many other processes that require the breaking of molecular bonds such as the breakdown of xenobiotics. Gene families involved in the neutralization of free radicals to maintain redox homeostasis are also involved in the breakdown of xenobiotics, these include CYPs and GSTs22,25. Expression of these genes is often induced by increases in ROS26. ROS can cause conformational changes that result in the release of a cofactor from a transcription factor allowing for nuclear localization of said transcription factor. One example is the release of the KEAP-1 protein from cap n collar transcription factor27.
In Ae. aegypti, the ortholog of the Drosophila melanogaster cap n’ collar gene is AAEL019563. The cap n’collar C isoform (CncC) is the invertebrate equivalent of the vertebrate nuclear factor erythroid 2-related factor (Nrf2), which is known to bind antioxidant response elements28. CncC binding to these promoters leads to the expression of the previously mentioned gene families under insecticide stress29. Importantly, CncC has been found to be constitutively active in resistant populations of Drosophila and Anopheles gambiae30,31. Studies have also found that this transcription factor is responsible for metabolic shifts in response to stress, which are illustrated by this data set as well32,33. Interestingly, in this dataset, AAEL019653 does not display a significant change until 24 hours after exposure, when it is downregulated (logFC=-1.04, FDR = 0.00019). It is unsurprising that this gene does not experience any upregulation throughout the time course as it is a constitutively expressed transcription factor which is released from its repressor upon increases in ROS within the cell.
Both lipid and carbohydrate metabolism were upregulated, along with acetyl-CoA biosynthesis. This is consistent with results in many other studies on the response to insecticide exposure10,23,34. Additional evidence supports the role of CncC in increased lipid and carbohydrate metabolism, as constitutive activation of this protein leads to a decrease in lipid stores in Drosophila and mouse cells showed increased glucose uptake in the same condition33,35. Acetyl-CoA is derived from the catabolism of proteins, lipids, and carbohydrates. Acetyl-CoA is a primary monitor of the metabolic state of an organism. When acetyl-CoA is plentiful, it is shuttled into lipid synthesis and while in a depleted state, it is transported into the mitochondria for ATP synthesis36. Other processes involved in ATP production include genes involved in NAD binding and FAD processes. Both of these molecules are high energy electron acceptors in the electron transport chain. This paired with ATP transmembrane transport overrepresentation clarifies the energetic stress the insect is under after exposure to a xenobiotic.
Heat Shock Proteins Exhibit the Largest Fold Changes Among Upregulated Genes Across Time
Eight of the top 20 genes with the largest fold change belong to the heat shock protein (HSP) family (Fig. 3). Each of the HSPs displays a similar patterns of transcript abundance over time, with an increase through 10 hours followed by a gradual decline. Acetone treated samples show a response in these gene as well, though expression levels return to baseline at 24 hours while permethrin treated samples still have elevated expression. HSPs are associated with response to environmental stressors and have been associated with insecticide exposure response in many insects37–42. HSPs have a wide array of functions, aiding in all aspects of protein processing while not acting as a part of the mature protein43. One of these functions is as chaperones to aid in initial protein folding or to repair damage following stress. The role these molecules play in insecticide response may be due to protein damage due to ROS. One study found protein oxidation products increase following inoculation with the larvicidal bacteria Bacillus thuringiensis kurstaki, however, direct evidence of these products has not been studied in adult mosquitoes following insecticide exposure44.
In a warming climate, the interactions between heat stress and insecticide resistance have the potential to pose major issues for vector control professionals45. There is even evidence that heat tolerance can protect from insecticide exposure and vice versa37,42. Continued research on the development of cross-resistance between heat and pyrethroids will be important for future vector control efforts.
Time Point Specific Pairwise Comparisons Identify 7 Detoxifying Genes Upregulated At All Time Points
To investigate time point specific detoxification signatures, we employed pairwise comparisons at each time point. 1245 (11.16%) were differentially expressed in at least 1 time point, 645 upregulated and 599 downregulated, using an FDR cut off of 0.05. Thirty-eight of the 645 upregulated genes were upregulated at all 3 time points (5.9%) while only 1 gene (a putative oxidase/peroxidase) was downregulated at all 3 time points (AAEL019639, FDR = 7.97x10− 9) (Fig. 4B).
Gene ontology analysis was used to explore functional characteristics of genes upregulated at all time points. Due to the small number of genes, the classic fisher algorithm was used which tests GO terms independently and is therefore less conservative than the TopGO default algorithm. Genes shared across time points fall into four primary categories: lipid metabolism, detoxifying enzymes, purine breakdown, and FAD binding (Fig. 4A). Among the detoxifying enzymes are CYP6M11, CYP4D38, esterase B1, CYP4C38, GPHX1, a microsomal GST and an ABC transporter. These are part of different stages of detoxification. It could be speculated that these detoxifying enzymes exhibit a more generalized function, aiding in the breakdown and excretion of various byproducts. Alternatively, the process of insecticide detoxification in insects might be so slow that they continue to sequester insecticides and break down entire molecules, even 24 hours after exposure. This could be an important consideration for vector control professionals who routinely perform insecticide resistance testing, as the time at which they are observing knockdown may not be an accurate depiction of resistance state.
Fatty Acid Beta Oxidation and ATP Production are Predominant Processes at 6-hours Post-Exposure
Six hours post exposure, 95 genes were upregulated while 15 genes were downregulated. Processes overrepresented at this time point are related to fatty acid beta oxidation and ATP production and transport. The most significantly upregulated gene is that coding for the protein carnitine O-palmitoyl transferase (whd) (AAEL005458, FDR = 2.36 x10− 5) which is involved in fatty acid beta oxidation46. Of the 95 upregulated genes, 34 (35.7%) are only upregulated at 6 hours, one of which is the transcription factor, Hnf4 (AAEL011323, FDR = 0.028). This factor is associated with a metabolic switch to oxidative phosphorylation in Drosophila and fatty acid beta oxidation in Ae. aegypti47,48. Additionally, it is likely that the withered gene (whd) is activated by Hnf4, as researchers found that the expression of whd significantly decreased upon RNAi knockdown of Hnf448.
Downregulated genes were more specific to the 6-hour time point, with 12 of the 15 only showing reduced abundance at 6 hours. However, there was no significant overrepresentation of function among this group based on a gene ontology analysis. An interesting gene within this group is a cytochrome P450, CYP9J22 (AAEL014619, FDR = 0.017), as the CYP9J family is implicated in insecticide detoxification. Because it is downregulated so early in the time course, it may be one of the first CYPs to respond to permethrin followed by a drastic downregulation, or it may just be an example of strain specific detoxifying gene responses7. The 6-hour time point shared 15 upregulated genes with the 10-hour time point, though this group did not contain any significantly overrepresented GO terms. However, among these genes is whd, indicating a continuation of increased fatty acid beta oxidation into 10 hours post exposure.
Carbohydrate Metabolism Elevates Ten Hours Post-Exposure and Sustains Through 24 Hours:
Gene expression changes at 10-hours post-exposure represent a switch in energy source as well as the peak of expression for genes displaying differential upregulation at all time points. At ten hours post exposure, 219 genes were upregulated, while 77 genes were downregulated. GO term overrepresentation among upregulated genes shows an increase in carbohydrate metabolism, with a continuation of fatty acid beta oxidation as well. This may indicate that mosquitoes have burned through much of their lipid stores by 10 hours post exposure, leaving them to function off of carbohydrates derived from sugar feeding. The most significant upregulated gene is that coding for glycine N-methyltransferase (GNMT) (AAEL012764, FDR = 6.31x10− 7). This gene has not been previously associated with insecticide response to our knowledge, however the role of the mouse ortholog in stress response has been investigated. Knockout of GNMT in mice resulted in a reduction in expression of detoxifying genes including CYPs, GSTs, catalase, and superoxide dismutase, as well as an increase in lipid peroxidation products49. Further study of the role of this gene in insecticide response may provide additional information on resistance mechanisms.
Eighty-one (37.0%) upregulated genes are only upregulated at 10 hours, though no significant GO term overrepresentation was found. Among these are 4 detoxifying genes previously associated with insecticide exposure: GSTD6, CYP9J9, CYP304B2, and CYP304C1, and 1 gene previously studied in response to insecticide exposure with negative results: GSTT129,50–53. Differing results in our study may be due to observation of expression changes over a longer span of time or merely strain specific differences. Forty-five of the 77 (58.4%) downregulated genes were only downregulated at 10 hours. GSTT4, a detoxifying gene previously linked to pyrethroid exposure, is among this group. Although it was found to be overexpressed in a study where GSTT1 was deemed insignificant53. This suggests that the role of individual GSTs in the detoxification response may vary between strains in their response to pyrethroid exposure.
Ten hours post exposure shares 85 upregulated genes with 24 hours post exposure, many of which are involved in oxidoreductase activity (among which are catalase, CYP9J6, CYP6BB2, and a CYP with no family assigned), carbohydrate metabolism, and organic acid metabolism. CYP9J6 and CYP6BB2 are well associated with pyrethroid response50,54,55. In some strains of Ae. aegypti, CYP6BB2 is increased in copy number, and has been specifically shown to aid in metabolism of permethrin56,57. Additionally, 9 HSP genes or genes with HSP-like binding regions are upregulated, consistent with results from the time course analysis.
Downregulated Genes at 24 Hours Post-Exposure are Tied to Energy Consumption:
Expression signatures 24-hours post-exposure were investigated to further assess broad differences creating the clustering divide observed in Fig. 1b. Twenty-four hours post exposure, 371 genes are upregulated while 476 genes are downregulated. The most significant upregulated gene was a cytochrome P450, possibly part of the CYP6B family, and experienced a 2-fold change (AAEL009018, FDR = 1.42x10− 7)58. This gene has been associated with both pyrethroid and bendiocarb resistance58,59. The specificity of its differential expression solely at 24 hours after exposure to permethrin in this strain is unknown.
The upregulated genes at this time point largely mirror the upregulated genes throughout the time course, however the downregulated genes are unique to 24 hours, with 445 of the 476 downregulated genes only downregulated at 24 hours. The GO hits fall within 6 categories: ATP binding, GTPase activity, protein processing, signal transduction, spermidine biosynthesis, and transcription (Fig. 5). Spermidine is associated with protection against the toxic effects of pyrethroids in zebrafish, though the downregulation of genes involved in its synthesis may indicate that the insect is no longer experiencing toxic effects at this time60. Rather, the insect is primarily breaking down byproducts of earlier metabolism of the pyrethroid molecules. Decreases in ATP binding, signal transduction, protein processing, GTPase activity, and transcription indicate an overall decrease in the molecular responses to pyrethroid exposure. This combination of activity suggests that the insect is ramping down their overall response to return to a basal physiological state.
Thirty-eight detoxifying enzymes exhibit an association with pyrethroid response
To further classify those xenobiotic response-associated genes that may not have exhibited a statistically significant change in expression, we conducted a weighted gene correlation network analysis utilizing the R package WGCNA .61 This analysis groups genes based on expression patterns over time as well as by their correlation with treatment. Thirty modules were identified with one highly correlated with treatment (Supp. Figure 3). This module contains 931 genes, 305 of which are differentially expressed using an FDR cutoff of 0.05 from the time course analysis (245 upregulated, 60 downregulated). Interestingly, no correlation with time was observed in this module. This suggests that the expression changes witnessed over this time course are attributed to the xenobiotic response, rather than merely to time-associated expression alterations. Forty-two genes within the oxidoreductase GO term are significant based on the time course analysis within this group (Fig. 6A). Oxidoreductase activity is mediated by many enzymes associated with detoxification functions.
Within the module most related to treatment are 25 cytochrome P450s (without filtering based on FDR). Eleven of these are significant when considering expression over time while 17 are significant in at least 1 time point (Fig. 6B). Expression patterns of these genes are inconsistent, but cluster into 5 groups based on these patterns. One CYP is significantly downregulated at 10- and 24-hours post exposure, CYP307A1. This gene is included in the Halloween gene cluster and is the ortholog of the Drosophila gene Spook. Spook is part of ecdysone synthesis, so this downregulation may have implications for the role of ecdysone metabolism in the insecticide stress response62,63.
We also found members of the CYP9J family is which is associated with pyrethroid resistance7, 64–67. Four genes from this family are present in the treatment associated module, 3 of which cluster together based on expression pattern (CYP9J6, CYP9J9, and CYP9J31). In prior work, researchers performed an extensive micro-array study on rhythmic gene expression in female Ae. aegypti heads with a special focus on the CYP9J family and found two groups of co-oscillating members of this family68. However, none of the genes explored in that study were differentially expressed in this study.
Nine GSTs/GPXs are present in the treatment-associated module, 5 of which are significant. Expression trends are more consistent within this group, with expression peaking at 10 hours post exposure (Fig. 6B). GSTX2 is associated with DDT resistance in South America, however researchers in Thailand performed in-vitro experiments with this protein derived from a Thai population and found no evidence of direct metabolism of DDT by this enzyme69,70. GPXH1 appears to be induced by xenobiotic exposure, however, overexpression of this gene is not associated with resistance71.
Four ABC transporters are in this group, 3 of which are upregulated. ABC transporters are responsible for transporting insecticide and byproducts of insecticide metabolism out of the cell to avoid accumulation of these toxic compounds72. Transcription of these are also initiated by CncC, and their role in xenobiotic detoxification has been validated by RNAi experiments in crop pests73.