Transcriptomic meta-signatures identified in Anopheles gambiae populations reveal previously undetected insecticide resistance mechanisms

Abstract

Increasing insecticide resistance in malaria-transmitting vectors represents a public health threat, but underlying mechanisms are poorly understood. Here, a data integration approach is used to analyse transcriptomic data from comparisons of insecticide resistant and susceptible Anopheles populations from disparate geographical regions across the African continent. An unbiased, integrated analysis of this data confirms previously described resistance candidates but also identifies multiple novel genes involving alternative resistance mechanisms, including sequestration, and transcription factors regulating multiple downstream effector genes, which are validated by gene silencing. The integrated datasets can be interrogated with a bespoke Shiny R script, deployed as an interactive web-based application, that maps the expression of resistance candidates and identifies co-regulated transcripts that may give clues to the function of novel resistance-associated genes.

Fig. 1

Distribution of available microarray datasets. Microarray datasets available from sub-Saharan Africa comparing insecticide resistant and susceptible: An. gambiae (Ag), An. coluzzii (Ac) or An. arabiensis (Aa). Resistance levels are characterised by the populations maximal recorded mortality in WHO discriminating dose assays: high resistance 0–33% mortality, moderate resistance 33–66%, low resistance 66–90% mortality, susceptible populations are those that consistently exhibited 100% mortality after exposure (See Supplementary Table 1). The insecticides that the populations have been exposed to prior to RNA extraction are represented by D DDT, P Pyrethroid, B Bendiocarb, NK Not known; unexposed mosquitoes have no corresponding letter. Figure created expressly for this manuscript by Manuela Bernardi

Fig. 2

qPCR analysis of top insecticide resistance candidates from multiple resistant populations. qPCR results of three 3–5 day old, unexposed, pyrethroid resistant Anopheles populations compared to the lab susceptible N’Gousso, three biological replicates and three technical replicates were used for each gene. Relative fold change (y), and resistant populations (x). Standard deviation bars are shown, with significance *p ≤ 0.05, **p ≤ 0.01, ***p ≤ 0.001 as calculated by ANOVA with Dunnett’s post hoc test

Fig. 3

RNAi phenotyping for potential candidates. RNAi was performed using two pyrethroid resistant colonies, Tiassalé and VK7. 3–5 day old females mosquitoes were injected with dsRNA from each of the respective transcripts. 72-h post injection, the mosquitoes were exposed to 0.05% WHO deltamethrin papers for 1 h and mortality scored 24 h later. a Mortality associated with dsRNA induced knockdown in Tiassalé mosquitoes. b dsRNAs showing significantly increased mortality in Tiassalé mosquitoes were then injected into a second resistant colony, VK7. c Knockdown levels for each dsRNA construct relative to GFP-injected controls (Tiassalé population only). Error bars represent standard deviation, three biological replicates and three technical replicates were used for each gene with each circle representing the mean value for each biological replicate in the qPCR data and significance is indicated by *p ≤ 0.05, **p ≤ 0.01, ***p ≤ 0.001, as computed by ANOVA with post hoc Tukey correction. Numbers on the bars represent number of mosquitoes tested under each condition

Fig. 4

RNAi phenotyping for potential transcription factor candidates. RNAi was performed using the pyrethroid resistant Tiassalé colony. Three to five-day-old females mosquitoes were injected with dsRNA from each of the respective transcripts. Seventy two-hours post injection, the mosquitoes were exposed to 0.05% WHO deltamethrin papers for 1 h and mortality scored 24 h later. a Deltamethrin induced mortality following dsRNA induced knockdown of the transcription factors Met and Dm. b Knockdown levels for each dsRNA construct relative to GFP-injected controls. Error bars represent standard deviation, three biological replicates and three associated technical replicates were used for each gene in the qPCR data and significance is indicated by *p ≤ 0.05, **p ≤ 0.01, ***p ≤ 0.001, as computed by ANOVA with post hoc Tukey correction. Numbers on the bars represent number of mosquitoes tested under each condition

Fig. 5

Graphical Outputs from IR-TEx App. These panels represent three graphical outputs from the IR-TEx App, each of which is available in downloadable tabular formats. These outputs result from a user inputted VectorBase ID, selected filtering criteria based on meta-data and finally, user inputted correlation value. Each panel here is AGAP001076-RA (CYP4G16). a Log₂ fold change of CYP4G16 (y) across each microarray experiment meeting pre-selected user criteria (x), here all data is selected. b Each point represents a dataset with significant differential expression of CYP4G16 and the associated approximate collection site. Green points show significant down regulation of CYP4G16, orange points show significant fold changes of 1-5 and red a fold change of >5 compared to susceptible controls. Map data source: Google Maps, 2018. c Log₂ fold change (y) across each microarray experiment meeting pre-selected user criteria (x) for each transcript showing a correlation of |r| ≥ 0.85 for CYP4G16

Fig. 6

Schematic of cuticular hydrocarbon synthesis in An. gambiae. Schematic showing the steps in cuticular hydrocarbon synthesis, adapted from Blomquist 1987 and Balabanidou et al., together with a list of putative members of each gene family in An. gambiae. All gene families are shown in blue and represent: FAS fatty acid synthases, elongases fatty acid elongases, reductases fatty acid reductases, desaturases fatty acid desaturases. General chemical classes at each step are shown in light red hexagons. Highlighted in green circles are the transcripts present in the correlation network of CYP4G16. The correlation network is significantly enriched for this pathway p_{(hypergeometric test)} = 1.44e⁻¹². Figure created expressly for this manuscript by Manuela Bernardi