Gene expression patterns and sequence polymorphisms associated with mosquito resistance to Bacillus thuringiensis israelensis toxins
© Després et al.; licensee BioMed Central Ltd. 2014
Received: 27 June 2014
Accepted: 16 October 2014
Published: 23 October 2014
Despite the intensive use of Bacillus thuringiensis israelensis (Bti) toxins for mosquito control, little is known about the long term effect of exposure to this cocktail of toxins on target mosquito populations. In contrast to the many cases of resistance to Bacillus thuringiensis Cry toxins observed in other insects, there is no evidence so far for Bti resistance evolution in field mosquito populations. High fitness costs measured in a Bti selected mosquito laboratory strain suggest that evolving resistance to Bti is costly. The aim of the present study was to identify transcription level and polymorphism variations associated with resistance to Bti toxins in the dengue vector Aedes aegypti. We used RNA sequencing (RNA-seq) for comparing a laboratory-selected strain showing elevated resistance to Bti toxins and its parental non-selected susceptible strain. As the resistant strain displayed two marked larval development phenotypes (slow and normal), each phenotype was analyzed separately in order to evidence potential links between resistance mechanisms and mosquito life-history traits.
A total of 12,458 genes were detected of which 844 were differentially transcribed between the resistant and susceptible strains. Polymorphism analysis revealed a total of 68,541 SNPs of which 12,571 SNPs exhibited more than 40% frequency difference between the resistant and susceptible strains, affecting 2,953 genes. Bti resistance is associated with changes in the transcription level of enzymes involved in detoxification and chitin metabolism. Among previously described Bti-toxin receptors, four alkaline phosphatases (ALPs) were differentially transcribed between resistant and susceptible larvae, and non-synonymous changes affected the protein sequence of one cadherin, six aminopeptidases (APNs) and four α-amylases. Other putative Cry receptors located in lipid rafts, such as flotillin and glycoside hydrolases, were under-transcribed and/or contained non-synonymous substitutions. Finally, immunity-related genes showed contrasted transcription and polymorphisms patterns between the two developmental resistant phenotypes, suggesting the existence of trade-offs between Bti-resistance, life-history traits and immunity.
The present study is the first to analyze the whole transcriptome of Bti-resistant mosquitoes by RNA-seq, shedding light on the importance of studying both transcription levels and sequence polymorphism variations to get a comprehensive view of insecticide resistance.
KeywordsRNA-seq Bacillus thuringiensis israelensis toxins Mosquito Bio-insecticide resistance Toxin receptors Lipid rafts Resistance costs Evolutionary trade-offs
Biological control strategies involving Bacillus thuringiensis israelensis (Bti) are increasingly used as alternatives to chemicals in order to limit the spreading of mosquito-borne pathogens transmission and nuisance. During sporulation, Bti produces a composite crystalline inclusion composed of four main toxins: Cry4Aa, Cry4Ba, Cry11Aa and Cyt1Aa. Upon larval ingestion, the crystalline inclusions dissolve in the alkaline pH of the insect’s gut, releasing Cry and Cyt protoxins that will be processed to active toxins by insect gut proteases . Cyt and Cry toxins are known to act in synergy to kill mosquito larvae, but the precise mechanisms of their interaction are not fully understood [2, 3]. Other Bacillus thuringiensis subspecies produce Cry toxins that are specific to their target organism, and most research has been done on Cry toxins mode of action in various insects (mostly in lepidopterans and coleopterans, but also in dipterans) and nematodes [4, 5]. The selectivity of Cry toxins is mainly due to the interaction with specific receptors present at the surface of midgut epithelium cells of the larva. The activated Cry toxins bind to specific protein receptors, insert into the membrane and form pores, resulting in midgut disruption, bacterial infection and death of the insect. In the recent years, several models have been proposed to describe the precise mode of action of Bt Cry toxins (reviewed in [6, 7]). The binding of Cry toxins can involve interaction with multiple specific receptors, including cadherin-like proteins but also aminopeptidases N (APNs) and alkaline phosphatases (ALPs) which are mostly located in lipid rafts (reviewed in ). In mosquitoes, besides cadherins, APNs and ALPs, α-amylases were also shown to bind to Cry4Ba and Cry11Aa toxins . As compared to what is known about the mode of action of Bt Cry toxins in lepidopterans and coleopterans, the understanding of Bti mode of action in mosquito lies far behind, mainly due to the lack of Bti-resistant mosquito populations to dissect the mechanisms of resistance. Indeed, the cocktail of toxins produced by Bti, and especially the presence of Cyt toxin, appear to hinder the evolution of resistance to Bti both in the laboratory and in the field . However, recent works conducted in our group revealed that Cyt1Aa toxin was strongly affected by the organic matter found in most mosquito breeding sites, inducing a decreased efficacy and a low-level persistence of Cyt toxins as compared to Cry toxins [10, 11]. This ultimately led to the modification of toxins proportions in field-persistent Bti (Cry toxins became the most abundant toxins while it is Cyt toxin in commercial Bti), as shown in leaf litters containing persistent Bti sampled in mosquito breeding sites several months after a treatment . This field-collected Bti was used to select a laboratory strain for 22 generations, resistant to Bti toxins (LiTOX strain ). Two phenotypes were observed in the LiTOX strain: larvae with a normal development time (5 days as the parental Bora-Bora strain) named LiTOX_N, and larvae with a slow development (about 10 days) named LiTOX_S. As compared to the susceptible parental Bora-Bora strain, LiTOX_N exhibited resistance ratio (RR50) of 30.2-fold, 13.7-fold and 6.3-fold for Cry4Aa, Cry4Ba and Cry11Aa respectively. LiTOX_S exhibited the same resistance ratio as LiTOX_N for Cry4Ba and Cry11Aa but was not resistant to Cry4Aa . In addition to an altered development, the resistant LiTOX strain exhibited lower female fecundity and lower egg survival, suggesting that Bti resistance is associated with high fitness costs . In the present study, we used the LiTOX strain to undertake a whole transcriptome analysis by RNA-seq of the two developmental phenotypes LiTOX_S and LiTOX_N versus the susceptible parental strain aiming at dissecting the molecular pathways to resistance to Bti Cry toxins at the whole organism level, and at understanding the trade-offs underlying Cry-toxin resistance and life history traits. Transcription level and polymorphism variations associated with each resistance phenotype were identified in regards of known and new putative mechanisms involved in Bti toxins resistance and the way they might interfere with mosquito life history traits is discussed.
Mapping and re-annotation results
A total of 197,175,220 short reads (75 bp) were sequenced, with an average of 33 million reads per sample. An average mapping rate of ~83% (164.66 million mapped reads) was reached. After filtering on base quality (150.86 million reads left) and on mapping quality, 129.32 million reads were retained for further analysis (Additional file 1: Table S1). A total of 12,942 transcripts showed a coverage >0.5 Reads Per Kilobase exon Model (RPKM) in at least one strain, including 2,717 transcripts showing a different structure as compared to genome annotation (21.8%) and 484 novel transcription events (NTEs) (3.7%) predicted based on spliced reads and distance from known transcripts (close or far, respectively). NTEs were not considered in further analyses. Because only 6,996 detected transcripts were functionally annotated in VectorBase (genome version AaegL2.1, vectorbase.org) we performed a Blast2GO analysis of all predicted peptides annotated as ‘conserved hypothetical protein’ or ‘hypothetical protein’ in Vectorbase against the protein database Swissprot (Blastp, E-value < 10-3, annotation cut-off = 55, GO weight = 5). This allowed re-annotation of 2,608 transcripts (Additional file 2: Table S2). We also re-annotated known mosquito Bti-toxin receptors based on bibliography [14, 15], ending up with about 77% of annotated transcripts in our dataset, including 80 putative Cry receptors (15 α-amylases, 14 ALPs, 12 cadherins and 39 APNs), and 40 trypsins/chymotrypsins potentially involved in Cry protoxin activation, and thereafter called ‘Bti candidates’. All the other transcripts were assigned to nine categories based on their annotation as follows: ‘Detoxification’, ‘Other enzymes’, ‘Cuticle’, ‘Immunity’, ‘Hormones and neurotransmitters signaling’, ‘Transcription factors’, ‘Intra or extracellular trafficking/chaperonins’, ‘Structure’ and ‘Unknown’.
Differential transcription analysis overview: Number (and proportion) of transcripts significantly over or under-transcribed in each LiTOX phenotype and altogether (LiTOX-N or LiTOX-S) as compared to the control parental strain
AC* test P value and FC > =2
Novel transcription event
Novel transcription event
SNPs detection and analysis
Non-synonymous changes in putative Cry toxin receptors, with their frequency and coverage in the susceptible and resistant phenotypes
Amino acid position
Amino acid change
A- > S
T- > I
M- > I
V- > I
A- > T
E- > A
T- > A
M- > T
D- > N
T- > M
E- > D
H- > Y
G- > R
K- > E
V- > M
A- > E
R- > K
V- > A
A- > S
N- > S
H- > L
T- > S
S- > N
I- > F
K- > R
R- > W
A- > V
A- > V
Two main mechanisms of resistance to Cry toxins have been described in insects so far: altered protoxin activation by gut proteases, and modification in transcription level and/or protein sequence of Cry receptors resulting in lower or failure in toxin binding [19–21]. The study of resistance to Cry toxins in insects has mostly been based on the study of a few candidate proteins first identified in lepidopterans, either toxin receptors or proteases, using proteomic approaches (2D-DIGE, ligand blotting, etc.). As compared to other Bacillus thuringiensis subspecies that typically produce a single Cry toxin, Bti produces a mixture of four different toxins, and resistance to Bti is likely to involve multiple distinct mechanisms, which may affect the physiology and metabolism of resistant insects. Altered metabolism in Bti resistant mosquitoes was revealed by fitness costs expressed at any mosquito life-stage (larval development time, fecundity, egg survival) . These costs can result from pleiotropic effects of resistant alleles on life-history traits (direct cost) or more indirectly from resource reallocation between different metabolic pathways.
Global changes in gene expression and metabolic trade-offs
The two LiTOX phenotypes share common global transcriptional responses: an over-transcription of genes involved in proteolytic activity, cell cytokinesis and wound healing, and an under-transcription of enzymes classically involved in the detoxification of chemicals such as chemical insecticides, plants allelochemicals and pollutants (cytochrome P450 monooxygenases, esterases, glutathion S-transferases, UDP-glycosyltransferases, etc.). Conversely, Bti toxins are large proteins and their degradation requires a proteolytic processing rather than the activity of detoxification enzymes. Therefore, in the absence of chemical selection pressure, the under-transcription of these enzymes might represent a metabolic trade-off in the LiTOX strain to compensate the increased expression of enzymes involved in protein degradation . The over-transcription of several proteases (endopeptidases, serine proteases) is in accordance with an observed increase in proteolytic activities of midgut enzymes of LiTOX larvae . This increased proteolytic activity in Bti resistant larvae could reflect a higher ability to degrade toxins. Despite this global trend, particular proteases such as the metalloproteinases AAEL011547, AAEL011550 and AAEL014514, and the chymotrypsin AAEL015105 or the trypsin AAEL006376, which are putatively involved in Cry protoxin activation, were under-transcribed in both LiTOX phenotypes. The lack of metalloproteinase activity in the LiTOX strain was evidenced by the analysis of enzymatic activities in the midgut of larvae using various specific enzyme inhibitors . Down-regulation of specific proteases might result in a lower and/or improper toxin activation which could also confer resistance. In lepidopterans, analysis of proteinase activities from gut extracts have shown the lack of trypsin-like and of chymotrypsin-like activity in Bt-resistant strains . Furthermore, the role of proteases in resistance to Bt toxin Cry1Ca1 was demonstrated in the caterpillar Spodoptera frugiperda where a serine-protease gene under-transcribed in the midgut of intoxicated larvae was further shown by RNAi-mediated knockdown to be involved in both reduced protoxin activation in the midgut and reduced susceptibility of caterpillars to toxins in bioassays .
Lipid rafts and Bti-toxins trafficking through the gut membrane
Lipid rafts are membrane domains enriched in glycosylphosphatidylinositol (GPI)-anchored proteins that have been implicated in transmembrane protein trafficking and were shown to play a role as portals of entry for many pathogens including viruses, bacteria and their toxins [25, 26]. Several studies suggested that lipid rafts might play a central role in Bt-toxins toxicity. For example, the presence of lipid rafts was shown to be essential for the toxicity of Cry1C on Sf9 Lepidoptera cells . However, the specific proteins associated with this mechanism have not been identified yet. A flotillin, a protein found in lipid rafts and described as a toxin transporter , was under-transcribed in both LiTOX phenotypes. Flotillin was shown to co-localize with Cry4Ba toxin in Aedes aegypti larval gut brush border membranes , and further investigation of its Cry-binding ability will contribute to characterize its role in Bti toxin resistance. Four glycoside-hydrolases were under-transcribed in both LiTOX phenotypes. Although never described as Cry toxin receptors, glycoside-hydrolases were shown to be associated with lipid rafts in the gut membrane of Aedes larvae . Lipid rafts appear to play a central role in pore-forming toxins insertion in the gut membrane by functioning as platforms to recruit GPI-anchored proteins including APNs, ALPs, and glycoside-hydrolases. The role of some GPI-anchored proteins (flotillin and glycoside-hydrolases) largely present in lipid rafts of Aedes larvae guts has been so far underexplored in studies aiming at understanding Cry toxins toxicity in insects, and the present study pinpoints them as new potential candidates in Aedes.
Differences between the two LiTOX phenotypes
Several gene families were differentially transcribed between the two developmental phenotypes: genes involved in cuticle/chitin metabolism, prophenoloxidases, and genes coding for enzymes involved in lipid metabolism.
In our dataset, the category of cuticle/chitin proteins contained the highest number of genes differentially transcribed (Figures 2 and 4b, Additional file 3: Table S3). Cuticle is an exoskeleton constituted of chitin that provides physical support and protects the insects from physical injuries, desiccation and infections. Cuticle is typically divided into hard cuticle, generally sclerotized and mechanically stiff, and soft cuticle, a more flexible region that allows locomotion and appendix movements . Chitin is always associated with chitin-binding proteins; although their exact function is to be investigated, it has been hypothesized that they could play important roles in cuticle metabolism and structure preservation . Nearly all the genes identified as differentially transcribed in this category contained the R&R chitin-binding motif and were associated with the soft cuticle (RR-1 cuticle proteins) or with the hard cuticle (RR-2) . The global under-transcription of these genes in the LiTOX_S could partly explain the slow development of this cohort, considering that a low chitin metabolism is believed to slow-down the molting process and increase the development duration. Nevertheless, it is intriguing why only genes encoding R&R-containing proteins were found differentially transcribed while genes encoding other chitin-binding proteins with peritrophin-A chitin-binding domains (ChtBD2 domain) and chitin metabolism enzymes (such as chitinases) were not differentially expressed. Also, it is challenging to know if this pattern is a side-effect consequence of the resistance phenotype or if it is part of the mechanisms of resistance developed by the insects against Cry toxins. A better understanding of the role of such proteins in the cuticle metabolism is mandatory to help deciphering their potential involvement in resistance phenotypes in insects.
Prophenoloxidases are involved in the melanization process, and as such they have a dual role: they contribute to cuticle formation (and larval development), and they participate in the immune response against pathogens by their encapsulation . In Ae. aegypti, both defensins and phenoloxidases were shown to participate to immune response against bacterial challenge . In our dataset, prophenoloxidases were under-transcribed in both LiTOX phenotypes, while defensins were under-transcribed only in LiTOX_N (Figure 2, Additional file 3: Table S3). Furthermore, a correlation was observed between age at pupation (development speed) and mosquito's ability to melanize Sephadex beads (immuno-capacity), suggesting that an increased immunity is costly and results in a slower development . The two developmental phenotypes observed in the resistant mosquito strain LiTOX might be the result of a trade-off between immunity and development.
When focusing on candidate genes, only few differences in transcription level and no significant difference in polymorphism were observed between the two LiTOX phenotypes (e.g. between LiTOX_N and LiTOX_S). The serine protease AAEL010139 was not expressed in LiTOX_N, the Cry4Aa-resistant cohort. Among the putative Cry-toxin receptors, four alkaline phosphatases (ALPs) were differentially transcribed between the two LiTOX phenotypes: ALP AAEL003286 was specifically under-transcribed in LiTOX_N, which may be linked to the high Cry4Aa resistance of this phenotype. Three other ALP were specifically under- (AAEL003309 and AAEL009077) or over-transcribed (AAEL000931) in LiTOX_S. Cry4Aa mode of action and mechanisms of resistance are the least investigated and understood among the Bti toxins. Functionally validating the involvement of these few candidate genes could help lifting the veil on the specific mechanisms associated with Cry4Aa resistance and provide a better understanding of the patterns of cross-resistance between Bti Cry toxins.
Polymorphisms rather than expression changes are found in putative Bti-toxins receptors
Among the 80 putative Cry receptors detected in larvae, only 5 were differentially transcribed in the LiTOX strain as compared to the susceptible strain. Among these, the α-amylase AAEL010540 was under-transcribed in both LiTOX phenotypes sharing resistance to Cry4Ba and Cry11Aa. In An. albimanus, an α-amylase (aamy1) has been described as a Bti toxin receptor, binding both Cry4Ba and Cry11Aa . Although this requires functional validation, this suggests that the Ae. aegypti α-amylase AAEL010540 (48% identity with aamy1) could bind Cry toxins. No cadherin or APN that were previously associated with Cry resistance in insects [34–37] or showed to bind Bti Cry toxins in mosquitoes [38–40] were differentially transcribed in the LiTOX strain, suggesting that if they are involved in Bti resistance in LiTOX, this is not related to an altered expression but rather to changes in protein sequence affecting their affinity for Cry toxins. Indeed, one cadherin, eight APNs and four α-amylases displayed differential SNPs leading to non-synonymous changes between the resistant and susceptible strains.
CAD-like proteins have been the most intensively studied putative Cry toxin receptor molecules in lepidopteran and coleopteran larvae [41–43]. It has been proposed that they act as the first receptors of Cry toxins, binding toxin monomers and facilitating further processing required for the pre-pore oligomer formation . Although twelve CAD-like proteins were detected in our dataset, none was differentially transcribed, and only one (AAEL018140) contained several differential SNPs (Table 2). This cadherin was previously shown by proteomic approaches to act as a Cry11Aa toxin receptor . It contains a N-terminal signal peptide, 11 cadherin repeats (CR1 to CR11), a membrane-proximal region (MPR) and a trans-membrane domain. The only differential non-synonymous SNP distinguishing Bti-resistant and susceptible strains (H160L) is located in the N-terminal region preceding the first cadherin repeat (CR1), while overlay assays and immune-blotting localized Cry11Aa and Cry4Aa toxin-binding regions in CR8-CR11 . More generally, Cry-toxin binding sites were usually described in the membrane-proximal region of insect cadherins, not in the N-terminal extracellular region . However, amino-acid substitution at any location can dramatically change the secondary structure of a protein; for instance, the three cadherin alleles conferring resistance to Bt in the pink bollworm Ectinophora gossypiella are not located within the Cry1Aa binding site but upstream . In this case, the deletion of several amino-acids upstream the toxin binding site was shown to alter the full-length cadherin and impede toxin binding. In the H160L mutation observed at relatively high frequency in the LiTOX strain (49-58%), a polar and positively charged amino acid (Histidine) is replaced by a non-charged hydrophobic amino acid (Leucine), which can result in a change in the spatial arrangement of the cadherin in the cell membrane. Further functional studies and in-silico protein modelling could allow investigating whether this change affects cadherin conformation and Cry toxin binding affinity. Altogether, our results suggest that cadherin might not have a central role in resistance to Bti in mosquitoes.
APNs and α-amylases
APNs and α-amylases are digestive enzymes that play a key role in larval nutrition and development, being involved in protein and glucose metabolism, respectively. APNs have been extensively studied as putative Cry toxin receptors in many insects and were classified into 8 classes based on their sequence identity [8, 14]. APNs from several classes were reported to bind with more or less affinity to the different Cry toxins in various insects, including AeAPN1 in Aedes and APN2 in Anopheles, but their precise role in Bti toxicity remains to be elucidated. Out of the 31 APNs detected in Ae. aegypti larvae, 12 contained differential SNPs, and non-synonymous differential changes affected the sequence of 8 of them (Table 2). In the APN2 AAEL005808 recently shown to be a functional receptor of Cry4Ba toxin , 4 non-synonymous changes were nearly fixed in the LiTOX strain (both phenotypes), suggesting ongoing selective sweep on this gene. One non-synonymous change affected the sequence of the APN1 AAEL012778 previously shown to bind Cry11Aa . Although not usually cited as a Cry-toxin receptor, an α-amylase was recently shown to bind Cry4Ba and Cry11Aa in Anopheles. Our dataset pinpoints five α-amylases potentially involved in Aedes resistance to Bti toxins: one α-amylase was under-transcribed only in LiTOX_N (resistant to all Cry toxins), and four other α-amylases presented sequence changes in both cohorts of the resistant LiTOX strain as compared to the susceptible strain. Because of their key role in insect nutrition, altered expression of these digestive enzymes as a resistance mechanism is likely to be costly in terms of larval development, and this might explain why the observed changes in these enzymes are mostly amino-acid changes rather than changes in expression levels. Indeed, change in the protein sequence might dramatically alter toxin binding site while having small effects on digestive efficiency. Further studies on the differential abilities of these α-amylases to bind Bti Cry toxins and to act as functional toxin receptors are required to determine their role in Bti toxicity. Also, molecular modelling of APN is necessary to identify the oligosaccharide that are involved in the binding of Cry toxins to verify if the SNPs described in this study can alter the Cry-APN interaction and partly explain the resistance phenotype observed .
The present study is the first to analyze the whole transcriptome of Bti-resistant larvae by deep sequencing, allowing studying change in gene expression level and sequence polymorphisms using the same dataset. It reveals dramatic modifications in the transcriptional profiles selected in resistant larvae, with a global transcriptional increase of genes coding for proteases and chitin-binding proteins, and a decrease in transcription of genes involved in detoxification and immunity. Whether these modifications are directly involved in Bti-toxin processing or more indirectly through complex metabolic compensations selected to limit the cost of resistance remains to be functionally investigated. Regarding candidate Bti-toxin receptors, our dataset relativizes the role of cadherin in resistance to Bti in mosquitoes, and highlights the importance of studying changes in the sequence rather than in transcription levels of APN and ALP, but also of other proteins involved in protein trafficking through the cell membrane such as flotillins and glucoside hydrolases.
Strains and Sample preparation
The Ae. aegypti laboratory strain Bora-Bora susceptible to all insecticides was selected using field-collected leaf-litter containing Bti toxins (LiTOX strain) as described in . At each generation, an average of 6,000 larvae was exposed to toxic leaf litter in order to obtain a mortality rate of 70%, so that about 1,800 adults constituted the next generation, limiting bottleneck effects in the LiTOX strain. Mosquito strains were reared in standard insectarium conditions (27°C, 14 h/10 h light/dark period, 80% relative humidity). The average generation turnover was 45 days and selection was carried out over 18 generations. Larvae were bred in the insectarium till they reached early 4th instar. LiTOX_S larvae were allowed to develop for 5 days more than other larvae. For each of the three phenotype (Bora-Bora, LiTOX_N and LiTOX_S) total mRNA was extracted from three pools of 60 4th-stage larvae (reared in different batches and originating from different egg-laying female pools) at the same time of the day, using the RNAqueous-4PCR kit (Applied Biosystems/Ambion, Austin, TX, USA), according to the manufacturer’s instructions. RNA quantity and quality was checked using bioanalyzer, and RNA from the 3 biological replicates was pooled in equal proportion to obtain one representative total RNA sample per phenotype (each made from 180 individuals). For each of the three phenotypes, two distinct cDNA libraries were prepared from the same pool of total mRNA and sent at the National Genotyping Center (Genoscope, France) and sequenced individually on the Illumina Genome Analyser II system (GAII) to assess technical variations.
Short read mapping and assembly on the reference genome
The Tophat algorithm (release 2.0.9, http://ccb.jhu.edu/software/tophat) was applied with defaults parameters to align all the reads onto the Aedes aegypti reference genome AaegL2.1 (http://ccb.jhu.edu/software/tophat) by taking into account, both already known (--no-novel-juncs) and novel ab initio splice exon-exon junctions .
Estimation of transcripts’ relative abundance and differential expression (DE) analysis
Bam files were then loaded into Genespring NGS Version 12 (Agilent) software. Reads were filtered on base quality (mean base quality >30 and <10 N per reads) and on mapping quality (alignment score > =97 and Mapping quality > = 250 and remove non-primary multiply mapped reads). Ae. aegypti genome was then re-annotated based on reads distribution and RPKM calculated for each known or new putative transcript using default parameters (min exon length percentile 10, min intro length percentile 10, max intron length percentile 90, min exon RPKM percentile 50, min gene RPKM percentile 50, min gene length percentile 10, min exon RPKM with respect to host gene percentage 75, minimum number of reads in exon 10). For each strain, RPKM correlation between technical replicates was calculated. As r2 were >0.9 for all strains (acceptable technical variation, see Additional file 7: Figure S1), technical replicates were pooled for further analyses, in order to reach a high coverage per transcript and per SNP position. Transcripts with >0.5 RPKM in at least one strain were retained for differential transcription analysis. Differential transcription level was tested for each transcript between each resistant phenotype (LiTOX_N and LiTOX_S) and the parental strain Bora-Bora using an Audic Clavery test designed to compare two cDNA libraries . This Bayesian method is based on the assumption that under the null hypothesis, read counts of the same gene in two libraries come from the same but unknown Poisson distribution. The posterior probability is obtained by Bayesian averaging (infinite mixture) of all possible Poisson distributions with mixing proportions equal to the posteriors under the flat prior. This probability is further adjusted for multiple testing (Benjamini and Hochberg correction). Despite this popular method has been validated by the statistician community , we recommend sequencing biological replicates rather than pooling them prior sequencing in future RNAseq analyses, as sequencing costs are now acceptable and the later approach allows a more robust identification of genes affected by variations across conditions. Transcripts showing an adjusted P- value <10-15 and a fold change ≥2 fold in either direction were considered differentially transcribed.
GO term analysis
The annotation terms from the GO ontology were retrieved from the VectorBase-UniprotKB links. For each resistant phenotype and each differential expression state (i.e. up or down regulated as compared to the Bora-Bora strain), GO terms significantly over-represented were determined using a corrected P-value threshold of 0.05.
SNP detection and analysis
Detection of polymorphisms was performed based on the 129.32 million reads passing quality filters with the following parameters: confidence score threshold =100, coverage >20 reads, base quality cut off =5, ignore locations within or next to homopolymer stretches >10 nucleotides. Among all detected polymorphism variations, only SNP substitutions were considered for differential polymorphism analyses. SNP allele frequencies were then computed between each LiTOX phenotype and the susceptible strain. In a previous study on the LiTOX strain, we have validated the reliability of mRNA sequencing pooled data for inferring population allelic frequencies on 269 SNP: the correlation between allelic frequencies obtained from pooled mRNA sequencing (180 pooled larvae) and individual genotypes (N = 28 individuals) obtained using a DNA Illumina GoldenGate array was highly significant (P <0.001, r =0.85) . Allele frequencies were considered as differential between a LiTOX phenotype and the susceptible strain (hereafter named as differential SNPs) if the following conditions were fulfilled: total read coverage at SNP position ≥30 and allelic frequency difference between both strains ≥40% in either direction. Genic effects of SNPs were computed by comparing SNP locations with reference genome annotation, and were defined as 5’UTR, synonymous, non-synonymous, 3’UTR and other (intronic or intergenic- i.e. close but not within gene boundaries).
Phylogenetic analysis and sequence alignment of R&R consensus domain from the differentially transcribed genes encoding chitin-binding proteins
Phylogenetic analysis (neighbor joining) was performed using MEGA 6.06 software on the chitin-binding domains identified using the CDS search software (https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi). Representatives from the two RR sub-groups from other insect species were added to support the tree (Anopheles gambiae: AnGCP2a-RR2, AAC05656; Bombyx mori: Bm-LCP17-RR1, FAA00504; BmLCP22-RR1, NP_001036828; BMEDG84A-RR2, BAA33195; BMWCP1A-RR2, BAB32475; Drosophila melanogaster: Dm-LCP1-RR1, NP_476619; Dm-EDG78-RR1, NP_524198; Dm-Gart-RR1, NP_476673; DMCcp84Aa-RR2, AAD19803; DMEDG84-RR2, P27780; Manduca sexta: MsLCP14-RR1, AAA29319; MsLCP14.6-RR1, Q94984; Tenebrio molitor: TM-LCP-A1A-RR2, P80681; TMACP20-RR2, P26967). Alignments were performed using the ClustalW software from PBIL Expasy tool (http://npsa-pbil.ibcp.fr/cgi-bin/npsa_automat.pl?page=/NPSA/npsa_clustalw.html).
Mosquito rearing was performed according to protocol CEEA-LR-13002 agreed by the French National Committee of ethics for animal experimentation.
RNAseq data are available at EBI/SRA/ArrayExpress (accession number E-MTAB-1635).
This work was funded by the Région Rhône-Alpes (CIBLE program) and the French Agence Nationale de la Recherche (ANR-08-CES-006-01 DIBBECO). We thank T. Gaude and S. Veyrenc for their technical assistance in mosquito rearing.
- Lacey LA: Bacillus thuringiensis serovariety israelensis and Bacillus sphaericus for mosquito control. J Am Mosq Control Assoc. 2007, 23 (2): 133-163.PubMedView ArticleGoogle Scholar
- Wirth MC, Park HW, Walton WE, Federici BA: Cyt1A of Bacillus thuringiensis delays evolution of resistance to Cry11A in the mosquito Culex quinquefasciatus. Appl Environ Microbiol. 2005, 71 (1): 185-189. 10.1128/AEM.71.1.185-189.2005.PubMed CentralPubMedView ArticleGoogle Scholar
- Perez C, Fernandez LE, Sun JG, Folch JL, Gill SS, Soberon M, Bravo A: 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-18308. 10.1073/pnas.0505494102.PubMed CentralPubMedView ArticleGoogle Scholar
- Ferré J, Escriche B, Bel Y, Vanrie J: Biochemistry and genetics of insect resistance to Bacillus-thuringiensis insecticidal crystal proteins. Fems Microbiol Lett. 1995, 132 (1–2): 1-7.View ArticleGoogle Scholar
- Griffitts JS, Haslam SM, Yang TL, Garczynski SF, Mulloy B, Morris H, Cremer PS, Dell A, Adang MJ, Aroian RV: Glycolipids as receptors for Bacillus thuringiensis crystal toxin. Science. 2005, 307 (5711): 922-925. 10.1126/science.1104444.PubMedView ArticleGoogle Scholar
- Soberon M, Gill SS, Bravo A: Signaling versus punching hole: How do Bacillus thuringiensis toxins kill insect midgut cells?. Cell Mol Life Sci. 2009, 66 (8): 1337-1349. 10.1007/s00018-008-8330-9.PubMedView ArticleGoogle Scholar
- Vachon V, Laprade R, Schwartz JL: Current models of the mode of action of Bacillus thuringiensis insecticidal crystal proteins: a critical review. J Invertebr Pathol. 2012, 111 (1): 1-12. 10.1016/j.jip.2012.05.001.PubMedView ArticleGoogle Scholar
- Pigott CR, Ellar DJ: Role of receptors in Bacillus thuringiensis crystal toxin activity. Microbiol Mol Biol Rev. 2007, 71 (2): 255-281. 10.1128/MMBR.00034-06.PubMed CentralPubMedView ArticleGoogle Scholar
- 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-757. 10.1111/j.1462-2920.2009.02117.x.PubMed CentralPubMedView ArticleGoogle Scholar
- Tetreau G, Stalinski R, Kersusan D, Veyrenc S, David JP, Reynaud S, Despres L: Decreased Toxicity of Bacillus thuringiensis subsp. israelensis to Mosquito Larvae after Contact with Leaf Litter. Appl Environ Microbiol. 2012, 78 (15): 5189-5195. 10.1128/AEM.00903-12.PubMed CentralPubMedView ArticleGoogle Scholar
- Tetreau G, Alessi M, Veyrenc S, Périgon S, David JP, Reynaud S, Despres L: Fate of Bacillus thuringiensis subsp. israelensis in the field: evidence for spore recycling and differential persistence of toxins in leaf litter. Appl Environ Microbiol. 2012, 78 (23): 8362-8367. 10.1128/AEM.02088-12.PubMed CentralPubMedView ArticleGoogle Scholar
- Paris M, Tetreau G, Laurent F, Lelu M, Despres L, David J-P: Persistence of Bacillus thuringiensis israelensis (Bti) in the environment induces resistance to multiple Bti toxins in mosquitoes. Pest Manag Sci. 2011, 67: 122-128. 10.1002/ps.2046.PubMedView ArticleGoogle Scholar
- Paris M, David JP, Despres L: Fitness costs of resistance to Bti toxins in the dengue vector Aedes aegypti. Ecotoxicology. 2011, 20: 1184-1194. 10.1007/s10646-011-0663-8.PubMedView ArticleGoogle Scholar
- Likitvivatanavong S, Chen JW, 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-2838. 10.1021/jf1036189.PubMed CentralPubMedView ArticleGoogle Scholar
- Chen JW, Aimanova KG, Fernandez LE, Bravo A, Soberon M, Gill SS: Aedes aegypti cadherin serves as a putative receptor of the Cry11Aa toxin from Bacillus thuringiensis subsp israelensis. Biochem J. 2009, 424: 191-200. 10.1042/BJ20090730.PubMed CentralPubMedView ArticleGoogle Scholar
- Rebers JE, Riddiford LM: Structure and expression of a Manduca sexta larval cuticle gene homologous to Drosophila cuticle genes. J Mol Biol. 1988, 203 (2): 411-423. 10.1016/0022-2836(88)90009-5.PubMedView ArticleGoogle Scholar
- Andersen SO: Amino acid sequence studies on endocuticular proteins from the desert locust, Schistocerca gregaria. Insect Biochem Mol Biol. 1998, 28 (5–6): 421-434.PubMedView ArticleGoogle Scholar
- Willis JH, Iconomidou V, Smith RF, Hamodrakas SJ: Cuticular proteins. Comprehensive Molecular Insect Science, Volume 4. Edited by: Gilbert L, Iatrou K, Gill S. 2005, Oxford: Elsevier, 79-110.View ArticleGoogle Scholar
- Jurat-Fuentes JL, Karumbaiah L, Jakka SRK, Ning CM, Liu CX, Wu KM, Jackson J, Gould F, Blanco C, Portilla M, Perera O, Adang M: Reduced Levels of Membrane-Bound Alkaline Phosphatase Are Common to Lepidopteran Strains Resistant to Cry Toxins from Bacillus thuringiensis. PLoS One. 2011, 6 (3): e17606-10.1371/journal.pone.0017606.PubMed CentralPubMedView ArticleGoogle Scholar
- Karumbaiah L, Oppert B, Jurat-Fuentes JL, Adang MJ: Analysis of midgut proteinases from Bacillus thuringiensis-susceptible and -resistant Heliothis virescens (Lepidoptera : Noctuidae). Comp Biochem Physiol B-Biochem Mol Biol. 2007, 146 (1): 139-146. 10.1016/j.cbpb.2006.10.104.PubMedView ArticleGoogle Scholar
- Rajagopal R, Sivakumar S, Agrawal N, Malhotra P, Bhatnagar RK: Silencing of midgut aminopeptidase N of Spodoptera litura by double-stranded RNA establishes its role as Bacillus thuringiensis toxin receptor. J Biol Chem. 2002, 277 (49): 46849-46851. 10.1074/jbc.C200523200.PubMedView ArticleGoogle Scholar
- Tetreau G, Bayyareddy K, Jones CM, Stalinski R, Riaz MA, Paris M, David JP, Adang MJ, Despres L: Larval Midgut Modifications Associated with Bti Resistance in the Yellow Fever Mosquito using Proteomic and Transcriptomic Approaches. BMC Genomics. 2012, 13: 248-10.1186/1471-2164-13-248.PubMed CentralPubMedView ArticleGoogle Scholar
- Tetreau G, Stalinski R, David JP, Despres L: Increase in larval gut proteolytic activities and Bti resistance in the dengue fever mosquito. Arch Insect Biochem Physiol. 2013, 82 (2): 71-83. 10.1002/arch.21076.PubMedView ArticleGoogle Scholar
- Rodriguez-Cabrera L, Trujillo-Bacallao D, Borras-Hidalgo O, Wright DJ, Ayra-Pardo C: RNAi-mediated knockdown of a Spodoptera frugiperda trypsin-like serine-protease gene reduces susceptibility to a Bacillus thuringiensis Cry1Ca1 protoxin. Environ Microbiol. 2010, 12 (11): 2894-2903. 10.1111/j.1462-2920.2010.02259.x.PubMedView ArticleGoogle Scholar
- Bayyareddy K, Zhu X, Orlando R, Adang MJ: Proteome Analysis of Cry4Ba Toxin-interacting Aedes aegypti Lipid Rafts using geLC-MS/MS. J Proteome Res. 2012, 11 (12): 5843-5855.PubMed CentralPubMedGoogle Scholar
- Rietveld A, Neutz S, Simons K, Eaton S: Association of sterol- and glycosylphosphatidylinositol-linked proteins with Drosophila raft lipid microdomains. J Biol Chem. 1999, 274 (17): 12049-12054. 10.1074/jbc.274.17.12049.PubMedView ArticleGoogle Scholar
- Avisar D, Segal M, Sneh B, Zilberstein A: Cell-cycle-dependent resistance to Bacillus thuringiensis Cry1C toxin in Sf9 cells. J Cell Sci. 2005, 118 (14): 3163-3171. 10.1242/jcs.02440.PubMedView ArticleGoogle Scholar
- Pust S, Dyve AB, Torgersen ML, van Deurs B, Sandvig K: Interplay between Toxin Transport and Flotillin Localization. Plos One. 2010, 5 (1): e8844-10.1371/journal.pone.0008844.PubMed CentralPubMedView ArticleGoogle Scholar
- Iconomidou VA, Willis JH, Hamodrakas SJ: Unique features of the structural model of ‘hard’ cuticle proteins: implications for chitin-protein interactions and cross-linking in cuticle. Insect Biochem Mol Biol. 2005, 35 (6): 553-560. 10.1016/j.ibmb.2005.01.017.PubMedView ArticleGoogle Scholar
- Jasrapuria S, Arakane Y, Osman G, Kramer KJ, Beeman RW, Muthukrishnan S: Genes encoding proteins with peritrophin A-type chitin-binding domains in Tribolium castaneum are grouped into three distinct families based on phylogeny, expression and function. Insect Biochem Mol Biol. 2010, 40 (3): 214-227. 10.1016/j.ibmb.2010.01.011.PubMedView ArticleGoogle Scholar
- Christensen BM, Li JY, Chen CC, Nappi AJ: Melanization immune responses in mosquito vectors. Trends Parasitol. 2005, 21 (4): 192-199. 10.1016/j.pt.2005.02.007.PubMedView ArticleGoogle Scholar
- Hillyer JF, Christensen BM: Mosquito phenoloxidase and defensin colocalize in melanization innate immune responses. J Histochem Cytochem. 2005, 53 (6): 689-698. 10.1369/jhc.4A6564.2005.PubMedView ArticleGoogle Scholar
- Koella JC, Boete C: A genetic correlation between age at pupation and melanization immune response of the yellow fever mosquito Aedes aegypti. Evolution. 2002, 56 (5): 1074-1079. 10.1111/j.0014-3820.2002.tb01419.x.PubMedView ArticleGoogle Scholar
- Xie RY, Zhuang MB, Ross LS, Gomez I, Oltean DI, Bravo A, Soberon M, Gill SS: Single amino acid mutations in the cadherin receptor from Heliothis virescens affect its toxin binding ability to Cry1A toxins. J Biol Chem. 2005, 280 (9): 8416-8425. 10.1074/jbc.M408403200.PubMedView ArticleGoogle Scholar
- Siqueira HAA, Gonzalez-Cabrera J, Ferré J, Flannagan R, Siegfried BD: Analyses of Cry1Ab binding in resistant and susceptible strains of the European corn borer, Ostrinia nubilalis (Hubner) (Lepidoptera: Crambidae). Appl Environ Microbiol. 2006, 72 (8): 5318-5324. 10.1128/AEM.00219-06.PubMed CentralPubMedView ArticleGoogle Scholar
- Sayed A, Nekl ER, Siqueira HAA, Wang HC, Ffrench-Constant RH, Bagley M, Siegfried BD: A novel cadherin-like gene from western corn rootworm, Diabrotica virgifera virgifera (Coleoptera: Chrysomelidae), larval midgut tissue. Insect Mol Biol. 2007, 16: 591-600.PubMedGoogle Scholar
- Peng DH, Xu XH, Ye WX, Yu ZN, Sun M: Helicoverpa armigera cadherin fragment enhances Cry1Ac insecticidal activity by facilitating toxin-oligomer formation. Appl Microbiol Biotechnol. 2010, 85 (4): 1033-1040. 10.1007/s00253-009-2142-1.PubMedView ArticleGoogle Scholar
- Likitvivatanavong S, Chen JW, Bravo A, Soberon M, Gill SS: Cadherin, Alkaline Phosphatase, and Aminopeptidase N as Receptors of Cry11Ba Toxin from Bacillus thuringiensis subsp. jegathesan in Aedes aegypti. Appl Environ Microbiol. 2011, 77 (1): 24-31. 10.1128/AEM.01852-10.PubMed CentralPubMedView ArticleGoogle Scholar
- Hua G, Zhang R, Abdullah MAF, Adang MJ: Anopheles gambiae cadherin AgCad1 binds the Cry4Ba toxin of Bacillus thuringiensis israelensis and a fragment of AgCad1 synergizes toxicity. Biochemistry. 2008, 47 (18): 5101-5110. 10.1021/bi7023578.PubMed CentralPubMedView ArticleGoogle Scholar
- Saengwiman S, Aroonkesorn A, Dedvisitsakul P, Sakdee S, Leetachewa S, Angsuthanasombat C, Pootanakit K: In vivo identification of Bacillus thuringiensis Cry4Ba toxin receptors by RNA interference knockdown of glycosylphosphatidylinositol-linked aminopeptidase N transcripts in Aedes aegypti larvae. Biochem Biophys Res Commun. 2011, 407 (4): 708-713. 10.1016/j.bbrc.2011.03.085.PubMedView ArticleGoogle Scholar
- Flannagan RD, Yu CG, Mathis JP, Meyer TE, Shi XM, Siqueira HAA, Siegfried BD: Identification, cloning and expression of a Cry1Ab cadherin receptor from European corn borer, Ostrinia nubilalis (Hubner) (Lepidoptera: Crambidae). Insect Biochem Mol Biol. 2005, 35 (1): 33-40. 10.1016/j.ibmb.2004.10.001.PubMedView ArticleGoogle Scholar
- Fabrick JA, Tabashnik BE: Binding of Bacillus thuringiensis toxin Cry1Ac to multiple sites of cadherin in pink bollworm. Insect Biochem Mol Biol. 2007, 37 (2): 97-106. 10.1016/j.ibmb.2006.10.010.PubMedView ArticleGoogle Scholar
- Fabrick J, Oppert C, Lorenzen MD, Morris K, Oppert B, Jurat-Fuentes JL: A Novel Tenebrio molitor Cadherin Is a Functional Receptor for Bacillus thuringiensis Cry3Aa Toxin. J Biol Chem. 2009, 284 (27): 18401-18410. 10.1074/jbc.M109.001651.PubMed CentralPubMedView ArticleGoogle Scholar
- Bravo A, Gomez I, Conde J, Munoz-Garay C, Sanchez J, Miranda R, Zhuang M, Gill SS, Soberon M: Oligomerization triggers binding of a Bacillus thuringiensis Cry1Ab pore-forming toxin to aminopeptidase N receptor leading to insertion into membrane microdomains. Biochimica Et Biophysica Acta-Biomembranes. 2004, 1667 (1): 38-46. 10.1016/j.bbamem.2004.08.013.View ArticleGoogle Scholar
- Morin S, Biggs RW, Sisterson MS, Shriver L, Ellers-Kirk C, Higginson D, Holley D, Gahan LJ, Heckel DG, Carriere Y, Dennehy TJ, Brown JK, Tabashnik BE: Three cadherin alleles associated with resistance to Bacillus thuringiensis in pink bollworm. Proc Natl Acad Sci U S A. 2003, 100 (9): 5004-5009. 10.1073/pnas.0831036100.PubMed CentralPubMedView ArticleGoogle Scholar
- Chen JW, Aimanova KG, Pan SQ, Gill SS: Identification and characterization of Aedes aegypti aminopeptidase N as a putative receptor of Bacillus thuringiensis Cry11A toxin. Insect Biochem Mol Biol. 2009, 39 (10): 688-696. 10.1016/j.ibmb.2009.08.003.PubMed CentralPubMedView ArticleGoogle Scholar
- Singh A, Sivaprasad C: Functional interpretation of APN receptor from M.sexta using a molecular model. Bioinformation. 2009, 3 (8): 321-325. 10.6026/97320630003321.PubMed CentralPubMedView ArticleGoogle Scholar
- Trapnell C, Pachter L, Salzberg SL: TopHat: discovering splice junctions with RNA-Seq. Bioinformatics. 2009, 25 (9): 1105-1111. 10.1093/bioinformatics/btp120.PubMed CentralPubMedView ArticleGoogle Scholar
- Audic S, Claverie JM: The significance of digital gene expression profiles. Genome Res. 1997, 7: 986-995.PubMedGoogle Scholar
- Tiňo P: Basic properties and information theory of Audic-Claverie statistic for analyzing cDNA arrays. BMC Bioinformatics. 2009, 10: 310-10.1186/1471-2105-10-310.PubMed CentralPubMedView ArticleGoogle Scholar
- Paris M, Marcombe S, Coissac E, Corbel V, David JP, Despres L: Investigating the genetics of Bti resistance using mRNA tag sequencing: application on laboratory strains and natural populations of the dengue vector Aedes aegypti. Evol Appl. 2013, 6 (7): 1012-1027.PubMed CentralPubMedGoogle Scholar
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. 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.