Gene Expression Profile of Bombyx mori Hemocyte under the Stress of Destruxin A

Abstract

Destruxin A (DA) is a cyclo-peptidic mycotoxin from the entomopathogenic fungus Metarhizium anisopliae. To uncover potential genes associated with its molecular mechanisms, a digital gene expression (DGE) profiling analysis was used to compare differentially expressed genes in the hemocytes of silkworm larvae treated with DA. Ten DGE libraries were constructed, sequenced, and assembled, and the unigenes with least 2.0-fold difference were further analyzed. The numbers of up-regulated genes were 10, 20, 18, 74 and 8, as well as the numbers of down-regulated genes were 0, 1, 8, 13 and 3 at 1, 4, 8, 12 and 24 h post treatment, respectively. Totally, the expression of 132 genes were significantly changed, among them, 1, 3 and 12 genes were continually up-regulated at 4, 3 and 2 different time points, respectively, while 1 gene was either up or down-regulated continually at 2 different time points. Furthermore, 68 genes were assigned to one or multiple gene ontology (GO) terms and 89 genes were assigned to specific Kyoto Encyclopedia of Genes and Genomes (KEGG) Orthology. In-depth analysis identified that these genes putatively involved in insecticide resistance, cell apoptosis, and innate immune defense. Finally, twenty differentially expressed genes were randomly chosen and validated by quantitative real-time PCR (qRT-PCR). Our studies provide insights into the toxic effect of this microbial insecticide on silkworm's hemocytes, and are helpful to better understanding of the molecular mechanisms of DA as a biological insecticide.

About the Authors:

Liang Gong

Contributed equally to this work with: Liang Gong, Xiurun Chen

Affiliation: College of Natural Resource and Environment, South China Agricultural University, Guangzhou, China

Xiurun Chen

Contributed equally to this work with: Liang Gong, Xiurun Chen

Affiliation: College of Natural Resource and Environment, South China Agricultural University, Guangzhou, China

Chenglan Liu

Affiliation: College of Natural Resource and Environment, South China Agricultural University, Guangzhou, China

Fengliang Jin

Affiliation: College of Natural Resource and Environment, South China Agricultural University, Guangzhou, China

Qiongbo Hu

* E-mail: [email protected]

Affiliation: College of Natural Resource and Environment, South China Agricultural University, Guangzhou, China

Introduction

Entomopathogenic fungi, such as Metarhizium anisopliae and Beauveria bassiana [1], [2], are very important natural factors for insect control. As a result, fungal insecticides are very attractive in pests IPM (Integrated Pest Management), particularly in areas of high chemical insecticide resistance [3], [4]. The entomopathogenic fungus M. anisopliae has been well studied, and commercially used to control termites [5], grasshoppers [6], thrips [7], whiteflies [8] and red spider mites [9]. Notably, Destruxins were firstly isolated from M. anisopliae, and they are important virulence factors which accelerate the death of infected insects [4], [10]–[11]. However, the molecular mechanisms of destruxins as insecticide have not been elucidated yet.

Destruxins are fungal secondary metabolites and bio-synthesized by non-ribosomal peptide synthases [12]. Chemically, destruxins are cyclic hexadepsipeptides composed of an α-hydroxy acid and five amino acid residues. So far, 39 destruxins analogues have been isolated from M. anisopliae and other fungal species [13]–[15]. The common analogues, Destruxin A (DA), Destruxins B (DB) and Destruxins (DE) exhibit a substantial insecticidal activities against many species of pests such as Plutella xylostella, Spodoptera litura, Manduca sexta and Pieris brassicae [16]–[18], etc.

There is limited information regarding the mechanisms of action of destruxins. DA suppresses the contractions of visceral muscles of Locusta migratoria, with an influx of extracellular Ca²⁺ [19]. Additionally, DA strongly inhibits the secretion rate of Rhodnius prolixus Malpighian tubules fluid without infecting on the production of intracellular ATP [18]. Other studies demonstrated that DB selectively inhibited the activity of V-type ATPase from different cells [20]–[21]. Some experiments have shown that destruxins are able to damage the innate immunity of insects. The function of the encapsulation and phagocytosis processes of insect hemocytes were found to be destroyed by destruxins [22]. In our previously survey, it was also found that DA could induce obviously morphologic alterations of silkworm's hemocytes in vivo, even at an extremely low dose [23]. Furthermore, another study showed that destruxinss may play a key role in suppressing the innate immune response of Drosophila melanogaster, by inhibiting the expression of antimicrobial peptides [24]. However, more studies are needed for better understanding of the host immunity response modulated by destruxins, which may be an important aspect of the mechanisms of destruxins acting on insects.

Over the past several years, next-generation sequencing (NGS) techniques, such as Illumina/Solexa, have been developed as high-throughput and low-cost sequencing platforms to investigate the de novo transcriptome, gene expression profiling, and detecting methylation patterns [25], [26]. In the present study, a digital gene expression (DGE) profiling analysis using Illumina sequencing technology was used to compare differentially expressed genes in the hemocytes of silkworm larvae with DA or mock treatment. Ten DGE libraries were constructed, sequenced, and assembled, and the unigenes with least 2.0-fold difference in expression were further assembled and analyzed. Finally, a number of the differentially expressed genes were confirmed by quantitative real-time PCR (qRT-PCR).

Results

Illumina sequencing and reads assembly

Ten DGE libraries of B. mori were sequenced including the DA-treated and control samples, which generated a number of row reads ranged from 7,041,039 to 7,652,389 for each of them. Following filtering the low quality reads, the total number of clean reads per library ranged from 7,007,499 to 7,607,918 million, and the percentage of clean reads in each library ranged from 99.31% to 99.56% (Figure 1). The alignment with reference transcriptome and reference genome showed that unique match ranged from 52.68% to 58.32%, and from 72.59% to 73.08%, respectively, which is the most critical subset of DGE libraries to identify a transcript precisely (Table 1). To evaluate if the number of detected genes increasing proportionally to total tags number, we performed the sequencing saturation analysis for each sample. The results showed that the number of detected genes was increasing as the number of reads was increasing. However, when the number of reads reached to 7.5 million, the growth rate of detected genes flattened, indicating that the number of detected genes tended to be saturated (Figure S1). We used the distribution of reads on the reference genes to evaluate the randomness. If the randomness is poor, reads preference to specific gene region will directly affect subsequent bioinformatics analysis. But our data (Figure S2) showed a even distribution with the dynamic range more than 9.560, which is a required ratio between the maximum and minimum expression level for RNA-Seq. To assess comparability of DGE data, we analyzed of the distributions of genes' coverage. The results are similar between the control and treated samples, indicating it is comparability (Figure S3).

[Figure omitted. See PDF.]

Figure 1. Quality assessment of reads in each library.

https://doi.org/10.1371/journal.pone.0096170.g001

[Figure omitted. See PDF.]

Table 1. Quality assessment of reads in each library.

https://doi.org/10.1371/journal.pone.0096170.t001

Functional annotation of differentially expressed genes

To check whether DA could result in significantly changes of gene expression in B.mori hemocytes, we identified and compared the differentially expressed genes between the DA-treated and control samples (Figure S4), which were further calculated by normalizing the number of unambiguous tags in each library to reads per kb per million reads (RPKM). The results revealed that many genes were significantly differentially expressed between the control and treated libraries. Furthermore, all the genes with the expression more than 2-fold changes were annotated by using Nr database, GO database and KEGG pathway database (Table S1). Among them, we found that the toxicity response of B.mori to DA was mainly associated with the insect innate immune response, xenobiotic detoxification and apoptosis. Interestingly, the profile of genes number affected by DA showed a curve type with a peak at the treatment of past 12 h (Figure 2).

[Figure omitted. See PDF.]

Figure 2. Screening of more than 2-fold differentially expressed genes in B.mori hemocytes after the treatment of DA in differential time intervals.

https://doi.org/10.1371/journal.pone.0096170.g002

At 1 h post-treatment, 10 genes were up-regulated with at least 2-fold changes (Figure 2). Among them, 5 genes of proteases and aminopeptidase (Bm_nscaf2674_066, Bm_nscaf2674_064, Bm_nscaf2674_063, Bm_nscaf2983_049 and Bm_nscaf2889_046) and 2 putative cuticle protein (Bm_nscaf2838_045 and Bm_nscaf2767_133) and the other 3 genes were annotated. It indicated that expression of digest-related genes was mainly changed in hemocytes in their early response to DA treatment. It is in accordance with the Gene Ontology annotation (Figure 3), which showed that genes involving in metabolic process and catalytic activity were much up-regulated at 1 h post treatment. Several signal pathways such as fatty acid metabolism and focal adhesion were also activated (Table S2).

[Figure omitted. See PDF.]

Figure 3. Summarizing of Gene Ontology functional classification of ten DGEs by comparing DA-treated and un-treated samples in differential time intervals, which described gene products in terms of their associated biological processes, cellular components and molecular functions.

https://doi.org/10.1371/journal.pone.0096170.g003

At 4 h post treatments, 20 and 1 genes were up-regulated and down-regulated, respectively (Figure 2). The highest up-regulated gene was Bm_nscaf2767_133 (9.4-fold), which was the only detected up-regulated gene at 4 time points (1 h, 4 h, 8 h and 12 h) and annotated as a putative cuticle protein. Whereas, Bm_nscaf2767_074 gene was the only detected down-regulated gene (-2.0-fold) and annotated as a transposon Ty3-G gap-Pol polyprotein in Clonorchis sinensis. Multiple signal transduction pathways were also activated at 4 h post treatment, such as MAPK (Mitogen-activated protein kinases) signaling pathway. Two genes (Bm_nscaf2136_098 and Bm_nscaf2888_274) which participate in MAPK pathway were annotated by KEGG (Table S2).

At 8 h post treatments, there were 18 up-regulated genes and 8 down-regulated genes (Figure 2). Bm_nscaf3058_128 was the highest up-regulated gene with 12.8-fold, which was annotated as peptidoglycan recognition protein (Table S1). Bm_nscaf2589_181 was down-regulated by -2.1-fold and annotated as E3 ubiquitin-protein ligase RNF31. These two genes play a key role in regulating immune and apoptopic signaling [27]. The other genes functions were annotated as well (Figure 2). Several signal pathways such as pancreatic secretion were also activated (Table S2).

At 12 h and 24 h post treatments, 74 and 8 genes were respectively up-regulated, as well as 13 and 3 genes were respectively down-regulated (Figure 2). Bm_nscaf2136_210 and Bm_nscaf2818_080 were the highest induced expression genes at both 12 h and 24 h post treatments, and were annotated as tropomyosin-2 isoform 2 and protease inhibitor 6, respectively. Lots of other genes were mainly found participating in developmental, multicellular organismal and signalling processes (Figure 3). ABC transporters and many other signal pathways were probably activated (Table S2 6).

Gene expression analysis and qRT-PCR validation

In order to confirm the integrity of the DGE data, qRT-PCR was carried out. Twenty genes were randomly selected including alkaliphilic serine protease, two putative cuticle protein CPH45, tubulointerstitial nephritis antigen precursor, serine protease precursor, Bmp-2 protein, fungal protease inhibitor F precursor and putative serotonin receptor and so on. The PCR products were sequenced, and further confirmed by blasting in NCBI database. The quantification data showed that the eighteen genes expression profile were consistent with the DGE results, although the ratios varied to some extent, but two genes (Bm_nscaf2838_045 and Fungal protease inhibitor F precursor) are not significantly expressed in the checked samples, although both revealed the same expression tendency as DGE results (Table 2).

[Figure omitted. See PDF.]

Table 2. Quantitative Real Time PCR (qRT-PCR) validation of DGE result.

https://doi.org/10.1371/journal.pone.0096170.t002

Discussion

Summary of the DGE results

Our DGE data firstly provided a comprehensive and global genes expression profile of B. mori hemocytes at different time points after DA treatment. The expression results with at least 2-fold changes showed that the numbers of up-regulated genes are 10, 20, 18, 74 and 8, as well as the numbers of down-regulated genes are 0, 1, 8, 13 and 3 after the treatment of 1 h, 4 h, 8 h, 12 h and 24 h, respectively. To sum up, the expression of 132 genes were significantly changed (in which 1, 2 and 12 genes were continually up-regulated at 4, 3 and 2 different time points, respectively, while 1 gene was either up or down-regulated continually at 2 different time points). Among them, 68 genes were assigned to one or multiple gene ontology (GO) terms and 89 genes were assigned to specific Kyoto Encyclopedia of Genes and Genomes (KEGG) Orthology. To sum up, 108 genes increased expression, while 25 genes decreased. The genes with significantly changed expression in this study were much less than a similarly study in P. xylostella, in which 1584 genes were observed with at least 2-fold changes, but comparing both, our data are only from hemocytes of B. mori, while the whole insect body of the fourth instar larvae were utilized in the DGE analysis of P. xylostella, which should be the reason why much more genes have been detected in P. xylostella [28].

Genes putatively involved in insect immune system

The immune system is generally divided into innate and adaptive immunity. Insect has only innate immunity system divided into humoral and cellular responses [29], [30]. The humoral response of innate immunity mainly includes three steps [31], [32], [33]: 1) identification of pathogen-associated molecular patterns (PAMPs) on pathogens by pattern recognition receptors (PRRs); 2) activation of the regulatory pathways; and 3) production of immune effectors to modulate cellular phagocytosis and molecular effectors such as antimicrobial peptides (AMPs). Peptidoglycan recognition proteins (PGRPs) are pattern recognition molecules that recognize bacteria and their unique cell wall component, peptidoglycan (PGN) [34]. In our study, four PGRPs (Bm_nscaf3058_131, Bm_nscaf3058_127, Bm_nscaf3058_130, Bm_nscaf3058_128) were significantly up-regulated. Specially, Bm_nscaf3058_128 was increased the expression of 12.78-fold after the treatment 8 h, and Bm_nscaf3058_127 was persistently up-regulated past the treatment of 8 h and 12 h with the changes of 2.06-fold and 2.93-fold, respectively (Table S1), indicating that B. mori innate immunity response to DA was activated through PGRP at a very early stage. DA also accelerated the expression of PGRPs of larvae of P. xylostella [28]. It suggests that PGRPs may be a common PRR that regulate the immunity of insects response to DA.

Insect PGRPs can activate the Toll signal pathway inducing production of antibacterial peptides (AMPs) such as cecropin and gloverin [35], [36]. In our study, B. mori Toll receptor was induced the expression of 2.39-fold after DA treatment of 8 h (Table S1). In consistence with our data, Toll was also up-regulated with 2-fold in P. xylostella larvae after injecting DA [28]. However, the scavenger receptor and C-type lectin were also up-regulated simultaneously in P. xylostella [28], while in our study, the expression changes of scavenger receptor and C-type lectin were not found. That may be caused by activating different immune cell responses to DA in different insect species or tissues, For example, in M. sexta hemocytes, Toll mRNA was significantly induced by Escherichia coli, Saccharomyces cerevisiae and Micrococcus lysodeikticus, but its transcript in the fat body was not induced by these microorganisms [37].

There were no significant changes of AMPs in our study, which may be due to either the threshold value of 2-fold is too high to find the difference, or indicating the evolutionary plasticity of B. mori by presenting novel proteins correlated with the response to DA [38]. But Pal et al (2007) found that DA could mediate a specific down-regulation of AMPs in Drosophila melanogaster, such as cecropins, attacin, metchnikowan, and diptericin [24], which led to a conclusion that DA has the potentiality to suppress the humoral immune response in D. melanogaster.

On the other hand, hemocytes are the primary mediators of cell-mediated immunity in insects including phagocytosis, nodulation, encapsulation, and melanization [39]. Phagocytosis is an ancient cellular process that plays an important role in host defense. Bm_nscaf2529_073 is homologue to D. melanogaster CG2765 acting as mediators of bacterial phagocytosis in cell line S2 [40]. In our study, this gene was significantly down-regulated. It suggests that gene Bm_nscaf2529_073 has an essential role of phagocytosis during the process of immunity of B. mori response to DA.

Genes putatively involved in insecticide detoxification and outside stress environment

The cytochrome P450 is a large and diverse group of enzymes that catalyze the oxidation of organic substances, such as xenobiotic [41]. Two cytochrome P450 genes, Bm_nscaf3005_56 and Bm_nscaf2827_03, were up-regulated with 4.2-fold and 3.1-fold after the treatment of 8 h, respectively (Table S1). However, in the case of P. xylostella, the expression of cytochrome P450 was not found to be changed, while Glutathione S-transferase gene was up-regulated after treated with DA[28].

Insect heat shock proteins (Hsp) are important molecular chaperone which were the first introduced by Tissieres et al. (1974) [42] due to its increased expression to the high temperature in D. melanogaster, and it also has been confirmed that Hsp involved in multiple physiological roles such as to increase lifespan, enhance stress resistance, and prevent apoptosis and neurodegenerative diseases [43], [44]. In our study, the expression of one 19.5 kDa Hsp and three 70 kDa Hsps were significantly decreased, but one 20.4 kDa Hsp was significantly increased, which provide a complex physiology role of this family protein in the process of B. mori response to DA.

Genes involved in early and persistently response to DA

In our experiment, numerous data of gene expression at five different time points (1 h, 4 h, 8 h, 12 and 24 h) after DA treatment were accomplished. However, we think that the genes involved in the early and continually response to DA should be paid more attention. B. mori cuticle proteins play essential roles in many physiological functions, during molting and metamorphosis [45]. Interestingly, two cuticle protein genes were found continual up-regulation in this research. Bm_nscaf2767_133 gene continually increased the expression with 4.62-fold, 9.36-fold, 9.69-fold and 4.7-fold in post-treatment 1, 4, 8 and 12 h, respectively, while Bm_nscaf2838_045 gene were up-regulated 6.8-fold in 1 h, 11.0-fold in 8 h and 4.7-fold in 12 h (Table S1). In consistence with our data, Zhang et al. (2010) [46] reported that the putative cuticle protein genes in adult of Leptinotarsa decemlineata could be highly induced by either insecticide azinphosmethyl (a type of organophosphate insecticide) 2–3 weeks after moulting, or dry environmental conditions. Consequently, we believe that insect cuticle proteins perhaps have other unknown functions.

Four trypsin genes (Bm_nscaf2674_064, Bm_nscaf2674_063, Bm_nscaf2983_049 and Bm_nscaf2674_066) were all significantly up-regulated at early or persistent time points. Trypin is a common serine protease involved in protein precursor cleavage used for signal transduction and cascade amplification, and can also activate specific defense mechanisms, such as complement activation, melanization, blood coagulation and antibacterial peptide synthesis [47]. So, we can hypothesize that DA stimulates silkworm hemocytes expressing immediately trypsins to start some signal transductions.

In conclusion

The most important finding in this study gives an overview of genes expression profile of B. mori hemocytes after the treatment of DA during different intervals, which can help to better understanding the immune response of B. mori hemocytes underling the stress of DA. But the results also highlight the physiological role of some proteins of B. mori response to DA, such as heat shock protein and cuticle protein. Taken them together, our findings are helpful to elucidate the molecular mechanisms of DA as a biologic insecticide.

Methods

Insects

The strain of Bombyx mori, P50, was provided by Dr. Ye Mingqiang (Sericulture&Agri-food Research Institute, Guangdong Academy Agricultural Science). The one-day old fifth instar larvae of silkworms were collected and used in this study.

Destruxin A and chemicals

Destruxins A (DA) was isolated and purified from strain MaQ10 of Metarhizium anisopliae as described [48], [49]. DA was identified and determined with purification rate of 95.7% by means of high performance liquid chromatography (HPLC) with a standard sample from Sigma-Aldrich Co. LLC (St. Louis, MO, USA) [23]. DA was diluted to 10000 µg/mL using dimethyl sulfoxide (DMSO, Sigma-Aldrich Co. LLC) as the stock solutions and stored at −20°C for further use.

Injection of DA into Silkworm Larvae

Firstly, DA stock was diluted with double distilled water to be 200 µg/mL. The injection dose of DA to each 5th instar larva was 10 µL. Before injection, the larvae of Bombyx mori were placed on ice for 5 min to paralyze. Before injection, the insect's surface were carefully washed with pure water, then, disinfected with a small amount of alcohol. Then, the needle tip was inserted in the soft part of prolegs of paralytic silkworm larvae and 10 µL DA working solution was carefully injected into silkworm larvae. Then, the injected larvae were maintained at 26°C. The control group was treated with 1% DMSO aqueous solution.

Collection of hemocytes

Hemocytes were collected 1, 4, 8, 12 and 24 h after injection, samples at each time point were dissected from 5 fifth instar larvae of B. mori. Before bleeding, the insect's surface were carefully washed with pure water, then, disinfected with a small amount of alcohol. When bleeding, the insect was placed on ice for a few minutes to paralyze, then, its hind leg was cut with a surgical scissors, and the blood was dripped into a centrifuge tube containing 650 µL anticoagulation buffer (KCl 69 mmol/L, NaCl 27 mmol/L, NaHCO₃ 2 mmol/L, dextrose 100 mmol/L, potassium citrate 30 mmol/L, citrate acid 26 mmol/L, Na₂-EDTA 10 mmol/L respectively, and pH 4.6, 420 mosm.) on ice. Then, the mixture of blood and buffer was centrifuged at the speed of 1000 r/min for 10 min at 4°C and the supernatant was abandoned. One ml RNAiso Plus (Trizol, Takara biotechnology (Dalian) CO., LTD) was added to the hemocytes and stored at −80°C for further use.

Gene expression profile sequencing

All collected hemocyte samples were used for RNA isolation. Total RNA from these samples was extracted using Trizol Total RNA Isolation Kit (Takara, Japan) according to manufacturer's protocol, followed by quality and quantity analysis using Nanodrop (Bio-Rad, USA) and 2100 Bioanalyzer (Agilent, USA). Digital gene expression (DGE) libraries were prepared by the Illumina gene expression sample prep kit (Illumina, San Diego, CA). Six micrograms of total RNA was used to isolate mRNAwith Oligo(dT) magnetic beads adsorption. The first and second cDNAs were synthesized and the bead-bound cDNAs were subsequently digested by Nla III, which recognizes and cuts on the CATG sites. The digested 3′cDNA fragments were purified with magnetic Beads, and connected to the Illumina adapter 1 at the sticky 5′-ends. The junction of adapter 1 and CATG site is the recognition site of Mme I, which cuts off the cDNA at a position of 17 bp downstream of the CATG site, thereby producing tags with adapter 1. After removing 3′-end fragments with magnetic beads precipitation, Illumina adapter 2 was ligated to the tags at the 3′-ends, generating tags with different adapters at both ends to form a tag library. The fragments are enriched by PCR amplification, 95 bp bands were purified by 6% PAGE gel electrophoresis. After denaturation, the single-stranded molecules were fixed onto the Illumina sequencing chip (flow cell) for sequencing by using the sequencing by synthesis (SBS) method. Each tunnel generated millions of raw 49 bp reads. Each time point provides two samples for sequencing as two biological replicates. The data sets are deposited at the NCBI Short Read Archive (http://www.ncbi.nlm.nih.gov/sra/) with an accession number for the treatment at each time point (Table 3).

[Figure omitted. See PDF.]

Table 3. NCBI SRA accession numbers for the treatments at each time point.

https://doi.org/10.1371/journal.pone.0096170.t003

Bioinformatics analysis of digital gene expression (DGE) tags

Before mapping the tags to the transcriptome database, the raw sequence data were filtered by removing low quality tags, such as tags with unknown sequences ‘N’, empty tags (reads with only adapter sequences but no tags), low complexity tags, and tags with a copy number of 1 (likely sequencing errors). A reference library containing all of the sequences of CATG plus 17 bases was created by searching the CATG sites in the transcriptome database. All clean tags were mapped to the reference library and allowed no more than one base mismatch. Clean tags, which were mapped to exactly one gene in the reference database, were designated as unambiguous tags for gene annotation. The number of unambiguous tags for each gene was calculated and normalized to RPKM (Reads Per Kb per Million reads) for the gene expression analysis [50], [51].

Screening and statistics analysis of differentially expressed genes

A rigorous algorithm was developed to identify differentially expressed genes (DEGs) in both treatment and control conditions. The total clean tag number of the sample 1 is N1, and total clean tag number of sample 2 is N2; gene A holds x tags in sample1 and y tags in sample2. The probability of gene expressed equally between two samples can be calculated with:

The false discovery rate (FDR) method was used to determine the P-value threshold in multiple testing [52]. FDR≤0.001 and the absolute value of log2^ratio ≥2 were used as the threshold to judge the significance of gene expression differences. For Gene Ontology (GO) enrichment analysis, we used the hypergeometric test to map all the differentially expressed genes to terms in the GO database, identifying GO terms significantly enriched for DEGs, and comparing them to the genome background.where N is the number of all genes with GO annotation; n is the number of DEGs in N; M is the number of all genes that are annotated to the certain GO terms; m is the number of DEGs in M. The calculated p-value goes through Bonferroni Correction, taking corrected p-value ≤0.05 as a threshold. GO terms fulfilling this condition are defined as significantly enriched GO terms in DEGs. This analysis is able to recognize the main biological functions that DEGs exercise.

The differentially expressed genes were also utilized in Kyoto Encyclopedia of Genes and Genomes (KEGG) ontology enrichment analyses to further understand their biological functions. KEGG Pathway analysis can identify significantly enriched metabolic pathways or signal transduction pathways in DEGs compared with the genome background. Pathways with Q-value≤0.05 were considered as significantly enriched in DEGs.

Validation of gene expression by quantitative real-time PCR

Quantitative real-time reverse transcription PCR (qRT-PCR) was performed to validate the Illumina sequencing data. Eight genes were randomly chosen for qRT-PCR analysis, which have the significantly differential expression in DGE data. The primers are shown in Table 2. Total RNA sample was extracted for the DGE experiment with three biological replicates. The RT-PCR was performed using the AMV RNA Kit (TaKaRa, Japan) according to the manufacturer's protocol. qRT-PCR reactions were performed in triplicate on a BioRad iQ5 real-time PCR detection system using 10 ng of cDNA, 0.2 µM of primers and SYBR Premix Ex TaqTM (TaKaRa) according to the manufacturer's protocol. Bombyx mori β-actin was used as a reference gene to standardize the level of other transcripts. The relative amounts of the transcripts were first normalized to the endogenous reference gene and then normalized to the gene expression level in the un-treated samples according to the 2^−ΔΔCt method [53], statistical analysis was performed using Sigma Plot 12.0 software based on t-test with p<0.05 representing significance.

Supporting Information

[Figure omitted. See PDF.]

Figure S1.

Sequence saturation analysis of each library.

https://doi.org/10.1371/journal.pone.0096170.s001

(TIF)

Figure S2.

Randomness assessment of each library.

https://doi.org/10.1371/journal.pone.0096170.s002

(TIF)

Figure S3.

Distribution of genes' coverage in each library.

https://doi.org/10.1371/journal.pone.0096170.s003

(TIF)

Figure S4.

Scattered plot of differential expression genes in DA-treated and un-treated samples in differential time intervals. RPKM means Reads Per Kb per Million read.

https://doi.org/10.1371/journal.pone.0096170.s004

(TIF)

Table S1.

Genes with more than 2-fold expression changes between the DA-treated and control samples are annotated by using Nr database, GO database and KEGG pathway.

https://doi.org/10.1371/journal.pone.0096170.s005

(DOC)

Table S2.

Details of KEGG pathway enrichment analysis of the genes with significant expression change.

https://doi.org/10.1371/journal.pone.0096170.s006

(DOC)

Author Contributions

Conceived and designed the experiments: QBH FLJ CLL. Performed the experiments: XRC LG. Analyzed the data: LG QBH. Contributed reagents/materials/analysis tools: XRC LG CLL FLJ QBH. Wrote the paper: LG QBH.

Citation: Gong L, Chen X, Liu C, Jin F, Hu Q (2014) Gene Expression Profile of Bombyx mori Hemocyte under the Stress of Destruxin A. PLoS ONE 9(5): e96170. https://doi.org/10.1371/journal.pone.0096170

References

1. Kanzok SM, Jacobs-Lorena M (2006) Entomopathogenic fungi as biological insecticides to control malaria. Trends in Parasitology 22: 49–51.

2. Gao Q, Jin K, Ying SH, Zhang Y, Xiao G, et al. (2011) Genome sequencing and comparative transcriptomics of the model entomopathogenic fungi Metarhizium anisopliaeand and M. acridum. PLoS Genetic 7: e1001264.

3. Enserink M (2005) Mosquito-killing fungi may join the battle against malaria. Science 308: 1531–1533.

4. Hu QB, An XC, Jin FL, Freed S, Ren SX (2009) Toxicities of destruxins against Bemisia tabaci and its natural enemy, Serangium japonicum. Toxicon 53: 115–121.

5. Hoe PK, Bong CFJ, Jugah K, Rajan A (2009) Evaluation of Metarhizium anisopliae var. anisopliae (Deuteromycotina: Hyphomycete) isolates and their effects on subterranean termite Coptotermes curvignathus (Isoptera: Rhinotermitidae). American Journal of Agricutural and Biological Sciences 4: 289–297.

6. Lomer CJ, Bateman RP, Johnson DL, Langewald J, Thomas M (2001) Biological control of locusts and grasshoppers. Annual Review of Entomology 46: 667–702.

7. Hunter DM, Milner RJ, Spurgin PA (2001) Aerial treatment of the Australian plague locust, Chortoicetes terminifera (Orthoptera: Acrididae) with Metarhizium anisopliae (Deuteromycotina: Hyphomycetes). Bulletin of Entomological Research 91: 93–99.

8. Tounou K, Agbokaa K, Poehlinga M, Raupacha K, Langewaldb J, et al. (2003) Evaluation of the entomopatogenic fungi Metarhizium anisopliae and Paecilomyces fumosoroseus (deuteromycotina: Hyphomycetes) for control of the green leafhopper Empoasca decipiens (homoptera: Cicadellidae) and potential side effects on the egg parasitoid Anagrus atomus (hymenoptera: Mymaridae). Biocontrol Science and Technology 13: 715–728.

9. Shi WB, Feng MG (2004) Lethal effect of Beauveria bassiana, Metarhizium anisopliae and Paecilomyces fumosoroseus on the eggs of Tetranychus cinnaharinus (acari: Tetranychidae) with a description of a mite egg bioassay system. Biological Control 30: 165–173.

10. Yi F, Zou C, Hu Q, Hu M (2012) The joint action of destruxins and botanical insecticides (rotenone, azadirachtin and paeonolum) against the cotton aphid, Aphis gossypii Glover. Molecules 17: 7533–7542.

11. Hu QB, Ren SX, An XC, Qian MH (2007) Insecticidal activity influence of destruxins on the pathogenicity of Paecilomyces javanicus against Spodoptera litura. Journal of Applied Entomology 131: 262–268.

12. Wang B, Kang Q, Lu Y, Bai L, Wang C (2012) Unveiling the biosynthetic puzzle of destruxins in Metarhizium species. Proceedings of the National Academy of Sciences of the United States of America 109: 1287–1292.

13. Morais RP, Lira SP, Seleghim MHR, Berlinck RGS (2010) A method for destruxin analysis by HPLC-PDA-ELSD-MS. Journal of the Brazilian Chemical Society 21: 2262–2271.

14. Che Y, Swenson DC, Gloer JB, Koster B, Malloch D (2001) Pseudodestruxins A and B: new cyclicdepsipeptides from the Coprophilous fungus Nigrosabulum globosum. Journal of Natural Prodoucts 64: 555–558.

15. Pedras MSC, Zaharia IL, Ward DE (2002) The destruxins: synthesis, biosynthesis, biotransformation, and biological activity. Phytochemistry 59: 579–96.

16. Amiri B, Ibrahim L, Butt TM (1999) Antifeedant properties of destruxins and their potential use with the entomogenous fungus Metarhizium anisopliae for Improved Control of Crucifer Pests. Biocontrol Science and Technology 9: 487–498.

17. Meng X, Xu XX, Hu JJ, Jin FL, Hu QB, et al. (2011) Toxicity and differential protein analysis following destruxin A treatment of Spodoptera litura (Lepidoptera: Noctuidae) SL-1 cells. Toxicon 58: 327–335.

18. Ruiz-Sanchez E, Orchard I, Lange AB (2010) Effects of the cyclopeptide mycotoxin destruxin A on the Malpighian tubules of Rhodnius prolixus (Stål). Toxicon 55: 1162–70.

19. Ruiz-Sanchez E, Lange AB, Orchard I (2010) Effects of the mycotoxin destruxin A on Locusta migratoria visceral muscles. Toxicon 56: 1043–1051.

20. Muroi M, Shiragami N, Takatsuki A (1994) Destruxin B, a specific and readily reversible inhibitor of vacuolar-type H+-translocating ATPase. Biochemical and Biophysical Research Communications 205: 1358–1365.

21. Nakagawa H, Takami A, Udagawa N, Sawae Y, Suda K, et al. (2003) Destruxins, cyclodepsipeptides, block the formation of actin rings and prominent clear zones and ruffled borders in osteoclasts. Bone 33: 443–455.

22. Vey A, Matha V, Dumas C (2002) Effects of the peptide mycotoxin destruxin E on insect haemocytes and on dynamics and efficiency of the multicellular immune reaction. Joural of Invertebrate Pathology 80: 177–187.

23. Fan JQ, Chen XR, Hu QB (2013) Effects of destruxin A on hemocytes morphology of Bombyx mori. Journal of Integrative Agriculture 12: 101–108.

24. Pal S, St Leger RJ, Wu LP (2007) Fungal peptide destruxin A plays a specific role in suppressing the innate immune response in Drosophila melanogaster. Journal of Biological Chemistry 282: 8969–8977.

25. Mardis ER (2008) Next-generation DNA sequencing methods. Annual Review of Genomics and Human Genetics 9: 387–402.

26. Rieber N, Zapatka M, Lasitschka B, Jones D, Northcott P, et al. (2013) Coverage bias and sensitivity of variant calling for four whole-genome sequencing technologies. PLoS One 8: e66621.

27. Tokunaga F, Nakagawa T, Nakahara M, Saeki Y, Taniguchi M, et al. (2011) SHARPIN is a component of the NF-κB-activating linear ubiquitin chain assembly complex. Nature 471: 633–636.

28. Han PF, Jin FL, Dong XL, Fan JQ, Qiu BL, et al. (2013) Transcript and protein profiling analysis of the destruxin a-induced response in larvae of Plutella xylostella. PLoS ONE 8: e60771 (2013)..

29. James RR, Xu J (2012) Mechanisms by which pesticides affect insect immunity. Journal of Invertebrate Pathology 109: 175–182.

30. Lavine MD, Strand MR (2002) Insect hemocytes and their role in immunity. Insect Biochemistry and Molecular Biology 32: 1295–1309.

31. Wang Q, Liu Y, He HJ, Zhao XF, Wang JX (2010) Immune responses of Helicoverpa armigera to different kinds of pathogens. BMC Immunology 11: 9.

32. Qian C, Cao X (2013) Regulation of Toll-like receptor signaling pathways in innate immune responses. Annals of the New York Academy Sciences 1283: 67–74.

33. Roman D (2004) Peptidoglycan recognition proteins (PGRPs). Molecular Immunology 40: 877–886.

34. Wang S, Beerntsen BT (2013) Insights into the different functions of multiple peptidoglycan recognition proteins in the immune response against bacteria in the mosquito, Armigeres subalbatus. Insect Biochemistry and Molecular Biology 43: 533–543.

35. Park JA, Kim Y (2012) Eicosanoid biosynthesis is activated via Toll, but not Imd signal pathway in response to fungal infection. Journal of Invertebrate Pathology 110: 382–388.

36. Zaidman-Rémy A, Poidevin M, Hervé M, Welchman DP, Paredes JC, et al. (2011) Drosophila immunity: analysis of PGRP-SB1 expression, enzymatic activity and function. PLoS One 6: e17231.

37. Ao JQ, Ling E, Yu XQ (2008) A Toll receptor from Manduca sexta is in response to Escherichia coli infection. Molecular Immunology 45: 543–52.

38. Schulenburg H, Kurz CL, Ewbank JJ (2004) Evolution of the innate immune system: the worm perspective. Immunological Reviews 198: 36–58.

39. Andreas V (2013) Evolutionary plasticity of insect immunity. Journal of Insect Physiology 59: 123–129.

40. Oppert B, Kramer KJ, Johnson D, Upton SJ, McGaughey WH (1996) Luminal proteinases from Plodia interpunctella and hydrolysis of Bacillus thuringiensis CryIAc protoxin. Insect Biochemistry and Molecular Biology 26: 571–584.

41. Liu J, Tawa GJ, Wallqvist A (2013) Identifying cytochrome P450 functional networks and their allosteric regulatory elements. PLoS ONE 8: e81980.

42. Tissieres A, Mitchell HK, Tracy UM (1974) Protein synthesis in salivary glands of Drosophila melanogaster: Relation to chromosome puffs. Journal of Molecular Biology 85: 389–398.

43. Liao PC, Lin HY, Yuh CH, Yu LK, Wang HD (2008) The effect of neuronal expression of heat shock proteins 26 and 27 on lifespan, neurodegeneration, andapoptosis in Drosophila. Biochemical and Biophysical Research Communications 376: 637–641.

44. Shakya B, Jyoti S, Naz F, Khan S, Afzal RM, et al. (2012) Effect of L-ascorbic acid on the hsp70 expression and tissue damage in the third instar larvae of transgenic Drosophila melanogaster (hsp70-lacZ) Bg (9). Toxicology International 19: 301–305.

45. Asano T, Taoka M, Shinkawa T, Yamauchi Y, Isobe T, et al. (2013) Identification of a cuticle protein with unique repeated motifs in the silkworm, Bombyx mori. Insect Biochemistry and Molecular Biology 43: 344–351.

46. Zhang J, Goyer C, Pelletier Y (2008) Environmental stresses induce the expression of putative glycine-rich insect cuticular protein genes in adult Leptinotarsa decemlineata (Say). Insect Molecular Biology 17: 209–216.

47. Jiang H, Kanost MR (2000) The clip-domain family of serine proteinases in arthropods. Insect Biochemistry and Molular Biology 30: 95–105.

48. Hu QB, Ren SX, Wu JH, Chang JM, Musa PD (2006) Investigation of destruxin A and B from Metarhizium strains in China, and the optimization of cultural conditions for the strain MaQ10. Toxicon 48: 491–498.

49. Pais M, Das BC, Ferron P (1981) Depsipeptides from Metarhizium anisopliae. Phytochemistry 20: 715–723.

50. Larsen PE, Collart FR (2012) BowStrap v1.0: Assigning statistical significance to expressed genes using short-read transcriptome data. BMC Research Notes 5: 275.

51. Zhang H, Wei L, Miao H, Zhang T, Wang C (2012) Development and validation of genic-SSR markers in sesame by RNA-seq. BMC Genomics 13: 316.

52. Kim KI, van de Wiel MA (2008) Effects of dependence in high-dimensional multiple testing problems. BMC Bioinformatics 9: 114.

53. Livak KJ, Schmittgen TD (2001) Analysis of relative gene expression data using real-time quantitative PCR and the 2–ΔΔCT method. Methods 25: 402–408.