Open Access

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

  • Pablo Emiliano Canton1,
  • Angeles Cancino-Rodezno2,
  • Sarjeet S. Gill3,
  • Mario Soberón1 and
  • Alejandra Bravo1Email author
BMC Genomics201516:1042

https://doi.org/10.1186/s12864-015-2240-7

Received: 7 September 2015

Accepted: 24 November 2015

Published: 9 December 2015

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.

Keywords

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

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)

Condition

 

EdgeR

DESeq2

Common

Only EdgeR

Only DESeq2

3 h

all

10

5

4

6

1

Up regulated

3

1

1

2

0

Down regulated

7

4

3

4

1

6 h

all

33

340

27

6

313

Up regulated

24

245

23

1

222

Down regulated

9

95

4

5

91

9 h

all

151

458

145

6

313

Up regulated

105

273

104

1

169

Down regulated

46

185

41

5

144

12 h

all

1059

1980

1038

21

942

Up regulated

623

1123

613

10

510

Down regulated

436

857

425

11

432

 

all

39

16

13

26

3

12 h no toxin

Up regulated

20

7

4

16

3

 

Down regulated

19

9

9

10

0

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

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

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

Cluster

GO.ID

Term

Annotated

Significant

Expected

ClassicFisher

A. Moderate down regulated

GO:0055114

oxidation-reduction process

769

75

29.4

2.10E-13

GO:0008272

sulfate transport

10

5

0.38

1.70E-05

GO:0006081

cellular aldehyde metabolic process

9

3

0.34

0.0039

GO:0032787

monocarboxylic acid metabolic process

70

8

2.68

0.0067

GO:0006749

glutathione metabolic process

14

3

0.54

0.0147

GO:0005978

glycogen biosynthetic process

6

2

0.23

0.0197

GO:1901700

response to oxygen-containing compound

6

2

0.23

0.0197

GO:0044270

cellular nitrogen compound catabolic process

385

17

14.72

0.0235

GO:0046700

heterocycle catabolic process

386

17

14.76

0.0264

GO:0006099

tricarboxylic acid cycle

18

3

0.69

0.0295

GO:0044262

cellular carbohydrate metabolic process

41

6

1.57

0.0318

GO:1901361

organic cyclic compound catabolic process

388

17

14.83

0.0326

GO:0006730

one-carbon metabolic process

8

2

0.31

0.035

GO:1901292

nucleoside phosphate catabolic process

357

13

13.65

0.0511

GO:0055085

transmembrane transport

576

34

22.02

0.0619

GO:0043648

dicarboxylic acid metabolic process

11

2

0.42

0.0638

GO:0007205

protein kinase C-activating G-protein coupled receptor signaling pathway

11

2

0.42

0.0638

GO:0044106

cellular amine metabolic process

20

3

0.76

0.0743

GO:0006979

response to oxidative stress

28

3

1.07

0.0897

B. Strong down regulated

GO:0006952

defense response

14

2

0.01

0.0021

GO:1902476

chloride transmembrane transport

6

1

0

0.0034

GO:0042742

defense response to bacterium

9

1

0.01

0.0051

C. Moderate up regulated

GO:0071702

organic substance transport

266

27

10.92

2.10E-05

GO:0008610

lipid biosynthetic process

103

13

4.23

0.00025

GO:0016192

vesicle-mediated transport

114

15

4.68

0.00091

GO:0043623

cellular protein complex assembly

32

7

1.31

0.00103

GO:0042398

cellular modified amino acid biosynthetic process

20

5

0.82

0.00105

GO:0007229

integrin-mediated signaling pathway

8

3

0.33

0.00329

GO:0006195

purine nucleotide catabolic process

341

24

14

0.00404

GO:0006955

immune response

17

4

0.7

0.00433

GO:0042742

defense response to bacterium

9

3

0.37

0.00478

GO:0032940

secretion by cell

18

4

0.74

0.00539

GO:0009207

purine ribonucleoside triphosphate catabolic process

340

24

13.96

0.00684

GO:0006777

Mo-molybdopterin cofactor biosynthetic process

11

3

0.45

0.00884

GO:0030029

actin filament-based process

35

6

1.44

0.00941

GO:0010466

negative regulation of peptidase activity

24

4

0.99

0.01544

GO:0006108

malate metabolic process

5

2

0.21

0.01547

GO:0007267

cell-cell signaling

14

3

0.57

0.0178

GO:0034613

cellular protein localization

136

11

5.58

0.01829

GO:0007264

small GTPase mediated signal transduction

189

12

7.76

0.02202

GO:0006470

protein dephosphorylation

54

6

2.22

0.02244

GO:0019751

polyol metabolic process

6

2

0.25

0.02258

GO:1901617

organic hydroxy compound biosynthetic process

6

2

0.25

0.02258

GO:0016051

carbohydrate biosynthetic process

26

4

1.07

0.04051

GO:0006457

protein folding

80

7

3.28

0.04524

GO:0007017

microtubule-based process

81

7

3.33

0.04786

GO:0046130

purine ribonucleoside catabolic process

340

24

13.96

0.04798

GO:0048285

organelle fission

34

4

1.4

0.04903

GO:0000165

MAPK cascade

9

2

0.37

0.04995

GO:0006270

DNA replication initiation

9

2

0.37

0.04995

GO:0033014

tetrapyrrole biosynthetic process

10

2

0.41

0.06079

GO:0046939

nucleotide phosphorylation

11

2

0.45

0.07233

GO:0016998

cell wall macromolecule catabolic process

11

2

0.45

0.07233

GO:0000902

cell morphogenesis

11

2

0.45

0.07233

GO:0044262

cellular carbohydrate metabolic process

41

4

1.68

0.08597

GO:0006879

cellular iron ion homeostasis

13

2

0.53

0.09725

GO:0006826

iron ion transport

13

2

0.53

0.09725

GO:0046854

phosphatidylinositol phosphorylation

13

2

0.53

0.09725

D. Strong up regulated

GO:0016042

lipid catabolic process

51

2

0.11

0.005

GO:0034314

Arp2/3 complex-mediated actin nucleation

7

1

0.01

0.015

GO:0070588

calcium ion transmembrane transport

17

1

0.04

0.035

GO:0032319

regulation of Rho GTPase activity

24

1

0.05

0.05

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

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

Condition

InteproID

Percentile

Description

6 h down regulated

IPR007588

1

Zync finger, FLYWCH-type

6 h up regulated

IPR002018

0.94

Carboxylesterase, typeB

IPR029058

1

Alpha/Beta hydrolase fold

9 h down regulated

IPR007087

0.9

Zync finger, C2H2

9 h up regulated

IPR013126

0.966

Heat shock protein 70 family

IPR029058

0.955

Alpha/Beta hydrolase fold

IPR002018

0.944

Carboxylesterase, typeB

IPR013783

0.933

Immunoglobulin-like fold

12 h down regulated

IPR007087

0.996

Zync finger, C2H2

IPR011009

0.977

Protein kinase-like domain

IPR013026

0.936

Tetratricopeptide repeat-containing domain

IPR013525

0.936

ABC-2 type transporter

12 h up regulated

IPR029058

1

Alpha/Beta hydrolase fold

IPR008978

0.997

HSP20-like chaperone

IPR027417

0.988

P-loop containing nucleoside triphoshpate hydrolase

IPR013126

0.969

Heat shock protein 70 family

IPR002048

0.969

EF-hand domain

IPR002018

0.958

Carboxylesterase, typeB

IPR001849

0.947

Pleckstrin homology domain

IPR015880

0.947

Zync finger C2H2-like

IPR002110

0.925

Ankyrin repeat

IPR001452

0.925

SH3 domain

IPR008985

0.925

Concanavalin A-like lectin/glucanases superfamily

IPR018502

0.925

Annexin repeat

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

Condition

Gene ID

Coverage

E-value

Identity

Homolog ID

Protein

Species

9 h Down regulated

AAEL013895

90 %

1e-32

53 %

XP_001868961.1

14.5 kDa salivary protein

Cq

 

AAEL009985

94 %

2e-59

65 %

XP_001868961.1

14.5 kDa salivary protein

Cq

 

AAEL000777

83 %

3e-15

60 %

ABU80079.1

CEC-A

Ag

 

AAEL004522

98 %

8e-29

65 %

ACA05605.1

gambicin

Am

9 h up regulated

AAEL007206

92 %

2e-22

48 %

XP_004526025.1

arhytmically expressed gene 5 protein-like

Cc

 

AAEL000811

91 %

1e-127

48 %

ETN67412.1

methusalah

Ad

 

AAEL002963

32 %

5e-09

57 %

XP_008549507.1

aras guanine nucleotide exchange factor L isoform X8

Md

12 h Down regulated

AAEL006228

98 %

9e-102

45 %

XP_001846856.1

PPARgamma constitutive conactivator 1

Cq

 

AAEL012622

98 %

5e-108

68 %

XP_001862944.1

ADAM metalloprotease

Cq

 

AAEL007780

96 %

2e-107

57 %

AAV90698.1

putative salivary secreted protein 275

Aa

 

AAEL009985

94 %

2e-59

65 %

XP_001868961.1

14.5 kDa salivary protein

Cq

 

AAEL006346

81 %

5e-62

55 %

XP_005177200.1

aBAG family molecular chaperone regulator 2-like isoform X1

Mud

 

AAEL003594

98 %

1e-169

61 %

XP_001862361.1

microtubule binding protein D-CLIP-190

Cq

 

AAEL013356

30 %

8e-62

62 %

XP_003401795.1

aeukaryotic initiation factor 4A-like isoform 2

Bi

12 h Up regulated

AAEL006493

98 %

7e-11

49 %

AAV90653.1

putative salivary secreted peptide

Aa

 

AAEL005874

100 %

4e-63

73 %

XP_003707525.1

aE3 ubiquitin-protein ligase AMFR-like

Mr

 

AAEL006493

98 %

7e-11

49 %

AAV90653.1

putative salivary secreted peptide

Aa

 

AAEL013220

75 %

3e-62

72 %

NP_648500.1

crimpled

Dm

 

AAEL011867

92 %

0

40 %

XP_003395464.1

amuscle M-line assembly protein unc-89 isoform X3

Apm

 

AAEL014310

86 %

8e-34

79 %

NP_648366.2

calcium-independent phospholipase A2 VIA, isoform B

Dm

 

AAEL007206

92 %

2e-22

48 %

XP_004526025.1

arhythmically expressed gene 5 protein

Cc

 

AAEL001576

94 %

1e-31

83 %

ETN61409.1

cysteine string protein

Ad

 

AAEL003567

99 %

6e-43

49 %

AAV90653.1

putative salivary secreted peptide

Aa

 

AAEL002578

97 %

1e-59

65 %

NP_649184.1

snakeskin, isoform A

Dm

 

AAEL000811

91 %

1e-127

48 %

ETN67412.1

methusalah

Ad

 

AAEL005428

98 %

2e-32

51 %

ABV44740.1

13.6 kDa midgut protein

Pp

 

AAEL002963

32 %

5e-09

57 %

XP_008549503.1

amyb-like protein J isoform X4

Md

Aa Aedes albopictus, Ad Anopheles darlingi, Ag Anopheles gambiae, Am Anopheles melas, Apm Apis mellifera, Bi Bombus impatiens, Cc Ceratitis capitata, Cq Culex quinquefasciatus, Dm Drosophila melanogaster, Md Microplitis demolitor, Mud Musca domestica, Mr Megachile roduntata, Pp Phlebotomus papatasi

aindicates protein classified as “PREDICTED”

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

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

Declarations

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.

Open AccessThis 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.

Authors’ Affiliations

(1)
Departamento de Microbiología, Instituto de Biotecnología, Universidad Nacional Autónoma de México
(2)
Facultad de Ciencias, Universidad Nacional Autónoma de México
(3)
Cell Biology and Neuroscience Department, University of California, Riverside

References

  1. WHO: Dengue and Severe Dengue, Fact sheet No 117 WHO. 2015 http://www.who.int/mediacentre/factsheets/fs117/en/
  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.View ArticlePubMedGoogle Scholar
  3. Devine G, Furlong M. Insecticide use: Contexts and ecological consequences. Agric Hum Values. 2007;24(3):281–306.View ArticleGoogle Scholar
  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 CentralView ArticlePubMedGoogle Scholar
  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 CentralView ArticlePubMedGoogle Scholar
  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 CentralView ArticlePubMedGoogle Scholar
  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 CentralView ArticlePubMedGoogle Scholar
  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.View ArticlePubMedGoogle Scholar
  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 CentralView ArticlePubMedGoogle Scholar
  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 CentralView ArticlePubMedGoogle Scholar
  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.View ArticlePubMedGoogle Scholar
  12. Mestre MB, Colombo MI. Autophagy and toxins: a matter of life or death. Curr Mol Med. 2013;13(2):241–51.View ArticlePubMedGoogle Scholar
  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 CentralView ArticlePubMedGoogle Scholar
  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.View ArticlePubMedGoogle Scholar
  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.View ArticlePubMedGoogle Scholar
  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 CentralView ArticlePubMedGoogle Scholar
  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 CentralView ArticlePubMedGoogle Scholar
  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 CentralView ArticlePubMedGoogle Scholar
  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.View ArticlePubMedGoogle Scholar
  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 CentralView ArticlePubMedGoogle Scholar
  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.View ArticlePubMedGoogle Scholar
  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 CentralView ArticlePubMedGoogle Scholar
  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 CentralView ArticlePubMedGoogle Scholar
  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 CentralView ArticlePubMedGoogle Scholar
  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 CentralView ArticlePubMedGoogle Scholar
  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.View ArticlePubMedGoogle Scholar
  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.View ArticlePubMedGoogle Scholar
  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 CentralView ArticlePubMedGoogle Scholar
  29. Love MI, Huber W, Anders S: Moderated estimation of fold change and dispersion for RNA-Seq data with DESeq2Genome Biol. 2014;15(12):550Google Scholar
  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 CentralView ArticlePubMedGoogle Scholar
  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 CentralView ArticlePubMedGoogle Scholar
  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 CentralView ArticlePubMedGoogle Scholar
  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 CentralView ArticlePubMedGoogle Scholar
  34. Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. 1990;215(3):403–10.View ArticlePubMedGoogle Scholar
  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 CentralView ArticlePubMedGoogle Scholar
  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 CentralView ArticlePubMedGoogle Scholar
  37. Nielsen E, Cheung AY, Ueda T. The regulatory RAB and ARF GTPases for vesicular trafficking. Plant Physiol. 2008;147(4):1516–26.PubMed CentralView ArticlePubMedGoogle Scholar
  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.View ArticlePubMedGoogle Scholar
  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 CentralView ArticlePubMedGoogle Scholar
  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.View ArticlePubMedGoogle Scholar
  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 CentralView ArticlePubMedGoogle Scholar
  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 CentralView ArticlePubMedGoogle Scholar
  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 CentralView ArticlePubMedGoogle Scholar
  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 CentralView ArticlePubMedGoogle Scholar
  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 CentralView ArticlePubMedGoogle Scholar
  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 CentralView ArticlePubMedGoogle Scholar
  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 CentralView ArticlePubMedGoogle Scholar
  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 CentralView ArticlePubMedGoogle Scholar
  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.View ArticlePubMedGoogle Scholar
  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 CentralView ArticlePubMedGoogle Scholar
  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 CentralView ArticlePubMedGoogle Scholar
  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 CentralView ArticlePubMedGoogle Scholar
  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.PubMedGoogle Scholar
  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 CentralView ArticlePubMedGoogle Scholar
  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 CentralView ArticlePubMedGoogle Scholar
  56. Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010;11(10):R106.PubMed CentralView ArticlePubMedGoogle Scholar
  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. Saldanha AJ. Java Treeview--extensible visualization of microarray data. Bioinformatics. 2004;20(17):3246–8.View ArticlePubMedGoogle Scholar
  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.View ArticlePubMedGoogle Scholar

Copyright

© Canton et al. 2015