Transcriptome Analysis in the Fat Body of Two Silkworm (Bombyx mori) Strains with Different Susceptibility to Fenvalerate

Sericulture is the main economic source of farmers in many developing countries such as China and India (Zhao et al., 2011). The total silk production of China accounts for 80% of the world’s total (Li et al., 2010). However, to date, pesticide poisoning and various diseases still cause great damage to the sericulture (Wang et al., 2013; Gu et al., 2014). Pyrethroid pesticide fenvalerate is an efficient and low toxic insecticide commonly used in mulberry fields in China (Shi et al., 2018). Pesticide pollution is becoming more and more serious due to the emergence of insecticide resistance (Tian et al., 2017). Silkworm is a non target insect that is extremely susceptible to insecticides in the ecological environment (Peng et al., 2011). Pesticide exposure may affect growth, reproduction and cocoon, which will bring huge economic losses to the sericulture (Wang et al., 2011).

As a main silk-producing insect, the silkworm is an economically important insect, and it is also used as a model insect for investigating the genome and molecular genetics of the Lepidoptera order (Zhao et al., 2011). The fat body is an important intermediate metabolic tissue in insects (Nath et al., 1997). As reflected by the transcriptome pattern observed from Drosophila melanogaster and B. mori, the fat body is a multifunctional tissue in insect (Tian et al., 2017; Hu et al., 2016). It plays a crucial role in nutrient synthesis, storage and substance metabolism, and participates in the metabolism and detoxification of toxic substanc (Tian et al., 2016). However, to date, the impact of fenvalerate exposure on the transcriptome pattern in the fat body of silkworms with different pesticide susceptibility has not been studied and reported. We tested 198 silkworm strains by bioassay for their fenvalerate tolerance and found that a notable difference exists in different strains (data unpublished). Of those silkworm strains tested, Lan5, a silkworm strain from Zhejiang province of China, is susceptible to fenvalerate and Mysore, a silkworm strain from India, is tolerant to fenvalerate.

High-throughput transcriptome sequencing is an efficient, accurate and economic technique to obtain new information about the expression of whole transcriptome (Surget-Groba et al., 2010; Wang et al., 2009). Using this technique, differentially expressed genes and critical target pathways in multiple samples can be analyzed conveniently and accurately (Qin et al., 2011). In this study, a transcriptome sequencing system was employed to examine differentially expressed genes in the fat body of two strains of silkworms with different susceptibility after exposure to fenvalerate. Tag annotation, gene ontology (GO) and KEGG pathway analyses were performed. The differentially expressed genes were verified using real-time PCR at the transcriptional level. These results can provide a reference for further research on the mechanism of silkworm’s tolerance to pesticides and the new control strategy of Lepidoptera pests.

According to the result of pre-experiment, fenvalerate solution (Sigma-Aldrich Trading Co, Ltd, USA) was diluted to a 0.02 mg/L concentration, which is less than LC50 for 24 h (Zhao et al., 2018). Mulberry leaves were immersed in the solutions for 10 s, and allowed to dry naturally, before being fed the fifth instar silkworms on the second day for three times a day.

Domesticated silkworm strain Lan5 (susceptible strain) and Mysore (tolerant strain) were provided by the National Center for Silkworm Genetic Resources Preservation of the Chinese Academy of Agricultural Sciences. They were reared at standard temperature under a photoperiod of 12 h of light and 12 h of dark. At 24 h after feeding with fenvalerate solution coated leaves, the larvae and the unfed control larvae were dissected. The fat bodies of ten larvae were collected and rinsed three times in diethylpyocarbonate (DEPC) treated double distilled H2O. All samples were stored at -70 °C until use. Each treatment was undertaken in triplicate.

For Illumina sequencing, total RNA was isolated from the fenvalerate-exposed fat body as well as the control fat body by using TRIzol reagent (Invitrogen, USA) according to the manufacturer’s protocol. To remove any residual DNA, samples were pretreated with RNase-free DNase I (TaKaRa) for 30 min at 37°C. The quality of RNA was assessed by formaldehyde agarose gel electrophoresis and was quantitated spectrophotometrically.

The RNA-seq libraries were prepared using TruSeqTM DNA Sample Preparation Kit-Set A (Illumina, San Diego, USA) following the manufacturer’s protocol (Illumina, San Diego, CA). Approximately 20 μg of total RNA from the fenvalerate-treated fat body or the wild fat body was used to isolate mRNA using Sera-mag Magnetic Oligo (dT) Beads (Illumina). To avoid cDNA synthesizing bias by priming, the purified mRNA was fragmented into small pieces (100-400 bp) using divalent cations at 94°C for 5 minutes. The double-stranded cDNA was first synthesized using the SuperScript Double-Stranded cDNA Synthesis kit (Invitrogen, Camarillo, CA) with random hexamer (N6) primers (Illumina). The libraries were sequenced using Illumina Highseq 2000 platform (Illumina, San Diego, USA) according to the manufacturer’s instructions.

The raw reads were filtered to remove low-quality sequences, the generated clean reads were mapped directly into deposited reference gene sequences generated from SilkDB v2.0 database (ftp://silkdb.org/pub/current/Gene/) using TopHat 2.0.6 software (Trapnell et al., 2010). The differentially expressed genes were identified at a Q value (corrected P value of T test with False discovery rate) of 0.05 (Wang et al., 2010). With regard to a gene which is only expressed in one sample, the differentially expressed gene was identified at a fold-change value of 1 (FC-test). Software Heatmap was used to plot the heatmap of differentially expressed genes based on the log10 of P value of two groups, ±10 was set as the extreme value.

The possible functions of all differentially expressed genes were determined using the gene ontology (GO) classification system (http://www.geneontology.org/). GO term with a corrected P value ≤ 0.05 was designated as significantly enriched term in differentially expressed genes. WEGO (Web Gene Ontology Annotation Plot) software was used for visualizing, comparing and plotting GO annotation results (Ye et al., 2006). Pathway enrichment analysis was carried out using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (http://www.genome.jp/kegg/). Pathway with a corrected P value ≤ 0.05 was designated as significantly enriched pathways in differentially expressed genes.

To confirm the RNA-Seq data, the quantitative real-time PCR was carried out with a real-time reverse transcription-PCR system (ABI 7300, Applied Biosystems, USA) with SYBR green. Total RNA was extracted from the fat body of both the treated and control larvae after 24 h using Trizol reagent (Takara). First-strand cDNA was synthesized with M-MLV Reverse transcriptase and an oligo (dT) primer. A total of 20 μl volume of real-time PCR reaction solution contained 0.2 µg cDNA, 5 pmol of each primer, and 10 μL of SYBR Green Real-time PCR Master Mix. PCR reactions were run in triplicates with 2 biological replicates with thermal cycling parameters at 95°C for 10 min followed by 40 cycles at 95°C for 30 sec, 60°C for 1 min according to the manufacturer’s protocol. The specific primers of 20 genes were designed based on the cDNA sequences (Table I). The housekeeping gene of B. mori actin3 was used as an internal control for the normalization. The relative expression level of these genes was estimated according to the 2-ΔΔCt method (Livak et al., 2001). All the samples were measured independently three times.

In our project, we sequenced 12 samples on Illumina HiSeq Platform in total and generated about 6.64 Gb per sample. The average genome mapping rate is 79.93% and the average gene mapping rate is 80.09%. 14,194 genes were identified in which 13,126 of them are known genes and 1,096 of them are novel genes. 11,666 novel transcipts were identified in which 5,776 of them are previously unknown splicing event for known genes, 1,096 of them are novel coding transcripts without any known features, and the remaining 4,794 are long noncoding RNA.

The differentially expressed genes in the infected fat body samples were identified by FDR≤0.001 and |log2 ratio|≥1. In the susceptible variety Lan5 DEG class, a total of 2363 genes were found with 490 up-regulated and 1873 down-regulated. In the tolerant variety Mysore DEG class, a total of 1611 genes were found with 746 up-regulated and 865 down-regulated.

Gene ontology (GO) assignments were used to understand the functions of these differentially expressed genes revealed by transcriptomics analysis. The DEGs were termed by GO ontology into three categories, namely cellular component, molecular function and biological process.

In the susceptible strain Lan5 DEG class, differentially expressed genes have been categorized into a total of 46 functional groups using the WEGO software. In each of the three main categories of the GO classification, “metabolic process”, “cellular process”, “cell”, “cell part”, “binding”, and “catalytic activity” terms are dominant. In almost all the terms, more down-regulated genes were observed than the up-regulated genes. Some genes related to the immune system process, behavior, locomotion, biological adhesion, cell killing, growth, cell junction and nucleic acid binding transcription factor activity were all down-regulated.

Meanwhile, in the tolerant strain Mysore DEG class, differentially expressed genes have also been categorized into a total of 46 functional groups using the WEGO software. In each of the three main categories of the GO classification, “metabolic process”, “cellular process”, “cell”, “cell part”, “binding” and “catalytic activity” terms are dominant. Interestingly, some genes related to the extracellular region, extracellular region part and cell junction were up regulated, whereas more down regulated genes were observed than the up regulated genes in other terms.

In order to further understand the biological functions of the differentially expressed genes, KEGG (http://www.genome.jp/kegg) ontology assignments were used to classify functional annotations of the identified genes. There are six branches for KEGG pathways: Cellular processes, environmental information processing, genetic information processing, human disease, metabolism, and organismal systems. There is a difference between up (down) regulated pathways observed in Mysore compared with Lan5. In Lan5, genes involved in metabolic pathways were the most significantly enriched. In Mysore, 12 metabolism pathways and 10 that were involved in organismal systems were identified.

20 co-regulated differentially expressed genes were verified by using qRT-PCR. The results show that 15 genes were concordant with the DGE libraries, but 3 genes from Lan5 (Cyp6b29, UGT40B4, Mor) and 2 genes from Mysore (Cyp4l6, GSTo4) were non-concordant. Although some differences were observed in the results of qRT-PCR and DGE analysis, the overall regulate trend of most genes was consistent, which indicated that the results of the DGE data were accurate.

Out of these determined genes, 6 genes were up-regulated, 6 genes were down-regulated and 8 genes exhibited opposite trends in the two DGE libraries, and in the 8 genes, 7 genes were down-regulated in Lan5 and up-regulated in Mysore. More of these kinds of genes were detoxification genes such as Cyp6b29, GSTs2 and CarE5A.

This study is the first report systematically analyzing the expression profile and regulation of genes in fat body of two silkworm strains with different susceptibility after exposure to fenvalerate. On a genome-wide scale, profiles of transcriptome provided information for deep insights into functional genomics by its ability to interpret the functional elements of the genome, biological pathways, and molecular mechanisms. RNA-Seq used in this study is a recently developed technology which is cost efficient and powerful for the rapid identification and analysis of majority part of the whole transcriptome (Snyder, 2009). In this study, RNA-Seq technology was employed to analyze the transcriptome change in fat body of two silkworm strains after feeding fenvalerate. GO and KEGG analysis were used to classify the functions of DEGs and identify the difference of pathways involved in the fenvalerate regulation.

Fenvalerate-exposed susceptible Lan5 and tolerant Mysore silkworms were subjected to DGE to perform a global analysis of their transcriptomes. The Lan5 silkworms yielded 2363 differentially expressed genes and the Mysore silkworms yielded 1611 differentially expressed genes. Among the 2363 differentially expressed genes in the Lan5 strain, the number of down-regulated genes (1873) was higher than the number of up-regulated genes (490), but more up-regulated genes (746) and less down-regulated genes (865) were detected in the Mysore strain. This difference in the two silkworm strains might be involved in different susceptibility or defense responses to fenvalerate, which is similar to the previous studies (Gao et al., 2014). However, the possible mechanism of the susceptibility differences between two silkworm strains and molecular information on these differentially expressed genes still needs to be clarified in our future research.

Many studies showed that the fat body is the major detoxification tissue in B. mori. There are some detoxification enzymes in fat body, such as cytochrome P450s, glutathione-S-transferase and carboxylesterase, which are able to hydrolyze insecticide and eliminate them from silkworm body through metabolism. Cytochrome P450s are a complex and ubiquitous superfamily of heme-containing enzymes that participate in metabolism of exogenous substrates (Li et al., 2014). There are 84 P450s genes in all and 78 of them are functional genes in the silkworm (Yamamoto et al., 2010). In this study, fenvalerate-induced up-regulation of many cytochrome P450, such as Cyp305b1, Cyp6ab4, Cyp9a21 and Cyp18a1, were observed in the fat body of both two silkworm strains, which suggested their role in detoxification. At the same time, the expression levels of Cyp4l6 and Cyp6b29 were up-regulated by fenvalerate in the fat body of tolerant Mysore, while down-regulation were observed in susceptible Lan5. On the other hand, GSTo4 and CarE5A showed the similar expression level changes in the fat body of two silkworm strains. We suppose that these four detoxification genes may be involved in the susceptibility differences between two silkworm strains.

As activators of the innate immune system in many organisms, HSPs are a family of proteins that help organisms protect from environmental-induced cellular damage. Wu et al. (2011) observed the up-regulation of several HSPs after BmCPV infection in silkworm. Gu et al. demonstrated that the transcription level of Hsp21.4 was found up-regulated in the fat body after exposure to phoxim for 24 h (Gu et al., 2015). In this study, the up-regulations of transcription levels of Hsp20.4 and Hsp20.8 were also found in the fat body of Mysore silkworm after exposure to fenvalerate for 24 h. This result suggested that HSPs may be functionally involved in protecting tissues against injury of fenvalerate.

In B. mori, 6 groups of antimicrobial peptides (AMPs) have been identified and characterized by a BLASTP search from silkworm genome database (Cheng et al., 2006). AMPs such as moricins, gloverins, and attacins play an important role in eliminating invades. In our study, the transcription levels of glv3, Mor and PGRP-S1 were observed down-regulated after exposure to fenvalerate in both two silkworm strains. Due to the down-regulation of these AMPs, the silkworm larvae may become more susceptible to infection of various microbes.

In conclusion, this is the first report that employed the DGE technique to examine the transcriptomes differences between two silkworm strains after exposure to fenvalerate. The current results can provide a reference for further research on the mechanism of silkworm’s tolerance to pesticides and the new control strategy of Lepidoptera pests.

This work was supported by the Foundation of Post Scientist in National Sericultural System (Grant No. CARS-22-ZJ0101), the Natural Science Foundation of Jiangsu Province (Grant No. BK20181228), and the Key Research & Development program of Zhenjiang (Grant No. NY2017017).

The authors have declared no conflict of interests.