Age and prior blood feeding of Anopheles gambiae influences their susceptibility and gene expression patterns to ivermectin-containing blood meals

Abstract

Background

Ivermectin has been proposed as a novel malaria transmission control tool based on its insecticidal properties and unique route of acquisition through human blood. To maximize ivermectin’s effect and identify potential resistance/tolerance mechanisms, it is important to understand its effect on mosquito physiology and potential to shift mosquito population age-structure. We therefore investigated ivermectin susceptibility and gene expression changes in several age groups of female Anopheles gambiae mosquitoes.

Methods

The effect of aging on ivermectin susceptibility was analyzed in three age groups (2, 6, and 14-days) of colonized female Anopheles gambiae mosquitoes using standard survivorship assays. Gene expression patterns were then analyzed by transcriptome sequencing on an Illumina HiSeq 2500 platform. RT-qPCR was used to validate transcriptional changes and also to examine expression in a different, colonized strain and in wild mosquitoes, both of which blood fed naturally on an ivermectin-treated person.

Results

Mosquitoes of different ages and blood meal history died at different frequencies after ingesting ivermectin. Mortality was lowest in 2-day old mosquitoes exposed on their first blood meal and highest in 6-day old mosquitoes exposed on their second blood meal. Twenty-four hours following ivermectin ingestion, 101 and 187 genes were differentially-expressed relative to control blood-fed, in 2 and 6-day groups, respectively. Transcription patterns of select genes were similar in membrane-fed, colonized, and naturally-fed wild vectors. Transcripts from several unexpected functional classes were highly up-regulated, including Niemann-Pick Type C (NPC) genes, peritrophic matrix-associated genes, and immune-response genes, and these exhibited different transcription patterns between age groups, which may explain the observed susceptibility differences. Niemann-Pick Type 2 genes were the most highly up-regulated transcripts after ivermectin ingestion (up to 160 fold) and comparing phylogeny to transcriptional patterns revealed that NPCs have rapidly evolved and separate members respond to either blood meals or to ivermectin.

Conclusion

We present evidence of increased ivermectin susceptibility in older An. gambiae mosquitoes that had previously bloodfed. Differential expression analysis suggests complex midgut interactions resulting from ivermectin ingestion that likely involve blood meal digestion physiological responses, midgut microflora, and innate immune responses. Thus, the transcription of certain gene families is consistently affected by ivermectin ingestion, and may provide important clues to ivermectin’s broad effects on malaria vectors. These findings contribute to the growing understanding of ivermectin’s potential as a transmission control tool.

Background

Anopheles gambiae is one of the most significant malaria vectors in sub-Saharan Africa. Current vector control methods predominantly rely on insecticide-based tools, including indoor residual spraying (IRS), and long-lasting insecticide-treated bed-nets (LLINs) [1]. The four classes of vector control insecticides: pyrethroids, organophosphates, organochlorines, and carbamates, have all been co-opted from agricultural use. So while highly successful, the combination of widespread agriculture and vector control insecticide applications has fostered genetic and behavioral resistance among An. gambiae and other malaria vector species [2]. Insecticide resistance is threatening the recent gains in controlling malaria worldwide, especially in Africa, highlighting the need for novel anti-vector strategies and application methods.

Ivermectin (IVM) is a semi-synthetic member of the avermectin drug family, a series of 16-membered lactone endectocides first isolated from Streptomyces avermetilis in 1975 [3]. It is commonly used in humans, livestock and pets to control endoparasites, including parasitic helminths such as Onchocerca volvulus [4, 5], Haemonchus contortus [6] and Dirofilaria immitus [7], and ectoparasites including scabies mites [8], bot flies and ticks [9, 10]. Due to IVM’s strong safety profile in humans [11] and the fact that adverse reactions are rare and generally mild [12, 13], it is distributed to entire communities in mass drug administrations (MDA) once or twice per year for control of onchocerciasis and lymphatic filariasis in Africa. Following from its activity against ectoparasite arthropods, IVM has been studied for its ability to cause mortality in multiple arthropod disease vectors that feed upon, but do not infest, their hosts such as tsetse flies, sandflies and mosquitoes [14–16]. IVM is particularly potent against Anopheline mosquitoes, including An. gambiae, causing significant mortality, delayed re-feeding, and reduced fertility when ingested via a blood meal [17–23]. MDA of IVM in West African villages was shown to significantly reduce An. gambiae survivorship and sporozoite rates in field-captured mosquitoes [22, 24, 25]. Clinical trials with repeated IVM MDA are underway to examine the potential to achieve sustained malaria control in communities.

As an anti-malaria vector tool, IVM has a unique administration route in that it is ingested through a blood meal rather than physical contact. In addition, IVM has a novel mode of action (MOA) for malaria vectors by acting as an agonist of glutamate-gated chloride channels [26], causing flaccid paralysis and eventual death. Due to its short pharmacokinetic persistence, many vectors will ingest a sub-lethal dose if they bite people more than 2 days post treatment [21]. To maximize its effectiveness, it is important to understand the gene expression and other molecular events in the vector following IVM ingestion to further characterize the drug’s effect on mosquito physiology, and the potential for mosquito resistance. Similar studies examining mosquitoes after exposure to traditional insecticides have mainly identified transcriptional changes to cuticular proteins, transporters, mitochondrial respiratory mechanisms, and detoxification processes. [27, 28]. Transcriptional responses to avermectin exposure have been studied in Pediculus humanus humanus [29] and Rhipicephalus microplus [30], with both studies observing elevated transcription of ATP-binding cassette (ABC) transporters. Drosophila have been selected for resistance against IVM and abamectin in two different studies, with target-site insensitivity attributed to IVM resistance [31] and increased expression of ABC transporters attributed to abamectin resistance [32].

Evaluation of a compound’s potential as a transmission control tool should involve characterization of the molecular events surrounding vector exposure. Specifically, we wanted to determine whether gene expression patterns in An. gambiae reflect up-regulation of classical insecticide detoxification genes in response to IVM or whether novel gene classes may be up or down regulated. Additionally, because of the relatively long extrinsic incubation period associated with malaria transmission and the lethal effect of IVM on vectors, it is important to understand IVM’s ability to affect mosquito population age structure and mosquito vectorial capacity. To this end, we first compared IVM susceptibility among different age groups of female An. gambiae. Next, RNA-seq technology was utilized to assess transcriptional responses following IVM ingestion within and between these age groups. Expression patterns were validated with quantitative PCR on individual samples and on naturally-exposed mosquitoes that blood fed on a treated person.

Methods

Mosquito colony

An. gambiae s.s (G3 strain) were reared at the Arthropod-Borne and Infectious Diseases Laboratory at Colorado State University (CSU). This colony was kept at a constant 28–29 °C and 70–80 % relative humidity on a 14:10 light:dark cycle. Larvae were hatched in room temperature water, provided ground TetraMin® fish food daily, and were allowed to pupate and emerge into an enclosed mesh-covered cage where they were provided water and 10 % sucrose ad libitum.

Mosquito bioassay and feeding schedule

Three age groups of mosquitoes were reared and assayed for IVM-induced mortality; 2 days-post-emergence (DPE), 6DPE and 14DPE. Each group was split into treatment and control groups, receiving an IVM or control blood meal at the specified age. Feeding was done using defibrinated calf blood (Colorado Serum Company, Denver, CO) containing 11.75 ng/ml IVM (10 mg/ml stock solubilized in dimethyl sulfoxide (DMSO) and diluted in phosphate buffered saline (PBS)). This concentration was chosen based on the human pharmacokinetic curve for IVM, corresponding approximately to the concentration found in human blood 36 hours post ingestion and representing a sub-lethal dose (<LC₅₀) to An. gambiae [21]. Control blood meals contained the same volume of DMSO alone in PBS. Blood feeding was performed using glass bell feeders (Lillie Glass Feeders, Smyrna, GA) with membranes of hog sausage casings, and held at 37 °C. Following feeding, mosquitoes were knocked down at 4 °C, sorted on chilled plates, and only fully engorged females were kept. For the 6DPE and 14DPE groups, mosquitoes were fed normal blood meals and allowed to oviposit as they aged to emulate natural conditions; only those that engorged on these normal blood meals were carried into later age classes to maintain intergroup similarities. Thus, for 2DPE groups, the IVM blood meal (or corresponding control) they ingested on day 2 was their first and only blood meal. 6DPE groups received one normal blood meal on day 2, allowed to oviposit on day 4, and then were given their IVM-containing blood meal (or corresponding control) on day 6. 14DPE groups received three normal blood meals on days 2, 6, and 10, allowed to oviposit 2 days after each, and then were given their IVM or control blood meal on day 14. Following treatment blood meals, mortality was recorded daily for five days. A total of five biological replicates were performed for each age group.

Statistical analysis of bioassays

Survival analysis of bioassays was performed on pooled replicates for each age (Log-rank (Mantel-Cox) test) to test for significant mortality from IVM. Comparison of mortality between age groups was analyzed using the Mantel-Haenszel hazard ratio with 95 % confidence intervals to account for control mortality (Prism 5; GraphPad Software, Inc.).

RNA-seq library preparation

Mosquito bioassay feeding schedules were repeated exactly as above using 2DPE and 6DPE mosquito groups fed on control or IVM-treated blood. Twenty-four hours after blood-feeding, 10 fully engorged females were flash-frozen and stored at −80 °C awaiting RNA extraction. Groups of 10 unfed 2DPE and 6DPE females were also collected and flash-frozen at the same time. Two biological replicates were performed at each age, resulting in 20 blood-fed individuals per treatment. RNA was isolated from each individual sample using a Trizol® (phenol-chloroform) extraction procedure and genomic DNA was removed using Turbo™ DNase (Ambion® Life Technologies, Grand Island, NY). Following isolation, RNA was analyzed with an Agilent 2100 Bioanalyzer to visually check for RNA degradation. Individual mosquito RNA was then quantified with a Qubit® 2.0 Fluorometer (Life Technologies, Grand Island, NY) and an equal amount of RNA (200 ng) from each individual was used to pool in groups of ten. Pooled RNA was then subjected to TruSeq Stranded mRNA processing (Illumina, Inc., San Diego, CA). This involved mRNA purification by poly-T oligo attached magnetic beads, fragmentation, cDNA synthesis, adapter ligation, and amplification, resulting in 12 barcoded cDNA libraries. Samples were sent to Beckman Coulter Genomics (Boston MA), where they were sequenced on an Illumina HiSeq 2500 to produce 338 million 2x100-bp paired-end reads. The total number of sequences obtained for each library can be found in Additional file 1: Table S1.

Sequencing data analysis

Data analysis was performed on a Linux platform using GMAP and GSNAP [33, 34] for indexing and alignment, respectively, and Trimmomatic [35] for quality control and removal of singletons. Reads were aligned to the AgamP3.22 genome build (VectorBase, http://www.vectorbase.org, Anopheles gambiae PEST, AgamP3.22) [36]. Low expression tags were filtered using a threshold of 100 counts per million (CPM). Statistical analysis was performed using the edgeR package [37] for R software [38]. Counts were first fit to a negative binomial generalized log-linear model (glmFit procedure) and the dispersion was estimated by plotting the Biological Coefficient of Variation (BCV) for each biological replicate (Additional file 1: Figure S1). Differential expression (DE) analysis was then performed using the glmLRT procedure. Significant DE genes were identified using the exact test for two groups of negative-binomial counts with α = 0.05 and the Benjamini-Hochberg false discovery rate p-value correction [39]. Overall, three blood feeding states were analyzed; non-blood-fed (NBF), blood-fed (BF), and IVM-blood-fed (IVM). DE contrasts are denoted in a X:Y format, whereby X is up- or down-regulated compared to Y. For example, a positive fold-change (FC) of 2.5 in the 2IVM:2BF contrast would indicate 2.5 up-regulation in the 2DPE IVM-blood fed group relative to the 2DPE control blood fed group. Likewise, a negative value would indicate down-regulation.

Gene ontology term associations

Gene Ontology (GO) terms for DE genes were retrieved using g:Profiler [40] with the “significant only” option un-checked so as to retrieve all GO terms. While GO terms fall into three different groups (biological process, molecular function, and cellular component), only biological process and molecular function were used for this analysis. Upon classifying GO terms into more general categories, the number of unique transcripts linked to each broad category was plotted for each contrast and direction of differential expression. Because a single gene may return multiple GO terms, transcripts were allowed to fall into multiple categories.

Phylogenetic analysis

To construct a bootstrapped phylogenetic tree, SEQBOOT in PHYLIP 3.69 [41] generated 1000 bootstrap copies. These were then analyzed by PROTDIST in PHYLIP to derive 1000 Dayhoff PAM (Point Accepted Mutation) matrices using the DCMut model [42] which is a version of the original PAM model [43]. The PAM model scales probabilities of change from one particular amino acid to another in terms of % change between two amino acid sequences. NEIGHBOR in PHYLIP then derived 1000 neighbor-joining trees [44] with An. christyi treated as the outgroup. CONSENSE in PHYLIP then derived a consensus tree and reported the number of times that each branch was recovered in the 1000 trees. TreeGraph_2.0.47-206_beta [45] was used to draw the consensus tree.

Real-time quantitative PCR confirmation

As an external validation of RNA-seq results, RT-qPCR primers were developed for the following five genes that showed significant differential expression following an IVM blood-meal: glutathione-S-transferase, delta class, member 3 (GSTd3, AGAP004382) ornithine decarboxylase (ODC, AGAP011806), ABC transporter, subfamily A, member 3 (ABCA3, AGAP006380), thioredoxin peroxidase (Tpx, AGAP011824), and cytochrome P450, family 9, subfamily J, member 5 (CYP9J5, AGAP012296). Validation was done only using 2DPE mosquitoes. Gene expression was normalized to ribosomal protein S7 (RPS7, AGAP010592) as an internal reference gene. Expression of this gene from RNA-seq was consistent (not DE) between IVM and BF groups, justifying its use as a reference. Sequences for all primer pairs can be found in Additional file 2: Table S2. Primers were designed using the web-based PrimerQuest tool [46]. Technical duplicates were performed for each sample-primer pair. Upon TurboDNase treatment, individual samples were pooled and normalized, using 4 ng of input RNA for the SensiFast™ SYBR & Fluorescein One-step kit (BioLine, London, UK) and the suggested 3-step cycling conditions. RNA samples for those labeled “G3” were derived from the same total RNA samples used for RNA-seq library preparation (An. gambiae G3 strain raised at CSU), representing two biological replicates, each consisting of three individual mosquito RNA extracts pooled together. To analyze these transcription patterns in different samples, gene expression was also assessed in mosquitoes originating from West Africa. An. gambiae s.s. (IRSS strain) were recently-colonized at the Institut de Recherche en Sciences de la Santé, in the city of Bobo Dioulasso, Burkina Faso. Additionally, wild Anopheles gambiae s.s. (fourth instar larvae and pupae) were captured from a rain water pool in the Kodeni district of Bobo Dioulasso, and reared into adults. Both the IRSS (colonized) and Kodeni (wild-type) mosquitoes were separated into treatment and control groups at 1DPE, and then at 2DPE the treatment groups were blood fed on the arm of a human who had taken a standard oral dose of IVM (150 μg/kg) 36-h prior, while the control group fed on a different, untreated human. This study was approved by the ethics committee of Colorado State University and written consent was obtained from each participating individual (IRB protocols 09-1148H and 12-3597H). PCR was performed on the Kodeni samples for species identification [47], differentiating between unique SNPs specific to An. gambiae and An. arabiensis. Using this method, only a total of four An. gambiae samples were identified, two from each treatment (IVM and control). Because of this limitation, qPCR results from these field samples represent pooled total RNA from two individuals as opposed to the three used for G3 and IRSS samples. Fold-changes in gene expression were calculated using the ΔΔCt method [48]. A Spearman rho correlation analysis (GraphPad Prism 6) was performed between the RNA-seq and RT-qPCR fold-change results using the same G3 RNA samples.

Results and discussion

Aging effects on IVM susceptibility

Significant mortality (p < 0.0001) was observed in the IVM-treated mosquitoes in each age group (Fig. 1). When comparing mortality in control blood fed mosquitoes between age groups, 6DPE mosquitoes had significantly higher survival than the 2DPE control group, likely reflecting the substantial cost of ingesting the first blood meal in life. They also survived significantly better than the 14DPE group, likely reflecting senescence caused by the accumulated effects of serial blood meals and ovipositions [49]. To most accurately compare the effect of IVM between age groups, the Mantel-Haenszel hazard ratio (HR) was calculated to assess the odds of death in treatment groups relative to their paired controls (Fig. 2). Comparing the 95 % confidence intervals for these HRs, mortality induced by this discriminating IVM concentration was lowest in 2DPE (HR = 5.47; 95 % CI: 4.20 – 7.14, N = 258), highest in 6DPE mosquitoes (HR = 27.89; 95 % CI: 19.84 – 39.02, N = 200), while the mortality was intermediate in 14DPE mosquitoes (HR = 11.49; 95 % CI: 7.54 – 17.52, N = 112).

Fig. 1

Effect of aging on mosquito survival following ingestion of 11.75 ng/mL IVM in calf blood. Mortality was recorded for five days following IVM-containing blood ingestion (and corresponding controls) in the three age groups analyzed. All IVM groups exhibited significant mortality compared with their controls (Log-rank Mantel-Cox test, p < 0.0001). Five biological replicates were performed at each age with the following total sample sizes: 2DPE (n = 258), 6DPE (n = 200), 14DPE (n = 112). Data presented as mean and S.E

Fig. 2

Mantel-Haenzel hazard ratio associated with IVM ingestion in three age groups. Mantel-Haenzel Hazard ratio computed by comparing each IVM treated group to its corresponding control, thus factoring in each age group’s control mortality. Error bars represent the 95 % confidence interval of the hazard ratio. Each age group showed significantly different hazard ratios as there was no crossover between 95 % CIs. 2DPE: HR = 5.47, 95 % CI = 4.20 – 7.14; 6DPE: HR = 27.89, 95 % CI = 19.84 – 39.02; 14DPE: HR = 11.49, 95 % CI = 7.54 – 17.52

Several studies have investigated the effect of mosquito aging on insecticide susceptibility and/or resistance profile maintenance. In the majority of studies, insecticide susceptibility increases with aging [50–53]. In a study with An. gambiae, Rajatileka et al. [50], demonstrated that blood feeding status had little to no impact on DDT, bendiocarb, and deltamethrin susceptibility, but that age was the more dominant factor contributing to insecticide susceptibility. In contrast, two Culex and Anopheles studies revealed increased insecticide tolerance after either single or multiple blood meals, as well as maintenance of resistant profiles [54, 55]. Another study investigating the effect of bloodfeeding on permethtrin tolerance in both a susceptible and resistant strain of An. funestus revealed that a previous blood meal had no effect on tolerance in the susceptible strain but significantly increased tolerance in the resistant strain [56]. In our study, we only tested a strain of An. gambiae that is susceptible to standard insecticides as well as to IVM, as we are not aware of any strains that display IVM resistance. Our data reveal that a previous blood meal and/or aging in An. gambiae significantly decreases tolerance to IVM in a susceptible strain, somewhat contrary to studies using classical insecticides.

Exposure to IVM is intrinsically linked to simultaneous blood meal ingestion and digestion. Our results uniquely show that rather than a steady increase in susceptibility with aging, it is highest in middle-aged mosquitoes that have already processed an earlier blood meal. We previously examined the IVM susceptibility between 2DPE and 8DPE An. gambiae, both ingesting their first blood meal in life. [21]. In contrast to our current data, those data identified no differences in susceptibility to sub-lethal IVM concentrations between age groups, suggesting that ingestion of a prior blood meal is the primary cause of the differential IVM susceptibility between young and old mosquitoes. Ingestion and processing of a blood meal is a taxing physiological commitment as the mosquito ingests up to three times its weight in blood within seconds to minutes [57]. Digesting this meal necessitates rapid diuresis as well as osmotic balancing, and it generates reactive oxygen species and the liberation of cytotoxic heme molecules [58]. The interaction of IVM with these and other blood meal associated molecular and physiological events is undoubtedly complex and not understood. The drug is known to be lipophilic and incorporates into cell membranes [59, 60], but its primary target, the glutamate-gated chloride channel, was shown to be present in An.gambiae nervous tissues of the thorax and head but absent from abdomen [26]. Given this lack of understanding, we analyzed differential gene expression in control and IVM-treated mosquitoes from the 2DPE and 6DPE groups.

RNA-seq – differential expression analysis

RNA-seq was performed on two biological replicates of 10 individuals from both the 2DPE and 6DPE groups, chosen because they exhibited the greatest susceptibility difference to a discriminating IVM concentration. Sequencing, alignment, and annotation identified 13,465 unique transcripts. Filtering of low-expression tags and normalization reduced this number to 2,307 which were then subjected to differential expression (DE) analysis. The complete list of DE genes in each contrast can be found in Additional files 3, 4, and 5.

Figures 3 and 4 summarize the results of DE analyses in the BF:NBF and IVM:BF contrasts. Transcripts were considered significantly differentially expressed if p < 0.05 and fold-change ≥ 2, indicated by red points on Fig. 3. Control blood feeding of both age groups (Figs. 3 and 4; BF:NBF) generally resulted in symmetrical changes, with similar numbers of up- and down-regulated transcripts. Consistent with previous microarray and sequencing studies [61–63], blood feeding alone caused dramatic DE, with blood digestion and vitellogenesis genes up-regulated up to 1000-fold and down-regulated transcripts more modestly changed (Figs. 3 and 4). One limitation of our study is that for technical and cost considerations, we only performed RNA-seq on two biological replicates for each age group, and we only analyzed a single time point post-blood meal digestion (24 h). Having three replicates may have made our analysis more precise, and it is clear that DE changes fluctuate significantly from the time of blood meal ingestion to completion of digestion [63, 64].

Fig. 3

Visualization of differential expression using volcano plots. Differential gene expression depicted as volcano plots for the BF:NBF and IVM:BF contrasts at both age points. Horizontal blue dotted lines indicate significance threshold (p < 0.05) whereas vertical lines indicate fold-change threshold (≥2-fold). Red points represent significantly differentially expressed transcripts by these criteria. Positive x-values represent up-regulation and negative x-values represent down-regulation

Fig. 4

Distribution of fold-changes and total number of DE genes associated with each contrast. Fold-changes were separated into 6 ranges (2–9, 10–19, 20–49, 50–99, 100–999, and >1000). The percentage of transcripts in each range was plotted for each contrast, separated by fold-change direction. The total number of DE transcripts in each contrast and direction is indicated numerically at the top of each bar

Considering age, it is notable that control blood feeding (BF:NBF) by 2DPE mosquitoes resulted in more than double the amount of DE transcripts than 6DPE mosquitoes (253 vs 111 up-regulated and 237 vs 94 down-regulated; Figs. 3 and 4). This is likely due to the fact that the 2DPE mosquitoes were digesting their first blood meal, while the 6DPE mosquitoes were digesting their second blood meal. Indeed, most An. gambiae have a ‘pre-gravid’ stage following emergence [65], whereby they utilize their first blood meal to replenish energy reserves and the second to complete vitellogenesis and mature their first egg batch. It may be that the significantly greater mortality in control mosquitoes that we observed in 2DPE mosquitoes relative to 6DPE mosquitoes (Fig. 1) is associated with this pre-gravid state and/or the higher number of DE genes in this group.

Analysis of DE between the 6DPE and 2DPE age groups within the same treatment categories (6IVM:2IVM, 6BF:2BF, and 6NBF:2NBF) (Additional file 6: Figure S2, Additional file 5) revealed few significantly DE genes, with most appearing in the NBF contrast. However, because these contrasts do not factor in the corresponding control groups, these analyses are less informative in studying the transcriptional response to IVM, and will not be discussed further. An alternative method to analyze the transcriptional response to IVM is to first compare the age-specific BF:NBF contrasts with the IVM:NBF contrasts, subtracting out the blood-feeding effect by focusing only on those DE genes unique to the IVM:NBF contrast. As expected, many DE genes are common between the IVM:BF and subtracted IVM:NBF contrasts within each age group (Additional files 7 and 8: Tables S6.4 and 7.4) (n = 22 in 2DPE, n = 59 in 6DPE) and have very similar fold-change values in both contrasts. However, there remained numerous DE genes unique to the subtracted IVM:NBF contrasts (Additional files 7 and 8: Tables S6.6 and 7.6) that may represent a synergistic interaction between the blood feeding and ivermectin effect. The fold-change values of the genes tended to be modest (−4 < FC < 4) and an examination of their gene descriptions does not identify any particular biological process.

In contrast to control blood feeding, IVM-containing blood meals resulted in asymmetrical transcriptional changes with the majority of significantly DE transcripts showing up-regulation at both ages (Figs. 3 and 4). Additionally, more genes were differentially-regulated in response to IVM in the 6DPE group relative to the 2DPE group, and the fold-changes were overall higher. Figure 5 shows a linear regression analysis comparing those significantly DE transcripts common to the IVM response in both 2-day and 6-day groups, while transcripts unique to one contrast are shown on their respective axes. This again revealed a trend towards higher up-regulation in 6DPE mosquitoes compared with 2DPE mosquitoes (linear regression, inset graph). In addition, four DE transcripts were down-regulated in 6DPE mosquitoes but up-regulated in 2DPE mosquitoes, as seen by their presence below the X-axis. The protein products of these transcripts are ATP-dependent RNA helicase DDX5/DBP2, tyrosine-3-monooxygenase, rap55, and maternal protein exuperantia. It is possible that one or more of these four transcripts are associated with the differential IVM susceptibility of 6DPE vs. 2DPE mosquitoes, but this is not obvious from their descriptions and all have only modest fold-change within each age group.

Fig. 5

Linear regression analysis comparing the gene expression after IVM ingestion between age groups. 62 transcripts were common to both contrasts while an additional 124 transcripts were unique to the 6-day contrast and 39 were unique to the 2-day contrast. Four transcripts showed a different direction of DE between the two contrasts while the remaining 58 showed the same direction of DE. The larger plot shows log₂ transformed fold-change values while the smaller plot shows untransformed data and the accompanying linear regression line and equation

Gene ontology annotations of DE genes in the IVM:BF contrasts in both age groups showed a similar category distribution pattern (Fig. 6), with the majority of transcripts related to: binding, metabolism, chromosome organization & nucleic acid processing, and enzymatic activity. Additionally, a significant number of DE genes were described as ‘Unknown/undefined’, which is attributed to the incomplete annotation of the An. gambiae genome. Lastly, specific to the 6IVM:6BF contrast, several ontology categories had a high-proportion of down regulated genes relative to up-regulated genes (cell cycle, regulation & homeostasis, transport, structural).

Fig. 6

Gene Ontology terms associated with DE genes following IVM ingestion. Gene ontology terms associated with differentially expressed genes 24 h following an IVM blood feed in both 2 and 6-day old mosquitoes. GO terms were gathered with g:Profiler and manually sorted into the 23 functional categories. Only biological process and molecular function terms were used, cellular component terms were removed. Individual transcripts were allowed to fall into multiple categories, based on the GO terms returned, to more accurately reflect the multi-functional aspect of many enzymes

RT-qPCR validation

We validated RNA-seq results of five DE transcripts by RT-qPCR (Fig. 7a). Expression values (log₂ FC) from RNA-seq were compared to those derived from RT-qPCR analysis on the same RNA-samples (G3) and from RNA derived from both colony (IRSS) and wild (Kodeni) mosquitoes naturally-exposed in Burkina Faso. As expected, transcript-specific fold-change in the same G3 strain RNA samples was highly concordant between the RNA-seq and RT-qPCR methods, which was corroborated in the correlation analysis (Fig. 7b). Fold-change values in the naturally-exposed mosquitoes from Burkina Faso (IRSS and Kodeni) appear to generally agree with RNA-seq, however the degree of up-regulation was generally less and lowest in the wild-type (Kodeni) mosquitoes. The discrepancy could partially be explained by the fact that we were unable to determine the concentration of IVM acquired from the blood of the human volunteer (the concentration range is inferred from published pharmacokinetic curves following a standard oral dose), so the concentration was likely slightly different from that fed to the G3 strain via membrane feeders.

Fig. 7

RT-qPCR validation of RNA-seq results. Five genes were selected for DE confirmation in the same RNA samples used for RNA-seq (G3), a laboratory raised colony from west-Africa (IRSS), and wild-caught mosquitoes from the Kodeni region of Burkina Faso (Kodeni). A: Log₂ fold-change comparison with standard deviation (does not apply to RNA-seq results), B: Correlation analysis between RNA-seq and RT-qPCR log₂ (fold-change) results from the same RNA samples (An. gambiae G3 strain). Spearman correlation coefficient is displayed (r). GSTd3: glutathione-S-transferase, delta class, member 3; ODC: ornithine decarboxylase; ABCA3: ABC transporter, subfamily A, member 3; Tpx: thioredoxin peroxidase; CYP9J5: cytochrome P450, family 9, subfamily J, member 5

Genes of interest

The following four groups of DE genes are highlighted based on their putative biological function and degree of differential expression that may relate to the observed age effect on IVM susceptibility, as well as increase our understanding of the broad physiological responses to IVM ingestion.

Detoxification genes

Several DE transcripts (Table 1) were identified for their involvement with either detoxification of insecticides or oxidative-stress response. These include enzymes classically related to insecticide resistance such as glutathione-S-transferase (GST), cytochrome P450s, and ABC transporters. Half of the transcripts were differentially expressed at similar fold-changes in both age contrasts. Two oxidative-stress response genes were significantly DE only in the 6DPE contrast, but thioredoxin peroxidase was highly up-regulated in both age contrasts. While both age contrasts showed up-regulation of a cytochrome P450 upon IVM ingestion, the exact CYP gene differed (CYP9J5 in 2DPE, CYP6Z1 in 6DPE).

Table 1 Detoxification-related transcripts differentially expressed after IVM ingestion

GSTs have been linked with resistance to all major insecticide classes, functioning to conjugate reduced glutathione to a target molecule to both neutralize electrophilic sites and render it more water soluble and easily excretable [66–69]. Cytochrome P450-mediated oxidation of insecticides has been widely studied as it relates to resistance [70–73], and increased ABC-transporter expression and activity in insects has been linked to detoxification of carbamates [74], organochlorines [75], organophosphates [76], and pyrethroids [77]. Similarly, ABC transporters have been studied for their potential involvement in the detoxification of IVM in Pediculus humanus [29], Culex pipiens [78], Rhipicephalus (Boophilus) microplus [30, 76], Caenorhabditis elegans [79], and Chironomus riparius [80]. The degree of up-regulation of the ABCA3 (AGAP006380) was largely the same between age groups and therefore would not explain the susceptibility differences observed. ABCG1, on the other hand, was down-regulated in the 6DPE group but not DE in the 2DPE comparison. Additionally, GSTd3 was up-regulated 3.6-fold in 2DPE but not DE in 6DPE. Taken together, there may be involvement of GSTs and ABC transporters in the response to IVM, however because fold-changes were generally mild, other factors are likely involved in explaining the differential age susceptibility.

Niemann-pick type C genes

The Niemann-Pick type C-2 (NPC2) family of genes are of interest because members were among the most highly up-regulated transcripts in our experiments (Table 2). One member (AGAP002851) is among the most highly DE genes in response to blood feeding within 2DPE mosquitoes (148 fold in our data set), it remains highly up-regulated in 6DPE mosquitoes in response to regular blood feeding (10 fold), and it is the highest up-regulated transcript related to aging and termination of the pre-gravid state (55 fold, 6NBF:2NBF, Table 1). Other members of this gene family responded more modestly to blood feeding in both age groups. Two NPC2 genes (AGAP002847, AGAP002848) were highly up-regulated in response to IVM ingestion at both ages, but the transcription of AGAP002847 in particular almost doubled in 6DPE vs. 2DPE mosquitoes.

Table 2 Niemann-Pick family members differentially expressed in multiple contrasts

In humans, these proteins are involved in cholesterol trafficking and homeostasis in the endosomal/lysosomal system [81, 82], and defects in either of the type C proteins, NPC1 or NPC2, cause two hereditary childhood neurodegenerative diseases (Niemann-Pick C diseases type 1 and type 2, respectively) [83]. While NPC1 is a large transmembrane-associated protein found in late endosomes [82], NPC2 is a small secretion signal-containing protein that is both secreted and found in lysomsomes [82, 84]. NPC2 contains a conserved domain called ML (MD-2 [myeloid differentiation factor-2]-related Lipid-recognition) that binds cholesterol and other lipids, such as fatty acids [85]. This ML domain is also a key domain in the human proteins MD-1 and MD-2 (myeloid differentiation factors), where it binds bacterial lipopolysaccharide (LPS) and mediates cellular innate immune responses to LPS through co-interactions with toll-like receptor 4 [86–88].

In insects, NPC2-related proteins appear to carry out a variety of functions related to their lipid-binding domains, including lipid adsorption, lipid trafficking and metabolism, pathogen recognition and immune-regulation. In insects they have been named either MD2-like, NPC2-like, or simply ML, but to remain consistent with their annotations in VectorBase [36], we refer to them here as NPC2 proteins and also reference some of their previous ML designations. These proteins have been described in the dipteran species Drosophila, Aedes, Culex, and Anopheles. In Drosophila, the eight NPC2 genes are transcribed variably in time and space within embryos and at least two of the genes (NPC2a and NPC2b) are necessary for sterol homeostasis and ecdysteroid biosynthesis [89]. Furthermore, three members of the Drosophila NPC2 protein family (NPC2a, NPC2e, and NPC2h) have been shown to bind LPS, lipid A, peptidoglycan and lipoteichoic acid, and play a role in immune signaling pathways involved in antimicrobial peptide production [90]. In Ae. aegypti, NPC1 may facilitate dengue virus entry while NPC2 proteins may play a role in modulating mosquito immunity to dengue virus [91]. The Anopheles NPC2 proteins (ML proteins) were extensively highlighted in Dong et al. [92] where the thirteen members were designated AgMDL1-13. They found that AGAP012352 (AgMDL1) and AGAP002857 (AgMDL2) were significantly transcribed in An. gambiae midguts upon P. falciparum ookinete invasion and that RNAi silencing of the former agonized oocyst development and growth of both gram negative and gram positive bacteria upon challenge.

We performed a phylogenetic analysis of NPC1 (Fig. 8) and NPC2 (Fig. 9) using current An. gambiae s. l. complex genome data from VectorBase [36] to try to better understand this diverse gene family and its myriad of functions. Anopheles christyi and An. epiroticus were included as outgroups in these analyses. As seen previously in Jupatanakul et al. [91], there are two distinct NPC1 proteins (NPC1a and NPC1b). Both NPC1a and NPC1b are monophyletic and have evolved during speciation within An. gambiae s.l., suggesting conserved functions, most likely related to basic sterol homeostasis and 20-hydroecdysone biosynthesis functions that have been ascribed to these proteins in D. melanogaster [89, 93, 94]. In our dataset, only An. gambiae NPC1a (AGAP008137) is modestly down-regulated and only in the 6IVM:6BF contrast (Table 2).

Fig. 8

Phylogenetic tree of the NPC1 gene family in Anopheles species. The Aedes aegypti Niemann-Pick C1a protein (AAEL003325) was used to identify the orthologues in Anopheles gambiae s.l., An. gambiae s.s., An. arabiensis, An. melas (note the duplication), An. merus, An. coluzzii, and An. quadrimaculatus. The Ae. aegypti Niemann-Pick C1b protein (AAEL009531) was used to identify the orthologues in Anopheles gambiae s.l. For both NPC1a and NPC1b, Anopheles christyi and An. epiroticus were included as outgroups. The amino acid alignment contained 1818 characters. Branches with >50 % boostrap scores are labelled

Fig. 9

Phylogenetic tree of the NPC2 gene family in Anopheles species. Aedes aegypti Niemann-Pick C2 proteins from four phylogentic groups (AAEL009557, AAEL001650, AAEL004120, AAEL006854) were used in BLASTp to identify the orthologues in Anopheles gambiae s.l. Anopheles christyi and An. epiroticus were included as outgroups. BLASTp returned a total of 108 orthologues in Anopheles species, separated into four groups based on the Ae. Aegypti orthologue. BLASTp and phylogenetic analysis of AAEL009557 returned forty proteins in six subgroups (group 1, blue font), AAEL001650 returned thirty proteins in five subgroups (group 2, gold font), AAEL004120 returned nine proteins in one subgroup (group 3, black font), and AAEL006854 returned twenty-nine proteins in four subgroups (group 4, pink font) The genome location of An. gambiae s.s orthologues are listed after each gene. Groups with >50 % bootstrap support are red. An. gambiae orthologues appear in a larger font with the IVM and blood meal-responsive transcripts in the largest font. Note the duplication appearing in the same physical location in group 4a

The NPC2 proteins from An. gambiae s.s [92], D. melanogaster [90] and Ae. aegypti [91] have been previously subjected to phylogenetic analysis. These latter two analyses placed these proteins in 4 broad groups. From this basis, BLASTp against Ae. aegypti representatives of the 4 previously-described groups (AAEL009557: group 1, AAEL001650: group 2, AAEL004120: group 3, AAEL006854: group 4) identified 108 NPC2-like proteins in An. gambiae s. l. and in the outgroup species An. christyi and An. epiroticus. Forty of these were similar to AAEL009557 (AaegML17, group 1) and phylogenetic analysis placed these into six subgroups (blue font in Fig. 9). Thirty proteins were similar to AAEL001650 (AaegML33, group 2) and phylogenetic analysis placed these into five subgroups (gold font in Fig. 9). Nine proteins were similar to AAEL004120 (AaegML1, group 3) and phylogenetic analysis placed these into one group (black font in Fig. 9). Lastly, 29 proteins were similar to AAEL006854 (AaegML13, group 4) and phylogenetic analysis placed these into four subgroups (pink font in Fig. 9). Some groups did not contain a representative from all 8 Anopheles species (e.g., subgroup 1a, band e; and subgroups 2a and b), but in general, each of the eight Anopheles species had one orthologous copy of each gene in each subgroup. The exceptions to this are in subgroup 1f, group 3 and subgroups 4a and 4b, where gene duplication and further evolution of the gene subgroup is apparent in some species. Indeed, eleven separate An. gambiae NPC2 genes are localized in chromosome 2R, nine of which are within region 13E. The phylogeny and gene localizations indicate ancient duplication event(s) within Anopheles followed by rapid diversification, and support data suggesting that NPC2 proteins possess a wide-range of lipid and/or sterol binding specificities and diverse functional characteristics, similar to the conclusions suggested by Huang et al. [89] Shi et al. [90], and Jupatanakul et al. [91] given the redundancy of NPC2 genetic structure in Drosophila and Aedes.

AGAP002851 is the main blood meal-responsive gene (Fig. 9, subgroup 4c, largest font), increasing especially in pre-gravid mosquitoes. Its Drosophila orthologue, NPC2e, has so far only been shown to be involved in innate immune signaling [90]. Consequently, the Anopheles protein may primarily be responding to the substantial bacterial bloom that occurs in the midgut lumen after blood meal ingestion [95], but it may also have a dual role in sterol homeostasis as has been shown for the Drosophila NPC2a protein. It is also important to note that insects cannot synthesize sterol from precursors, but must obtain it from food [96]. The IVM-responsive NPC2 genes (AGAP002848 and AGAP002847) appear to be evolutionarily distinct from the blood meal-responsive NPC2 (group 2 versus group 4, Fig. 9) and were placed in two different subgroups (2a and 2b), suggesting they carry out different functions. The high up-regulation of separate NPC proteins following either a normal blood meal or IVM blood meal paired with the evidence for very diverse binding and function of these proteins in insects suggests direct involvement with the An. gambiae response to IVM. Furthermore, the significantly higher upregulation of AGAP002847 in the 6DPE contrast compared to the 2DPE contrast may influence the differential age susceptibility we observed. These data highlight the need for future studies to investigate this potential connection and further characterize these proteins.

Peritrophic matrix associated genes

The mosquito peritrophic matrix (PM) is a complex structure of chitin-binding and chitin-associated proteins [97, 98]. It lines the luminal side of the midgut and functions to compartmentalize the digestive process in the lumen and protect the midgut epithelium from abrasion during digestion. Our experiment only allows us to follow transcript changes that affect the thick, Type 1 PM in An. gambiae females that is secreted upon a blood meal rather than the Type 2 PM constitutively produced in larval guts. The primary Type 1 PM proteins of An. gambiae are pre-translated in midgut cell secretory vesicles and released into the lumen upon midgut distension following a blood meal [99, 100], and their transcriptional profiles follow a pattern of high steady-state transcription in NBF and post BF mosquitoes, with only a significant drop at 3 h post BF [64]. Comparison of our dataset with work done by Dinglasan et al. [97] revealed DE of a high number of putative peritrophic matrix-associated genes upon IVM ingestion in both age contrasts (Table 3), suggesting that PM gene transcription is responding directly to IVM treatment. In general the 2DPE contrast showed both an increased number of significantly DE transcripts as well as higher fold-changes than the 6DPE contrast, but interestingly, both age contrasts had DE genes unique to their own age group. Additionally, gene ontology analysis revealed down-regulation of chitin processing transcripts upon blood-feeding (Additional file 9: Figure S3) yet up-regulation of this category after an IVM blood meal (Fig. 6).

Table 3 Peritrophic matrix-associated transcripts differentially expressed after IVM ingestion

IVM is known to significantly decrease the rate of blood meal digestion [21, 101], even though the glutamate-gated chloride channel target is not expressed in the An. gambiae midgut [26]. The significant DE of PM-related transcripts in both age contrasts may be associated with this reduced blood meal processing. There is also good evidence that a primary function of the PM is to bind the abundant heme released from red blood cell digestion, which is toxic and can disrupt cell membranes [102]. The observed increased PM gene transcription may be in response to the prolonged presence of heme in IVM-treated females. Following this hypothesis, it is notable that the less IVM-susceptible age group (2DPE) shows overall higher DE transcription among PM-associated genes in response to IVM than the more IVM-susceptible age group (6DPE), perhaps contributing to their increased tolerance to IVM toxicity.

Immune response genes

Table 4 lists significantly DE genes previously-linked with mosquito innate immune responses against microorganisms. The older age contrast showed significant DE among many immune-related transcripts. Most were modestly up-regulated, but the expression of two proteases was changed by more than 15-fold (AGAP010968 & AGAP007142). On the other hand, only four immune function genes were significantly DE in the 2DPE contrast.

Table 4 Immunity-related transcripts differentially expressed after IVM ingestion

The up-regulation of so many immune-responsive elements specific to the 6DPE contrast, which have taken a previous blood meal, may suggest that IVM is further influencing the mosquito-microbial interactions and/or innate immune memory from the previous blood meal [103, 104]. It will be important in future studies to examine the impact that IVM may have on the microflora of surviving mosquitoes, as well as how this interaction might influence the previously observed anti-sporogonic effects of the drug [105].

More complex interactions between IVM and the significant DE genes mentioned in all tables above are likely, and some genes have clear linkages across the groups we have highlighted. For example, in the case of immune response genes, it is clear that NPC2 proteins serve as pattern recognition receptors for bacteria lipid moieties. The dramatic up-regulation of some NPC2 transcripts in response to IVM treatment may connect to the up-regulation of downstream immune response pathways (Table 4). Similarly, the PM directly interacts with the microflora bloom that occurs in the gut of blood fed mosquitoes, and Dinglasan et al. [97] previously isolated two NPC2 proteins (AGAP002851, AGAP002848) from dissected An. gambiae peritrophic matrices. Also, the PM protein ICHIT (AGAP006432) is immune-responsive following both bacterial and Plasmodium invasion [106], but it also has protein-protein aggregation properties, potentially acting as “glue” between chitin-binding and non-chitin binding proteins in the PM [97]. Lastly, detoxification responses connect to cholesterol processing by NPC proteins and to blood digestion and the immune response. It is known that human ABC transporters are essential for sterol homeostasis by efflux of cholesterol and other lipids from cells [107], and anti-oxidant enzymes like the highly up-regulated thioredoxin peroxidase (AGAP011824) are important for redox homeostasis in the gut that may be disrupted in IVM-treated mosquitoes through either prolonged blood meal digestion and/or microflora dysregulation.

Conclusions

We have shown here that older An. gambiae mosquitoes that have ingested previous blood meals are more susceptible to IVM compared to young mosquitoes which have not ingested a prior blood meal. Following from this difference, DE analysis showed that IVM in a blood meal induces significant transcription of more than 100 genes, most of which are up-regulated. Gene expression analyses revealed that IVM in a blood meal mostly induced transcription of non-canonical genes, including PM-associated genes, immune response genes and NPC genes. In regards to the differential susceptibility observed between 2- and 6-day old mosquitoes, several gene expression pattern differences may play a role in producing this phenotype including: a) the NPC gene AGAP002847, which changed from 80.6 upregulated in the 2DPE contrast to 159.5 fold upregulated in the 6DPE contrast, b) the PM-associated genes in which expression was generally higher in the younger, less susceptible age group, suggesting a difference in midgut structure and/or makeup, and c) the immune-related genes, many of which were significantly upregulated in 6DPE mosquitoes only, which suggests an interaction between IVM and immune responses and may be contributing to their increased mortality after IVM ingestion. In contrast, transcriptional evidence for involvement of traditional xenobiotic detoxification mechanisms was limited to the down-regulation of an ABC transporter in 6-day mosquitoes and up-regulation of a GST in 2-day mosquitoes. In addition to new insights into the physiological and molecular effects of IVM on An. gambiae, some unexpected findings have added to the knowledge of the under-studied insect NPC proteins and their interesting evolution and diverse functions. While traditional detoxification mechanisms may be party involved, overall the data suggest complex midgut interactions resulting from IVM ingestion that likely involves blood meal digestion physiological responses, midgut microflora, and innate immune responses. These findings have increased our understanding of IVM’s effects on malaria vectors and provide starting points for future studies looking more closely at mosquito-IVM interactions.

Supporting information

The data set supporting the results of this article is available in the NCBI Sequence Read Archive (SRA) repository. Accession #: SRP060547.

Additionally, expression data from this study is available on the VectorBase Expression Browser under the title “Ivermectin-containing blood meals at different ages” (Seaman et al., 2015) (https://www.vectorbase.org/expression-browser/).