Skip to main content
Download PDF

A differential transcriptional profile by Culex quinquefasciatus larvae resistant to Lysinibacillus sphaericus IAB59 highlights genes and pathways associated with the resistance phenotype

Parasites & Vectors volume 12, Article number: 407 (2019) Cite this article

Abstract

Background

The study of the mechanisms by which larvae of the Culex quinquefasciatus mosquito survive exposure to the entomopathogen Lysinibacillus sphaericus has benefited substantially from the generation of laboratory-selected colonies resistant to this bacterium. One such colony, RIAB59, was selected after regular long-term exposure of larvae to the L. sphaericus IAB59 strain. This strain is characterized by its ability to produce the well known Binary (Bin) toxin, and the recently characterized Cry48Aa/Cry49Aa toxin, able to kill Bin-resistant larvae. Resistance to Bin is associated with the depletion of its receptor, Cqm1 α-glucosidase, from the larvae midgut. This study aimed to identify novel molecules and pathways associated with survival of the RIAB59 larvae and the resistance phenotype.

Methods

A transcriptomic approach and bioinformatic tools were used to compare the profiles derived from the midguts of larvae resistant and susceptible to L. sphaericus IAB59.

Results

The RNA-seq profiles identified 1355 differentially expressed genes (DEGs), with 673 down- and 682 upregulated transcripts. One of the most downregulated DEGs was cqm1, which validates the approach. Other strongly downregulated mRNAs encode the enzyme pantetheinase, apolipoprotein D, lipases, heat-shock proteins and a number of lesser known and hypothetical polypeptides. Among the upregulated DEGs, the top most encodes a peroxisomal enzyme involved in lipid metabolism, while others encode enzymes associated with juvenile hormone synthesis, ion channels, DNA binding proteins and defense polypeptides. Further analyses confirmed a strong downregulation of several enzymes involved in lipid catabolism while the assignment of DEGs into metabolic pathways highlighted the upregulation of those related to DNA synthesis and maintenance, confirmed by their clustering into related protein networks. Several other pathways were also identified with mixed profiles of down- and upregulated transcripts. Quantitative RT-PCR confirmed the changes in levels seen for selected mRNAs.

Conclusions

Our transcriptome-wide dataset revealed that the RIAB59 colony, found to be substantially more resistant to Bin than to the Cry48Aa/Cry49Aa toxin, developed a differential expression profile as well as metabolic features co-selected during the long-term adaptation to IAB59 and that are most likely linked to Bin resistance.

Background

Culex quinquefasciatus is one of the most abundant mosquito species associated with urban areas, particularly those which are characterized by precarious sanitation. This antropophilic species is a known vector of human pathogens and it has been the target of control programmes worldwide. Biolarvicides based on the entomopathogenic bacterium Lysinibacillus sphaericus are some of the most effective agents used to control Cx. quinquefasciatus populations [1]. Binary (Bin) toxin, its major insecticidal factor, is highly active against larvae from this mosquito species and displays long persistence in its breeding sites [2]. However, as most insecticides that have been used in the field, its effectiveness faces the emergence of resistant mosquito populations [3]. Therefore, the major challenge for L. sphaericus utilization is to overcome the resistance that can be selected against the Bin toxin, the active ingredient of the commercially available biolarvicides.

Bin is synthesized as a heterodimeric protoxin and its two subunits, BinA (~43 kDa) and BinB (~51 kDa), act in synergy to reach optimal toxicity in equimolar amounts [4]. Its larvicidal activity against Cx. quinquefasciatus larvae requires protoxin ingestion, midgut processing of protoxin into toxin by serine-proteases [5] and specific binding to midgut-bound receptors [6]. These receptors were identified as Cqm1 glycosylphosphatidylinositol (GPI) anchored α-glucosidases [7, 8], maltases that may play roles in the metabolism of carbohydrates [9, 10]. Binding of the Bin toxin to the Cqm1 receptors allows it to form pores and penetrate the membranes of midgut-cells [11,12,13]. Larval intoxication by Bin induces cell autophagy and cytoplasmatic vacuolization, one of the most prominent cytopathological alterations recorded in the midgut epithelium of treated individuals [14,15,16]. Mitochondrial damages have also been recorded as a marked effect of Bin toxin and these have been recently demonstrated to be involved with the activation of the intrinsic apoptosis pathway [17].

Resistance to the Bin toxin has been recorded in several occasions [18,19,20,21,22,23] and the major mechanism responsible for that is the failure of this toxin to bind to the midgut-cells, mostly due to the lack of Cqm1 receptors [19, 20, 24]. Further studies defined that mutations in the cqm1 gene can disrupt the expression of full-length GPI anchored proteins and cause high resistance ratios (> 5000-fold) [7, 25,26,27,28,29]. Laboratory-selected colonies highly resistant to the Bin toxin have been used as models to investigate resistance and to screen for other bacterial strains and toxins that could be active against these larvae and overcome Bin resistance [21, 23]. In this context, the L. sphaericus IAB59 strain was thus found to be toxic to Bin-resistant larvae [21, 30] and able to produce the Bin toxin as well as Cry48Aa/Cry49Aa, another binary toxin responsible for the larvicidal effect against those larvae [31, 32]. Previous studies have indicated that the Cry48Aa/Cry49Aa toxin was able to interact with receptors from the midgut of susceptible [33] and also with Bin resistant larvae lacking the Cqm1 receptor [34]. Specific receptors involved in this interaction with Cry48Aa/Cry49Aa are being studied and potential candidates have been recently tentatively identified [35]. The L. sphaericus IAB59 strain was also employed to select a laboratory resistant Cx. quinquefasciatus colony (RIAB59) and, after 72 generations, a high resistance level was achieved to the IAB59 crystals-spores used to treat the larvae throughout the selection process [21, 36]. During the selection procedure, a high resistance specific to the Bin toxin arose much earlier [21, 37] than the resistance to the whole set of IAB59 toxins [36]. In fact, the resistance level of the RIAB59 larvae to Cry48Aa/Cry49Aa only was much lower than that recorded for the Bin toxin [34].

The high level of resistance of the Cx. quinquefasciatus RIAB59 larvae to the Bin toxin [36] was previously characterized and shown to be associated with the cqm1REC allele. A frame shift change within this allele leads to the loss of the GPI-anchored Bin receptors in the midgut epithelium [36, 38]. Individuals from the RIAB59 colony are homozygous for the cqm1REC and the resistance is recessively inherited [36]. The disruption of Cqm1 expression also represents the lack of a maltase that can presumably lead to an impact on the digestion of carbohydrates by the insect larvae. Nevertheless this colony has been successfully maintained in the laboratory under intermittent selection with L. sphaericus IAB59 strain for more than 200 generations. Assessment of biological parameters as fecundity, fertility and pupae weight in RIAB59 individuals did not reveal significant differences in comparison to a susceptible reference colony [36, 38]. These results suggest that those individuals developed mechanisms to balance the potential costs associated to the resistance phenotype. In this context we employed here a differential expression analysis at the transcriptome scale comparing the RIAB59 larvae with individuals from a susceptible colony.

Methods

Culex quinquefasciatus colonies

CqSLab is a laboratory reference colony susceptible to insecticidal compounds while RIAB59 is a resistant laboratory colony selected after continuous exposure to the L. sphaericus IAB59 strain as described in previous studies [36, 38]. The IAB59 strain produces Binary 1 toxin, named Bin1 [39], as well as Cry48Aa/Cry49Aa toxins [31]. Bin1 has a similar larvicidal activity and capacity to bind to midgut receptors as Binary 2, produced by other L. sphaericus strains, such as 2362, 1593 and C3–41 [37]. The strong ratio of resistance (RR) to the Bin toxin (RR > 5000) achieved by the RIAB59 colony is known to be due to the individuals being homozygous to the cqm1REC allele [38], therefore not being able to express the midgut bound Cqm1 receptor [7]. Since the resistance to L. sphaericus IAB59 has been achieved [36] this colony has been exposed to the bacterium every five generations, to periodically verify the status of in vivo resistance. In this study larvae from generation F194 were used, with the last exposure to L. sphaericus IAB59 having taken place in generation F189. The larvae samples collected for this analysis from the two colonies were thus not directly exposed to the L. sphaericus IAB59. Both susceptible and resistant colonies have been maintained in the insectarium of the Instituto Aggeu Magalhães (IAM)-FIOCRUZ for more than five years under controlled conditions of 26 ± 1 °C, 70 % relative humidity and a 14h:10h (light:dark) photoperiod. Larvae were reared in dechlorinated water and fed with cat food (Friskies®, Nestlé Purina,  Ribeirão Preto, Brazil). Adults were maintained on a 10% sugar solution and females were also artificially fed with defibrinated rabbit blood.

RNA extraction, mRNA library construction and sequencing

Total RNA from three pools of 20 larvae midguts, from each colony, was extracted using an RNeasy Mini Kit (Qiagen, Hilden, Germany), following the manufacturer’s instructions. RNA purity and concentration were assessed using a NanoDrop 2000™ spectrophotometer (Thermo Fisher Scientific, Waltham, USA) and Qubit 2.0 Fluorometer (Thermo Fisher Scientific). The RNA integrity was evaluated by agarose gel electrophoresis. Paired-end libraries were prepared from total RNA using the TruSeq Stranded mRNA Library Prep kit (Illumina, San Diego, USA), following standard procedures, and sequenced using a MiSeq™ Reagent Kit V3 (Illumina, San Diego, USA; 150 cycles) on an Illumina MiSeq Sequencer of the IAM-FIOCRUZ.

RNA-seq data analysis

The quality of the sequenced reads was checked applying the FastQC tool (www.bioinformatics.babraham.ac.uk/projects/fastqc/) and we observed that, on average, all bases had Phred scores higher than 30, hence removal of low quality reads was not needed. Each library was then mapped against the genome assembly of the Cx. quinquefasciatus Johannesburg strain, CpipJ2 (file: culex-quinquefasciatus-johannesburgscaffoldscpipj2.fa). This was downloaded from the VectorBase database (https://www.vectorbase.org/), applying STAR aligner v.2.5.3 [40] with default parameters, except for the quantMode option which generates a GeneCounts file considering stranded libraries. The data on gene counts were next loaded to the R environment and all replicates organized in a matrix, which was transformed to log2 scale using the function rlogTransformation from the DESeq2 package [41]. This matrix was used as input to the prcomp function in order to perform the principal components analysis (PCA). The function plot was then applied to plot the first and second principal components calculated from each sample. The R package DESeq2 was used to perform the differential expression analysis, considering only genes represented by at least five reads for all three biological replicates, in at least one condition (resistant or susceptible). Genes with absolute values of log2 fold changes equal or greater than 1 and with FDR corrected P-values lower than 0.05 were selected for further analyses. It is worth noting that the real fold change value can be obtained by calculating the log2 fold change value [e.g. log2 fold change(sensitive/resistant) = 4 can be converted as 2|4| = 16; this means a gene expression rate 16 times higher/lower]. The normalized read counts for each gene analyzed in all samples was next used as input to the function heatmap.2 from the gplots R package in order to creat a heatmap visualization from the expression data. In addition, the same data was also used to create a MA plot using the function plot, where the differentially expressed genes were highlighted. In order to further explore and refine the annotation of several DEGs tagged as hypothetical proteins in the Cx. quinquefasciatus genome, three additional annotation steps using Blastx alignments with default parameters were performed. mRNA nucleotide sequences from all DEGs were used as queries against the most recent versions of the following databases: (i) UniProtKB curated database; (ii) non-redundant database from NCBI retaining only the top hit for the annotation; and (iii) non-redundant database from NCBI using Blast2GO [42] retaining the top five hits for annotation. The identity of some genes was also obtained by manual annotation (Additional file 1: Table S1). All databases were downloaded up to July 12th 2018.

The STRINGdb R package [43] was used to assign Kyoto Encyclopedia of Genes and Genomes (KEGG) [44] pathway terms to the Cx. quinquefasciatus differentially expressed genes (DEGs) (http://www.vectorbase.org). The Pathview package [45] was then applied to map the DEGs and their fold change measures to the enriched KEGG pathways. In addition, KEGG terms were employed to perform an enrichment functional analysis of those genes using the function get_enrichment from the STRINGdb R package. Here, it is worth to mention the background frequency of KEGG pathways for the Cx. quinquefasciatus genome comes from the STRING database (https://string-db.org/). Furthermore, the DEGs were also visualized into protein networks present at STRING database. The protein networks built from the lists of DEGs were analyzed using the software Cytoscape v.3.5.1 (http://www.cytoscape.org) along with the plug-in AutoAnnotate v.1.2 (http://www.baderlab.org/Software/AutoAnotate). Here, the networks were clustered using the algorithm Community Cluster (GLay) since it does not demand a weighted network. The cluster annotations were assigned based on the three most frequent words present in the Vectorbase annotation of the protein members from each cluster.

qRT-PCR validation of RNA-seq data

The total RNA used in this assay was the same as described above (RNA extraction, mRNA library construction and sequencing). The reactions were performed using a QuantiTect® SYBR Green RT-PCR® Kit one step (Qiagen), following the manufacturer’s instructions and using specific primers for each gene selected (Additional file 2: Table S2). 18S was used as an endogenous control gene [46]. The samples were analyzed in a QuantStudio® 5 System (Thermo Fisher Scientific) and the relative quantification was performed using Applied Biosystems™ Analysis software, Relative Quantification Analysis Module v.3.3 [47]. Means and standard errors from three biological replicates from each colony were compared using Studentʼs t-test in GraphPad Prism v.5.0.0 for Windows (GraphPad Software Inc., San Diego, USA), considering a P-value < 0.05 statistically significant.

Results

Overview of the RIAB59 resistant colony and the sequencing results

The Cx. quinquefasciatus RIAB59 colony selected for this study has been previously shown to display a strong and stable resistance ratio (RR) to the Bin toxin (RR > 5000-fold) [36, 38]. The specific resistance of this colony to the Cry48Aa/Cry49Aa toxin was also previously reported [34], but it was first re-evaluated here and quantified prior to the sequencing effort. The diagnostic bioassays showed that the RR between the LC50 of Cry48Aa/Cry49Aa for the RIAB59 larvae and that for the susceptible individuals was approximately 15-fold (data not shown). Early fourth-instar larvae from the resistant RIAB59 colony (without recent exposure to the IAB59 bacterium), as well as from a susceptible Cx. quinquefasciatus colony from which it was derived (CqSLab), were then harvested and had their midguts dissected and subjected to RNA extraction followed by whole transcriptome shotgun sequencing (RNA-seq). A total of 19,288,276 reads were sequenced and mapped against the Cx. quinquefasciatus genome. This analysis yielded 7045 genes (data not shown) that fulfilled threshold parameters of being represented by at least five reads from each of three biological replicates for either the resistant or the susceptible colonies. These correspond to roughly 36% of the total gene repertoire of the annotated Cx. quinquefasciatus genome. Information on the sequencing data and libraries is presented in Table 1.

Table 1 Statistics for the sequenced RNA-seq libraries of Culex quinquefasciatus larvae from a Lysinibacillus sphaericus RIAB59 resistant colony and a susceptible one
Full size table

To compare the overall profile of gene expression by the larvae from the two colonies, a principal components analysis (PCA) was first carried out with the sets of annotated genes from them. Two distinct profiles of gene expression consistently reproduced with the biological replicates used were observed (Fig. 1a). Next, a heatmap was also generated comparing their overall gene expression and distinct profiles were confirmed between the two colonies which were reproduced by the biological replicates (Fig. 1b). These results are consistent with the resistance to IAB59 being associated with differential expression of a significant number of genes. It is important to notice, however, that the differential expression revealed by this study was not an induced response derived from larvae challenged with L. sphaericus IAB59 but a constitutive profile which is likely to result from an adaptation to the selection process.

Fig. 1
figure 1

Gene expression profile generated by RNA-seq comparing the Culex quinquefasciatus larvae from a Lysinibacillus sphaericus resistant (R) and susceptible (S) colonies. a Principal components analysis (PCA) of the gene expression profile. b Heatmap for upregulated (green) and downregulated (red) genes in biological replicates from resistant (R2, R3, R4) and susceptible larvae (S5, S7, S8). c MA plot of differentially expressed genes. Upregulated and downregulated genes are marked in green and red, respectively

Full size image

Landscape analysis of differentially expressed transcripts in the RIAB59 colony

To individually assess the involvement of differentially expressed genes (DEGs) in the Cx. quinquefasciatus resistance to the IAB59 toxins (mainly Bin), a compilation of the DEGs found comparing the resistant and susceptible colonies was then carried out. The overall distribution of the genes found after sequencing and annotation according to the log2 fold changes in number of reads observed between the resistant and susceptible colonies is represented in Fig. 1c. A total of 1355 DEGs with absolute values of log2 fold change ≥ 1 and P-value ≤ 0.05 were identified (Additional file 1: Table S1), with 673 of those genes downregulated in the RIAB59 colony while 682 were upregulated. Most of the DEGs identified varied within a range of two- to four-fold in the number of reads (log2 fold changes between 1 and 2), although several were found with much greater changes in expression (log2 fold changes ≥ 3). After an effort to identify the DEGs using the automatic annotation available at the VectorBase database as well as blast searches against Uniprot and GenBank databases, roughly 90% of the DEGs were identified while only 158 out of 1355 DEGs were found to encode hypothethical proteins. From the set of identified DEGs, one of the greatest scores of downregulation for the RIAB59 colony, with a log2 fold change of 4.65 (25-fold decrease) (P = 9.14E-304), is the cqm1 gene (CPIJ013173; annotated as “Neutral & basic aminoacid transport protein rBAT” in VectorBase), directly responsible for the high resistance phenotype to the Bin toxin. Considering its role as the Binary toxin receptor and the fact that its absence from the larvae midgut has been previously linked to resistance to L. sphaericus, the gene encoding the Cqm1 α-glucosidase is a robust marker for this phenotype. Here, the identification of its mRNA within the most downregulated transcripts from the resistant colony confirms not only the resistance phenotype to Bin but also the validity of the methods applied in this study.

Identification of the top differentially expressed genes (DEGs)

From the 1355 DEGs, those genes showing log2 fold changes greater than 3 and significant P-values (< 0.05) are highlighted in Tables 2 and 3. For clarity and simplification, when multiple genes were found coding for identical or nearly identical proteins, the IDs for all genes are indicated but only the values for the gene with the highest fold change are shown. Among the downregulated DEGs in the RIAB59 colony (Table 2), the gene encoding the Cqm1 α-glucosidase is found as the fourth most downregulated. To assess whether other changes in mRNA expression could help pinpoint specific proteins as potential receptors for the Cry48Aa/Cry49Aa toxin, a search for possible toxin ligands whose expression differed significantly between the susceptible and resistant colonies was first carried out. None of the candidates indicated by a previous study [35] showed differential expression except for two: the vanin 1/pantetheinase precursor, with the highest level of downregulation seen; and the apolipoprotein D transcripts (Table 2). A pantetheinase, an enzyme that acts in the catalysis of D-pantetheine into pantothenate (vitamin B5) [48], has been recently identified as one of several ligands of the Cry48Aa/Cry49Aa toxin [35], while apolipoprotein D has previously been found as a ligand for the Cry8Ea toxin in a coleopteran [49]. However, these proteins seem be rather involved in other metabolic features of the resistant larvae than being receptors, as will be further discussed. Likewise, other downregulated DEGs listed in Table 2 code for proteins that are likely not involved with specific binding by the toxin, but nevertheless may be relevant for the resistance phenotype. An example is the transcript coding for an ankyrin, a transmembrane protein belonging to pathways involved in pathogen recognition, activation of defense responses, macroautophagy and membrane transport [50,51,52,53]. Two chitin binding proteins showed substantial upregulation (CPIJ014193) and downregulation (CPIJ014195) and, although the specific role of these proteins remains unclear, they are known to participate in bacteria-host interactions during pathogen infection [54, 55]. Other strongly downregulated genes include two lipases, a fatty acid hydroxylase, heat-shock proteins, cytochrome P450 and glutathione S-transferase homologues. Some of these and several other genes found to be under strong downregulation code for proteins that are localized to the mitochondria with at least one of those involved in apoptosis, a pathway that might also be associated with yet another gene found within this list, coding for a p53 and DNA damage-regulated polypeptide. Various genes coding for hypothetical polypeptides and proteins less well known are also listed in Table 2, but their roles need further investigation.

Table 2 Top most downregulated genes in the Culex quinquefasciatus larvae from the RIAB59 resistant colony, revealed by the RNA-seq analysis
Full size table
Table 3 Top most upregulated genes in the Culex quinquefasciatus larvae from the RIAB59 resistant colony, revealed by the RNA-seq analysis
Full size table

Among the upregulated mRNAs in the resistant colony (Table 3), the top most is the one encoding 2-hydroxyacyl-CoA lyase 1, an enzyme involved in lipid metabolism within the peroxisome [56, 57]. The identification among the top five upregulated transcripts of two genes encoding the enzymes farnesol dehydrogenase and farnesoic acid 0-methyl transferase, both implicated in the biosynthesis of the juvenile hormone in mosquitoes [58, 59], is another relevant finding in this profile. Upregulated DEGs found in the resistant colony also include genes associated with chromatin and DNA related processes (see CPIJ019303, CPIJ015649 and CPIJ010133). Other mRNAs listed in Table 3 include those encoding ion channels, a hypothetical protein with chitin binding domain and lesser known proteins such as the ficolins, a group of oligomeric lectins with roles in innate immunity in vertebrates [60, 61], and hexamerins, hexameric serum proteins whose levels have been shown to be altered upon infection in insects [62]. Transcripts encoding mitochondrial proteins and related to apoptosis are also found upregulated. Overall the differences between the most down- or upregulated genes are consistent with likely substantial changes in expression and abundance of individual proteins encoded by these genes in the resistant larvae.

Categorization of differentially expressed transcripts

Considering the identification of several genes associated with lipid metabolism among the most downregulated transcripts in the RIAB59 resistant colony, a more detailed search for DEGs whose transcripts encoded related proteins was carried out. In addition to those genes listed in Table 2, several other transcripts encoding proteins associated with lipid catabolism were also seen to be downregulated with values of log2 fold changes between 1.9 (3.7-fold decrease) and 3 (8-fold decrease) (Additional file 1: Table S1). These downregulated transcripts include not only lipases, generally involved with the digestion of triacylglycerol molecules (CPIJ004230, CPIJ004226, CPIJ001036 and so on), but also enzymes which are predicted to be part of fatty acid and lipid degradation pathways, both in the mitochondria as well as in the peroxisome (CPIJ011600, CPIJ003870, CPIJ004138).

Next, to enable the identification of further metabolic pathways that can be associated with resistance to IAB59 and allow a wider understanding of the biological processes involved, the whole set of DEGs was assigned to defined KEGG metabolic pathways using a streamlined approach. Amongst the 1355 DEGs (Additional file 1: Table S1), 325 were then assigned to 123 KEGG pathways of which 36 showed terms specifically enriched for the downregulated DEGs from the resistant colony, of which four are involved in lipid metabolism (Table 4). These include pathways associated with the metabolism of sphingolipid (cqu00600), glycerophospholipid (cqu00564), ketone bodies (cqu00072) and other lipids (cqu00565). Two other pathways are linked with organelles associated with catabolism, phagosome (cqu04145) and peroxisome (cqu04146), with their downregulation being consistent with the previously described data. Several metabolic pathways are also indentified having multiple downregulated DEGs, including pathways involved with the metabolism of carbohydrates, nucleotides and amino acids. Regarding the upregulated genes from the resistant larvae, these were also assigned to 36 KEGG enriched pathways (Table 5). Signal transduction pathways, such as the mTOR (cqu04150) dependent pathway, and cellular processes such as autophagy (cqu04140) and the ubiquitin mediated proteolysis (cqu04120) are examples of biological processes that might be relevant for the resistant larvae due to their differential upregulation. It is also noteworthy that several pathways related to DNA replication and repair were enriched and mostly associated with upregulated transcripts, such as nucleotide excision repair (cqu03420), mismatch repair (cqu03430), homologous recombination (cqu03440) and non-homologous end-joining (cqu03450).

Table 4 KEGG pathways (36) with enriched downregulated terms in Culex quinquefasciatus larvae from the RIAB59 resistant colony
Full size table
Table 5 KEGG pathways (36) with enriched upregulated terms in Culex quinquefasciatus larvae from the RIAB59 resistant colony
Full size table

We also opted to examine selected KEGG pathways in order to have a clearer idea of the impact of the upregulated or downregulated DEGs. For this analysis the DEGs were individually assigned to high quality KEGG pathways annotated for Cx. quinquefasciatus, providing a more detailed analysis of the effects upon each pathway and including DEGs which might have been missed by the first approach. For some pathways, however, genes were detected coding for proteins that are involved either in the activation or inhibition of processes within the same pathway and the final outcomes on the regulation of those pathways are difficult to determine, although the data clearly indicate that these processes are under strong modulation. The mTOR pathway, for instance, showed eleven upregulated and three downregulated genes in the resistant colony (Fig. 2a). The autophagy pathway, in contrast, showed 12 upregulated and one downregulated gene with a more consistent upregulated profile (Fig. 2b). Other examples include pathways associated with DNA replication and repair, where several upregulated transcripts related to these pathways are found (Additional file 1: Table S1), although not included in Table 3 due to smaller values of log2 fold changes. These include transcripts encoding the DNA polymerase delta catalytic subunit (CPIJ018287), the DNA repair protein Rad50 (CPIJ004115) and the replication factor C/RFC (CPIJ017892), all found represented in two or more of the KEGG maps related to DNA metabolism (Additional file 3: Figure S1). Overall the KEGG analysis proved to be a usefull resource in order to produce a functional view of the differential expression profile of the resistant larvae, with analysis of the up- and downregulated genes being more informative. However, it is important to notice that relevant DEGs can be absent in these pathways. One example is the cqm1 gene, of great functional relevance in this study but is missing from the currently available KEGG maps.

Fig. 2
figure 2

KEGG signaling pathways of mTOR (a) and autophagy (b) with differentially expressed genes in Culex quinquefasciatus larvae from a Lysinibacillus sphaericus resistant colony (RIAB59). Upregulated and downregulated genes are represented as green and red rectangles, respectively. Non-colored rectangles denote proteins which were not sequenced, not differentially expressed or do not have homologs annotated in the current Cx. quinquefasciatus genome

Full size image

STRING evaluation of protein-protein interactions among DEGs

In order to analyze our findings at a different functional level, the DEGs were also analyzed with the STRING tool which provides a functional protein-protein interaction view of the data. The DEGs genes formed distinct clusters named according to their representation in the DEGs gene annotation (Additional file 4: Table S3). A set of downregulated genes (239) grouped in eight clusters (Fig. 3a), of which 169 formed a single cluster designed as “Wd-repeat protein spliceosome”, with a marked and concentrated network of molecular interactions. This major cluster is formed by genes encoding heat-shock proteins and ribonucleoproteins and displays interactions with all other clusters including the “cytochrome p450 71b36”, the second most represented one (28 genes), as well as six other smaller clusters (6–8 genes) which include the “catalyse phosphorolytic purine” were the vanin 1/pantetheinase gene is found. Considering a set of upregulated 226 genes in the resistant larvae, the analysis grouped them into ten clusters, some of which are in agreement with the pathways identified by the KEGG analysis (Fig. 3b). The “protease factor DNA” cluster includes the highest gene number of DEGs (69) and several of them are related to DNA metabolism, such as replication, repair, and chromosomes maintenance. Similarly to “DNA protease factor” the “protein f-box wd40” cluster group genes related with DNA metabolism and include some coding for proteins with cell cycle roles. Other clusters may also be related to specific metabolic processes enhanced in the resistant larvae, some of which were also implied by the analyses of KEGG pathways. Noteworthy are several of the genes found as part of the fourth cluster, “part electron fe-s”, which code for protein kinases responsible for cell signaling and protein regulation, implying changes in cell signaling and again the cell cycle as important for the resistance phenotype. Another interesting finding from the STRING data is the identification of clusters of completely distinct cytochrome p450 genes among both downregulated (cytochrome P450 71 b3) and upregulated (P450 4d8 cytochrome) DEGs. These may indicate marked changes in the expression of defense proteins by the resistant larvae.

Fig. 3
figure 3

STRING clusters which display the functional protein-protein interaction of differentially expressed genes of Culex quinquefasciatus larvae from a Lysinibacillus sphaericus resistant colony (RIAB59). a Downregulated genes. b Upregulated genes

Full size image

Validation of differential gene expression

To validate the RNA-seq data and the analyses derived from them, five DEGs were selected to have their relative expression profile evaluated by quantitative real time RT-PCR, comparing the resistant RIAB59 and the susceptible colonies. The selection considered both down- and upregulated DEGs in the resistant larvae, their potential biological relevance and also included genes with different ranges of log2 fold changes. The cqm1 gene (CPIJ013173) was first chosen, since it is the internal marker for the resistant phenotype to the Bin toxin. Indeed, the evaluation of its relative expression levels by the qRT-PCR analysis showed a marked reduction in the resistant larvae (log2 fold change of 4.6 or a 25-fold decrease), as expected and confirming the RNA-seq results (Fig. 4a). The next two targets were chosen from the repertoire of downregulated genes, the vanin 1 pantetheinase (CPIJ017593) and caspase 3 (CPIJ012580). Pantetheinase showed the highest values of downregulation (log2 fold change of 7.1 or a 128-fold decrease) while the caspase 3 gene was selected due to the fact that it has an important role as a marker of apoptosis, despite the lower score of downregulation seen by us (log2 fold change of 2.1 or 4-fold decrease). Their relative transcription quantification performed by qRT-PCR confirmed reductions in levels of both mRNAs (log2 fold change of 2.0 or a 4.1-fold decrease for the pantetheinase and log2 fold change of 0.5 with a 1.4-fold decrease for the caspase transcript), with the reduction in levels of the pantetheinase mRNA being not as strong, for reasons unknown (Fig. 4b). For the last two genes, in addition to the scores of log2 fold changes, another relevant criteria for selection was their upregulation and participation in several KEGG pathways. The chosen Rac serine/threonine kinase gene (CPIJ007754) participates in five pathways including FoxO, mTOR and autophagy, while the DNA polymerase-delta gene (CPIJ018287) is a component of eight pathways related to DNA recombination/repair. For these genes the values of log2 fold changes according to the RNA-seq analysis for the resistant larvae (from Additional File 1: Table S1) were 1.35 (2.5-fold increase) and 2.4 (5.3-fold increase), respectively, with equivalent small but significant (P < 0.05) increases in expression, confirmed by the qRT-PCR (Fig. 4c). Overall these analyses then confirm that the differential expression patterns found for these selected genes with the qRT-PCR are in agreement with the respective profile revealed by RNA-seq.

Fig. 4
figure 4

Relative expression levels of Culex quinquefasciatus genes from a Lysinibacillus sphaericus resistant larvae (R) compared to a susceptible ones (S), by quantitative real-time PCR. a Cqm1 gene. b Caspase-3 (Casp) and pantetheinase (PTT) genes. c DNA polymerase delta catalytic subunit (Pol-delta) and Rac serine/threonine kinase (Rac-Ser-Threo) genes. Gene expression levels are relative to those from the endogenous control gene 18S used for normalization. Means and standard errors were obtained from three biological replicates. *P < 0.0001, **P < 0.008, ***P < 0.005, ****P < 0.05

Full size image

Discussion

Resistance of Cx. quinquefasciatus to the Bin toxin from L. sphaericus is one of a few examples where resistance can be directly linked to mutations in a gene coding for the toxin receptor, in this case the Cqm1 maltase [7, 25,26,27,28,29]. Indeed, the expression pattern seen here for the cqm1 transcripts in resistant larvae, in agreement with the previous biological and molecular data [26, 36, 38], was a remarkable marker of Bin resistance in view of the substantial downregulation detected for this mRNA. The selection of resistance, however, possibly involves further changes and this study has identified several other genes and pathways whose expression are substantially altered in the resistant larvae. These are most likely associated with the strong resistance phenotype to Bin (RR > 3000-fold) rather than a response to the Cry48Aa/Cry49Aa toxin whose profile of low-resistance profile (RR~15-fold) detected in this study did not allow a reliable analysis. Nevertheless it is not possible to rule out that some differential responses exhibited by RIAB59 larvae could be related to their low-resistance to the Cry48Aa/Cry49Aa toxins. These changes could not be necessarily related to receptors, since Bin and Cry48Aa/Cry49Aa toxins do not depend on the same binding sites [34, 63], but eventually to other pathways specifically associated to their mode of action.

While the key role of Cqm1 as the Cx. quinquefasciatus receptor for the Bin toxin is known, candidate midgut ligands of the Cry48Aa/Cry49Aa toxin have only recently been identified, but have not yet been functionally validated as specific receptors. Those putative ligands include aminopeptidases/metalloproteases, alkaline phosphatases and maltases, some of them orthologs to molecules previously identified as receptors for other insecticidal Cry toxins [35]. In this study, however, none of these candidates showed any significant changes in expression at the transcriptional level. This might be linked to the low resistance level to the Cry48Aa/Cry49Aa toxin detected here for the RIAB59 larvae, as described above, although mutations changing binding sites would not be detected by our approach nor would changes in genes inducible by exposure to the toxins. A possible exception, nevertheless, was the vanin 1/pantetheinase, since this protein was found among the previously identified candidate and the downregulation of its gene detected here is higher than that observed for the cqm1 gene. At this stage, however, pantetheinase’s downregulation might rather reflect changes in the lipid metabolism associated to the resistance to the Bin toxin. This is due to the fact that vitamin B5, the product of the vanin 1/pantetheinase, is an important cofactor for lipid biogenesis and degradation [48]. In addition, an RNA-seq study using larvae resistant to the 2362 strain that produces Bin toxin only showed a similar decrease in levels of pantetheinase mRNAs (our unpublished data), which strongly suggests its association with the Bin toxin resistance. Another protein with possible alternative roles in the resistance phenotype is apolipoprotein D (ApoD), a carrier of lipophilic molecules that has been shown to play a critical role in anti-oxidation and anti-apoptosis [64]. ApoD is a member of the lipocalin protein family, consisting of a large and diverse group of extracellular proteins with several functions in cellular regulation and homeostasis [65]. ApoD has also recently been shown to have important functions associated with stress response, longevity regulation and control of metabolism [64].

The downregulation of genes related to lipid catabolism in the resistant colony corroborates a previous finding which detected a remarkable accumulation of lipid inclusions in the midgut epithelial cells from non-treated fourth instar larvae from another Cx. quinquefasciatus colony (R2362) resistant to the Bin toxin [15]. Lipid synthesis in mosquitos has been observed as a response to bacterial pore-forming toxins [66, 67] and recently the accumulation of lipid droplets observed in Ae. aegypti challenged with viruses and bacteria was also clearly correlated with its involvement in mosquito immunity [68]. An intriguing observation is the strong upregulation of the gene encoding peroxisomal 2-hydroxyacyl-CoA lyase, presumably also involved in lipid degradation. In animal cells, fatty acid oxidation can occur in both mitochondria and peroxisomes [69] with these two ubiquitous organelles being metabolically connected, so that any alteration in mitochondrial function may induce changes in peroxisomal physiology [69,70,71]. It would then not be surprising if a chronic mitochondrial dysfunction induced as an adaptive host response to the action of pore-forming toxins might lead to a compensatory change in the peroxisomes as a consequence of the resistance event.

Upregulation of genes involved in DNA metabolism was one of the most marked differential expression feature exhibited by the resistant colony, and indicates that the resistant colony developed an expression profile to respond to DNA damage. The mode of action of bacterial pore forming toxins also involves apoptosis and direct evidence has been specifically provided for Bin [17] and Cry toxins [72]. The intrinsic apoptotic pathway is specifically activated in response to mitochondria damage, resulting in the cytochrome C release from these organelles and caspases activation that leads to other events, including DNA degradation [72]. A higher expression of caspase 3 in susceptible larvae treated with Bin toxin was an effect caused by this toxin to the cells [73]. The action of Bin toxin induces autophagy in the midgut cells and subsequent cytoplasmic vacuolization as it has been well documented in larvae subjected to Bin treatment [16, 17]. In resistant larvae several genes from the autophagy pathway (cqu04140) show differential regulation and these may be related to the lower capacity by the Bin toxin to act upon the resistant cells and provoke autophagy.

The identification of several genes and pathways whose expression are substantially changed in the resistant colony provides a framework for understanding some of the mechanisms involved in the resistance to L. sphaericus toxins, mainly Bin, as well as more generall mechanisms required by the larvae to survive the exposure to the bacterium. Indeed, a broader profile of features associated with resistance of such larvae, other than the major gene already shown to be required for resistance to the Bin toxin, is essential for the understanding the strategies used by these insects for their survival and for their ability to overcome the toxic effects. This aspect is of particular importance for the management of insecticide resistance in the field.

Conclusions

The dataset provided here shows that the resistant RIAB59 larvae is associated with a gene expression profile distinct from susceptible larvae. This is correlated not only with lack of the receptor which primarily confers the Bin resistance status, but also with metabolic pathways that might enhance the survival of the larvae to the effects mainly caused by the Bin toxin from the IAB59 strain. DEGs in pathways associated with lipid catabolism, DNA metabolism, apoptosis and immune response were specifically regulated and are likely to be associated with the response to L. sphaericus and with the stable resistance phenotype. This study also highlights other adaptative features that are likely to be crucial for the capacity of these individuals to successfully maintain this phenotype on a long-term basis.

Availability of data and materials

The transcriptome dataset is available (release date 2019-09-30) at https://www.ncbi.nlm.nih.gov/sra/PRJNA551311.

Abbreviations

Bin:

Binary

Cqm1:

Culex quinquefasciatus maltase 1

DEG:

differentially expressed genes

GPI:

glycosylphosphatidylinositol

KEGG:

Kyoto Encyclopedia of Genes and Genomes

RIAB59:

Culex quinquefasciatus strain resistant to Lysinibacillus sphaericus strain IAB59

RR:

ratio of resistance

References

  1. Lacey LA, Grzywacz D, Shapiro-Ilan DI, Frutos R, Brownbridge M, Goettel MS. Insect pathogens as biological control agents: back to the future. J Invertebr Pathol. 2015;132:1–41.

    Article  CAS  Google Scholar 

  2. Lacey LA. Bacillus thuringiensis serovariety israelensis and Bacillus sphaericus for mosquito control. J Am Mosq Control Assoc. 2007;23:133–63.

    Article  CAS  Google Scholar 

  3. Silva Filha MHNL, Berry C, Regis L. Lysinibacillus sphaericus: toxins and mode of action, applications for mosquito control and resistance management. In: Dhadialla TS, Gill SS, editors. Insect midgut and insecticidal proteins. Oxford: Academic Press; 2014. p. 89–176.

    Chapter  Google Scholar 

  4. Berry C. The bacterium, Lysinibacillus sphaericus, as an insect pathogen. J Invertebr Pathol. 2012;109:1–10.

    Article  Google Scholar 

  5. Broadwell AH, Baumann P. Proteolysis in the gut of mosquito larvae results in further activation of the Bacillus sphaericus toxin. Appl Environ Microbiol. 1987;53:1333–7.

    CAS  PubMed  PubMed Central  Google Scholar 

  6. Nielsen-LeRoux C, Charles J-F. Binding of Bacillus sphaericus binary toxin to a specific receptor on midgut brush-border membranes from mosquito larvae. Eur J Biochem. 1992;210:585–90.

    Article  CAS  Google Scholar 

  7. Romão TP, Chalegre KDM, Key S, Ayres CFJ, Oliveira CMF, De-Melo-Neto OP, et al. A second independent resistance mechanism to Bacillus sphaericus binary toxin targets its alpha-glucosidase receptor in Culex quinquefasciatus. FEBS J. 2006;273:1556–68.

    Article  Google Scholar 

  8. Silva-Filha MH, Nielsen-LeRoux C, Charles J-F. Identification of the receptor for Bacillus sphaericus crystal toxin in the brush border membrane of the mosquito Culex pipiens (Diptera: Culicidae). Insect Biochem Mol Biol. 1999;29:711–21.

    Article  CAS  Google Scholar 

  9. Nascimento NAD, Ferreira LM, Romão TP, Correia DMDC, Vasconcelos CRDS, Rezende AM, et al. N-glycosylation influences the catalytic activity of mosquito α-glucosidases associated with susceptibility or refractoriness to Lysinibacillus sphaericus. Insect Biochem Mol Biol. 2017;81:62–71.

    Article  Google Scholar 

  10. Sharma M, Lakshmi A, Gupta GD, Kumar V. Mosquito-larvicidal binary toxin receptor protein (Cqm1): crystallization and X-ray crystallographic analysis. Acta Crystallogr Sect F Struct Biol Commun. 2018;74:571–7.

    Article  CAS  Google Scholar 

  11. Boonserm P, Moonsom S, Boonchoy C, Promdonkoy B, Parthasarathy K, Torres J. Association of the components of the binary toxin from Bacillus sphaericus in solution and with model lipid bilayers. Biochem Biophys Res Commun. 2006;342:1273–8.

    Article  CAS  Google Scholar 

  12. Lekakarn H, Promdonkoy B, Boonserm P. Interaction of Lysinibacillus sphaericus binary toxin with mosquito larval gut cells: binding and internalization. J Invertebr Pathol. 2015;132:125–31.

    Article  CAS  Google Scholar 

  13. Pauchet Y, Luton F, Castella C, Charles J-F, Romey G, Pauron D. Effects of a mosquitocidal toxin on a mammalian epithelial cell line expressing its target receptor. Cell Microbiol. 2005;7:1335–44.

    Article  CAS  Google Scholar 

  14. Charles JF. Ultrastructural midgut events in Culicidae larvae fed with Bacillus sphaericus 2297 spore/crystal complex. Ann Inst Pasteur Microbiol. 1987;138:471–84.

    Article  CAS  Google Scholar 

  15. de Melo JV, Vasconcelos RHT, Furtado AF, Peixoto CA, Silva-Filha MHNL. Ultrastructural analysis of midgut cells from Culex quinquefasciatus (Diptera:Culicidae) larvae resistant to Bacillus sphaericus. Micron. 2008;39:1342–50.

    Article  Google Scholar 

  16. Opota O, Gauthier NC, Doye A, Berry C, Gounon P, Lemichez E, et al. Bacillus sphaericus binary toxin elicits host cell autophagy as a response to intoxication. PLoS ONE. 2011;6:e14682.

    Article  CAS  Google Scholar 

  17. Tangsongcharoen C, Chomanee N, Promdonkoy B, Boonserm P. Lysinibacillus sphaericus binary toxin induces apoptosis in susceptible Culex quinquefasciatus larvae. J Invertebr Pathol. 2015;128:57–63.

    Article  CAS  Google Scholar 

  18. Chevillon C, Bernard C, Marquine M, Pasteur N. Resistance to Bacillus sphaericus in Culex pipiens (Diptera: Culicidae): interaction between recessive mutants and evolution in southern France. J Med Entomol. 2001;38:657–64.

    Article  CAS  Google Scholar 

  19. Nielsen-LeRoux C, Charles JF, Thiéry I, Georghiou GP. Resistance in a laboratory population of Culex quinquefasciatus (Diptera: Culicidae) to Bacillus sphaericus binary toxin is due to a change in the receptor on midgut brush-border membranes. Eur J Biochem. 1995;228:206–10.

    Article  CAS  Google Scholar 

  20. Nielsen-LeRoux C, Pasteur N, Me Pr Tre JR, Charles J, Sheikh H, Chevillon C. High resistance to Bacillus sphaericus binary toxin in Culex pipiens (Diptera: Culicidae): the complex situation of west mediterranean countries. J Med Entomol. 2002;39:729–35.

    Article  CAS  Google Scholar 

  21. Pei G, Oliveira CMF, Yuan Z, Nielsen-LeRoux C, Silva-Filha MH, Yan J, et al. A strain of Bacillus sphaericus causes slower development of resistance in Culex quinquefasciatus. Appl Environ Microbiol. 2002;68:3003–9.

    Article  CAS  Google Scholar 

  22. Rao DR, Mani TR, Rajendran R, Joseph AS, Gajanana A, Reuben R. Development of a high level of resistance to Bacillus sphaericus in a field population of Culex quinquefasciatus from Kochi, India. J Am Mosq Control Assoc. 1995;11:1–5.

    CAS  PubMed  Google Scholar 

  23. Wirth MC, Georghiou GP, Malik JI, Abro GH. Laboratory selection for resistance to Bacillus sphaericus in Culex quinquefasciatus (Diptera: Culicidae) from California, USA. J Med Entomol. 2000;37:534–40.

    Article  CAS  Google Scholar 

  24. Oliveira CMF, Silva-Filha MH, Nielsen-LeRoux C, Pei G, Yuan Z, Regis L. Inheritance and mechanism of resistance to Bacillus sphaericus in Culex quinquefasciatus (Diptera: Culicidae) from China and Brazil. J Med Entomol. 2004;41:58–64.

    Article  CAS  Google Scholar 

  25. Chalegre KDM, Romão TP, Tavares DA, Santos EM, Ferreira LM, Oliveira CMF, et al. Novel mutations associated with resistance to Bacillus sphaericus in a polymorphic region of the Culex quinquefasciatus cqm1 gene. Appl Environ Microbiol. 2012;78:6321–6.

    Article  CAS  Google Scholar 

  26. Chalegre KDM, Tavares DA, Romão TP, de Menezes HSG, Nascimento NA, Oliveira CMF, et al. Co-selection and replacement of resistance alleles to Lysinibacillus sphaericus in a Culex quinquefasciatus colony. FEBS J. 2015;282:3592–602.

    Article  CAS  Google Scholar 

  27. Darboux I, Charles JF, Pauchet Y, Warot S, Pauron D. Transposon-mediated resistance to Bacillus sphaericus in a field-evolved population of Culex pipiens (Diptera: Culicidae). Cell Microbiol. 2007;9:2022–9.

    Article  CAS  Google Scholar 

  28. Darboux I, Pauchet Y, Castella C, Silva-Filha MH, Nielsen-LeRoux C, Charles J-F, et al. Loss of the membrane anchor of the target receptor is a mechanism of bioinsecticide resistance. Proc Natl Acad Sci USA. 2002;99:5830–5.

    Article  CAS  Google Scholar 

  29. Guo QY, Cai QX, Yan JP, Hu X-M, Zheng DS, Yuan ZM. Single nucleotide deletion of cqm1 gene results in the development of resistance to Bacillus sphaericus in Culex quinquefasciatus. J Insect Physiol. 2013;59:967–73.

    Article  CAS  Google Scholar 

  30. Nielsen-LeRoux C, Rao DR, Murphy JR, Carron A, Mani TR, Hamon S, et al. Various levels of cross-resistance to Bacillus sphaericus strains in Culex pipiens (Diptera: Culicidae) colonies resistant to B. sphaericus strain 2362. Appl Environ Microbiol. 2001;67:5049–54.

    Article  CAS  Google Scholar 

  31. Jones GW, Nielsen-LeRoux C, Yang Y, Yuan Z, Dumas VF, Gomes Monnerat R, et al. A new Cry toxin with a unique two-component dependency from Bacillus sphaericus. FASEB J. 2007;21:4112–20.

    Article  CAS  Google Scholar 

  32. Jones GW, Wirth MC, Monnerat RG, Berry C. The Cry48Aa-Cry49Aa binary toxin from Bacillus sphaericus exhibits highly restricted target specificity. Environ Microbiol. 2008;10:2418–24.

    Article  CAS  Google Scholar 

  33. Guo QY, Hu XM, Cai QX, Yan JP, Yuan ZM. Interaction of Lysinibacillus sphaericus Cry48Aa/Cry49Aa toxin with midgut brush-border membrane fractions from Culex quinquefasciatus larvae. Insect Mol Biol. 2016;25:163–70.

    Article  CAS  Google Scholar 

  34. de Melo JV, Jones GW, Berry C, Vasconcelos RHT, Oliveira CMF, Furtado AF, et al. Cytopathological effects of Bacillus sphaericus Cry48Aa/Cry49Aa toxin on binary toxin-susceptible and -resistant Culex quinquefasciatus larvae. Appl Environ Microbiol. 2009;75:4782–9.

    Article  Google Scholar 

  35. Rezende TMT, Romão TP, Batista M, Berry C, Adang MJ, Silva-Filha MHNL. Identification of Cry48Aa/Cry49Aa toxin ligands in the midgut of Culex quinquefasciatus larvae. Insect Biochem Mol Biol. 2017;88:63–70.

    Article  CAS  Google Scholar 

  36. Amorim LB, Oliveira CMF, Rios EM, Regis L, Silva-Filha MHNL. Development of Culex quinquefasciatus resistance to Bacillus sphaericus strain IAB59 needs long term selection pressure. Biol Control. 2007;42:155–60.

    Article  Google Scholar 

  37. Silva-Filha MHNL, Oliveira CMF, Regis L, Yuan Z, Rico CM, Nielsen-LeRoux C. Two Bacillus sphaericus binary toxins share the midgut receptor binding site: implications for resistance of Culex pipiens complex (Diptera: Culicidae) larvae. FEMS Microbiol Lett. 2004;241:185–91.

    Article  CAS  Google Scholar 

  38. Amorim LB, de Barros RA, Chalegre KDM, Oliveira CMF, Narcisa Regis L, Silva-Filha MHNL. Stability of Culex quinquefasciatus resistance to Bacillus sphaericus evaluated by molecular tools. Insect Biochem Mol Biol. 2010;40:311–6.

    Article  CAS  Google Scholar 

  39. Priest FG, Ebdrup L, Zahner V, Carter PE. Distribution and characterization of mosquitocidal toxin genes in some strains of Bacillus sphaericus. Appl Environ Microbiol. 1997;63:1195–8.

    CAS  PubMed  PubMed Central  Google Scholar 

  40. Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29:15–21.

    Article  CAS  Google Scholar 

  41. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550.

    Article  Google Scholar 

  42. Götz S, García-Gómez JM, Terol J, Williams TD, Nagaraj SH, Nueda MJ, et al. High-throughput functional annotation and data mining with the Blast2GO suite. Nucleic Acids Res. 2008;36:3420–35.

    Article  Google Scholar 

  43. Franceschini A, Szklarczyk D, Frankild S, Kuhn M, Simonovic M, Roth A, et al. STRING v9.1: protein-protein interaction networks, with increased coverage and integration. Nucleic Acids Res. 2013;41:D808–15.

    Article  CAS  Google Scholar 

  44. Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28:27–30.

    Article  CAS  Google Scholar 

  45. Luo W, Brouwer C. Pathview: an R/Bioconductor package for pathway-based data integration and visualization. Bioinformatics. 2013;29:1830–1.

    Article  CAS  Google Scholar 

  46. Liu N, Li T, Reid WR, Yang T, Zhang L. Multiple Cytochrome P450 Genes: their constitutive overexpression and permethrin induction in insecticide resistant mosquitoes, Culex quinquefasciatus. PLoS ONE. 2011;6:e23403.

    Article  CAS  Google Scholar 

  47. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-delta delta C(T)) method. Methods. 2001;25:402–8.

    Article  CAS  Google Scholar 

  48. Kaskow BJ, Proffit JM, Blangero J, Moses EK, Abraham LJ. Diverse biological activities of the vascular non-inflammatory molecules—the vanin pantetheinases. Biochem Biophys Res Commun. 2012;417:653–8.

    Article  CAS  Google Scholar 

  49. Shu C, Tan S, Yin J, Soberón M, Bravo A, Liu C, et al. Assembling of Holotrichia parallela (dark black chafer) midgut tissue transcriptome and identification of midgut proteins that bind to Cry8Ea toxin from Bacillus thuringiensis. Appl Microbiol Biotechnol. 2015;99:7209–18.

    Article  CAS  Google Scholar 

  50. Bennett V, Baines AJ. Spectrin and ankyrin-based pathways: metazoan inventions for integrating cells into tissues. Physiol Rev. 2001;81:1353–92.

    Article  CAS  Google Scholar 

  51. Dubreuil RR. Functional links between membrane transport and the spectrin cytoskeleton. J Membr Biol. 2006;211:151–61.

    Article  CAS  Google Scholar 

  52. Popelka H, Klionsky DJ. Structural basis for extremely strong binding affinity of giant ankyrins to LC3/GABARAP and its application in the inhibition of autophagy. Autophagy. 2018;14:1847–9.

    Article  Google Scholar 

  53. Yang Y, Zhang Y, Ding P, Johnson K, Li X, Zhang Y. The ankyrin-repeat transmembrane protein BDA1 functions downstream of the receptor-like protein SNC2 to regulate plant immunity. Plant Physiol. 2012;159:1857–65.

    Article  CAS  Google Scholar 

  54. Wang P, Granados RR. Molecular structure of the peritrophic membrane (PM): identification of potential PM target sites for insect control. Arch Insect Biochem Physiol. 2001;47:110–8.

    Article  CAS  Google Scholar 

  55. Tetreau G, Dittmer NT, Cao X, Agrawal S, Chen Y, Muthukrishnan S, et al. Analysis of chitin-binding proteins from Manduca sexta provides new insights into evolution of peritrophin A-type chitin-binding domains in insects. Insect Biochem Mol Biol. 2015;62:127–41.

    Article  CAS  Google Scholar 

  56. Casteels M, Sniekers M, Fraccascia P, Mannaerts GP, van Veldhoven PP. The role of 2-hydroxyacyl-CoA lyase, a thiamin pyrophosphate-dependent enzyme, in the peroxisomal metabolism of 3-methyl-branched fatty acids and 2-hydroxy straight-chain fatty acids: figure 1. Biochem Soc Trans. 2007;35:876–80.

    Article  CAS  Google Scholar 

  57. Jenkins B, de Schryver E, Van Veldhoven PP, Koulman A. Peroxisomal 2-hydroxyacyl-CoA lyase is involved in endogenous biosynthesis of heptadecanoic ccid. Molecules. 2017;22:1718.

    Article  Google Scholar 

  58. Mayoral JG, Nouzova M, Navare A, Noriega FG. NADP + -dependent farnesol dehydrogenase, a corpora allata enzyme involved in juvenile hormone synthesis. Proc Natl Acad Sci USA. 2009;106:21091–6.

    Article  CAS  Google Scholar 

  59. Vannini L, Ciolfi S, Spinsanti G, Panti C, Frati F, Dallai R. The putative-farnesoic acid O-methyl transferase (FAMeT) gene of Ceratitis capitata: characterization and pre-imaginal life expression. Arch Insect Biochem Physiol. 2010;73:106–17.

    CAS  PubMed  Google Scholar 

  60. Ma YJ, Lee BL, Garred P. An overview of the synergy and crosstalk between pentraxins and collectins/ficolins: their functional relevance in complement activation. Exp Mol Med. 2017;49:e320.

    Article  CAS  Google Scholar 

  61. Wesener DA, Dugan A, Kiessling LL. Recognition of microbial glycans by soluble human lectins. Curr Opin Struct Biol. 2017;44:168–78.

    Article  CAS  Google Scholar 

  62. Burmester T, Scheller K. Ligands and receptors: common theme in insect storage protein transport. Naturwissenschaften. 1999;86:468–74.

    Article  CAS  Google Scholar 

  63. Rezende TMT, Romão TP, Batista M, Berry C, Adang MJ, Silva-Filha MHNL. Identification of Cry48Aa/Cry49Aa toxin ligands in the midgut of Culex quinquefasciatus larvae. Insect Biochem Mol Biol. 2017;88:63–70.

    Article  CAS  Google Scholar 

  64. Zhou Y, Wang L, Li R, Liu M, Li X, Su H, et al. Secreted glycoprotein BmApoD1 plays a critical role in anti-oxidation and anti-apoptosis in Bombyx mori. Biochem Biophys Res Commun. 2018;495:839–45.

    Article  CAS  Google Scholar 

  65. Flower DR. The lipocalin protein family: a role in cell regulation. FEBS Lett. 1994;354:7–11.

    Article  CAS  Google Scholar 

  66. 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:1135–45.

    Article  CAS  Google Scholar 

  67. Canton PE, Cancino-Rodezno A, Gill SS, Soberón M, Bravo A. Transcriptional cellular responses in midgut tissue of Aedes aegypti larvae following intoxication with Cry11Aa toxin from Bacillus thuringiensis. BMC Genomics. 2015;16:1042.

    Article  Google Scholar 

  68. Barletta ABF, Alves LR, Silva MCLN, Sim S, Dimopoulos G, Liechocki S, et al. Emerging role of lipid droplets in Aedes aegypti immune response against bacteria and Dengue virus. Sci Rep. 2016;6:19928.

    Article  CAS  Google Scholar 

  69. Fransen M, Lismont C, Walton P. The peroxisome-mitochondria connection: how and why? Int J Mol Sci. 2017;18:1126.

    Article  Google Scholar 

  70. Le Borgne F, Demarquoy J. Interaction between peroxisomes and mitochondria in fatty acid metabolism. Open J Mol Integr Physiol. 2012;2:27–33.

    Article  Google Scholar 

  71. Demarquoy J. Crosstalk between mitochondria and peroxisomes. World J Biol Chem. 2015;6:301.

    Article  Google Scholar 

  72. Portugal L, Muñóz-Garay C, de Castro DL, Soberón M, Bravo A. Toxicity of Cry1A toxins from Bacillus thuringiensis to CF1 cells does not involve activation of adenylate cyclase/PKA signaling pathway. Insect Biochem Mol Biol. 2017;80:21–31.

    Article  CAS  Google Scholar 

  73. Tangsongcharoen C, Jupatanakul N, Promdonkoy B, Dimopoulos G, Boonserm P. Molecular analysis of Culex quinquefasciatus larvae responses to Lysinibacillus sphaericus bin toxin. PLoS ONE. 2017;12:e0175473.

    Article  Google Scholar 

Download references

Acknowledgements

The authors thank the Programme for Technological Development in Tools for Health-PDTIS-FIOCRUZ for use of its facilities and the insectarium from IAM-FIOCRUZ for the technical support.

Funding

This work was supported by Conselho Nacional de Pesquisa-Brazil (CNPq Grant 458492/2014-0) and Fundação de Amparo à Ciência e Tecnologia de Pernambuco-Brazil (FACEPE Grants APQ-1659-2.01/15, APQ-1616-2.13/15 and IBPG-1364-2.13/13).

Author information

Authors and Affiliations

  1. Instituto Aggeu Magalhães-FIOCRUZ, Av. Moraes Rego s/n Cidade Universitária, Recife, PE, 50740-465, Brazil

    Tatiana Maria Teodoro Rezende, Antonio Mauro Rezende, Gabriel Luz Wallau, Crhisllane Rafaele Santos Vasconcelos, Osvaldo Pompílio de-Melo-Neto, Maria Helena Neves Lobo Silva-Filha & Tatiany Patrícia Romão

Authors
  1. Tatiana Maria Teodoro Rezende
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Antonio Mauro Rezende
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Gabriel Luz Wallau
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Crhisllane Rafaele Santos Vasconcelos
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Osvaldo Pompílio de-Melo-Neto
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Maria Helena Neves Lobo Silva-Filha
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Tatiany Patrícia Romão
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

TMTR: acquisition, analysis and interpretation of data, manuscript writing and revision. AMR: data analyses and manuscript writing. GLW: data analyses, interpretation of the data and manuscript writing. CRSV: data analysis and manuscript writing. OPDMN: data analysis, manuscript writing and intellectual contributions. MHNLSF: conception of the study, data analysis and manuscript writing. TPR: conception of the study, definition of goals, methods and colony to be studied, experimental execution support, data analysis, manuscript writing, revision and financial support. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Maria Helena Neves Lobo Silva-Filha or Tatiany Patrícia Romão.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional files

Additional file 1: Table S1.

Differentially expressed genes of Culex quinquefasciatus larvae from a Lysinibacillus sphaericus resistant colony (RIAB59) compared to a susceptible colony.

Additional file 2: Table S2.

Primers used to perform qRT-PCR reactions to evaluate the expression of Culex quinquefasciatus genes.

Additional file 3: Figure S1.

Examples of KEGG pathways that displayed mostly upregulated genes related to DNA synthesis and maintenance in Culex quinquefasciatus larvae from a Lysinibacillus sphaericus resistant colony compared to a susceptible one. a Nucleotide excision repair (cqu03420). b Mismatch repair (cqu03430). c Homologous recombination (cqu03440). d Non-homologous end-joining (cqu03450).

Additional file 4: Table S3.

Description of STRING clusters found for downregulated and upregulated genes of Culex quinquefasciatus larvae from a Lysinibacillus sphaericus resistant colony (RIAB59) compared to a susceptible colony.

Rights and permissions

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

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Rezende, T.M.T., Rezende, A.M., Luz Wallau, G. et al. A differential transcriptional profile by Culex quinquefasciatus larvae resistant to Lysinibacillus sphaericus IAB59 highlights genes and pathways associated with the resistance phenotype. Parasites Vectors 12, 407 (2019). https://doi.org/10.1186/s13071-019-3661-y

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s13071-019-3661-y

Keywords

Parasites & Vectors

ISSN: 1756-3305

Contact us