Gene 549 (2014) 85–96

Contents lists available at ScienceDirect

Gene journal homepage: www.elsevier.com/locate/gene

Analysis of midgut gene expression profiles from different silkworm varieties after exposure to high temperature Qing Rong Li, Yang Xiao, Fu Quan Wu, Ming Qiang Ye, Guo Qing Luo, Dong Xu Xing, Li Li, Qiong Yang ⁎ The Sericulture and Agri-Food Research Institute, Guangdong Academy of Agricultural Sciences, GuangZhou, China

a r t i c l e

i n f o

Article history: Received 21 April 2014 Received in revised form 2 July 2014 Accepted 17 July 2014 Available online 18 July 2014 Keywords: Silkworm Midgut Digital gene expression High temperature

a b s t r a c t The silkworm is a poikilothermic animal, whose growth and development is significantly influenced by environmental temperature. To identify genes and metabolic pathways involved in the heat-stress response, digital gene expression analysis was performed on the midgut of the thermotolerant silkworm variety ‘932’ and thermosensitive variety ‘HY’ after exposure to high temperature (932T and HYT). Deep sequencing yielded 6,211,484, 5,898,028, 5,870,395 and 6,088,303 reads for the 932, 932T, HY and HYT samples, respectively. The annotated genes associated with these tags numbered 4357, 4378, 4296 and 4658 for the 932, 932T, HY and HYT samples, respecti'vely. In the HY-vs-932, 932-vs-932T, and HY-vs-HYT comparisons, 561, 316 and 281 differentially expressed genes were identified, which could be assigned to 179, 140 and 123 biological pathways, respectively. It was found that some of the biological pathways, which included oxidative phosphorylation, related to glucose and lipid metabolism, are greatly affected by high temperature and may lead to a decrease in the ingestion of fresh mulberry. When subjected to an early period of continuous heat stress, HSP genes, such as HSP19.9, HSP23.7, HSP40-3, HSP70, HSP90 and HSP70 binding protein, are up-regulated but then reduced after 24 h and the thermotolerant ‘932’ strain has higher levels of mRNA of some HSPs, except HSP70, than the thermosensitive variety during continuous high temperature treatment. It is suggested that HSPs and the levels of their expression may play important roles in the resistance to high temperature stress among silkworm varieties. This study has generated important reference tools that can be used to further analyze the mechanisms that underlie thermotolerance differences among silkworm varieties. © 2014 Elsevier B.V. All rights reserved.

1. Introduction High temperatures can affect cellular structure and metabolism (Gibbs, 2002; Liu et al., 2006; Walter et al., 1990). For insects, environmental tolerance is an important factor for determining whether a population can survive within a particular ecological environment. Silkworms are poikilothermic, and thus environmental temperature influences their growth and development—especially in tropical and subtropical regions. We have observed that silkworm varieties associated with different regions, voltinisms, or systems exhibit extremely variable thermotolerances (Li QR, unpublished). Why these silkworm varieties have different thermotolerances, however, remains unclear. Stress-response systems, such as heat shock proteins (HSPs) and the antioxidase system, can improve organismal thermotolerance

Abbreviations: DGE, digital gene expression; DEGs, differentially expressed genes; RIN, RNA integrity number; TPM, transcripts per million clean tags; qRT-PCR, quantitative realtime reverse transcription PCR; GO, gene ontology. ⁎ Corresponding author. E-mail address: [email protected] (Q. Yang).

http://dx.doi.org/10.1016/j.gene.2014.07.050 0378-1119/© 2014 Elsevier B.V. All rights reserved.

(Pirkkala et al., 2001; Sorensen and Loeschcke, 2007; Stanley and Fenton, 2000). HSPs that are induced by high temperature may be involved in maintaining the integrity of cellular proteins (Denlinger and Yocum, 1998). Based on molecular weight, HSPs can be divided into five families, which include HSP100, HSP90, HSP70, HSP60, and sHSP (only sHSP is not conserved) (Kim et al., 1998; Li et al., 2009; Waters and Rioflorido, 2007; Waters et al., 2008). Among poly-, bi-, and univoltine silkworm strains, heat shock responses vary considerably and levels of thermotolerance increase as larval development proceeds (Manjunatha et al., 2010). However, it is still unclear which genes are transcriptionally regulated by high temperature and what mechanisms underlie these regulatory processes. Global transcriptome analysis could help characterize regulatory mechanisms that underlie high-temperature tolerance in the silkworm. Genome-wide transcriptional profiling is an important and powerful tool that has been used in silkworm to identify transcriptional programs underlying developmental processes (e.g., the silkworm life cycle) and defense responses to pathogens (Huang et al., 2009; Liu et al., 2009). Serial analysis of gene expression and microarray analysis have allowed us to visualize global changes in transcript abundance associated with spatial, temporal, or conditional differences (Huang et al., 2009; Liu et al.,

86

Q.R. Li et al. / Gene 549 (2014) 85–96

2009). Li et al. subjected silkworms (both parental strains and hybrids) to high-temperature stress and measured changes to the proteome and changes to protein phosphorylation levels within the posterior silk gland (Li et al., 2012). These results indicated that high temperatures alter the expression of proteins related to both the heat response and silk synthesis in silkworm. Moreover, proteins that were differentially expressed between the hybrids and parental lines likely explain the observed heterosis. Recently developed RNA deep-sequencing (RNA-seq) technologies, such as digital gene expression (DGE) (Wang et al., 2009), have substantially altered our understanding of the eukaryotic transcriptome in terms of extent and complexity. Here we analyzed two silkworm varieties, ‘932’ and ‘HY’, which have been selected based on temperature-tolerance tests involving N200 available silkworm varieties. The thermotolerant variety, ‘932’, is a bivoltine race that contains the China system. Its female parent is the polyvoltine race ‘Jiubaihai’, which can resist high temperature and high humidity. Its male parent is the bivoltine race ‘7302’. Its filial generation, named ‘932’, was bred in high-temperature and high-humidity conditions, endowing it with the ability to resist and adapt to relatively high temperatures (35–40 °C) and relative humidity (95%). In contrast, the Sericultural Research Institute of the Chinese Agricultural Academy cultivated the ‘HY’ variety to exhibit exceptional cocoon-silk traits (e.g., high yields and good quality). ‘HY’ is a bivoltine race that contains the Japan system. Its female parent is the high-yield race ‘782’, and its male parent is the high-yield and good-quality race ‘758’. Its filial generation, named ‘HY’, was bred in appropriate environmental conditions, so ‘HY’ has relatively weak resistance and adaptability to conditions of high temperature and humidity. To analyze transcriptional differences between thermotolerant and thermosensitive races, RNA-seq was used to analyze gene expression within the midgut of ‘932’ and ‘HY’ following exposure to high temperature. Bioinformatic analyses were then used to identify differentially expressed genes (DEGs) and differences in metabolic or function pathways. Through this study we hoped to identify genes and biological pathways that underlie differences in high-temperature tolerance among silkworm varieties and establish a theoretical platform for further investigating transgenic, high-resistance breeding.

2.3. Isolation of total RNA Total RNA was isolated using the Trizol Reagent (Invitrogen, Carlsbad, CA, USA). Each sample consisted of larval midgut from at least three silkworms. To remove residual genomic DNA, RNA samples were incubated with 10 U DNase (Roche Applied Science) at 37 °C for 20 min. RNA was then purified using RNeasy columns (Qiagen). RNA integrity was then confirmed using an Agilent 2100 Bioanalyzer. Samples had an average RIN (RNA integrity number) value of 9.3 (minimum RIN value was 8.9). A minimum of 6 μg total RNA was used for Illumina sequencing (Illumina Inc., San Diego, CA, USA). 2.4. Construction and evaluation of DGE-tag profiles

2. Materials and methods

DGE libraries for four samples (from 932T, 932, HYT, and HY) were constructed using Illumina gene expression sample-preparation kits. Briefly, oligo (dT)-coated magnetic beads were used to capture mRNAs from the four total-RNA samples. First- and second-strand cDNAs were synthesized, and bead-bound cDNA was digested with NlaIII. The 3′-cDNA fragments attached to the oligo (dT)-beads were then ligated to the Illumina GEX NlaIII adapter 1, which contains an MmeI recognition site. Samples were incubated with MmeI, which cuts 17 bp downstream of the recognition site, thereby leaving an adapter-1 tag. The 3′ fragments were then released from the beads via magnetic-bead capture, and an Illumina GEX adapter 2 was introduced at the MmeI site. The resulting adapter-ligated cDNA tags were then amplified using 15 cycles of PCR. The 105-bp PCR fragments were purified via 6% polyacrylamide Tris–borate–EDTA gel and subjected to Illumina/Solexa sequencing. The four tag libraries were constructed using in situ amplification and were subsequently deep-sequenced using an Illumina Genome Analyzer. Sequencing quality was evaluated and the data were summarized using the Illumina/Solexa pipeline software. Saturation of the four libraries was also analyzed. For the raw data, adaptor sequences were filtered, and the types of clean tags were represented as distinct clean tags. Subsequently, we classified clean tags and distinct clean tags according to their copy number within the library and calculated their percentage of the total clean and distinct tags. The raw data were deposited in the GEO database under submission number GSE55990.

2.1. Experimental animals

2.5. Annotation and analysis of DGE profiles

Two silkworm varieties were selected for analysis, namely ‘932’ (thermotolerant) and ‘HY’ (thermosensitive). Both are bivoltine.

For annotation, all tags were mapped to the reference sequence, which included silkworm unigenes from the silkworm Genomics Database (http://silkworm.genomics.org.cn/). Expression of each gene was estimated from the frequency of clean tags and then normalized to the number of transcripts per million clean tags (TPM) ('t Hoen et al., 2008), which is the standard method used for DGE analysis (Morrissy et al., 2009). TPM indicates reads per kilobase of transcript per million sequenced reads. For differential expression analysis, fold changes were assessed using the log2 ratio (TPM-T/TPM-932, TPM-HYT⁄TPMHY) after expression abundances were normalized to TPM. After assignment, gene ontology (GO) and pathway annotation and enrichment were analyzed using the NCBI COG (http://www.ncbi.nlm.nih. gov/COG) (Tatusov et al., 2003), the GO Database (http://www. geneontology.org/) (Gene Ontology Consortium, 2008), and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database (http://www.genome.jp/kegg/) (Kanehisa, 2002). The DGE analysis data (Gene ID and counts in four samples) were submitted to the GEO database under submission number GSE55990.

2.2. High-temperature treatment and sampling 2.2.1. Sampling for digital gene expression profile sequencing For both ‘932’ and ‘HY’, 50 healthy fifth-instar larvae were placed into a cultivation cabinet at 35 ± 0.5 °C and 75 ± 2% relative humidity. Control ‘932’ and ‘HY’ animals were placed at 27 ± 0.5 °C and 75 ± 2% relative humidity. After ~18 h the midgut was dissected from each silkworm, mulberry leaf digestate was removed from the enteric cavity, and samples were frozen in liquid nitrogen for storage. Samples were designated 932T (‘932’ after exposure to high temperature), 932 (‘932’ after control treatment), HYT (‘HY’ after exposure to high temperature), and HY (‘HY’ after control treatment). 2.2.2. Sampling to determine the mRNA expression levels of several HSPs, after heat treatment One hundred healthy fifth-instar larvae were placed into the same heat and control cultivation conditions as described in Section 2.2.1, for both the ‘932’ and ‘HY’ varieties. After 6, 12, 18, 24, 30 and 36 h, the midgut was dissected from each silkworm and frozen in liquid nitrogen for storage.

2.6. Real-time PCR Initially, to validate the DGEs obtained from Solexa sequencing, 12 genes: including chemosensory protein 11 (CSP, BGIBMGA004041-TA), juvenile hormone diol kinase (JHK, BGIBMGA008813-TA), glucan binding

Q.R. Li et al. / Gene 549 (2014) 85–96

87

Table 1 Primers used for qRT-PCR. Gene

Gene number

Primer

Primer sequence

Gene

Gene number

Primer

Primer sequence

CoADE

BGIBMGA011307-TA

BGIBMGA008814-TA

BGIBMGA001399-TA

CBD

BGIBMGA009067-TA

EIF2A

BGIBMGA007571-TA

SOD

BGIBMGA001307-TA

GST

BGIBMGA011820-TA

ANK

BGIBMGA001794-TA

GBP

BGIBMGA000353-TA

CaBP

BGIBMGA004115-TA

PSA

BGIBMGA012485-TA

CSP

BGIBMGA004041-TA

HSP19.9

BGIBMGA004540-TA

HSP23.7

BGIBMGA004515-TA

HSP40-3

BGIBMGA003153-TA

HSP90

BGIBMGA004612-TA

HSP70BP

BGIBMGA013070-TA

GCCGTTGGAGCAAAGACT CCCTCATGGCGAGTTCAT TTGAAACCAGGAAGGGAGC AGGGTCGGTAAGAATGAACTC CGGAGCAGACCTTCAAAC GGTTGCCGTGTAGTGGTG TGAGAACGCCCGATACGA CCGCCGAGAAATCTTGTG GACTTAGGGTTTCCGTTCC TTTGGCTTCGTGATGTCC TTGAATCAGGCTCTACCGA AGTGAGGCTTTCTTTGTTGC GCAAACACGAGGAGAAGAAAGA TGGGTGCGATTACAGACAACA CAAATCTAGGTTGCGTGGAG TGCGAGTAAATTGTGGGTGT TGGGAGGGTTTGCTATGTTT GGCTTCGCCTTGCTCTGATT

JHK

FBA

F R F R F R F R F R F R F R F R F R

ActA3

NM_001126254

F R F R F R F R F R F R F R F R F R

ATTAGCCAAAGCCCGAGGA CCCATAGCGACAGCCATTCT GTTGCCTTGTTGGCGACT CATTATTAGTTTCCCATTGAC CATGGTGGTCCCAGTTCTGC CCCAAGTCATCAGGGTCAGC TTCAGGCAATACGAACAAGC GCAAGCCGTGACCAGAGTGT TCCTGCCCTACCTTGAGACG GAATCATAGCGGGTTCAGTTT CTGTTCGCTGCTGTCTACTGT TGTTTCTGGGCTGGAGTAC GTTCACGACTACTTCCGACCTT TATCCGTGCTCGTCCTGTTT CAGGTGGCGAACTCCTCATT TCCTTCCTTTGTGACCGATA CTCCCTCGAGAAGTCCTACGAACT GGATGCCGCACGATTCCATAC

protein (GBP, BGIBMGA011608-TA), superoxide dismutase (SOD, BGIBMGA001307-TA), eukaryotic translation initiation factor 2A (EIF2A, BGIBMGA007571-TA), glutathione S-transferase omega 1 (GST, BGIBMGA011820-TA), Ankyrin repeat domain (ANK, BGIBMGA001794TA), calcium-binding protein (CaBP, BGIBMGA004115-TA), acyl-coenzyme A dehydrogenase (CoADE, BGIBMGA011307-TA), phosphoserine aminotransferase (PSA, BGIBMGA012485-TA), fructose-1,6-bisphosphatase (FBA, BGIBMGA001399-TA) and chitin binding domain 3 protein (CBD, BGIBMGA009067-TA), were subjected to quantitative real-time reverse transcription PCR (qRT-PCR) analysis using a Roche LightCycler480. The mRNA expression levels of some heat stress related proteins, which were differentially expressed between the heat treatment and control comparison, were then analyzed by qRT-PCR at different times after heat treatment. These proteins included HSP19.9 (BGIBMGA004540-TA), HSP23.7 (BGIBMGA004515-TA), HSP40 homolog 3 (BGIBMGA003153TA), HSP70 (BGIBMGA007950-TA), HSP90 (BGIBMGA004612-TA), and HSP70 binding protein (BGIBMGA013070-TA). The BmactinA3 gene (NM_001126254.1) was used as the endogenous control (Franzetti et al., 2012; Mounier and Prudhomme, 1986; Zhang et al., 2009). The cDNA was synthesized using 2 μg of

total RNA. PCR primers were designed using Primer Premier 5 software (Premier/Canada) and are listed in Table 1. The qRT-PCR analysis was performed by the SYBR-GREEN1 fluorescent relative quantitative approach, and a thermal denaturing step was used to generate a melting curve for the verification of amplification specificity. Three independent replicates were assessed from each data set and statistical analysis was performed using the 2−ΔΔCT method (Livak and Schmittgen, 2001).

3. Results 3.1. Tag identification and saturation analysis Using Solexa sequencing, the number of sequenced tags for the 932, 932T, HY, and HYT samples totaled 6,211,484, 5,898,028, 5,870,395, and 6,088,303, respectively. After filtering the adaptor sequences and removing low-quality tags, 6,139,297, 5,828,615, 5,799,074, and 6,005,746 tags remained. The number of distinct clean tags for 932, 932T, HY, and HYT totaled 50,153, 49,051, 47,463, and 55,780, respectively (Table 2). Distinct clean tags were then grouped based on copy number within the library

Table 2 Categorization and abundance of tags for the 932, 932T, HY, and HYT samples. Clean tags: tags that remain after filtering low-quality tags from the raw data. Distinct tags: different kinds of tags. Unambiguous tags: clean tags that remain after the removal of tags that mapped to reference sequences from multiple genes.

Raw data Clean tags All tag mapping to gene

Unambiguous tag mapping to gene

All tag-mapped genes Unambiguous tag-mapped genes Mapping to genome

Unknown tag

Summary

932

932T

HY

HYT

Total Distinct tag Total number Distinct tag number Total number Total % of clean tag Distinct tag number Distinct tag % of clean tag Total number Total % of clean tag Distinct tag number Distinct tag % of clean tag Number % of ref genes number % of ref genes Total number Total % of clean tag Distinct tag number Distinct tag % of clean tag Total number Total % of clean tag Distinct tag number Distinct tag % of clean tag

6,211,484 116,988 6,139,297 50,153 2,423,430 39.47% 11,031 21.99% 1,246,300 20.30% 10,436 20.81% 4357 29.80% 4173 28.54% 2,525,525 41.14% 25,848 51.54% 1,190,342 19.39% 13,274 26.47%

5,898,028 116,457 5,828,615 49,051 1,898,826 32.58% 10,823 22.06% 982,577 16.86% 10,215 20.83% 4378 29.94% 4182 28.60% 2,991,884 51.33% 24,577 50.10% 937,905 16.09% 13,651 27.83%

5,870,395 113,350 5,799,074 47,463 1,220,423 21.05% 10,003 21.08% 744,278 12.83% 9497 20.01% 4296 29.38% 4077 27.88% 3,513,487 60.59% 24,323 51.25% 1,065,164 18.37% 13,137 27.68%

6,088,303 133,417 6,005,746 55,780 1,773,412 29.53% 11,564 20.73% 993,226 16.54% 11,033 19.78% 4658 31.85% 4450 30.43% 2,878,578 47.93% 28,180 50.52% 1,353,756 22.54% 16,036 28.75%

88

Q.R. Li et al. / Gene 549 (2014) 85–96

Fig. 1. Distinct clean tag copy numbers for the 932, 932T, HY, and HYT samples. Copy numbers of more than 50% of the distinct clean tags in each sample are 2–5, followed by 6–10 copies accounting for about 14% and only 4–5% of the distinct clean tag copy numbers reach 50–100.

Fig. 2. The gene number mapped by all clean tags and unambiguous clean tags. The four graphs denote 932, 932T, HY and HYT respectively. The percentage of genes identified (y axis) increases as the total number of tags (x axis) increases.

Q.R. Li et al. / Gene 549 (2014) 85–96

89

(Fig. 1), which revealed that N55% of the distinct clean tags had a copy number of 2–5. This was true for each sample. Sequencing saturation is achieved when unique tags are no longer detected as the total number of tags continues to rise. To estimate whether the sequencing depth was sufficient to represent the entire transcriptome, sequencing saturation was analyzed for each of the libraries. Genes mapped by all clean tags and unambiguous clean tags increased with the total number of tags, but when sequencing counts surpassed 3 million, the number of detected genes was saturated (Fig. 2).

two varieties reared under normal conditions. Exceptions were HSP23.7 (BGIBMGA004515-TA) and HSP25.4 (BGIBMGA005781-TA) which showed higher and lower levels, respectively, in the thermotolerant variety, ‘932’, compared with the thermosensitive variety, ‘HY’. However some HSPs, such as HSP25.4 (BGIBMGA005781-TA), HSP19.9 (BGIBMGA004540-TA), HSP90 (BGIBMGA004612-TA), HSP23.7 (BGIBMGA004515-TA), HSP70-3 (BGIBMGA007950-TA) and HSP70 binding protein (BGIBMGA013070-TA), were up-regulated in the 932-vs-932T comparison, and showed higher levels of mRNA expression in the 932T compared with the HYT (Fig. 5).

3.2. Analysis of tag mapping

3.4. Comparison and analysis of gene ontology

A reference gene database (http://silkworm.genomics.org.cn/) was used to map tags. For distinct clean tags, 19.8–20.8% mapped unambiguously to the unigene database, 50.1%–51.5% mapped to the silkworm genome database, and 26.5%–28.8% did not map to the unigene virtual-tag database. Between 27.9% and 30.4% of the reference genes were unambiguously associated with a tag (Table 1). Mapping distributions associated with the four tag databases are shown in Fig. 3. We also found that in each of the four samples N 6% of the distinct clean tags perfectly matched the antisense strand.

After all tags were identified and mapped, about 4000 genes were unambiguously mapped (Table 1). GO analysis was used to functionally classify the four libraries. We analyzed differences in GO functional categories between the two control groups (HY-vs-932) and between experimental and control groups (932-vs-932T and HY-vs-HYT). Under normal conditions (27 ± 0.5 °C, and 75 ± 2% relative humidity), the 932 and HY control samples yielded 4296 and 4357 genes, respectively (Table 1), which could be categorized into 49 functional groups (Fig. 6A). Within the three main GO-classification categories (biological process, cellular component, and molecular function), there were 25, 14, and 10 functional groups, respectively. We also analyzed GO functional categories of DEGs associated with the HY-vs-932 comparison and found that 44 of these 49 GO functional categories contained DEGs (Fig. 6B). Genes categorized as biological adhesion (GO:0032091), cell junction (GO:0006810), and pigmentation (GO:0032401) were expressed at lower levels in 932 compared with HY. Genes associated with multi-organism process (GO:0042742), rhythmic process (GO:0031667), synapse (GO:0007269), antioxidant activity (GO:0010035), electron carrier activity (GO:0051234), and receptor activity (GO:0009755, GO:0007186) were expressed at higher levels in 932 compared with HY.

3.3. Screening for differentially expressed genes Methods for detecting DEGs within the DGE profile (Audic and Claverie, 1997) were used to calculate a false discovery rate of b 0.001 and a differential multiple of N2. Results showed that 316, 165 and 216 genes were up-regulated and 245, 151 and 65 genes were downregulated in the HY-vs-932, 932-vs-932 and HY-vs-HYT comparisons, respectively (Fig. 4). Worth mentioning is the mRNA expression changes of HSPs in the experiment. No significant differences were detected in the mRNA expression of most of the heat shock proteins between the

Fig. 3. Mapping pie chart of distinct clean tags for 932, 932T, HY, HYT samples. “Perfect match (PM)(Sense)” indicates that the tag perfectly matched the sense strand, “1 tag- N 1 gene” indicates that 1 tag matched 1 gene, “1 tag- N n gene” indicates that one tag matched n genes (n N 1), “1 Mismatch (MM)(Sense)” indicates the tag contained 1 nucleotide that did not match the sense strand, “PM(AntiSense)” indicates that the tag perfectly matched the antisense strand, “1 MM(AntiSense)” indicates that the tag contained 1 nucleotide that did not match the antisense strand, “PM Genome 1 tag- N 1 position” indicates that one tag matched one position within the genome; “PM Genome 1 tag- N n position” indicates that one tag matched n positions within the genome (n N 1), “1 MM Genome” indicates that the tag contained 1 nucleotide that did not match the genome, “Unknown Tag” indicates that the tag did not match a gene or a position within the genome.

90

Q.R. Li et al. / Gene 549 (2014) 85–96

400

Number of Genes

350

369

Down-regulated

316

300 250

Up-regulated

245

264 216

200

165

151

150 100

65

50 0 HY vs 932

HYT vs 932T

932 vs 932T

HY vs HYT

Fig. 4. DEGs associated with the comparisons HY-vs-932, 932-vs-932T, and HY-vs-HYT. HY-vs-932 is the comparison between 932G and HY controls (27 °C), 932-vs-932T is the comparison between 932 exposed to control conditions (27 °C, 932) and 932 exposed to high-temperature stress (35 °C, 932T), HY-vs-HYT is the comparison between HY exposed to control conditions (27 °C, HY) and HY exposed to high-temperature stress (35 °C, HYT), HYT-vs-932T is the comparison between HYT and 932T. The red column indicates the number of upregulated genes, and the green column indicates the number of down-regulated genes.

To compare the effects of high temperature on gene expression in thermotolerant and thermosensitive varieties, we analyzed GO functional enrichment for DEGs associated with 932-vs-932T (Fig. 7A) and HY-vs-HYT (Fig. 7B) comparisons. We identified 316 and 281 DEGs in the 932-vs-932T and HY-vs-HYT comparisons (Fig. 4), which were enriched for 43 and 39 functional groups, respectively (Fig. 7). In 932vs-932T, DEGs categorized as immune system process (GO:0007017), multi-organism process (GO:0042742), reproduction (GO:0003006, GO:0007269), rhythmic process (GO:0031667, GO:0048511), synapse (GO:0007269), electron carrier activity (GO:0009055), and enzyme regulator activity (GO:0030234, GO:0004866) were down-regulated by high temperature, whereas DEGs categorized as positive regulation of biological process (GO:0006289 etc.) were up-regulated. Other GO terms were associated with DEGs that were both up- and downregulated (Fig. 7A). In HY-vs-HYT, DEGs categorized as biological adhesion (GO:0031589, GO:0007162), cell proliferation (GO:0008283), multi-organism process (GO:0019725), reproduction (GO:0048610), membrane-enclosed lumen (GO:0016591), antioxidant activity (GO:0016209), enzyme regulator activity (GO:0004866), and receptor activity (GO:0004930) were up-regulated by high temperature, whereas DEGs categorized as nucleic acid binding transcription factor activity (GO:0001071) were down-regulated. Other GO terms were associated with DEGs that were both up- and down-regulated (Fig. 7B).

3.5. Pathway enrichment analysis of DEGs In comparisons between thermotolerant and thermosensitive control samples (HY-vs-932), several biological pathways including glycerolipid metabolism, oxidative phosphorylation, protein digestion and absorption, steroid biosynthesis, and fat digestion and absorption were significantly enriched for DEGs (Table 3) and mRNA levels of most DEGs in those pathways were up-regulated compared to the thermosensitive variety (Figs. 8A, B, C, E, F). In comparisons between groups subjected to high temperature and control groups, several biological pathways including oxidative phosphorylation, protein digestion and absorption, glycolysis, glutathione metabolism, tyrosine and tryptophan biosynthesis, citrate cycle etc., were significantly enriched for DEGs associated with 932-vs-932T (Table 4); and mRNA levels of most DEGs in pathways of oxidative phosphorylation, glycolysis and citrate cycle showed down-regulated compared to the control groups (Figs. 8D, E, F). Similarly, biological pathways including protein digestion and absorption, glycerolipid and glycerophospholipid metabolism, steroid biosynthesis, ribosome biogenesis in eukaryotes, antigen processing and presentation, vitamin digestion and absorption, and PPAR signaling pathway were significantly enriched for DEGs associated with HY-vs-HYT (Table 5). 3.6. Confirmation of DEGs using qRT-PCR To confirm the subset of DEGs, 12 genes were selected for qRT-PCR analysis. In the 932-vs-932T, HY-vs-HYT comparisons, respectively, relative expression levels of 12 genes mostly showed a significant change in the same direction between the qRT-PCR and RNA-seq methodologies (Fig. 9). Only three genes in the HY-vs-HYT comparison, acylcoenzyme A dehydrogenase (CoADE) (Fig. 9E), phosphoserine aminotransferase (PSA) (Fig. 9J) and fructose-1,6-bisphosphatase (FBA) (Fig. 9K), exhibited a conflicting significant change in the opposite direction, between qRT-PCR and RNA-seq data. This result indicated a substantial conformity between qRT-PCR and RNA-seq data. The inconsistent verification of a few genes may be due to technology-dependent difference, and this discrepancy between qRT-PCR and the DGE data had also been reported in the previous studies (Hegedus et al., 2009). 3.7. Gene expression analysis of heat stress related proteins

Fig. 5. Heatmap showing expression levels difference of some HSPs among four samples, 932, 932T, HY and HYT. The color of a point in the figure represents the gene level (reference to the value of TPM). The four columns represent samples HY, HYT, 932 and 932T. Lines represent different genes. It can be seen from the figure that gene expression levels are highest in 932T.

To further analyze the expression changes and differences shown by HSPs under heat stress, the relative mRNA expression levels of HSP19.9, HSP23.7, HSP40-3, HSP70, HSP90 and HSP70 binding protein were analyzed by qRT-PCR. The mRNA levels from the intestines were assessed at different times after high temperature treatment and the difference

Q.R. Li et al. / Gene 549 (2014) 85–96

91

Fig. 6. GO terms associated with the HY-vs-932 comparison. GO terms associated with total mapped genes (HY and 932) (A) and DEGs (HY-vs-932) (B) are shown. Based on sequence identity, N4000 total mapped genes for 932 and HY (A) and 561 DEGs in HY-vs-932 (B) were categorized into 49 and 44 functional groups, respectively.

in expression was analyzed between the heat tolerant and heat sensitive varieties. The results showed that mRNA levels of all HSPs were higher in the heat tolerant variety under early continuous high temperature processing (approximately before 24 h) than in the normal control at the same stage of development. However, levels were lower than in the normal control, at the same developmental stage, after 24 h (Figs. 10A, C, E, G, I, K). It was found that, with the exception of HSP70, mRNA levels of HSPs in the midgut of the thermotolerant variety were significantly higher than those in the thermosensitive variety during continuous high temperature treatment (Figs. 10B, D, F, H, J, L). 4. Discussion and conclusion Digital gene expression techniques, such as next generation highthroughput deep sequencing and bioinformatic analysis, are able to detect changes in gene expression in organisms maintained under different conditions. Hegedus et al. (2009) obtained a sequencing depth of greater than five million tags per zebrafish sample, using the Solexa/Illumina DGE system, and over 70% of all tags mapped to the

transcript database (Hegedus et al., 2009). The midgut is an important organ for digestion and the absorption of nutrients. It is also the first natural barrier against external pathogenic microorganisms. Here, to investigate transcriptional changes associated with the high-temperature response of different silkworm races, we analyzed a thermotolerant race, ‘932’, and a relatively thermosensitive race, ‘HY’. These strains were subjected to 27 °C or 35 °C for 18 h, and global gene expression was assayed in the midgut. The conditions were selected on the basis of our preliminary laboratory results to create the continuous heat stress of the experiment (Li et al., 2013). Deep sequencing resulted in greater than six million tags per sample and the number of annotated genes was 4300. Interestingly, 26.47–28.75% of all distinct tags were unknown, implying that these tags may represent unassembled sequences or new transcripts. High temperature seriously affects some biological pathways (Downer and Kallapur, 1981; Malik and Reddy, 2008; Sohal et al., 1985). It has been shown that many insect metabolic pathways, such as lipid metabolism and lipid peroxidation, are severely affected by high-temperature stress (Downer and Kallapur, 1981). In addition,

92

Q.R. Li et al. / Gene 549 (2014) 85–96

Fig. 7. GO terms associated with DEGs identified after exposure to high temperature. GO terms associated with DEGs identified in the 932-vs-932T (A) and HY-vs-HYT (B) comparisons are shown. (A) A total of 316 DEGs were categorized into 43 functional groups. (B) A total of 281 DEGs were categorized into 39 functional groups.

high temperatures may accelerate the rate of carbohydrate consumption in insects. Ivanovic et al. observed a decrease in midgut amylolytic and proteolytic activities at lower temperatures; when Morimus funereus was reared at 30 °C, stored glycogen was consumed within three days (Ivanovic, 1991). We performed a comparative analysis of

metabolic pathways enriched by DEGs in the HY-vs-932, 932-vs-932T and HY-vs-HYT comparisons. Firstly, in the HY-vs-932 comparison, DEGs were primarily involved in the following metabolic pathways: glycerolipid metabolism, oxidative phosphorylation, protein digestion and absorption, steroid biosynthesis, and fat digestion and absorption

Table 3 Significantly enriched pathways within DEGs associated with HY-vs-932. #

Pathway

DEGs with pathway annotation (433)

All genes with pathway annotation (8473)

P value

Pathway ID

1 2 3 4 5 6 7 8 9 10

Metabolic pathways Pancreatic secretion Glycerolipid metabolism Oxidative phosphorylation Protein digestion and absorption Parkinson's disease Steroid biosynthesis Fat digestion and absorption Alpha-Linolenic acid metabolism Drug metabolism-other enzymes

131 (30.25%) 42 (9.7%) 25 (5.77%) 20 (4.62%) 26 (6%) 18 (4.16%) 13 (3%) 20 (4.62%) 12 (2.77%) 17 (3.93%)

1387 (16.37%) 286 (3.38%) 160 (1.89%) 125 (1.48%) 197 (2.33%) 122 (1.44%) 71 (0.84%) 154 (1.82%) 68 (0.8%) 123 (1.45%)

1.052100e −13 3.857217e −10 4.783716e −07 4.668776e −06 7.366638e −06 4.313358e −05 5.201517e −05 0.0001050476 0.0001457576 0.0001602123

ko01100 ko04972 ko00561 ko00190 ko04974 ko05012 ko00100 ko04975 ko00592 ko00983

Q.R. Li et al. / Gene 549 (2014) 85–96

93

Fig. 8. Heatmap showing the expression level difference of some biological pathways. The color of a point in the figure represents the gene level (reference to the value of TPM). Columns represent different samples and lines represent different genes. In comparison ‘932’ with ‘HY’, expression levels of some different expression genes (DEGs) of Glycerolipid metabolism, oxidative phosphorylation and steroid biosynthesis were shown on A, B and C, respectively; D showed that of oxidative phosphorylation in comparison 932 with 932T; some DEGs levels of glycolysis and citrate cycle among four samples, 932, 932T, HY and HYT, were shown on E and F, respectively.

(Table 3). Most DEGs were expressed at higher levels in the thermotolerant race, ‘932’ (Figs. 8A, B, C). These data may indicate that the activity of the above mentioned pathways are kept at a higher level in the thermotolerant ‘932’ variety than in the heat sensitive ‘HY’ variety, under normal conditions. Secondly, in the 932-vs-932T comparison, the DEGs were primarily involved in some significantly different metabolic pathways, such as oxidative phosphorylation (Fig. 8D), glycolysis

(Fig. 8F), glutathione metabolism, and tricarboxylic acid cycle (Fig. 8E), were down-regulated by continuous high temperature stress. Previous studies and observations (unpublished) have shown that continuous high temperature treatment can significantly reduce the amount of fresh mulberry ingested by the silkworm and body size is significantly smaller than that of individuals reared under normal conditions. Continuous high temperature may impact the above mentioned

Table 4 Significantly enriched pathways within DEGs associated with 932-vs-932T. #

Pathway

DEGs with pathway annotation (243)

All genes with pathway annotation (8473)

P value

Pathway ID

1 2 3 4 5 6 7 8 9 10

Metabolic pathways Oxidative phosphorylation Protein digestion and absorption Glycolysis/gluconeogenesis Rheumatoid arthritis Glutathione metabolism Parkinson's disease Cysteine and methionine metabolism Collecting duct acid secretion Huntington's disease

71 (29.22%) 16 (6.58%) 16 (6.58%) 8 (3.29%) 7 (2.88%) 8 (3.29%) 11 (4.53%) 6 (2.47%) 5 (2.06%) 14 (5.76%)

1387 (16.37%) 125 (1.48%) 197 (2.33%) 56 (0.66%) 49 (0.58%) 64 (0.76%) 122 (1.44%) 39 (0.46%) 31 (0.37%) 214 (2.53%)

2.65617e−07 4.906006e−07 0.0001641922 0.0001761520 0.0004495068 0.0004506031 0.0007299597 0.0007711817 0.001718148 0.003363982

ko0100 ko00190 ko04974 ko00010 ko05323 ko00480 ko05012 ko00270 ko04966 ko05016

94

Q.R. Li et al. / Gene 549 (2014) 85–96

Table 5 Significantly enriched pathways within DEGs associated with HY-vs-HYT. #

Pathway

DEGs with pathway annotation (205)

All genes with pathway annotation (8473)

P value

Pathway ID

1 2 3 4 5 6 7 8 9 10

Pancreatic secretion Protein digestion and absorption Alpha-linolenic acid metabolism Glycerolipid metabolism Glycerophospholipid metabolism Steroid biosynthesis Ribosome biogenesis in eukaryotes Antigen processing and presentation Vitamin digestion and absorption PPAR signaling pathway

23 (11.22%) 17 (8.29%) 7 (3.41%) 11 (5.37%) 9 (4.39%) 6 (2.93%) 7 (3.41%) 4 (1.95%) 6 (2.93%) 6 (2.93%)

286 (3.38%) 197 (2.33%) 68 (0.8%) 160 (1.89%) 139 (1.64%) 71 (0.84%) 95 (1.12%) 34 (0.4%) 79 (0.93%) 90 (1.06%)

3.684626e−07 5.20151e−06 0.001216536 0.001695617 0.006496887 0.007232946 0.00802847 0.00874273 0.01198170 0.02160844

ko04972 ko04974 ko00592 ko00561 ko00564 ko00100 ko03008 ko04612 ko04977 ko03320

pathways, to adversely affect their normal role in growth and development of the silkworm. Thirdly, biological pathways, which included protein digestion and absorption, glycerolipid and glycerophospholipid metabolism, PPAR signaling pathway, and vitamin digestion and absorption, were enriched for DEGs associated with the HY-vs-HYT comparison. Interestingly, the antigen processing and presentation pathway was also associated with HY-vs-HYT. This pathway is involved in immune function, bacterial and viral infection (Brodsky, 1992; Chain et al., 1986) and thus may help confer high-temperature resistance to the temperature-sensitive silkworm, HY. After being subjected to high temperature, the insect's ability to resist adverse environmental stimuli is significantly reduced. This can lead to susceptibility to pathogen infection, which activates the immune-related response. The ability to resist

stress may be weaker in the heat sensitive variety, ‘HY’, than in the thermotolerant variety, ‘932’. In addition, our results showed that significant differences in the level of mRNA expression of some HSPs exist between the thermotolerant variety and the thermosensitive variety, in both the normal t and high temperature group. The results of the qRT-PCR analysis showed that mRNA levels of all detected HSPs and HSP70 binding protein are up-regulated during early continuous high temperature processing (approximately before 24 h), but are lower than in the normal temperature control after 24 h. The role of HSPs in thermotolerance has been examined in plants, insects and other animals. In yeast cells, the major protein responsible for thermotolerance is hsp 100 but in Drosophila cells hsp 70 appears to be the primary protein involved

Fig. 9. Quantitative RT-PCR validation of tag-mapped genes. In the graph, columns are the result of real-time PCR; lines represent the trends of digital gene expression profiles based on the transcripts per million clean tags; ‘932’ and ‘HY’ in the x-axis represent the sample of different varieties.

Q.R. Li et al. / Gene 549 (2014) 85–96

100

150

100

14000 12000 10000 6000

4000

HSP23.7

16000

CK T

932-HSP23.7

16000

Relative expression level (%)

Relative expression level (%)

Relative expression level (%)

HSP19.9

200

150

18000

18000

250 CK T

932-HSP19.9

14000

Relative expression level (%)

250

200

95

932T HYT

12000 10000 8000 6000 4000

2000

2000

Time

0 18h

24h

30h

36h

6h

12h

18h

A

30h

400

200

Time

0 24h

30h

400

200

6h

12h

24h

30h

1500 1000

Time 6h

12h

18h

I

2000

1000

24h

30h

36h

4000

2000

12h

18h

24h

30h

Time 0

36h

6h

12h

1500 1000 500

Time 12h

18h

18h

24h

30h

24h

30h

36h

H

932-HSP70BP

25

HSP70BP

2000

6h

932T HYT

1000

25

0

36h

3000

Time 6h

500 0

30h

HSP70

932T HYT

2500

24h

6000

3000

36h

Relative expression level (%)

2000

18h

D

0

HSP90

3500

12h

G

3000

2500

6h

5000

F Relative expression level (%)

Relative expression level (%)

18h

Time

0

36h

4000

Time

36h

932-HSP90

30h

932-HSP70

600

0

3500

24h

5000

E 3000

18h

932T HYT

Relative expression level (%)

Relative expression level (%)

Relative expression level (%)

600

18h

12h

C

HSP40-3

12h

6h

800

CK T

6h

Time

0

36h

B

932-HSP40-3

800

24h

Relative expression level (%)

12h

36h

20

Relative expression level (%)

6h

Time

0

15

10

5

Time

0 6h

12h

J

18h

24h

K

30h

36h

20

15

10

5

Time 0 6h

12h

18h

24h

30h

36h

L

Fig. 10. qRT-PCR analysis of HSP19.9, HSP23.7, HSP40-3, HSP70, HSP90 and HSP70 binding protein during continuous high temperature stress in the midgut of the silkworm thermotolerant race ‘932’ (A, C, E, G, I and K) and between the thermotolerant ‘932’ and sensitive ‘HY’ strains (B, D, F, H, J and L, respectively). In the qRT-PCR analysis, relative expression levels of HSP19.9, HSP23.7, HSP40-3, HSP70, HSP90 and HSP70 binding protein were calculated after normalizing against BmActinA3 expression as an internal control, using the same RNA samples. Each data point represents the mean ± SE (n = 3).

(Parsell et al., 1993). Carmel et al. found a high positive correlation between the levels of Hsp40 expression and induced thermotolerance in fruit flies, which implies a significant contribution of the Hsp40 gene to thermoadaptation (Carmel et al., 2011). It was found that increasing the magnitude and duration of pretreatment increased Hsp70 concentrations and improved the tolerance to more severe stress in Drosophila melanogaster (Krebs and Feder, 1998). In this study, it has been found that the thermotolerant race, ‘932’, has higher mRNA levels of some HSPs, apart from HSP70, than the thermosensitive variety during high temperature continuous treatment. It can be suggested that expression levels of HSPs, such as HSP19.9, HSP23.7, HSP40-3 and HSP90, may play important roles in the different resistance to high temperature stress among silkworm varieties. Manjunatha et al. reported that heat shock responses vary considerably among poly-, bi-, and univoltine silkworm strains, and levels of thermotolerance increase as larval development proceeds (Manjunatha et al., 2010). Calabria et al. reported that flies that carry the warm-climate chromosome arrangement O(3 + 4) have higher basal levels of Hsp70 protein than their cold-climate O(st) counterparts. The O(3 + 4) carriers are also more heat tolerant (Calabria et al., 2012). In conclusion, in this study, we found that continuous high temperature stress has a detrimental influence on some metabolic or functional pathways, such as glucose and lipid metabolism and oxidative phosphorylation, which may be related to the normal growth and development of the silkworm. Expression levels of some HSPs may play important roles in thermotolerance among different silkworm varieties. More functional studies and analyses are still needed.

Acknowledgements This study was supported by grants from the National Natural Science Foundation of China (NSFC) (No. 31201856) to QR Li, and the Presidential Foundation of the Guangdong Academy of Agricultural Sciences (201214) to QR Li.

References Audic, S.,Claverie, J.M., 1997. The significance of digital gene expression profiles. Genome Res. 7, 986–995. Brodsky, F.M., 1992. Antigen processing and presentation: close encounters in the endocytic pathway. Trends Cell Biol. 2, 109–115. Calabria, G., Dolgova, O., Rego, C., Castañeda, L.E., Rezende, E.L., Balanyà, J., Pascual, M., Sørensen, J.G.,Loeschcke, V.,Santos, M., 2012. Hsp70 protein levels and thermotolerance in Drosophila subobscura: a reassessment of the thermal co-adaptation hypothesis. J. Evol. Biol. 25 (4), 691–700. Carmel, J., Rashkovetsky, E., Nevo, E., Korol, A., 2011. Differential expression of small heat shock protein genes Hsp23 and Hsp40, and heat shock gene Hsr-omega in fruit flies (Drosophila melanogaster) along a microclimatic gradient. J. Hered. 102 (5), 593–603. Chain, B.M., Kay, P.M., Feldmann, M., 1986. The cellular pathway of antigen presentation: biochemical and functional analysis of antigen processing in dendritic cells and macrophages. Immunology 58, 271–276. Denlinger, D.L., Yocum, G.D., 1998. Physiology of heat sensitivity. In: Hallman, G.J., Denlinger, D.L. (Eds.), Temperature Sensitivity in Insects and Application in Integrated Pest Management, 29. Westview Press, Oxford. Downer, R.G., Kallapur, V.L., 1981. Temperature induced changes in lipid composition and transition temperature of flight muscle mitochondria of Schistocerca gregaria. J. Therm. Biol. 6, 189–194. Franzetti, E., Huang, Z.J., Shi, Y.X., Xie, K., Deng, X.J., Li, J.P., Li, Q.R., Yang, W.Y., Zeng, W.N., Casartelli, M., Deng, H.M., Cappellozza, S., Grimaldi, A., Xia, Q., Tettamanti, G., Cao, Y.,

96

Q.R. Li et al. / Gene 549 (2014) 85–96

Feng, Q., 2012. Autophagy precedes apoptosis during the remodeling of silkworm larval midgut. Apoptosis 17 (3), 305–324. Gene Ontology Consortium, 2008. The Gene Ontology project in 2008. Nucleic Acids Res. 36, D440–D444. Gibbs, A.G., 2002. Lipid melting and cuticular permeability: new insights into an old problem. J. Insect Physiol. 48, 391–400. Hegedus, Z., Zakrzewska, A., Agoston, V.C., Ordas, A., Rácz, P., Mink, M., Spaink, H.P., Meijer, A.H., 2009. Deep sequencing of the zebrafish transcriptome response to mycobacterium infection. Mol. Immunol. 46, 2918–2930. Huang, L.L.,Cheng, T.C.,Xu, P.Z.,Cheng, D.J.,Fang, T.,Xia, Q.Y., 2009. A genome-wide survey for host response of silkworm, Bombyx mori during pathogen Bacillus bombyseptieus infection. PLoS One 4, e8098. Ivanovic, J., 1991. Metabolic response to stressors. In: Ivanovic, J., Jankovic Hladni, M. (Eds. ), Hormones and Metabolism in Insect Stress. CRC Press, pp. 27–67. Kanehisa, M., 2002. The KEGG database. Novartis Found. Symp. 247, 91–101 (discussion 101–103, 119–128, 244–152). Kim, K.K., Kim, R., Kim, S.H., 1998. Crystal structure of a small heat-shock protein. Nature 394, 595–599. Krebs, R.A.,Feder, M.E., 1998. Hsp70 and larval thermotolerance in Drosophila melanogaster: how much is enough and when is more too much? J. Insect Physiol. 44, 1091–1101. Li, Z.W.,Li, X.,Yu, Y.Q.,Xiang, Z.H.,Kishino, H.,Zhang, Z., 2009. The small heat shock protein (sHSP) genes in the silkworm, Bombyx mori, and comparative analysis with other insect sHSP genes. BMC Evol. Biol. 9, 215–229. Li, J.S., Ye, L.P., Lan, T.Y., Yu, M.L., Liang, J.S., Zhong, B.X., 2012. Comparative proteomic and phosphoproteomic analysis of the silkworm (Bombyx mori) posterior silk gland under high temperature treatment. Mol. Biol. Rep. 39, 8447–8456. Li, Q.R., Xiao, Y., Wu, F.Q., Zheng, Q., Yang, Q., Luo, G.Q., Ye, M.Q., Xing, D.X., Li, L., Tang, C.M., Wang, Z.J.,Dai, F.W., 2013. Effects of high temperature and humidity on ultrastructure of intestinal epithelial cells in different silkworm races. Sci. Seric. 39 (5), 65–69. Liu, H.T., Huang, W.D., Pan, Q.H., Weng, F.H., Zhan, J.C., Liu, Y., Wan, S.B., Liu, Y.Y., 2006. Contributions of PIP2-specific-phospholipase C and free salicylic acid to heat acclimation-induced thermotolerance in pea leaves. J. Plant Physiol. 163, 405–416. Liu, S.P.,Zhang, L.,Li, Q.B.,Zhao, P.,Duan, J.,Cheng, D.J.,Xiang, Z.H.,Xia, Q.Y., 2009. MicroRNA expression profiling during the life cycle of the silkworm (Bombyx mori). BMC Genomics 10, 455. Livak, K.J.,Schmittgen, T.D., 2001. Analysis of relative gene expression data using real-time quantitative PCR and the 2(−Delta Delta C (T)) method. Methods 25, 402–408. Malik, F.A., Reddy, Y.S., 2008. Effect of high temperature on haemolymph sugar levels in three selected silkworm races. Acta Entomol. Sin. 51 (11), 1113–1120.

Manjunatha, H.B., Rajesh, R.K., Aparna, H.S., 2010. Silkworm thermal biology: a review of heat shock response, heat shock proteins and heat acclimation in the domesticated silkworm, Bombyx mori. J. Insect Sci. 10, 204. Morrissy, A.S., Morin, R.D., Delaney, A., Zeng, T., McDonald, H., Jones, S., Zhao, Y., Hirst, M., Marra, M.A., 2009. Next generation tag sequencing for cancer gene expression profiling. Genome Res. 19, 1825–1835. Mounier, N., Prudhomme, J.C., 1986. Isolation of actin genes in Bombyx mori: the coding sequence of a cytoplasmic actin gene expressed in the silk gland is interrupted by a single intron in an unusual position. Biochimie 68 (9), 1053–1061. Parsell, D.A.,Taulien, J.,Lindquist, S., 1993. The role of heat-shock proteins in thermotolerance. Philos. Trans. R. Soc. Lond. B Biol. Sci. 339 (1289), 279–285 (29). Pirkkala, L.,Nykanen, P.,Sistonen, L., 2001. Roles of the heat shock transcription factors in regulation of the heat shock response and beyond. FASEB J. 15, 1118–1131. Sohal, R.S., Muller, A., Koletzko, B., 1985. Effect of age and ambient temperature on n-pentane production in adult housefly, Musca domestica. Mech. Ageing Dev. 29 (3), 317–326. Sorensen, J.G., Loeschcke, V., 2007. Studying stress responses in the post-genomic era: its ecological and evolutionary role. J. Biosci. 32, 447–456. Stanley, K., Fenton, B., 2000. A member of the HSP60 gene family from the peach potato aphid, Myzus persicae (Sulzer). Insect Mol. Biol. 9, 211–215. 't Hoen, P.A., Ariyurek, Y., Thygesen, H.H., Vreugdenhil, E., Vossen, R.H., de Menezes, R.X., Boer, J.M., van Ommen, G.J., den Dunnen, J.T., 2008. Deep sequencing-based expression analysis shows major advances in robustness, resolution and inter-lab portability over five microarray platforms. Nucleic Acids Res. 36, e141. Tatusov, R.L.,Fedorova, N.D.,Jackson, J.D.,Jacobs, A.R.,Kiryutin, B.,Koonin, E.V.,Krylov, D.M., Mazumder, R.,Mekhedov, S.L.,Nikolskaya, A.N., et al., 2003. The COG database: an updated version includes eukaryotes. BMC Bioinformatics 4, 41. Walter, M.F.,Petersen, N.S., Biessmann, H., 1990. Heat shock causes the collapse of the intermediate filament cytoskeleton in Drosophila embryos. Dev. Genet. 11, 270–279. Wang, Z.,Gerstein, M.,Snyder, M., 2009. RNA-Seq: a revolutionary tool for transcriptomics. Nat. Rev. Genet. 10, 57–63. Waters, E.R., Rioflorido, I., 2007. Evolutionary analysis of the small heat shock proteins in five complete algal genomes. J. Mol. Evol. 65, 162–174. Waters, E.R., Aevermann, B.D., Sanders-Reed, Z., 2008. Comparative analysis of the small heat shock proteins in three angiosperm genomes. Cell Stress Chaperones 13, 127–142. Zhang, X.,Hu, Z.Y.,Li, W.F.,Li, Q.R.,Deng, X.J.,Yang, W.Y.,Cao, Y.,Zhou, C.Z., 2009. Systematic cloning and analysis of autophagy-related genes from the silkworm, Bombyx mori. BMC Mol. Biol. 27 (10), 50.