The Japanese pine sawyer, Monochamus alternatus Hope (Coleoptera: Cerambycidae), is a major pest of coniferous forests. Bacillus thuringiensis Berliner (Bacillales: Bacillaceae) (Bt) is an ecosystem-friendly agent for the control of pest insects. To identify the genes that are potentially associated with Bt toxicology, differentially expressed gene (DGE) libraries of M. alternatus after 20 h of exposure to a sublethal concentration of Bt were prepared from dissected midgut tissue. The transcripts were sequenced using the Illumina sequencing platform. After quality checks, the total numbers of clean reads were 29.20 and 41.53 million in the control and Bt-exposed libraries, respectively. We found that 947 and 8,895 unique genes were significantly up- and down-regulated, respectively. The annotation of the non-redundant sequences using the Gene Ontology database identified 3,837 sequences involved in cellular components, 3,190 sequences involved in molecular functions, and 6,930 sequences involved in biological processes, and showed that the genes were distributed across 50 categories. The Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis indicated that the top 20 enriched pathways included metabolism of xenobiotics by the cytochrome P450 pathway and steroid hormone biosynthesis. Seventy-four DEGs were enriched in the metabolism of xenobiotics by the cytochrome P450 pathway, including the DEGs encoding Cyp1A1 (cytochrome P450, family 1, subfamily A, polypeptide 1), Cyp1A2 (cytochrome P450, family 1, subfamily A, polypeptide 2), Cyp3A (cytochrome P450, family 3, subfamily A), aldehyde dehydrogenase (NAD[P] ) and ADH1_7 (alcohol dehydrogenase 1/7), all of which were up-regulated. However, the DEGs encoding cathepsin L and CBR1 (carbonyl reductase 1) were down-regulated. Therefore, the effects of Bt on M. alternatus are complex and time dependent. This study is a first step toward achieving an understanding of the profile of the secondary targets of Bt in M. alternatus, and the results may provide insights for further exploration of the functions of the target genes in detoxification and resistance to Bt.
The Japanese pine sawyer, Monochamus alternatus Hope (Coleoptera: Cerambycidae), is a major pest of coniferous forests, especially pines (Pinus species; Pinales: Pinaceae), and is also the key vector of the exotic pinewood nematode, Bursaphelenchus xylophilus (Steiner et Buhrer) (Nematoda: Aphelenchoididae) in eastern Asia (Kobayashi et al. 1984). Chemical insecticides have been used to control insects like the pestiferous longhorned beetles for many years. However, pesticide resistance, environmental pollution, and inaccessible larvae within the wood of trunks and branches have largely prevented successful long-horned beetle control in tree plantations (D'amico et al. 2004; Chen et al. 2005; Schloss et al. 2006). It is vital that the disruption of the ecosystem through the misuse of chemical pesticides be avoided; therefore, a safe, efficient, and eco-friendly pest control strategy is required (Srinivasan 2008). Consequently, the biopesticide Bacillus thuringiensis (Berliner) (Bacillales: Bacillaceae) (Bt) can be an alternative and more ecosystem-friendly option for the integrated management and efficient control of longhorned beetles and for the development of integrated pest management tools.
Bacillus thuringiensis is a Gram-positive bacterium that produces parasporal crystals composed of polypeptide protoxins during sporulation. To exert their toxic effect after ingestion by susceptible insects, these Cry proteins are proteolytically activated by gut proteases, and the active toxins bind to receptors located on the brush border midgut epithelium. Breakdown of the gut epithelium through a pore formation mechanism on the target membrane or an alternative cell death process involving the adenylyl cyclase/PKA signaling pathway have been proposed as responsible for Cry cytotoxic effects (Zhang et al. 2006; Bravo et al. 2007).
Cry toxins are toxic to susceptible insect pests yet comparatively harmless to plants and vertebrates, making their application a more environmentally safe alternative compared with conventional pesticides (Bravo et al. 2004; Chen et al. 2005; Vaughn et al. 2005; Christouo et al. 2006; Rausell et al. 2007). As these toxins affect only a narrow range of target species, leaving vertebrates and most beneficial invertebrates unaffected, they are currently being used worldwide as biopesticides formulated in commercial sprays or introduced in transgenic plants (Loseva et al. 2002; Shelton et al. 2002; Bravo & Soberon 2008; Rodrigo-Simon et al. 2008; Sanahuja et al. 2011).
The Cry3 class of toxins exhibits specific toxicity toward a number of coleopteran larvae, in both diet bioassays and field trials (Genissel et al. 2003; Chen et al. 2005; Rausell et al. 2007; Park et al. 2009). Bacillus thuringiensis Cry3Aa fused to a cellulase-binding peptide shows increased toxicity against longhorned beetles. Cry3Aa is the only Bt coleopteran active toxin that has been tested against a broad range of Coleopteran species (23 species, 60% of which were susceptible) (van Frankenhuyzen 2009).
Studying the short-term transcriptome responses of insects to xenobiotics could lead to the discovery of novel molecular mechanisms contributing to insecticide detoxification and tolerance (David et al. 2010). In the present study, we studied transcription profiles in M. alternatus upon oral infection with the entomopathogen B. thuringiensis. We compared gene expression changes between oral infection and control, using an Illumina next generation sequencing approach (RNA-Seq). RNA-Seq is a powerful tool that enables quantification of transcript levels on a genome-wide scale (Wang et al. 2009b). Detailed gene expression profile knowledge of B. thuringiensis interaction with M. alternatus should facilitate the development of more effective Bt based products against M. alternatus and other pest Coleoptera.
Materials and Methods
INSECTS
Monochamus alternatus larvae were originally obtained from dead, infested Pinus massoniana Lamb. (Pinales: Pinaceae) trees. The larvae were reared on artificial diets (Xu et al. 2009) in complete darkness at room temperature with a 55% relative humidity until they reached the 3rd instar larval.
INSECTICIDE TREATMENTS
Commercial biopesticide formulations based on B. thuringiensis var. tenebrionis purchased from Xirui Biological Technology Co., Ltd., Taiyuan, China, were used in the study. Based on previous studies, newly molted M. alternatus 3rd instar larvae were exposed to Bt insecticides orally at the dose of 10 billion spores/ mL. One mL Bt insecticide was compounded into 9 mL of the artificial diet. After incorporation of the insecticide into the diet, the larvae were transferred onto the treated diet in individual glass vials. Then, the glass vials were transferred to the growth chamber under the above-mentioned conditions. Only live larvae with no obvious signs of intoxication were selected for transcript abundance measurements. Distilled water was used as the control treatment agent. The mortality did not exceed 10% in the control group.
RNA EXTRACTION AND PURIFICATION
At the end of the 20 h Bt exposure, midguts of 3rd instar M. alternatus larvae fed control diet or diet contaminated with Bt were dissected, and extraneous tissue attached to midguts was carefully removed and pooled with tissues from the remainder of the body. Dissected midguts were placed in RNA/ater(ambion.com) overnight at 4 °C or snap-frozen in liquid nitrogen. Three replicates of 30 midguts each for both treatment and control were performed. The total RNA from the midgut samples was extracted using TRIzol reagent (Life Technologies, Carlsbad, California) following the manufacturer's instructions. The quality and concentration of the RNA samples were examined by denatured agarose gel electrophoresis and spectrophotometer analysis. The high-quality total RNA was further purified using the RNeasy mini kit (Qiagen GmBH, Germany) and RNase-Free DNase set (Qiagen GmBH, Germany).
CONSTRUCTION OF THE cDNA LIBRARY AND ILLUMINA SEQUENCING
The RNA of 90 midgut samples in 3 replicates used for sequencing was first treated with DNase I to degrade any possible DNA contamination. The mRNA was then enriched using oligo (dT) magnetic beads and was fragmented into short fragments. The first-strand cDNA was then synthesized using random hexamer-primer. Buffer, dNTPs, RNase H, and DNA polymerase I were added to synthesize the second strand. The double-strand cDNA was purified with magnetic beads. End reparation and 3C-end single nucleotide A (adenine) addition were then performed, and the sequencing adaptors were then ligated to the fragments. The fragments were enriched by PCR amplification. An Agilent 2100 Bioanalyzer and an ABI StepOnePlus Real-Time PCR system were used to qualify and quantify the library. The library products from 90 samples were sequenced using an IlluminaHiSeq™ 2000 instrument. In total, 41,529,187 paired reads most of which were 150 bp in length were obtained.
TRANSRIPTOME DE NOVO ASSEMBLY
The raw reads were filtered by removing reads with adaptors, reads with a percentage of unknown nucleotides greater than 5%, and the low-quality reads in which the percentage of low-quality bases (base quality <10) was greater than 20%. De novo transcriptome assembly was performed using the short reads assembler Trinity (Ding et al. 2003). Trinity combines 3 independent software modules: Inchworm, Chrysalis, and Butterfly, applied sequentially to process large volumes of RNA-Seq reads. Trinity partitions the sequence data into many individual de Bruijn graphs, each representing the transcriptional complexity at a given gene or locus, and then processes each graph independently to extract full-length splicing isoforms and to tease apart transcripts derived from paralogous genes. The sequences produced by Trinity are called unigenes. When multiple samples from the same species are sequenced, sequence splicing and redundancy removing can be further processed in the unigenes from the assembly of each sample to acquire the non-redundant ones with the sequence as long as possible.
SCREENING OF DIFFERENTIALLY EXPRESSED GENES (DEGs)
Referring to “The significance of digital gene expression profiles” (Audic & Claverie 1997), we developed a strict algorithm to identify differentially expressed genes (DEGs) between 2 samples. The total clean tag number of sample 1 is N1 and total clean tag number of sample 2 is N2 ; gene A holds x tags in sample 1 and y tags in sample 2. The probability of gene A being expressed equally in the 2 samples can be calculated using the following equation:
The P-value corresponds to the differential gene expression test. The False Discovery Rate (FDR) is a method for determining the threshold P-value in multiple tests. We assumed that we selected R differentially expressed genes and that S of the R genes really showed differential expression and the other V genes were false positives. If we decide that the error ratio Q = V / R must remain below a cutoff (e.g., 1%), we should preset the FDR to a number that is not greater than 0.01 (Benjamini & Yekutieli 2001). We used “FDR < 0.001” and the “absolute value of log2Ratio ≥ 1” as the threshold to judge the significance of the gene expression difference.
GENE ONTOLOGY ENRICHMENT ANALYSIS OF DEGs
Gene Ontology (GO) Enrichment Analysis provides all of the GO terms that are significantly enriched in DEGs and filters the DEGs that correspond to biological functions. This method first maps all of the DEGs to GO terms in the database ( http://www.geneontology.org/), calculating the gene numbers for each term, and then uses a hypergeometric test to determine the GO terms that are significantly enriched with DEGs. The calculated P-value is subjected to the Bonferroni correction, using the corrected P-value <0.01 as a threshold. The GO terms fulfilling this condition are defined as significantly enriched GO terms. This analysis is able to recognize the main biological functions of the DEGs (Gene Ontology Enrichment Analysis Software Toolkit GOEAST).
With the Nr annotation, we used the Blast2GO program (Conesa & Gotz 2008) to obtain the GO annotations of the DEGs. Blast2GO has been cited by other researchers more than 150 times and is a widely recognized GO annotation software. After obtaining the GO annotations for the DEGs, we used the WEGO software (Ye et al. 2006) for the GO functional classification of the DEGs and to understand the distribution of the gene functions of M. alternatus from the macro level.
KEGG PATHWAY ENRICHMENT ANALYSIS OF DEGs
The molecular pathway analysis was conducted by using Pathfinder (release 63.0) based on the Kyoto Encyclopedia of Genesand Genomes (KEGG) database ( http://www.genome.jp/kegg/). Pathway enrichment analysis identifies metabolic pathways or signal transduction pathways that are significantly enriched in DEGs. Pathways with a P-value <0.01 are considered significantly enriched.
VERIFICATION OF DEGS WITH RT-qPCR
We verified the DEG results by reverse transcriptase quantitative polymerase chain reaction (RT-qPCR) using the AACT method. Thirteen genes were randomly selected for RT-qPCR verification. Cytoplasmic actin of M. alternatus was used as the endogenous control and was normalized to be in equal quantity in each cDNA treatment. The primers for the RT-qPCR target genes were designed online using GenScript Real-Time PCR Primer Design ( https://www.genscript.com/ssl-bin/app/primer). The melting curves for each template-primer pair were examined for non-specific amplification. The RT-qPCR data were acquired on a LightCycler Real-Time PCR instrument using SYBR Premix ExTaq™ (TaKaRa, Japan) and universal thermocycler conditions according to the manufacturer's protocol (Roche Diagnostics).
Briefly, the PCR assay was conducted in a 20 µL final volume with a PCR mix containing 2 µL cDNA at a concentration of 4 ng/µL (i.e., a quantity of cDNA corresponding to 20 ng of total RNA), 10 µL SYBR Premix Ex Taq™, 0.4 µL PCR forward primer, 0.4 µL PCR reverse primer, and 7.2 µL dH20. The amplification scheme was 95 °C for 30 s, 95 °C for 5 s, and 60 °C for 20 s. The dissociation curve method was applied according to the manufacturer's protocol (60 °C to 95 °C) to ensure the presence of a single specific PCR product. The annealing temperature was 60°C, and RT-qPCR was performed using 40 cycles. For each cDNA, 3 RT-qPCRs were performed, and standard curves were generated. To minimize the positional effect of variations in the block temperature on the PCR results, the PCRs for the target and reference genes were placed on the cycler in a randomized block design. The threshold cycle (CT) and the relative expression levels were calculated using the Light-Cycler480 1.5 software (Roche Diagnostics). For this analysis, 3 technical replicates were performed for each of the 3 biological replicates, 30 midgut samples were contained in each biological replicate, and they were the same samples that were sequenced with RNAseq.
Results
QUALITY ASSESSMENT OF THE READS
We sequenced 2 libraries of M. alternatus: control (SRA submission number: SRP026138) and Bt-exposed (SRA submission number: SRP030482). The total numbers of clean reads were 29.20 and 41.53 million, respectively, after removing the low-quality reads. The clean reads comprised more than 98.8% of the raw reads in each library, indicating the high quality of the sequencing.
GENE EXPRESSION PROFILES AFTER BT EXPOSURE
To identify the DEGs after exposure to Bt, the differences in gene expression were analyzed by pair-wise comparisons of control and Bt-exposed M. alternatus. We found that 947 and 8,895 unique genes were significantly up- and down-regulated, respectively. The unigenes detected with at least 2-fold differences in the 2 libraries are shown in Fig. 1 (FDR < 0.001). The red and green dots represent transcripts that are found in higher and lower abundance, respectively, by more than 2-fold. The blue dots represent transcripts that differed by less than 2-fold between the 2 libraries, which were designated as “no difference in expression.”
GENE ONTOLOGY FUNCTIONAL CLASSIFICATION OF DEGs
We obtained the expression patterns of the DEGs annotated to their given GO terms. Several GO categories were enriched in each term. The annotation of the non-redundant sequences using the GO database yielded better results, identifying 3,837 sequences as cellular components, 3,190 sequences involved in molecular functions, and 6,930 sequences involved in biological process, and showed that the genes were distributed among 50 categories, including immune system process, developmental process, metabolic process, reproductive process, response to stimulus, synapse, antioxidant activity, enzyme regulator activity, and receptor activity. Cellular process and metabolic process were the 2 largest groups in the biological process category, containing 1,247 and 1,131 unigenes, respectively. Similarly, catalytic activity (1,485 unigenes) and binding terms (1,246 unigenes) were dominant in the molecular function category, and the most abundant GO cellular components were cell and cell part (970 unigenes). Furthermore, the most highly up- and down-regulated genes belonged to carboxylesterase activity and hydrolase activity, respectively.
Fig. 1.
Comparison of the gene expression levels between the control (CK) and Bt-exposed Monochamus alternatus. To compare the gene expression levels between the 2 libraries, each library was normalized to 1 million tags. The x-axis represents log10 of the reads per kb per million reads (RPKM) of the control sample, and the y-axis indicates log10 of the RPKM of the treated sample. The expression level of each gene is included in the volcano plot. The red dots represent transcripts that are more prevalent in the Bt-treated library, the green dots show those present at a lower frequency in the Bt-treated library, and the blue dots indicate transcripts that did not change significantly. The parameters “FDR < 0.001” and “absolute value of log2(Treated/Control) ≥ 1” were used as the thresholds to judge the significance of the gene expression difference.
KEGG PATHWAY ENRICHMENT ANALYSIS OF DEGs
The KEGG pathway enrichment analysis of DEGs showed that these targets were enriched in 44 pathways with P-values less than 0.01. We ranked those pathways according to the P-value of each enriched pathway in ascending order and selected the top 20 pathways (Table 1) for further mechanistic analyses. Among the top 20 pathways enriched, we found that the metabolism of xenobiotics by the cytochrome P450 pathway ranked 2nd and that steroid hormone biosynthesis ranked 9th.
CONFIRMATION OF THE GENE EXPRESSION RESULTS BY RT-qPCR
To validate the Illumina expression profiles, we randomly selected 13 genes that showed differential expression for further analysis by RT-qPCR. Actin, which was constitutively expressed in our previous work, was chosen as a reference gene for data normalization. The graph in Fig. 2 demonstrates a regression in which the x-axis is Illumina fold change and the y-axis is RT-qPCR results. The regression equation is y = 0.231X + 0.34. The correlation coefficient and P-value for the regression are 0.809 and 0, respectively.
Table 1.
Top 20 pathways of enriched DEGs with P-values less than 0.01.
Discussion
Bt strains are usually used for controlling coleopteran larvae but not adults (Qi et al. 2011).The mode of action of Bt toxins generally involves multiple steps. However, the signaling pathway model has recently been challenged due to lack of supporting experimental data (Vachon et al. 2012). Although the exact mode of action of Bt toxins is not completely understood, it is clear that a number of genes expressed in the insect gut are involved in Bt toxicity (Soberón et al. 2010; Vachon et al. 2012). In this communication, we introduced preliminary results from our analysis of the differentially expressed genes in surviving M. alternatus larvae after exposure to Bt at a sublethal level compared with unexposed larvae. We found that the effects of Bt on M. alternatus at 20 h post-exposure are complex. The analysis of the host responses to Bt infection may contribute to an enhanced understanding of the mechanistic underpinnings of the degree of toxicological response and immunological adaptation to Bt.
Genes usually interact with each other to play roles in certain biological functions. Pathway-based analysis helps to further understand the biological functions of genes. KEGG (Kyoto Encyclopedia of Genes and Genomes) is the major public pathway-related database (Kanehisa et al. 2008). In addition, the pathways associated with these unique candidate targets yielded new insights that will lead to a better understanding of their functions and relations.
Further research showed that there were 74 DEGs enriched in the metabolism of xenobiotics by the cytochrome P450 pathway, including the DEGs encoding Cyp1A1 (cytochrome P450, family 1, subfamily A, polypeptide 1), Cyp1A2 (cytochrome P450, family 1, subfamily A, polypeptide 2), Cyp3A (cytochrome P450, family 3, subfamily A), aldehyde dehydrogenase (NAD[P]+) and ADH1_7 (alcohol dehydrogenase 1/7), all of which were up-regulated. However, the DEGs encoding cathepsin L and CBR1 (carbonyl reductase 1) were down-regulated.
Cytochrome P450s (CYPs) are a diverse class of enzymes with versatile functions including metabolic detoxification and resistance (Scott & Wen 2001). An enrichment among up-regulated genes including CYP genes was observed in the recently collected population of Tribolium castaneum (Herbst) (Coleoptera: Tenebrionidae) exposed to Bt (Behrens et al. 2014). We also found that several CYPs were induced. The observed up-regulation of stress and CYP genes may thus function to alleviate immunity-related self-damage.
In many invertebrate groups, cathepsin L proteases have been identified as major components of the gut digestive enzymes. It has also been demonstrated that they participate in other functions, such as immunological processes and tissue remodeling during insect metamorphosis (Matsumoto et al. 1995; Tryselius & Hultmark 1997; Hashmi et al. 2002; Baum et al. 2007; Wang et al. 2009a). The inhibition of cathepsin L protein activity by a protease inhibitor decreases the development and reproduction rate in Myzus persicae (Sulzer) (Aphidomorpha: Aphididae) (Cristofoletti et al. 2003). We found that the cathepsin L gene of M. alternatus was down-regulated. Depletion of intestinal proteases is considered able to induce (and deregulate) the synthesis of other proteases, which could result in toxic effects altering insect development (Sapountzis et al. 2014). We inferred that Bt may act as a toxin in M. alternatus by disrupting its development. The cathepsin L gene could be knocked out by RNA interference to observe if the development of M. alternatus would be delayed.
Another pathway that was closely connected with Bt is steroid hormone biosynthesis. There were 53 DEGs enriched in this pathway, and these included the up-regulated DEGs encoding UDP glucuronosyltransferase (UGT) and the down-regulated DEGs encoding 17β-hydroxysteroid dehydrogenases, 3-oxo-5-beta-steroid 4-dehydrogenase, 3-oxo-5-alpha-steroid 4-dehydrogenase 1, and estronesulfotransferase. These genes could possibly be related to ecdysonedependent tissue repair.
UDP glucuronosyltransferases are responsible for the process of glucuronidation, a major part of phase II metabolism. Glucuronosyltransferases are a very essential phase II detoxifying enzyme system involved in the metabolism of chemicals (Bock 2006). The reaction catalyzed by the UGT enzyme involves the addition of a glucuronic acid moiety to xenobiotics and a major pathway for foreign chemical removal for most drugs, toxins, and endogenous substances. The level of glucuronosyltransferase was up-regulated in this study. Increased levels of detoxification enzymes such as glucuronosyltransferase are the likely one of the mechanisms for chemoprevention of Bt. UGT upregulation can also be a bio-marker of pathogen- or immunity-related oxidative stress.
The 17β-hydroxysteroid dehydrogenases (17β-HSDs), among the most important enzymes involved in both the endocrine and the immune system, constitute a class of enzymes involved in the last step of steroid synthesis (Labrie et al. 1997). They also may be an indicator of ecdysone-dependent tissue repair or remodeling. Although a number of 17β-HSDs have been identified and characterized in vertebrates, the information from invertebrates is scarce. Most steps of the steroidogenic and steroid metabolic pathways described for vertebrates have been demonstrated to occur in invertebrates (Janer & Porte 2007). Some endocrine disruptor chemicals induce changes in these metabolic pathways and might lead to changes in steroid levels. After exposure to endocrine-disrupting chemicals (bisphenol A or 2,2′,4,4′-tetrabromodiphenyl ether), the expression of 17β-HSD transcripts was down-regulated in digestive glands in Mytilus galloprovincialis Lamarck (Mytiloida: Mytilidae) (Zhang et al. 2014). Similarly, in the present analysis, the mRNA expression level of 17β-HSD decreased significantly at Bt exposure. The alteration of gene transcription of 17β-HSDs after exposure to Bt may hint to the function of steroid metabolism and steroidogenic process in M. alternatus.
Overall, studying the short-term transcriptome responses of insects to xenobiotics could lead to the discovery of novel molecular mechanisms contributing to insecticide detoxification and tolerance. The transcriptional profile of M. alternatus' response to a sublethal insecticide dose not only provides novel insights into the molecular mechanisms of Bt toxicity, but also identifies potential biomarkers that may be explored in future studies to assess the insecticide's secondary effects and improve its efficacy and safety. Although much of this analytical work is descriptive, the output is rich in identifying many new gene products that may be associated with fundamental aspects of M. alternatus in toxicology and biology. Analysis of genes that demonstrated differential expression in our studies will assist in a further understanding of networks linking Bt and its potential targets. Identifying the molecular mechanisms of Bt toxicity is necessary for its future use for controlling M. alternatus and other coleopteran insects. We plan to evaluate the hypothetical function of several specific genes highlighted in the current study.
Acknowledgments
This research was supported by the National Natural Science Foundation of China (Grant No. 31470653).
