Transcriptomic Analysis of Aedes aegypti in Response to Mosquitocidal Bacillus thuringiensis LLP29 Toxin

Globally, Aedes aegypti is one of the most dangerous mosquitoes that plays a crucial role as a vector for human diseases, such as yellow fever, dengue, and chikungunya. To identify (1) transcriptomic basis of midgut (2) key genes that are involved in the toxicity process by a comparative transcriptomic analysis between the control and Bacillus thuringiensis (Bt) toxin (LLP29 proteins)-treated groups. Next-generation sequencing technology was used to sequence the midgut transcriptome of A. aegypti. A total of 17130 unigenes, including 574 new unigenes, were identified containing 16358 (95.49%) unigenes that were functionally annotated. According to differentially expressed gene (DEG) analysis, 557 DEGs were annotated, including 226 upregulated and 231 downregulated unigenes in the Bt toxin-treated group. A total of 442 DEGs were functionally annotated; among these, 33 were specific to multidrug resistance, 6 were immune-system-related (Lectin, Defensin, Lysozyme), 28 were related to putative proteases, 7 were lipase-related, 8 were related to phosphatases, and 30 were related to other transporters. In addition, the relative expression of 28 DEGs was further confirmed through quantitative real time polymerase chain reaction. The results provide a transcriptomic basis for the identification and functional authentication of DEGs in A. aegypti.

and has been found to be highly toxic to mosquitoes. It is a promising Bt strain to control mosquitoes because of its high toxicity and short life cycle 19,20 . The Cry toxins are mainly involved in (1) protoxin solubilization, (2) protoxin proteolytic activation by specific proteases, (3) interaction between putative receptors and active toxins, (4) oligomerization of toxin, and (5) insertion to epithelial cells. Thus, forming pores in the plasma membrane of midgut cell, and eventually leading to cell death 16,17 . The use of Cry toxins against insect pests depends on the interaction between Cry toxins and other receptors, such as alkaline phosphatase (ALP), aminopeptidase-N (APN), and ATP-binding cassette (ABC) transporters 16,[21][22][23][24] . Therefore, it is important to understand the interaction of Cry toxins with other midgut proteins.
As toxin binding is partially dependent on receptor glycosylation 25 , proteins that interact with carbohydrate-recognition domains (CRD) might affect toxicity. Furthermore, the transcriptomic sequencing is a useful method to determine the composition and function of DEG genes of many complex insects, including Drosophila melanogaster 26 , A. aegypti and Anopheles gambiae 27 , Bombyx mori 28 , and Maruca vitrata 29 . In the present study, next-generation sequencing technology was used to sequence the transcriptome of A. aegypti exposed to LLP29 protein, and successfully built a transcript database; In addition, putative insecticide resistance transcript, immune-system related transcript, stressed-related transcript, detoxification-related molecules, and selected DEGs were confirmed through quantitative real time polymerase chain reaction (qPCR). Furthermore, the present study provides valuable information that can be used to develop new genetics-based strategies and novel molecular tools to control A. aegypti.

Results and Discussion
Sequence analysis and functional annotation. LLP29 is highly toxic to many mosquitoes, such as A.
aegypti, A. albopictus and Culex quinquefasciatus 19,20 . Meanwhile, because of fewer plasmids and short life cycle, it is good for reducing the cost of Bt application and has potential as a good biocontrol agent in mosquitoes 20 . In order to further understand the transcriptomic basis of midgut and its main functional genes after Bti LLP29 infection, transcriptomic analysis of A. aegypti in response to LLP29 was carried out in the present research, other than only single Cry toxin, such as Cry11Aa 30 . After complete sequencing and quality check, a total of 9.32 Gb of clean data were obtained. The Q30 base percentage of each sample was above 85.77%, which showed reliable base identification and limited base error detection. We generated 40,678,664 and 33,418,822 raw reads from the control and Bacillus thuringiensis (Bt) toxin (LLP29 proteins)-treated groups and 29,434,409, 24,268,858 reads were mapped to reference genome of A. aegypti. A total of 17130 unigene sequences were identified in the midgut of A. aegypti (Table S1). These unigenes were compared by BLAST with those in Swiss-Prot 31 , GO 32 , NR 33 , COG 34 , and KEGG databases 35 . A total of 16358 (95.49%) unigenes were functionally annotated by BLAST searches against all databases (Table 1). Furthermore, with the development of sequencing technique, 574 new unigenes that not present in any published annotation were identified in the transcriptome of A. aegypti, including 391 unigenes that were functionally annotated through different databases (Table 2). They will contribute further to the transcriptome of A. aegypti and bioinformatic approaches, and greatly benefit the knowledge of Bt mechanism on immune response 36 .
Furthermore, all unigenes were functionally classified based on the KEGG pathway analysis. A total of 6811 unigenes were annotated from the KEGG pathway database, and were distributed into 258 KEGG pathways (Table S4). The highly abundant KEGG pathways were RNA transport (ko03013) with 168 unigenes, ribosome (ko03010) with 162 unigenes, purine metabolism (ko00230) with 158 unigenes, protein processing in the endoplasmic reticulum (ko04141) with 151 unigenes, spliceosome (ko03040) with 136 unigenes. Ninety nine unigenes were related to three different xenobiotic-metabolism pathways related to cytochrome p450 (ko00980, ko00982, and ko00983) and other 28 unigenes were related to the innate immune system, such as toll-like receptor (TLR) pathways (ko04620), and reactive to Gram-positive bacterial and fungal infections 39 . The other 32 putative unigenes were related to the TNF pathway (ko04668), while the IMD signaling pathway was considered a homolog of TNF pathway based on the records of the KEGG map and responsive on Gram-negative bacteria 39,40 . Further, the JAK/STAT pathway was found to be related to 26 unigenes, and the JAK/STAT pathway (ko04630) is involved in multiple developmental processes and immune responses 40 , such as inhibition of microbial infection and repair of septic damage 40,41 . Differential gene expression. To identify specific target genes that affect different biological processes, the gene expression profiles were compared to examine changes in gene activity among the control and Bti LLP29 protein-treated groups. To identify significant differences in gene expression, FPKM with fold change >2 and FDR <0.01 were used as the thresholds. Based on this parameter, 457 unigenes were found to be differentially expressed in the midgut of A. aegypti (Table S5). Among the differentially expressed genes, 226 were upregulated and 231 were downregulated in response to Bti toxin. The number of upregulated unigenes was nearly similar to the number of downregulated unigenes ( Fig. 4; Table S5). In this study, several homologous genes potentially implicated in the immune response to LLP29 toxin were found to be highly differentially expressed and similar results were also obtained in other previous studies in response to insecticides or single toxin, such as Cry11Aa 30,42-44 .  Table S7). To identify the COG function, 165 DEGs were identified, which were further classified into 21 main categories basing on the COG function classification. The amino acid transport and metabolism group (16.69%) is the highest class, followed by the general function prediction (16.2%) (Fig. 6). According to the KEGG pathway analysis, 172 differential genes were identified that are involved in various biological pathways of A. aegypti ( Fig. 7; Table S8). Further, these 172 DEGs were annotated from the KEGG pathway database, and the most enriched unigenes were mapped into 50 different KEGG pathways. They were further categorized into 5 classes: metabolism was the most represented class with 61 DEGs, including 26 different metabolism pathways which might play a role in LLP29 toxin resistance. Among them Alanine, aspartate and glutamate metabolism (ko00400), followed by Amino sugar and nucleotide sugar metabolism (ko00520), Glycolysis/Gluconeogenesis (ko00010), Pyruvate metabolism (ko00620), Biosynthesis of amino acids (ko01230), and Carbon metabolism (ko01200), a highly enriched pathways. This was followed by the organismal systems class with 27 DEGs, including 10 different system pathways, the most enriched was the Insulin signaling pathway (ko04911), followed by PPAR signaling pathway (ko03320), Estrogen signaling pathway (ko04915), Antigen processing and presentation (ko04612), Salivary secretion (ko04970) and GABAergic synapse pathway (ko04727). In environmental information processing class with 25 DEGs, including 9 different signal transduction pathways, the most enriched pathway was the PI3K-Akt signaling pathway (ko04151), followed by AMPK signaling pathway (ko04152) and MAPK signaling pathway (ko04010). However, signal transduction pathways also been reported to have an important role in response to Cry11Aa toxin 45,46 . Within genetic information processing class with 10 DEGs which included two different pathways, the protein processing in endoplasmic reticulum (ko04141) was the most highly enriched pathway and the other pathway was aminoacyl-tRNA biosynthesis (ko00970). The cellular processing class had 10 DEGs, including three pathways, focal adhesion (ko04510), peroxisome (ko04146), and phagosome (ko04145) (Fig. 7).
Specific protein types related to multidrug resistance and immune system. Based on the BLASTx results and searches for target genes related to protein families with specific multidrug resistance/defense, 68 multidrug resistance protein genes were identified in the midgut transcriptome of A. aegypti. The specific domain related to multidrug resistance was confirmed using Pfam and Interpro databases 47,48 . Among these unigenes, 11 ATP-binding cassettes (ABC) transporter unigenes were differentially expressed, including 5 unigenes that were upregulated and 6 unigenes that were downregulated. However, most of the (ABC) transporter unigenes were reported downregulated in the intoxication with Cry11Aa toxin to A. aegypti 30,44 . Furthermore, the highly expressed genes were also searched, and the results revealed 4 ABC transporter unigenes whose FPKM values were >100 (Table S9). Except for Cry11Aa, there are more mosquitocdial toxins harbored in LLP29, such as Cry4, Cry10, Cyt1 and Cyt 2, etc. With different toxins infection, the complex immune response of ABC transporter unigenes needs further study.
The cytochrome P450 is another multidrug resistance protein involved in many physiological functions, such as resistance to insecticides and hormone metabolism 49,50 . The association of cytochrome P450 with propoxur, dichlorvos, and pyrethroid resulted in insecticide resistance in C. pipiens mosquito 51 . A total of 175 cytochrome  44 ,and intoxication of Cry11Aa toxin to A. aegypti larvaes 30 , and these monooxygenases play an important role in the degradation of insecticides 52 . Furthermore, they also respond to Cry1Ab protoxin in Lepidopterans, C. fumiferana and M. sexta 53 . Eighteen highly expressed unigenes of cytochrome P450 with FPKM value >100 were identified in the present study (Table S9).
The cytosolic matrix glutathione S-transferases (GSTs) have been reported as important detoxification-related protein in the metabolism of insecticides 54 . In insects, the cytosolic matrix GSTs are classified into six groups (delta, epsilon, omega, sigma, theta, and zeta) according to the sequence similarity, immunoreactivity, specificity of substrate, and sensitivity to inhibitors 55,56 . Meanwhile, delta and epsilon classes are very unique in insects 57 . In the present study, 26 GST unigenes were identified. Among them, two GST unigenes were differentially expressed, including one upregulated and one downregulated while the same gene was found up regulated in response of Cry11Aa toxin to A. aegypti larvaes 30 . However, it is important to notice that we found similar results as coded by other previous studies 42, 44 , while six unigenes were highly expressed with a FPKM value > 100 (Table S9). The identification GST unigenes in the midgut transcriptome of A. aegypti is very important. Glutathione S-transferase plays a major role in the detoxification of secondary metabolites and insecticides 58,59 , thus providing high resistance to insecticides 60 .
Studies have reported that insect carboxylesterases are involved in xenobiotic metabolism, and have been classified into 13 clades 61 . In the present study, 46 putative carboxylesterases/ juvenile hormone esterases were also identified. Among these, 9 unigenes were differentially expressed, including 4 upregulated and 5 downregulated DEGs while 5 unigenes were highly expressed (FPKM value > 100) (Table S9). However most of the unigenes related to carboxylesterases were downregulated at the exposure of LLP29 toxin at 14 hours of treatment. It is notable that most of the carboxylesterases related unigenes were also down regulated in response to insecticides or toxins in other reported studies 30,43,44 .
Lectins are a class of carbohydrate binding proteins ubiquitously expressed in plants, animals, bacteria, and viruses 62 . The lectin is a major constituent of humoral immune system of insects, and is involved in the process of self/non-self-recognition. Interestingly, lectin proteins can recognize specific carbohydrate structures, and can bind to the cell surface sugars to agglutinate the cells 63 . Some insect lectins recognize polysaccharide chains on the surface of pathogens, and probably are involved in self-defense. In the present study, 47 lectin genes were identified in the midgut of A. aegypti, including G-specific type and C-type lectin genes. Only one differentially expressed lectin unigene, which was upregulated against the Bti LLP29 protein, was identified. Five highly expressed lectin unigenes with FPKM value > 100 (Table S9) were also identified. Further, defensin proteins have been reported to exhibit a strong innate immune response to a wide range of pathogenic bacteria 64 . In the present study, four defensin-related unigenes were identified; however, only one defensin-related unigene was differentially expressed and was downregulated in response to the Bti LLP29 protein (Table S9). Lysozyme is another antimicrobial enzyme that takes part in innate immunity. It catalyzes the hydrolysis of specific 1,4-beta linkages between N-acetylmuramic acid and N-acetyl-D glucosamine residues in peptidoglycan, thus lysing the bacterial cell 65 . In the present study, 9 unigenes related to lysozyme, including one upregulated and two downregulated DEGs (Table S9). These putative immune response unigenes were further validated by RT-qPCR (Fig. 8).
In the present study, the digestion-related genes, including various protease groups, were examined. The results revealed that 28 differentially expressed unigenes were putative proteases, including serine protease, trypsin, chymotrypsin, venom serine protease, carboxypeptidase, and other metalloproteinases. The serine proteases are involved in Cry protoxin activation via proteolytic elimination of peptide fragments 66 . In addition, 45 highly expressed unigenes related to a protease with FPKM value of >100 were identified (Table S10). Earlier studies have reported that chitinases catalyze the hydrolysis of chitin. As chitin is present in the cell wall of many green algae and fungi, and in the exoskeleton of numerous crustacean and insects 67 , chitinases can be used in the biocontrol of plant pathogenic fungi and insects 68 . In the present study, 4 DEGs related to chitinases, including 3 upregulated unigenes, were identified in the midgut of Bti toxin-treated A. aegypti. Five chitinase-related unigenes highly expressed with FPKM value > 100 were identified (Table S11). The chitinases related genes were also reported in previous studies in response to insecticides or toxins [42][43][44] . However, no chitinase related gene was reported in the single Cry11Aa toxin to A. aegypti larvaes 30 .
In the present study, 308 transporter protein-related unigenes were identified in the midgut of A. aegypti treated with Bti LLP29 protein. These included (1) 30 differentially expressed unigenes-5 sugar transporters, 7 sodium-dependent nutrient transporters, 3 amino acid transporters, 3 zinc transporters, 3 synaptic vesicle proteins, 2 mitochondrial carriers, and 7 different transporter-like genes; and (2) 13 highly expressed unigenes with FPKM value > 100 (Table S11). The genes related to lipid and fatty acid metabolism, binding, and transport have been reported to play a key role in the movement and homeostasis of lipid-related substances in the midgut 53,73 . However, based on BLASTx and Pfam search, no lipophorins were identified in the midgut by some studies 74,75 . In the present study, one triacylglycerol lipase 76 was identified. The hemolymph of some insects contains lipoproteins that bind to xenobiotics during microbial infection 75 . Therefore, the expression of lipid-related genes might be a defense mechanism in response to Bt LLP29 toxin. The combination of lipid-related and highly expressed transporter-related unigenes suggests the role of midgut in the transportation of lipid-bound proteins and peptides. qPCR analysis of selected different types of related genes. To verify the differential expression values of the transcriptomic data, qPCR was used in this study. A total of 28 genes were randomly selected from different families, including multidrug resistant genes: ABC transporter (3); cytochrome P450 (3); carboxylesterase (1); juvenile hormone (1); GST (1); lectin (1); defensin (1); lysozyme (1); proteases including serine proteinase (1), trypsin (1), carboxypeptidase (1), metalloproteinase (1), and venom serine protease; lipase (2); chitinase (1); alkaline phosphatase (1); inositol polyphosphate (1); sodium transporter (2); sugar transporter (2); zinc transporter protein (1); and ammonium transporter (1). Among these, majority of the genes presented relatively similar transcriptomic and qPCR data in terms of expression (Fig. 8). While the expression of some genes, such as gene7939, gene9072, gene9074, and gene9635, was 2-3 fold higher in the transcriptomic data when compared with that in the qPCR data (Fig. 8).

Materials and Methods
Sample preparation. A. aegypti were reared in an environmentally controlled room at 28 °C and 85% relative humidity (RH) with a photoperiod of 14:10 h light/dark. The Bt strain LLP29 was isolated from leaves of Mognolia denudate Desr. (Magnoliales: Magnoliaceae) 19,20 and were grown in nutrient broth sporulation medium containing erythromycin with final concentration of 50 ug/ml at 30 o C temperature for three days, following cell autolysis, spores and inclusion protein was harvested and washed three times with distilled water and 1 M NaCl. After resuspending the inclusions were purified according to the previous report 77 . The total protein (Bti toxin LLP29 protein) of concentration 0.255 μg/mL was fed to fourth instar A. aegypti larvae (100 larvae/cup) in 50 ml of dechlorinated water for 14 hours. For negative control, the larvae were not treated. After exposure the severely affected and dead larvae were discarded, and the survived larvae were used for RNA extraction. Both samples were used in triplicates. Midguts were dissected from both samples carefully under stereomicroscope (VHX-5000, KEYENCE), and placed in RNA-hold (Transgen biotech.) and stored at −80 degrees until further use.
Illumina sequencing and cDNA library. Total RNA was isolated from the midgut of the fourth-instar larvae using Trizol reagent (Invitrogen, USA). The isolated total RNA was precipitated with RNase free DNase I to remove genomic DNA using Amp Grade (Invitrogen, USA), following the instructions of the manufacturer. The quality of RNA was checked by agarose gel electrophoresis, and also by using the NanoDrop ND-1000 spectrophotometer (Thermo Scientific, USA). To prepare RNA samples, 1.5 μg RNA was used per sample as initial material. Sequence library was created through NEBNext1 Ultra ™ RNA Library Prep Kit for Illumina1 (NEB, USA), according to the instruction from the manufacturer. Messenger RNA was purified from the total RNA using poly-T oligo-attached magnetic beads. Fragmentation was carried out using RNA fragmentation kit, followed by the synthesis of first and second strand complementary DNAs using M-MuLV reverse transcriptase (RNase H-), DNA polymerase I, and RNaseH. The complementary DNA fragments were end-repaired and ligated to a NEB Next Adaptor after adenylation of the 3′ ends, and then purified with the AMPure XP system (Beckman Coulter, Beverly, USA) to create a cDNA library. The library quality was assessed on the Agilent Bioanalyzer 2100 system. After clustering, the products were sequenced on an Illumina Hiseq. 2500 platform to generate 100 bp paired-end reads (Biomarker Technologies Co., Ltd., Beijing, China). The Illumina GA Pipeline (version 1.3) was used to filter the raw data to eliminate the adapter sequences, low-quality reads and empty reads.

Bioinformatic analyses.
Reads were compared and mapped to reference genome of A. aegypti Liverpool strain (GCA_000004015.2) by using Bowtie and topHAT2 software 78 . RSEM software was used to convert the number of reads mapping to each unigenes as FPKM (Fragments Per Kilobase of transcript per Million mapped reads) and estimated the expression level of corresponding unigenes 79,80 . Furthermore, EB Seq software was used for differential expression analysis 81 . Fragments per kilobase of exon model per million mapped reads (FPKM) of fold change >2 and false discovery rate (FDR) of <0.01 were used as the thresholds to identify significant differences in gene expression. All unigenes were compared with the Swiss-Prot 31 , Gene Ontology (GO) 32 , Non Redundant (NR) 33 , Clusters of Ortholgous Groups (COG) 34 , Kyoto Encyclopedia of Genes and Genomes (KEGG) 35 , and other databases using the BLAST software 82 . The E-value of the Basic Local Alignment Search Tool (BLAST) parameter was set at 1−5. The amino acid sequence of candidate unigene-related proteins was investigated through Pfam database to gain annotated information 47 .
Quantitative real-time PCR validation. The selected DEGs were confirmed by (q-PCR). Total RNA from each sample was extracted as described earlier. Reverse transcription was carried out using the PrimeScript RT reagent kit with gDNA Eraser (TOYOBO). Quantitative real-time PCR was conducted on the CFX96 Real-Time System (Bio-Rad, Hercules, CA, USA) by SybrGreen method with Premix Ex Taq II (Takara, Kyoto, Japan). The qPCR cycling condition: 95 °C for 30 s, followed by 40 cycles at 95 °C for 30 s, and 60 °C for 35 s was used to confirm selected genes using specific primer and 40 s rRNA as control and efficiency values were detected at the range of 94-99% (Table S12). The qPCR product specificity was analyzed using the melting curve. The relative gene expression values were calculated by the 2-ΔΔCt method 83 .

Conclusions
LLP29 is a highly mosquitocidal Bti strain with independent intellectual property that harbors fewer plasmids and then in shorter life cycle 19,20 . In order to further know Bt mechanism and provide the theory base of mosquito biocontrol application, LLP29 was used as a novel Bti candidate, and the transcriptomic basis of midgut as well as the main functional genes of A. aegypti after Bt infection was searched by transcriptome Illumina resequencing in the present study. In the present study, the midgut of control and LLp29 toxin-treated A. aegypti mosquitoes was sequenced by Illumina sequencing. A total of 17130, including 16358 unigenes were functionally annotated by BLAST searches against various databases. Furthermore, 574 new unigenes, including 391 functionally annotated unigenes, were identified. A total of 557 DEGs, including 226 upregulated and 231 downregulated unigenes, were annotated by comparing the Bti toxin-treated and control groups. Furthermore, 442 DEGs were functionally annotated. Furthermore, specific DEGs related to multidrug resistance, immune response, stressed-related, detoxification, lipid/fatty acid, phosphatases, and other transporter proteins were identified. In addition, 28 randomly selected DEGs were further confirmed through qPCR. Majority of genes exhibited similar transcriptomic and qPCR data in terms of expression. The results of the present study are more valuable than the ones only immunity to Cry11Aa toxin in Bti and it will be useful for further research on interaction between Cry toxins and other proteins. Due to the presence of many putative proteins and these midgut proteins may compete with Cry toxin to bind with receptors and protect different receptors from toxin binding, thus altering the toxicity of Cry toxins. Further research is needed to identify the other midgut proteins that can interfere with interaction between Cry toxins and different receptors, and also to identify polysaccharides or chemicals that counterpart to neutralize the effect of these midgut proteins to increase the toxicity of Cry toxins for mosquito control.
Data accessibility. The transcriptome raw data has been submitted to the National Center for Biotechnology Information (NCBI). Read data can be accessed directly in Sequence Read Archive (SRA) with accession number SRP150934.