Skip to main content

Transcriptional cellular responses in midgut tissue of Aedes aegypti larvae following intoxication with Cry11Aa toxin from Bacillus thuringiensis

Abstract

Background

Although much is known about the mechanism of action of Bacillus thuringiensis Cry toxins, the target tissue cellular responses to toxin activity is less understood. Previous transcriptomic studies indicated that significant changes in gene expression occurred during intoxication. However, most of these studies were done in organisms without a sequenced and annotated reference genome. A reference genome and transcriptome is available for the mosquito Aedes aegypti, and its importance as a disease vector has positioned its biological control as a primary health concern. Through RNA sequencing we sought to determine the transcriptional changes observed during intoxication by Cry11Aa in A. aegypti and to analyze possible defense and recovery mechanisms engaged after toxin ingestion.

Results

In this work the changes in the transcriptome of 4th instar A. aegypti larvae exposed to Cry11Aa toxin for 0, 3, 6, 9, and 12 h were analyzed. A total of 1060 differentially expressed genes after toxin ingestion were identified with two bioconductoR packages: DESeq2 and EdgeR. The most important transcriptional changes were observed after 9 or 12 h of toxin exposure. GO enrichment analysis of molecular function and biological process were performed as well as Interpro protein functional domains and pBLAST analyses. Up regulated processes include vesicular trafficking, small GTPase signaling, MAPK pathways, and lipid metabolism. In contrast, down regulated functions are related to transmembrane transport, detoxification mechanisms, cell proliferation and metabolism enzymes. Validation with RT-qPCR showed large agreement with Cry11Aa intoxication since these changes were not observed with untreated larvae or larvae treated with non-toxic Cry11Aa mutants, indicating that a fully functional pore forming Cry toxin is required for the observed transcriptional responses.

Conclusions

This study presents the first transcriptome of Cry intoxication response in a fully sequenced insect, and reveals possible conserved cellular processes that enable larvae to contend with Cry intoxication in the disease vector A. aegypti. We found some similarities of the mosquito responses to Cry11Aa toxin with previously observed responses to other Cry toxins in different insect orders and in nematodes suggesting a conserved response to pore forming toxins. Surprisingly some of these responses also correlate with transcriptional changes observed in Bti-resistant and Cry11Aa-resistant mosquito larvae.

Background

Vector-borne diseases remain some of human kind’s biggest health challenges that impose a significant toll on afflicted populations. The World Health Organization (WHO) estimates more than 500 million cases of malaria worldwide, transmitted by Anopheles mosquitoes, while Aedes mosquitoes transmit different viruses that cause hundreds of millions of annual infections of yellow fever, dengue, and chikungunya fever [1]. About 40 % of the human population is at risk since they live in Aedes inhabited locations. Importantly effective vaccines or antiviral drugs for most of these diseases do not exist, except for yellow fever. Thus efforts to control mosquito populations remain an important strategy to reduce infection rates. Biological control through the use of insecticidal proteins produced by the bacterium Bacillus thuringiensis (Bt) represents a promising alternative to control mosquitoes due to the emergence of mosquito populations that are resistant to chemical insecticides [1]. In addition, the recalcitrant nature of chemical insecticides and their non-target effects pose environmental and health disadvantages [2, 3].

Bt bacteria produce a wide variety of toxins, of which the pore-forming 3-domain Cry and the Cyt toxin families have been widely studied, and used in mosquito control due to their highly specific insecticidal activities [4]. In particular, the Bt subsp. israelensis (Bti) strain has been used in the control of A. aegypti. Bti mainly produces a combination of Cry11Aa, Cry4Ba, Cry4Aa and Cyt1Aa toxins as crystalline inclusions during the sporulation phase. Although Bti has been used for more than four decades in vector control programs worldwide, resistant populations have not been observed in the field. This lack of resistance is usually explained by the combined action of Bti Cry and Cyt toxins which are synergistic towards mosquito larvae [5, 6]. Also, it was shown that Cyt1Aa overcomes lab-generated resistance to individual Cry toxins [7].

The mode of action of Cry toxins has been characterized mainly in lepidopteran insect species. However, conservation of the three-dimensional structure among Cry toxins specific to different insect orders as well as conservation of cell receptors and other biochemical evidence suggest a common mechanism of action of Cry toxins in different insect orders [4]. Briefly, susceptible larvae ingest the toxin crystals, from which protoxins are solubilized and processed by midgut proteases. Activated toxin monomers bind to GPI-membrane anchored proteins, such as aminopeptidases (APN) or alkaline phosphatases (ALP). Subsequently, an additional binding event occurs with a cadherin protein, which favors further processing of monomers and induces conformational changes that generate an oligomeric toxic structure. This oligomer gains higher affinity to GPI-anchored receptors and subsequently inserts into the membrane, opening a non-selective pore [8]. Osmotic shock mediated by this pore is thought to be the main cause of cell death [8]. For the mosquitocidal Cry toxins produced by Bti, similar Cry-receptor molecules have been identified in different mosquito species like cadherins, APN and ALP indicating a conserved mode of action in mosquitoes (reviewed in [9]). In addition an α-amylase was shown to bind Cry4Ba and Cry11Aa toxins and proposed as toxin receptor in Anopheles albimanus [10].

Although much is known about the molecular interactions of Cry toxins with its membrane receptors, much less is known about the cellular responses after Cry toxin action. Bacterial pore-forming toxins have been studied in a number of systems including mammalian-cells, where toxin-receptor binding or damage to the membrane elicits a diverse array of responses depending on the toxin dose, cell type, and duration of exposure. Mechanisms such as apoptosis and pyroptosis have been observed, as well as lipid synthesis, membrane repair, endo/exocytosis and activation of autophagy [1115]. Shedding of membrane proteins [16] and membrane vesicle blebbing are other observed effects [13]. It has been reported that phosphorylation of MAP kinases is frequently observed during the responses elicited in different cells to several pore-forming toxins. Specifically JNK kinase was activated by streptolysin-O [17] and Cry5B [18]; ERK 1/2 by anthralysin-O [16]; and p38 MAPkinase by streptolysin-O [17], anthralysin-O [16], α-toxin [14, 19], pneumolysin [20, 21], β-toxin [22], aerolysin [23], Cry5B [18], Cry3Aa [24], Cry1Ab and Cry11Aa toxins [25].

The targets of these MAP kinases include transcriptional factors implying that changes in gene expression can be expected during cellular responses to damage induced by a pore-forming toxin. Previous studies have indicated that this is indeed the case for Cry toxins in insects and nematodes. For example, transcript changes in response to Cry toxins have been reported in lepidopteran larvae, such as Spodoptera frugiperda [26], Choristoneura fumiferana [27], and Lymatria dispar [28], as well as in the coleopteran Tenebrio molitor [24]. However, limitations on the availability of a reference genome restricted the amount of information obtained, and phylogenetic distances between these insects may influence the response mechanisms to Cry toxins. Therefore, any new information about transcriptome changes after Cry toxin exposure in larvae brings us closer to unraveling a shared response in insects against these important toxins. In this paper we present a transcriptomic study of the response of the mosquito A. aegypti in the midgut tissue after intoxication with Cry11Aa.

Results

High-throughput Illumina sequencing

Bti produces several mosquitocidal toxins. For this study, we chose Cry11Aa since it shows high toxicity to A. aegypti larvae with an LC50 value of 450.9 ng/ml (confidence limits, 350.5–557.1 ng/ml). Also, because it was previously shown that Cry11Aa affects MAPK p38 phosphorylation during intoxication of A. aegypti larvae [25]. In this work we analyzed the transcriptome of 4th instar A. aegypti larvae exposed to Cry11Aa toxin (500 ng/ml) for different times. Control larvae at the shorter and longer incubation times in the absence of toxin were also included in our analyses. The Illumina reads obtained were checked for overrepresented sequences or artifacts, but none were found. All samples produced around 20–60 million reads with good quality. Reads were aligned to the A. aegypti transcripts as indicated in methods. Additional file 1: Table S1 shows the alignment rates for all samples. Above 65 % of all reads were aligned to A. aegypti transcripts uniquely. In addition 5 % to 12 % of total reads aligned to A. aegypti genome but not to the transcriptome, which could represent new genes or new splicing alternatives that have not been characterized before (data not shown).

Differential expression analysis

We used two different algorithms, DESeq2 and EdgeR [29, 30], to identify genes with significant differential expression. Table 1 shows that DESeq2 reported more genes as differentially expressed genes (DEG) than EdgeR for every exposure time except at 3 h. Also, most of the genes detected by EdgeR were also detected by DESeq2, with the subset of genes defined as DEG only by EdgeR being small overall. We observed that as the exposure time to Cry11Aa increases, the number of DEGs also increases, where the number of positively regulated genes is larger than the number of negatively regulated ones.

Table 1 Number of genes reported in differential expression analysis by DESeq2 (P-adj < 0.05) and EdgeR (FDR < 0.05)

Across all timepoints, a total of unique 1060 DEG common to both packages were used for further analyses (Additional file 2: Table S2). As control we analyzed the DEG of non-toxin exposed larvae after incubation for 0 or 12 h by RNAseq in triplicate. In contrast to the 1038 common DEG seen at 12 h of Cry11Aa exposure, only 13 DEG were found in the control non-toxin exposed larvae after 12 h, with nine genes down regulated and four up regulated genes (Additional file 3: Table S3). Of these genes found in control larvae, only 6 genes are also found in the 1060 DEGs found for all toxin exposure conditions. Two of them, AAEL003841 and AAEL003816, correspond to a defensin antimicrobial peptide and a close homolog, suggesting that this process is not related to the response to Cry11Aa toxicity.

Biological Process and Molecular Function GO enrichment analysis

We performed a GO enrichment analysis of biological process annotations on the subset of DEG that were common in the DESeq2 and EdgeR results. Figure 1 represents these overrepresented terms for each toxin exposure time. Some related GO terms are represented as a single term for simplicity. It is not surprising to find almost no enriched terms at 3 h of Cry11Aa exposure since few genes were detected as differentially expressed under this condition. However, there is a clear increase in this number as exposure time increases, as is expected from the number of DEG for each condition. As exposure time to Cry11Aa increases there is a progressive enrichment of terms related to cellular remodeling and component shuttling in positively regulated genes: cytoskeleton, endocytosis, lipid metabolism, cell-wall metabolism and vesicular transport.

Fig. 1
figure1

Enrichment analysis for biological process GO terms. For each experimental condition, terms reported as enriched by TopGO on differentially expressed genes (as described in Methods) are colored green if present in down regulated genes, red if found in up regulated genes, and purple if found in both sets of genes. Some related GO terms are shown as one for simplicity

In addition, our data suggest that small GTPase and PIP3 signaling contribute to a general response to the membrane damage induced by Cry11Aa. It was previously reported that signaling through Rho-GTPase participates in the regulation of vesicular traffic [31], indicating that this response could have a role in the removal and repair of damaged cell membranes after Cry toxin intoxication. Protein integrity is also conserved through the up regulation of chaperones Hsp70 and Hsp20, and this is reflected in the enrichment of the terms for protein folding, complex assembly and post-translational modification.

Among the terms that were down regulated we found some transmembrane transport functions such as those involved in sulfate and carbohydrate transport. Other terms that were highly repressed include those related to catabolic pathways, cell replication, and PKC signaling. The GO term for immune response was initially negatively regulated, but as intoxication proceeds, immune response showed a positive regulation. Importantly our studies were done with purified Cry11Aa-crystals by sucrose gradient centrifugation, limiting the number of bacteria or bacterial wall components present in the crystal samples used for intoxication. These data may represent a survival response of larvae that attempt to defend against possible future bacterial intrusions.

To compare all conditions at the same time instead of through individual enrichment, we performed a clustering analysis on the log2 fold values of the 1060 common DEG as described in Methods, and determined those genes with an overall similar expression profile throughout the time of toxin exposure. Four different clusters could be observed (Fig. 2). Clusters A and B show groups of genes with down regulation patterns. Genes with strong down regulation (when difference between maximal and minimal expression values are over 2 fold) are grouped in cluster B while moderately down regulated genes (when the difference between maximal and minimal expression values are equal or lower than 2 fold) are in cluster A. Up regulated genes were grouped in clusters C and D. Cluster C showed those with moderate up regulation and cluster D those with strong up regulation. Additional file 2: Table S2 lists all genes that are present in each cluster. The only exception was AAEL008767 gene, corresponding to a serine protease, which showed a unique expression profile with a peak induction at 3 h. With these separate clusters we performed again a GO enrichment analysis of biological process terms as described above and results are presented in Table 2. In general the GO enrichment results of these clusters resemble the GO terms shown in Fig. 1. We observed moderate down regulation of terms involved in oxidation processes, TCA enzymes, transmembrane transport, aldehyde, monocarboxylic acid and glutathione metabolisms, catabolic process of nitrogen compounds, heterocycle and nucleosides, and PKC signaling. The strongest down regulation at 9 h were found in bacterial defense terms. The cluster with moderate up regulation is enriched with terms of vesicular transport, lipid biosynthetic process, protein complex formation and protein folding, cellular morphogenesis, cytoskeleton binding, initiation of DNA replication, integrin and cell-cell signaling, as well as MAPK and small GTPase signaling. The strong up regulation cluster was enriched for actin nucleation, lipid catabolism, calcium ion transport, and Rho GTPase regulation (Table 2).

Fig. 2
figure2

Clustering of differentially expressed genes by expression profile. Log2 fold expression profiles from 0 to 12 h of Cry11Aa exposure for all genes determined as differentially expressed in at least one experimental condition and clustered using DEseq2 and Cluster 3.0 as described in Methods. We considered strong response when difference between maximal and minimal expression values are over 2 fold and moderate response when difference between maximal and minimal expression values are equal or lower than 2 fold. Cluster a, represents genes with moderate down regulation; cluster b, strong down regulation; cluster c, moderate up regulation and cluster d, strong up regulation

Table 2 TopGO enrichment of biological process terms in differentially expressed gene clusters

Of the 17,478 predicted genes for A. aegypti, only about a third are annotated with a biological process GO term. Therefore, there are many DEG for which a role in biological process could not be assigned in our GO analysis shown in Fig. 1. We performed a similar enrichment analysis but using a different GO ontology, that was focused to the molecular function instead of the biological process. Figure 3 shows this GO molecular function enrichment analysis. At 3 h the only enriched term for both up and down regulated genes is for serine-type endopeptidase activity. This could be an artifact due to the small number of DEG detected for this condition. An alternate explanation is that the serine protease response may be down-regulated when the insects stop feeding as a response to intoxication, and other serine proteases may be up-regulated due to a starvation response. From 6 h onwards, up regulated genes are enriched in terms related to phospholipid interactions, such as lipid binding or phospholipase A2 activity. We also observed enrichment of carboxylyase activity in both up and down regulated genes at 6 h, but this activity is only further enriched at 9 and 12 h of toxin exposure. Here again, MAPK activity and terms related to both ARF and Rho GTPase activity are enriched in up regulated genes from 6 h and forward, supporting a possible role these two signaling pathways in the response to Cry intoxication. Probably related to the activity of these proteins, nucleotide binding and phosphotransferase terms are also enriched at 12 h. The scavenger receptor term in these exposure times may point to the recycling of membrane proteins. In the down regulated genes at 9 and 12 h we observed enrichment of terms related to ion binding and transport as well as terms related to oxidation-reduction processes, although other different terms involved in oxidation-reduction processes were observed in the group of up regulated genes. Interestingly, the ALP and APN activity terms were down regulated at 6 and 12 h, respectively. Proteins in both these families function as receptors for Cry toxins [4]. However, these genes that are down regulated do not correspond to the ALP or APN proteins that have been shown to bind Cry11Aa toxin in mosquito epithelial cells [32].

Fig. 3
figure3

Molecular functions of genes differentially expressed. For each experimental condition, genes were categorized into different categories as reported by TopGO and enrichment on differentially expressed genes (as described in Methods) are colored green indicating down regulated genes, red indicating up regulated genes, and purple if found in both sets of genes

Identification of differential expressed genes that were not annotated

We found 676 DEG that were not annotated with GO biological process terms. To extend our analysis, these 676 genes were analyzed to obtain their Interpro protein functional domains [33] as described in Methods. Table 3 shows the representative domains that are distributed in the different sets of genes, and Additional file 4: Table S4 shows a list with all domains found with frequency ≥ 90th percentile (0.9). For the up regulated genes at 6 h of toxin exposure we detected protein domains related to carboxylases, and these are also up regulated at 9 and 12 h of Cry11Aa intoxication. At these times we also found further evidence of genes that support the results of up regulated genes identified by GO enrichment of biological functions that were presented in Fig. 1 such as chaperones involved in protein folding and complex assembly, annexin involved in cytoskeleton and EF-Hand domains involved in Ca2+ binding. Regarding DEG that were down regulated we found domains of zinc fingers, type-C lectins, and ATPse domains of replication proteins. ABC transporter domains were also found to be down regulated.

Table 3 Representative Interpro domains identified on differential expressed genes without biological process GO annotation

After performing the Interpro analysis we realized that 154 DEG still did not have any identifiable protein domain. These 154 genes were further analyzed by performing pBLAST analysis on the NR database [34], though few proteins were found to have significant homology to other gene targets. Table 4 shows the identified genes. We found homologs of phospholipase A2 family, which supports the importance of activating this process in response to Cry11Aa intoxication. Again, a component of Ras GTPase signaling such as the homolog of Ras guanine exchange factor was identified in genes induced after Cry11Aa exposure. Among the down regulated genes we found some antimicrobial peptides such as CEC-A and gambicin, and also an ADAM metalloprotease.

Table 4 Homologs to differentially expressed genes without biological process GO annotation found by pBLAST analysis

Validation through RT-qPCR

To corroborate the RNAseq data, we selected 14 genes from the common DEG to validate their expression by RT-qPCR at the same time-points analyzed by RNAseq. Among these genes, two genes related to MAPK signaling: JNK (AAEL008622) and Ets domain containing protein (AAEL006533) were analyzed. The closest homolog for AAEL006533 is Elk-1 transcription factor, a known downstream target of MAPK-p38. Therefore its positive regulation may confirm the activation of this pathway, as the transcript for p38 itself does not show differential expression, which is consistent with previous observations [25]. According to Vectorbase annotation A. aegypti has to two isoforms of JNK, where AAEL008622 shows higher homology to human JNK2. We also selected one gene for phospholipase A2 activity (AAEL001523), two genes for vesicular traffic such as lipoma preferred partner and sorting nexin (AAEL007704, AAEL005655), two metalloproteases (AAEL002661, AAEL005992), three genes for transporter proteins related to transport such as monocarboxylate, sucrose and sugar transporter proteins (AAEL008347, AAEL003633, AAEL012903), a CoA-Ligase (AAEL014664) and an ALP (AAEL003905). We also analyzed the gene for the antimicrobial peptide defensin (AAEL003841) and a hypothetical protein with high homology to defensin (AAEL003816) that were also observed in control larvae.

Transcript abundance ratios were determined at all times after toxin exposure. As controls we analyzed the transcript abundance ratios of control larvae without toxin exposure for different times or after exposure to the non-lethal Cry11Aa mutants E97A and V142D. The mutant E97A is not toxic since it is unable to form toxin oligomers required for pore formation [35], while the mutant V142D is able to oligomerize but it does not insert into the lipid bilayer of the cell membrane [36]. Additional file 5: Figure S1 shows the relative expression of all 14 genes at different times of the four conditions analyzed (after Cry11Aa toxin ingestion, control without-toxin, or with non-toxic mutants Cry11Aa-E97A and Cry11AaV142D). In order to compare data of Additional file 5: Figure S1 with our results of transcription expression by RNAseq a correlation coefficient for each target gene between the RT-qPCR and their corresponding RNAseq data is presented in Fig. 4. Most target genes show strong positive correlation with the data of the biological replicate of the LC50 curve, but not strong positive correlation with data of either non-exposed larvae or non-lethal Cry11Aa exposed larvae. Overall correlation of this set of genes is 79 % with the Cry11Aa LC50 biological replicate, which is more than the closest value of 43 % with Cry11Aa V142D data. These results indicate that our RNAseq data and observations have good experimental validation. Results also indicate that the expression values observed are more likely a response to the successful formation of toxin pores in the membranes of the intestinal epithelial cells than to the presence of the toxin, and that changes observed are not merely due to the growth or circadian regulation of genes in the mosquito larvae during the time of the experiment.

Fig. 4
figure4

RNAseq and RT-qPCR correlation. Transcript abundance ratios were determined for each gene by SYBR Green RT-qPCR on a biological replicate of a 12 h Cry11Aa exposure curve, a control curve of non-toxin exposed A. aegypti larvae, or larvae exposed to non-lethal Cry11Aa mutants. Spearman correlation was determined between corresponding RNAseq values and log2 transformed abundance ratios. Strong positive correlations (above 0.5) are colored in white, weak positive correlation (between 0.3 and 0.5) are colored in light gray, and low, negative or no correlation values are colored in dark gray. Genes are indicated to be up regulated or down regulated as determined by RNAseq analysis. The percentage of genes with strong positive correlation between each curve and RNAseq data is shown

As an additional validation analysis of RNAseq data we determined the changes in expression of JNK protein by western blot analysis (Additional file 6: Figure S2). RNAseq data showed that JNK (AAEL008622) was up regulated after 9 and 12 h of Cry11Aa exposure. Our analysis indicated a direct correlation of the expression of JNK protein with JNK (AAEL008622) gene expression after 9 and 12 h of Cry11Aa exposure (Additional file 6: Figure S2).

Discussion

Here we show the analysis of the transcriptional response of A. aegypti midgut epithelial cells after the first 12 h of Cry11Aa exposure. With intoxication, cellular responses are likely diverse in an attempt by the host to avert toxicity. We took care not to include severely intoxicated or dead larvae in the sequenced samples. That way we could say that the larvae analyzed were mildy intoxicated and maybe able to defend from the toxin up to the point when they were sampled. We found that the longer the exposure the more genes were differentially expressed. This is consistent with previous microarray data obtained (Gill SS, data not published), where high doses of toxin generated a widespread repression of genes while an LD50 or lower showed a general profile favoring positive differential expression. Also in the coleopteran pest Tenebrio molitor the timepoint where the most differential gene expression was observed was after 12 h incubation with Cry3Aa toxin [24].

Due to the experimental bioassays conditions we cannot be sure that all individual larvae consume Cry toxin at the same rate and in the same amount. It is possible that in the shorter exposure times not all larvae had fed Cry11Aa crystals, which could explain why the overall transcriptional response is milder than at longer exposure times. However, at 9 and 12 h of Cry11Aa toxin exposure we saw a strong cellular response and a similar profile of enriched biological functions. The control larvae at 0 and 12 h in absence of toxin administration showed a limited number of DEG, indicating that the responses observed here correspond to transcriptional changes due to Cry11Aa toxin action.

Up-regulated responses

Our results indicate that exposure to an LC50 of Cry11Aa induces the expression of genes related to vesicular traffic, protein folding, lipid metabolism and cell membrane metabolism. Chaperone, phospholipase A2 and carboxylesterase activities were also found among the induced functions at longer intoxication times. Small GTPases of the Ras superfamily are known regulators of cell junctions, cytoskeleton and vesicular traffic in cells [37, 38]. In particular, they play key roles in maintaining epithelial homeostasis [38], as observed in mammalian intestinal epithelium [31] and also insect epithelium [39]. Here we observed that components of the ARF and Rab family are induced, as well as many nucleotide exchangers and regulators of Rho signaling, especially at 12 h of exposure. These data, coupled with the induction of genes related to actin nucleation and cell junction establishment, point to an active process of epithelial tissue homeostasis. We also found homologs of Drosophila melanogaster crimpled and snakeskin genes induced at 12 h; both have a role in maintenance of cell junctions [40, 41]. It is important to mention that these genes were shown to be expressed in the gut of D. melanogaster [40, 41]. The gene for Rab-11 was induced in A. aegypti after 12 h of toxin ingestion which is consistent with previous reports that showed that Rab-5 and Rab-11 had an important role in the response to Cry5B in C. elegans [42]. Rab-5 and Rab-11 of C. elegans are thought to regulate a pathway through which Cry5B exposed cells may remove and recycle damaged membranes, since it was shown that absence of either of these two proteins results in a hypersensitive phenotype to Cry toxin exposure [42]. Endocytosis is also utilized in the clearance of membranes damaged by toxin pores of S. aureus α-toxin [15] and S. pyogenes streptolysin-O [43]. Thus, while the specific genes involved may vary the overall cellular responses detected in this work are similar to the responses observed in other organism to different pore forming toxins.

Small GTPase signaling networks involve MAPK proteins and, since genes related to these transduction modules are induced, it could partly explain the cellular response to Cry toxins provided by JNK and p38 [31]. It was previously reported the participation of the p38 MAPK signaling pathway in the response and defense against Cry11Aa toxin in A. aegypti [18, 25].

JNK and p38 have been shown to regulate a great number of genes in response to Cry5B in C. elegans [18]. The p38 MAPK pathway may be inducing expression of Ets domain containing transcription factors, which in turn could be the reason for an up regulation of metalloproteases [44]. Our data also suggest the involvement of signaling via PIP3 in the larval defense response. Arachidonic acid, the product of phospholipase A2 activity, has been reported to stimulate cell adhesion via a novel MAPK p38-RhoA pathway [45]. Both arachidonic acid and PIP3 are messenger metabolites derived from membrane phospholipids. These data may indicate a link between some enriched functional processes detected in our analysis, and the potential role that they may play in defending A. aegypti from the effects of pore forming toxins. However, more precise biochemical experiments are required to evaluate the cross-talk between these signaling pathways.

Down-regulation response

A. aegypti’s defense response involves the down-regulation of several processes. Among these we identify genes related detoxification proteins such a cytochrome 450 and ABC transporters. Most of these down-regulated processes were observed at 12 h of toxin exposure. It is notable that similar expression patterns were observed in Bti tolerant A. aegypti strains [46], in the Cry11Aa resistant A. aegypti [47], as well as in chemical insecticide resistant strains of A. aegypti [48]. It was suggested that this down-regulation may be a metabolic tradeoff for the energy spent on other defense processes [46, 48]. In the absence of pathogen-associated molecular patterns (PAMPs) from bacterial components to stimulate innate immunity in our purified Cry11Aa crystals, our results suggest that the cells divert metabolism to membrane remodeling and tissue junction sealing, rather than functions for which no stimuli is detected. This energy saving measures could explain the initial down-regulation of cell proliferation, although this process seems to be induced at the last time-point analyzed, perhaps to substitute irreparably damaged cells after wound-healing. Our analysis also revealed an overall repression of transmembrane transport, especially that concerning sulfate metabolism and carbohydrates, which seems inconsistent with attempts to reestablish solute gradients after the disruption produced by Cry11Aa pores at the cell membrane. However, the cell appears to shift first to restore membrane integrity and repress any energy consuming activities not related to this process.

As mentioned above we found genes with ABC transporter domains down-regulated DEG. This is relevant since some Cry1Ac resistance phenotypes in different lepidopteran species, such as Plutella xylostella, Bombyx mori, Heliothis virescens and Helicoverpa armigera larvae are linked to mutations in the ABCC2 transporters [4951], and negative differential expression of genes involved in detoxification activity, including ABC transporters and cytochrome 450, has been found in mosquitoes tolerant to Bti [46]. Many detoxification processes require energy consumption to pump offending agents out of the cell. We show here that when membrane is damage due to the pore formation of the toxin, the cell responds by resealing the cell membrane and removing damaged parts of the membrane instead of activation of transmembrane transports and detoxification. A possible explanation is that down regulation of detoxification-related genes could be an energy saving process as previously suggested for Bti- and Cry11Aa- resistant Ae. aegypti mosquitoes that also showed under-transcription of enzymes classically involved in the detoxification [4648].

Comparison with mosquito resistant to Bti toxins and with Cry toxin intoxications in different insect orders

It is interesting to mention that the convergence of some processes from the transitory response to Cry11Aa intoxication with that of a genetically established resistance to Bti toxins [4648] may suggests that resistance involves the selection of mutations that could confer a metabolic state similar to a constitutive toxin defense response. Awareness of the pathways utilized by midgut cells to contend with toxin damage may allow us to better predict possible resistance mechanisms that could occur by field application of Bt formulations.

Importantly, there is a great deal of overlap in genome wide transcription changes and functional enrichment of A. aegypti reported here with that found for intoxication of T. molitor with Cry3Aa [24]. For example, both analyses found that toxin treated larvae show a repression of enzymes related to antioxidant processes, as well as functions associated with DNA replication and cell division. Carbohydrate and tricarboxylic acid enzymes are also repressed in both organisms. Signaling through MAPK and small GTPases are induced in both organisms upon intoxication, as are functions related to maintaining cell integrity and cytoskeleton complexes. All of these data suggest that insects from different orders respond similarly to Cry toxin action, supporting a conserved defense response to Cry toxins in these insect species.

Implications in the Cry11Aa mechanism of action

Finally, it appears that only fully functional pore forming Cry11Aa toxins is able to elicit the full range of transcriptional changes observed in this work. Neither Cry11Aa-mutant toxins tested in RT-qPCR assays showed high positive correlation with RNAseq data, where Cry11Aa E97A, which is unable to form a stable oligomeric pre-pore structure, was most different. The higher correlation observed with Cry11Aa V142D for some genes is similar to that observed for non-toxin exposed larvae, but could also be due to non-complete loss of function. The oligomeric structure formed by this toxin may not form a pore stable enough to cause a lethal phenotype, but that is nonetheless sufficient to be sensed by the cellular machinery.

Our data indicated that binding to the cell membrane was not enough to provoke a response, involving an intracellular cascade activation of a receptor coupled protein G that activates adenylate cyclase and PKA activities as was proposed in Zhang et al., 2006 [52]. First, our results with the mutant toxins indicated that full pore formation is required to elicit the entire response program, at least for the genes tested. We observed activation of PIP3 response and small monomeric GTPases families which are found in the cytosol and have a molecular weight of about 21 kDa, in contrast to the proposed participation of a heterotrimeric G proteins, specifically from G s family that participates in activation of adenylate cyclase [52]. Finally, we did not observed changes in PKA expression.

Conclusions

The burden imposed on global human health by mosquito borne diseases is tremendous. Biological control of mosquito populations through the use of Bt toxins is an important part of strategies to diminish this burden. However, it is important to consider emergence of resistant insect vector populations due to widespread use of Bti preparations. Although Bti has been used for decades worldwide seemingly without such resistance events surfacing, the risk remains. In this work we took a glimpse at the physiological response that the A. aegypti midgut mounts upon action of a Cry toxin, and how the epithelial cells therein change their gene expression to defend and restore tissue integrity. This knowledge contributes to the field of Bt toxins in several ways. First, it complements and extends earlier transcriptomic studies performed on other insects exposed to Cry that may allow us to determine differences and similarities in response between distinct insect species. Second, it improves our understanding of the overall molecular interactions and the mode of action of Cry toxins with target cells, building upon the in vitro biochemical data with data of in vivo effects. Lastly, we now possess information that may be helpful to improve our use of Bti toxins for the biological control of mosquitoes. Knowing what are the defense response pathways employed by mosquito larvae to survive the activity of Bt toxins may potentially help in the design of better strategies or formulations that could limit the defense response in mosquito An alternative could be the use of inhibitors of pathways that may weaken the mosquito’s defense mechanisms. Another alternative is the release of mosquitoes with knock-out or mutant alleles in proteins of the defense pathways. It is clear that these ideas are just speculations that require much further work before their field application to minimize the impact of vector-borne diseases on human populations.

Methods

Production, purification of Cry11Aa toxins and insect bioassay

Bt strains were grown for three days on HCT medium supplemented with 10 μg/ml of erythromycin as previously described [5] for the production of crystals of Cry11Aa or its mutants E97A and V142D [35, 36]. Crystals were obtained by washing the cultures with 0.3 M NaCl/0.01 M EDTA and purified by discontinuous sucrose gradients [53]. The medium lethal concentration (LC50) of wild-type Cry11Aa pure crystals was determined after 24 h exposure of 4th instar A. aegypti larvae to different doses in a response curve in triplicate as previously reported [5] and calculated by ProBit analysis (POLO-Plus, LeOra Software).

Intoxication and tissue dissection

For each time of each biological replicate of the intoxication assay, a dose of 500 ng Cry11Aa pure crystals per ml was administered to three cups containing each ten 4th-instar A. aegypti larvae in 100 ml of water, these experiments were performed four times. After 3, 6, 9 and 12 h of incubation with Cry11Aa crystals the midgut tissue was dissected from the larvae, taking care to remove the food bolus. Control conditions of non-intoxicated larvae at time point zero and 12 h without toxin were also analyzed in triplicate. Dead larvae or larvae that were severely affected by intoxication were discarded. Each biological replicate was conducted at the same time of the day to minimize gene expression changes due to circadian rhythms. The dissected midguts were placed in RNAlater (Qiagen) for subsequent processing and stored at −70 °C per the manufacturer’s instructions.

RNA extraction and sequencing

Strength of differential expression analysis increases with the number of replicates used to estimate the expression levels, for this reason four biological replicates were sequenced on an Illumina platform. RNAlater buffer was removed from dissected midguts and the tissues homogenized with 27G syringes in TRIzol reagent (Ambion, Life Technologies). Cellular debris was removed by centrifugation at 12,000 xg for 10 min at 4 °C, and the supernatant transferred to new tubes. RNA was extracted following the manufacturer’s instructions. Total RNA was dissolved in DEPC treated water and quantified by Nanodrop 1000 (Thermo Scientific). RNA was treated with DNAse I (Thermo Scientific), and samples were cleaned with RNeasy columns (QIAGEN) following the kit’s protocol. The integrity and concentration of eluted RNA were determined by Bioanalyzer 2100 (Agilent Technologies). For sequencing, mRNA libraries were prepared from 2 μg of total RNA using standard Illumina protocols using the TruSeq RNA Sample Preparation Kit. The mRNA libraries were sequenced in 72 bp pair-ended format on a Genome Analyzer IIx (Illumina) at the Unidad Universitaria de Secuenciación Masiva y Bioinformática (UUSMB) of the Instituto de Biotecnología-UNAM facilities, or in 100 bp pair-ended format on a Illumina Hi-Seq 2000 at Laboratorio Nacional de Genómica para la Biodiversidad (LANGEBIO) of Centro de investigación y de Estudios Avanzados del Instituto Politécnico Nacional (CINVESTAV-IPN) facilities. Quality control of reads was performed at UUSMB facilities. Read quality control include data for GC content, for overrepresented reads (for example, contaminating ribosomal RNA reads). Quality per base and more parameters were obtained using the FastQC utility. A specific read is flagged as overrepresented if that identical sequence appears above the 0.1 % of total reads obtained in the sequencing of an individual sample. Longer reads were trimmed with FastX-Trimmer to 72 pb to remove adapter and lower quality bases at the ends.

Differential expression analysis

RNAseq reads were aligned to A. aegypti transcripts v3.2 (Vectorbase) [54] using TopHat 2 v 2.0.12 [55] in transcriptome mode and base features file v 3.2, allowing for 2 mismatches. For each sample, a count table for all genes was generated with HTSeq Count v0.5.4p5 [56]. Aligned reads where filtered to remove reads with multiple alignments or ambiguous assignment by executing the script in “intersection non-empty” mode. Only reads with unique alignments were further analyzed. A count table for all samples of four biological replicates was compiled and used for differential expression analysis by the BioconductoR packages DEseq2 v 1.4.5 [29] and EdgeR v 3.6.8 [30], both in generalized linear model mode that included both time of toxin exposure and biological replicate of the experiment as possible sources of variation to control for differences among individuals dissected on different days. Results for DEG were obtained from comparisons of each intoxication time to non-intoxicated control. A gene was considered significant if the adjusted p-value was < 0.05 for DESeq2 or if False Discovery Rate was < 0.05 for EdgeR. For each intoxication time DEG reported by both tools were recovered by a custom R script.

Functional analysis

GO enrichment analysis of biological process annotations was performed with the BioconductoR package TopGO v 2.16.0 [57] using the “Weight” algorithm. GO annotations for all genes were retrieved from Vectorbase, formatted, and used to detect enrichment of GO terms in up or down regulated DEG from each intoxication condition on genes common to both differential expression analysis tools. We used the “Weight” algorithm because it yields more informative GO terms due to the way it processes the GO ontology hierarchy, and allowed us to obtain information of more particular metabolic processes. A p-value of 0.1 was used as cutoff for enriched terms. GO annotations terms with less than five annotated genes were excluded from analysis to reduce statistical artifacts. Genes without a GO biological process term were further analyzed, and their predicted Interpro functional domains were retrieved from Vectorbase. Frequency for domains was counted with a custom Perl script and those above the 90th percentile were kept. Finally, we conducted pBLAST queries [34] against NR database with genes for which no Interpro domain could be retrieved. Only matches having more that 40 % of query identity and an e-value of < e-8 were considered significant.

Clustering analysis

Through DEseq2 a variance stabilized transformation was performed on all log2 fold changes for all samples. These data were filtered to retain only those genes that were determined to have significant differential expression in at least one treatment and were part of the genes common to DESeq2 and EdgeR analysis. Data for each gene were summarized as mean of replicates per treatment. Clustering of individual DEGs that followed the same expression pattern over the 4 time points was performed with Cluster 3.0 (Stanford University), centering values by mean and using complete linkage and euclidian distances as parameters. Java Treeview v1.1.6r4 [58] was used to visualize results and to recover clusters of genes with common expression profiles throughout the intoxication experiment. GO enrichment analysis was performed within these clusters as before.

RT-qPCR target selection and validation

Targets for validation by RT-qPCR were selected from DEG. These targets were selected according to three criteria: high statistical significance of differential expression, absolute log2 fold change larger than 2, and statistical significance in more than one tested condition. Some target genes do not cover all criteria but were nonetheless selected due to their very high significance or log2 fold change. Primers were designed to amplify products of 100–180 bp. All primer pairs were designed in the same exon with exception of two pairs of primers for AAEL003633 and AAEL008192 genes, which were located in two exons separated by an intron. For cDNA template synthesis from RNA samples we used Superscript III Reverse Transcriptase (Invitrogen) following the manufacturer’s instructions. For all primers, a control amplification was performed with the total RNA in absence of reverse transcriptase, in these controls the signal was absent or extremely low indicating that contamination of our total RNA samples was negligible.

Dynamic range for signal detection was determined by SYBR-green (Thermo Scientific) RT-qPCR amplification of some selected genes on cDNA prepared using 5000, 500, 50, 5 and 0.5 μg of total RNA from A. aegypti midgut on Illumina Eco equipment. To determine primer amplification efficiency, 1.5 μg of total RNA from each time of a LC50 intoxication replicate were mixed for cDNA synthesis. This template was used for SYBR green RT-qPCR amplification in 10-fold serial dilutions for each primer pair.

RT-qPCR was performed (Roche Light Cycler 480) using a Cry11Aa LC50 intoxication biological replicate that was not sequenced. The RT-qPCR data for all 14 selected genes presented in Additional file 5: Figure S1 were obtained for 0, 3, 6, 9 and 12 h after administration of Cry11Aa. In the case of Cry11Aa Cry11Aa E97A, and Cry11Aa V142D mutants the RT-qPCR data were obtained after 0, 9, and 12 h since these times were the ones that showed most important significant changes in expression. The controls of no-toxin exposed larvae were also done at 0, 9 and 12 h. Primer efficiency correction was used in ΔΔCq relative quantitation analysis. The ribosomal protein S3 gene was used for reference and normalization. Transcript abundance ratios between conditions were log2 transformed and the Spearman correlation determined with corresponding RNAseq log2 fold data for each gene.

Western blot

Midguts were dissected from 4th instar A. aegypti larvae and pools of 50 entire midguts in 100 μl of rehydration buffer (7 M urea, 2 M thiourea, 4 % CHAPS, 40 mM dithiothreitol, 0.5 % pharmalyte or IPG buffer [GE Life Sciences], 0.002 % bromophenol blue, 2.5 ml) containing protease inhibitors (Complete, Roche Diagnostics) and were kept at −80 °C until processed. Midgut pools were homogenized with a motorized pellet pestle (Motor Sigma-Aldrich Z359971.1EA) on ice. Ten μg of midgut protein was fractionated by 10 % SDS-PAGE and transferred to nitrocellulose membrane. Western blot was performed as described by Aranda et al., (1996) [59]. JNK protein was detected after incubation of 1 h with polyclonal anti-JNK (1/5,000; JNK2 N-18: sc-827 Santacruz Biotech. Inc). Visualization was performed with goat anti-rabbit secondary antibody coupled to horseradish peroxidase (Amersham) (1/10,000; 1 h), followed by Super Signal chemiluminescence Luminol substrate (Pierce), according to the manufacturer’s instructions.

Data accessibility

Raw RNA-seq data and gene expression abundance measurements for this study have been deposited at Gene Expression Omnibus under accession code GSE74785. Read data can accessed directly in Sequence Read Archive with accession code SRP065731.

Abbreviations

ALP:

alkaline phosphatases

APN:

aminopeptidases

Bt :

Bacillus thuringiensis

Bti :

Bacillus thuringiensis subsp. israelensis

DEG:

differentially expressed genes

WHO:

World Health Organization

References

  1. 1.

    WHO: Dengue and Severe Dengue, Fact sheet No 117 WHO. 2015 http://www.who.int/mediacentre/factsheets/fs117/en/

  2. 2.

    Jirakanjanakit N, Saengtharatip S, Rongnoparut P, Duchon S, Bellec C, Yoksan S. Trend of temephos resistance in Aedes (Stegomyia) mosquitoes in Thailand during 2003–2005. Environ Entomol. 2007;36(3):506–11.

    Article  PubMed  Google Scholar 

  3. 3.

    Devine G, Furlong M. Insecticide use: Contexts and ecological consequences. Agric Hum Values. 2007;24(3):281–306.

    Article  Google Scholar 

  4. 4.

    Bravo A, Likitvivatanavong S, Gill SS, Soberon M. Bacillus thuringiensis: A story of a successful bioinsecticide. Insect Biochem Mol Biol. 2011;41(7):423–31.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  5. 5.

    Perez C, Fernandez LE, Sun J, Folch JL, Gill SS, Soberon M, et al. Bacillus thuringiensis subsp. israelensis Cyt1Aa synergizes Cry11Aa toxin by functioning as a membrane-bound receptor. Proc Natl Acad Sci U S A. 2005;102(51):18303–8.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  6. 6.

    Canton PE, Zanicthe Reyes EZ, de Escudero R, Bravo A, Soberon M. Binding of Bacillus thuringiensis subsp. israelensis Cry4Ba to Cyt1Aa has an important role in synergism. Peptides. 2011;32(3):595–600.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  7. 7.

    Wirth MC, Georghiou GP, Federici BA. CytA enables CryIV endotoxins of Bacillus thuringiensis to overcome high levels of CryIV resistance in the mosquito, Culex quinquefasciatus. Proc Natl Acad Sci U S A. 1997;94(20):10536–40.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  8. 8.

    Pardo-Lopez L, Soberon M, Bravo A. Bacillus thuringiensis insecticidal three-domain Cry toxins: mode of action, insect resistance and consequences for crop protection. FEMS Microbiol Rev. 2013;37(1):3–22.

    Article  CAS  PubMed  Google Scholar 

  9. 9.

    Likitvivatanavong S, Chen J, Evans AM, Bravo A, Soberon M, Gill SS. Multiple receptors as targets of Cry toxins in mosquitoes. J Agric Food Chem. 2011;59(7):2829–38.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  10. 10.

    Fernandez-Luna MT, Lanz-Mendoza H, Gill SS, Bravo A, Soberon M, Miranda-Rios J. An alpha-amylase is a novel receptor for Bacillus thuringiensis ssp. israelensis Cry4Ba and Cry11Aa toxins in the malaria vector mosquito Anopheles albimanus (Diptera: Culicidae). Environ Microbiol. 2010;12(3):746–57.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  11. 11.

    Gurcel L, Abrami L, Girardin S, Tschopp J, van der Goot FG. Caspase-1 activation of lipid metabolic pathways in response to bacterial pore-forming toxins promotes cell survival. Cell. 2006;126(6):1135–45.

    Article  CAS  PubMed  Google Scholar 

  12. 12.

    Mestre MB, Colombo MI. Autophagy and toxins: a matter of life or death. Curr Mol Med. 2013;13(2):241–51.

    Article  CAS  PubMed  Google Scholar 

  13. 13.

    Keyel PA, Loultcheva L, Roth R, Salter RD, Watkins SC, Yokoyama WM, et al. Streptolysin O clearance through sequestration into blebs that bud passively from the plasma membrane. J Cell Sci. 2011;124(Pt 14):2414–23.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  14. 14.

    Husmann M, Dersch K, Bobkiewicz W, Beckmann E, Veerachato G, Bhakdi S. Differential role of p38 mitogen activated protein kinase for cellular recovery from attack by pore-forming S. aureus alpha-toxin or streptolysin O. Biochem Biophys Res Commun. 2006;344(4):1128–34.

    Article  CAS  PubMed  Google Scholar 

  15. 15.

    Husmann M, Beckmann E, Boller K, Kloft N, Tenzer S, Bobkiewicz W, et al. Elimination of a bacterial pore-forming toxin by sequential endocytosis and exocytosis. FEBS Lett. 2009;583(2):337–44.

    Article  CAS  PubMed  Google Scholar 

  16. 16.

    Popova TG, Millis B, Bradburne C, Nazarenko S, Bailey C, Chandhoke V, et al. Acceleration of epithelial cell syndecan-1 shedding by anthrax hemolytic virulence factors. BMC Microbiol. 2006;6:8.

    PubMed Central  Article  PubMed  Google Scholar 

  17. 17.

    Stassen M, Muller C, Richter C, Neudorfl C, Hultner L, Bhakdi S, et al. The streptococcal exotoxin streptolysin O activates mast cells to produce tumor necrosis factor alpha by p38 mitogen-activated protein kinase- and protein kinase C-dependent pathways. Infect Immun. 2003;71(11):6171–7.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  18. 18.

    Kao CY, Los FC, Huffman DL, Wachi S, Kloft N, Husmann M, et al. Global functional analyses of cellular responses to pore-forming toxins. PLoS Pathog. 2011;7(3):e1001314.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  19. 19.

    Kloft N, Busch T, Neukirch C, Weis S, Boukhallouk F, Bobkiewicz W, et al. Pore-forming toxins activate MAPK p38 by causing loss of cellular potassium. Biochem Biophys Res Commun. 2009;385(4):503–6.

    Article  CAS  PubMed  Google Scholar 

  20. 20.

    Ratner AJ, Hippe KR, Aguilar JL, Bender MH, Nelson AL, Weiser JN. Epithelial cells are sensitive detectors of bacterial pore-forming toxins. J Biol Chem. 2006;281(18):12994–8.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  21. 21.

    Stringaris AK, Geisenhainer J, Bergmann F, Balshusemann C, Lee U, Zysk G, et al. Neurotoxicity of pneumolysin, a major pneumococcal virulence factor, involves calcium influx and depends on activation of p38 mitogen-activated protein kinase. Neurobiol Dis. 2002;11(3):355–68.

    Article  CAS  PubMed  Google Scholar 

  22. 22.

    Nagahama M, Shibutani M, Seike S, Yonezaki M, Takagishi T, Oda M, et al. The p38 MAPK and JNK pathways protect host cells against Clostridium perfringens beta-toxin. Infect Immun. 2013;81(10):3703–8.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  23. 23.

    Huffman DL, Abrami L, Sasik R, Corbeil J, van der Goot FG, Aroian RV. Mitogen-activated protein kinase pathways defend against bacterial pore-forming toxins. Proc Natl Acad Sci U S A. 2004;101(30):10995–1000.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  24. 24.

    Oppert B, Dowd SE, Bouffard P, Li L, Conesa A, Lorenzen MD, et al. Transcriptome profiling of the intoxication response of Tenebrio molitor larvae to Bacillus thuringiensis Cry3Aa protoxin. PLoS One. 2012;7(4):e34624.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  25. 25.

    Cancino-Rodezno A, Alexander C, Villasenor R, Pacheco S, Porta H, Pauchet Y, et al. The mitogen-activated protein kinase p38 is involved in insect defense against Cry toxins from Bacillus thuringiensis. Insect Biochem Mol Biol. 2010;40(1):58–63.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  26. 26.

    Rodriguez-Cabrera L, Trujillo-Bacallao D, Borras-Hidalgo O, Wright DJ, Ayra-Pardo C. Molecular characterization of Spodoptera frugiperda-Bacillus thuringiensis Cry1Ca toxin interaction. Toxicon. 2008;51(4):681–92.

    Article  CAS  PubMed  Google Scholar 

  27. 27.

    Meunier L, Prefontaine G, Van Munster M, Brousseau R, Masson L. Transcriptional response of Choristoneura fumiferana to sublethal exposure of Cry1Ab protoxin from Bacillus thuringiensis. Insect Mol Biol. 2006;15(4):475–83.

    Article  CAS  PubMed  Google Scholar 

  28. 28.

    Sparks ME, Blackburn MB, Kuhar D, Gundersen-Rindal DE. Transcriptome of the Lymantria dispar (gypsy moth) larval midgut in response to infection by Bacillus thuringiensis. PLoS One. 2013;8(5):e61190.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  29. 29.

    Love MI, Huber W, Anders S: Moderated estimation of fold change and dispersion for RNA-Seq data with DESeq2Genome Biol. 2014;15(12):550

  30. 30.

    Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26(1):139–40.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  31. 31.

    Citalan-Madrid AF, Garcia-Ponce A, Vargas-Robles H, Betanzos A, Schnoor M. Small GTPases of the Ras superfamily regulate intestinal epithelial homeostasis and barrier function via common and unique mechanisms. Tissue Barriers. 2013;1(5):e26938.

    PubMed Central  Article  PubMed  Google Scholar 

  32. 32.

    Likitvivatanavong S, Chen J, Bravo A, Soberón M, Gill SS. Cadherin, alkaline phosphatase and aminopeptidase N as receptors of Cry11Ba toxin from Bacillus thuringiensis jegathesan in Aedes aegypti. Appl Environ Microbiol. 2011;77(1):24–31.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  33. 33.

    Mitchell A, Chang HY, Daugherrty L, Fraser M, Hunter S, Lopez R, et al. The InterPro protein families database: the classification resource after 15 years. Nucleic Acids Res. 2015;43(Database issue):D213–21.

    PubMed Central  Article  PubMed  Google Scholar 

  34. 34.

    Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. 1990;215(3):403–10.

    Article  CAS  PubMed  Google Scholar 

  35. 35.

    Munoz-Garay C, Rodriguez-Almazan C, Aguilar JN, Portugal L, Gomez I, Saab-Rincon G, et al. Oligomerization of Cry11Aa from Bacillus thuringiensis has an important role in toxicity against Aedes aegypti. Appl Environ Microbiol. 2009;75(23):7548–50.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  36. 36.

    Carmona D, Rodriguez-Almazan C, Munoz-Garay C, Portugal L, Perez C, de Maagd RA, et al. Dominant negative phenotype of Bacillus thuringiensis Cry1Ab, Cry11Aa and Cry4Ba mutants suggest hetero-oligomer formation among different Cry toxins. PLoS One. 2011;6(5):e19952.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  37. 37.

    Nielsen E, Cheung AY, Ueda T. The regulatory RAB and ARF GTPases for vesicular trafficking. Plant Physiol. 2008;147(4):1516–26.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  38. 38.

    Citi S, Guerrera D, Spadaro D, Shah J. Epithelial junctions and Rho family GTPases: the zonular signalosome. Small GTPases. 2014;5(4):1–15.

    Article  PubMed  Google Scholar 

  39. 39.

    Warner SJ, Longmore GD. Distinct functions for Rho1 in maintaining adherens junctions and apical tension in remodeling epithelia. J Cell Biol. 2009;185(6):1111–25.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  40. 40.

    Yanagihashi Y, Usui T, Izumi Y, Yonemura S, Sumida M, Tsukita S, et al. Snakeskin, a membrane protein associated with smooth septate junctions, is required for intestinal barrier function in Drosophila. J Cell Sci. 2012;125(Pt 8):1980–90.

    Article  CAS  PubMed  Google Scholar 

  41. 41.

    Nilton A, Oshima K, Zare F, Byri S, Nannmark U, Nyberg KG, et al. Crooked, coiled and crimpled are three Ly6-like proteins required for proper localization of septate junction components. Development. 2010;137(14):2427–37.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  42. 42.

    Los FC, Kao CY, Smitham J, McDonald KL, Ha C, Peixoto CA, et al. RAB-5- and RAB-11-dependent vesicle-trafficking pathways are required for plasma membrane repair after attack by bacterial pore-forming toxin. Cell Host Microbe. 2011;9(2):147–57.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  43. 43.

    Corrotte M, Fernandes MC, Tam C, Andrews NW. Toxin pores endocytosed during plasma membrane repair traffic into the lumen of MVBs for degradation. Traffic. 2012;13(3):483–94.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  44. 44.

    Ghosh S, Basu M, Roy SS. ETS-1 protein regulates vascular endothelial growth factor-induced matrix metalloproteinase-9 and matrix metalloproteinase-13 expression in human ovarian carcinoma cell line SKOV-3. J Biol Chem. 2012;287(18):15001–15.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  45. 45.

    Garcia MC, Ray DM, Lackford B, Rubino M, Olden K, Roberts JD. Arachidonic acid stimulates cell adhesion through a novel p38 MAPK-RhoA signaling pathway that involves heat shock protein 27. J Biol Chem. 2009;284(31):20936–45.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  46. 46.

    Despres L, Stalinski R, Tetreau G, Paris M, Bonin A, Navratil V, et al. Gene expression patterns and sequence polymorphisms associated with mosquito resistance to Bacillus thuringiensis israelensis toxins. BMC Genomics. 2014;15:926.

    PubMed Central  Article  PubMed  Google Scholar 

  47. 47.

    Lee S-B, Aimanova KG, Gill SS. Alkaline phosphatase and aminopeptidase are altered in a Cry11Aa resistant strain of Aedes aegypti. Insect Biochem Mol Biol. 2014;54:112–21.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  48. 48.

    Despres L, Stalinski R, Faucon F, Navratil V, Viari A, Paris M, et al. Chemical and biological insecticides select distinct gene expression patterns in Aedes aegypti mosquito. Biol Lett. 2014;10(12):20140716.

    PubMed Central  Article  PubMed  Google Scholar 

  49. 49.

    Guo Z, Kang S, Zhu X, Xia J, Wu Q, Wang S, et al. Down-regulation of a novel ABC transporter gene (Pxwhite) is associated with Cry1Ac resistance in the diamondback moth, Plutella xylostella (L.). Insect Biochem Mol Biol. 2015;59:30–40.

    Article  CAS  PubMed  Google Scholar 

  50. 50.

    Baxter SW, Badenes-Perez FR, Morrison A, Vogel H, Crickmore N, Kain W, et al. Parallel evolution of Bacillus thuringiensis toxin resistance in lepidoptera. Genetics. 2011;189(2):675–9.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  51. 51.

    Gahan LJ, Pauchet Y, Vogel H, Heckel DG. An ABC transporter mutation is correlated with insect resistance to Bacillus thuringiensis Cry1Ac toxin. PLoS Genet. 2010;6(12):e1001248.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  52. 52.

    Zhang X, Candas M, Griko NB, Taussig R, Bulla LA. A mechanism of cell death involving an adenylyl cyclase/PKA signaling pathway is induced by the Cry1Ab toxin of Bacillus thuringiensis. Proc Natl Acad Sci USA. 2006;103:9897–902.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  53. 53.

    Thomas WE, Ellar DJ. Bacillus thuringiensis var israelensis crystal delta-endotoxin: effects on insect and mammalian cells in vitro and in vivo. J Cell Sci. 1983;60:181–97.

    CAS  PubMed  Google Scholar 

  54. 54.

    Giraldo-Calderon GI, Emrich SJ, MacCallum RM, Maslen G, Dialynas E, Topalis P, et al. VectorBase: an updated bioinformatics resource for invertebrate vectors and other organisms related with human diseases. Nucleic Acids Res. 2015;43(Database issue):D707–713.

    PubMed Central  Article  PubMed  Google Scholar 

  55. 55.

    Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg SL. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 2013;14(4):R36.

    PubMed Central  Article  PubMed  Google Scholar 

  56. 56.

    Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010;11(10):R106.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  57. 57.

    Alexa A, Rahnenfuhrer J: topGO: Enrichment analysis for Gene Ontology. In., 2.6.0 edn; 2010: R package version 2.23.0 https://www.bioconductor.org/packages/3.3/bioc/html/topGO.html.

  58. 58.

    Saldanha AJ. Java Treeview--extensible visualization of microarray data. Bioinformatics. 2004;20(17):3246–8.

    Article  CAS  PubMed  Google Scholar 

  59. 59.

    Aranda E, Sanchez J, Peferoen M, Guereca L, Bravo A. Interactions of Bacillus thuringiensis crystal proteins with the midgut epithelial cells of Spodoptera frugiperda (Lepidoptera: Noctuidae). J Invertebr Pathol. 1996;68(3):203–12.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgments

We thank Lizbeth Cabrera, and Ricardo Grande, Alejandro Sánchez, Verónica Jiménez and Jêrome Verleyen of the UUSMB of the Instituto de Biotecnología-UNAM, and LANGEBIO of CINVESTAV facilities for their technical assistance. Research was funded in part through grant NIH 2R01 AI066014, CONACYT 0221200, PEC acknowledges CONACyT for PhD fellowship.

Author information

Affiliations

Authors

Corresponding author

Correspondence to Alejandra Bravo.

Additional information

Competing interests

Authors declare that they have no competing interests.

Authors’ contributions

PEC performed most of the experimental and bioinformatics analysis of this study. ACR performed larvae dissection, protein extraction and western blotting. PEC, AB, MS and SSG participate in data analysis, wrote and corrected the manuscript. PEC, MS and AB conceived and designed the project. All authors read and approved the final manuscript.

Additional files

Additional file 1: Table S1.

TopHat2 alignment rates. Number of total reads is rounded to closest million. (DOCX 73 kb)

Additional file 2: Table S2.

Complete list of 1060 differentially expressed genes across all experimental conditions (3-12 h of Cry11Aa toxin administration) and gene description (if available), by expression profile cluster. (DOCX 325 kb)

Additional file 3: Table S3.

List of differentially expressed genes in control (non-toxin exposed larvae) RNAseq data at 12 h. (DOCX 68 kb)

Additional file 4: Table S4.

List of all Interpro domains found in differential expressed genes without biological process GO annotation. (DOCX 92 kb)

Additional file 5: Figure S1.

RNAseq and log2 transformed RT-qPCR fold changes. Data is presented for RNAseq and RT-qPCR at 3, 6, 9, and 12 h of an LC50 Cry11Aa treatment. Also shown are RT-qPCR values for 9 and 12 h of unexposed larvae or Cry11Aa non-toxic mutants treatments. All log2 fold changes are referred to control larvae at the start of respective treatment. Panel A and B show genes with high correlation between RNAseq data and RT-qPCR of Cry11Aa. Panel C shows the three genes with low or no positive correlation between these data. (ZIP 779 kb)

Additional file 6: Figure S2.

Characterization of JNK expression by Western Blot analysis using α-JNK antibody. Presence of JNK was determined on protein extracts of midguts of A. aegypti larvae exposed to Cry11Aa for 9 or 12 h (lanes 4 and 5, respectively). Proteins were separated by SDS-PAGE and transferred to PVDF membranes for Western Blotting with human αJNK2 antibody. Non-toxin exposed larvae dissected at corresponding times were used as controls (lanes 1 to 3). (TIF 1042 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Canton, P.E., Cancino-Rodezno, A., Gill, S.S. et al. Transcriptional cellular responses in midgut tissue of Aedes aegypti larvae following intoxication with Cry11Aa toxin from Bacillus thuringiensis . BMC Genomics 16, 1042 (2015). https://doi.org/10.1186/s12864-015-2240-7

Download citation

Keywords

  • Aedes aegypti
  • Bacillus thuringiensis subsp israelensis
  • RNA-seq
  • biological control
  • defense response
  • Cry toxins