Transcriptome analysis of hsp18.3 functions and regulatory systems using RNA-sequencing in the red flour beetle, Tribolium castaneum

1. lntroduction

Heat shock proteins (HSPs) are a multi-gene family coding for abundant molecular chaperones that are typically expressed and synthesized in cells when they are exposed to stress and that were originally identified in Drosophila melanogaster (Ritossa and Vonborstel 1964). Usually,HSPs can be divided into five families, including HSP100,HSP90, HSP70, HSP60, and small HSPs (sHSPs), based on their molecular weight and the homology of their amino acid sequences (Lindquist and Craig 1988; Csermely et al. 1998). All of these families are highly evolutionarily conserved, except for the sHSP family, which is more diverse than other four (Li Z et al. 2009). Nevertheless,there are still some common characteristics among sHSPs in terms of distribution, structure and function. It has been reported that sHSPs are abundant and ubiquitous in almost all organisms ranging from prokaryotes to eukaryotes and including humans (Waters and Rio florido 2007; Aevermann and Waters 2008). The molecular weights of this family range from 12 to 42 kDa (Kim et al. 1998; Waters et al.2008). Additionally, these sHSP sequences have a relatively conserved α-crystallin domain that comprises approximately 80–100 amino acid residues located near the C-terminal region (Fu et al. 2006). There is also a conserved β-sheet sandwich in a sHSP secondary structure that is responsible for dimer formation (van Montfort R et al. 2001; van Montfort R L et al. 2001). This multimeric structure formed by sHSPs has the function of a molecular chaperone, which binds to the proteins and prevents them from thermal denaturation(Lee et al. 1995; Nakamoto and Vigh 2007).

In addition to helping the unfolding proteins maintain their correct states, binding to denatured proteins and preventing irreversible protein aggregation during high temperature stress (Haslbeck et al. 2005; Chen and Zhang 2015), sHSPs can also develop the protection function in response to other stresses, such as cold, drought, oxidation and hypertonic stress, even heavy metals (Kegel et al. 1996; Waters et al. 2008). Concurrently, the results showed that sHSPs were associated with multiple cellular activities, including organization of the cytoskeletal integrity, maintaining cellular communication and signaling transduction as well as modulating membrane lipid polymorphism (Arrigo 2000;Tsvetkova et al. 2002; Haslbeck et al. 2005; Bakthisaran et al. 2015). Furthermore, these sHSPs also participated in many other physiological processes, including cell cycle regulation and cell differentiation (Liang et al. 2007;Bakthisaran et al. 2015), interfering with apoptosis and defending against diseases, and enhancing the immune response (Arrigo 2000; Altincicek et al. 2008).

Recent studies have showed that an increasing number of sHSPs can be found in different insect species and have various functions. Knockdown of hsp22.6 remarkably reduced temperature tolerance in Apis cerana cerana, and Acchsp22.6 was induced by abiotic and biotic stresses,including insecticides (cyhalothrin and pyridaben) as well as hormone (20-hydroxyecdysone) and bacteria(Ascosphaera apis) treatments (Zhang et al. 2014).Similarly, hsp27 transcript levels were dramatically upregulated in response to temperature changes and chemical exposure in both Chironomus riparius and Lucilia cuprina(Martinez-Paz et al. 2014; Singh and Tiwari 2016). A significant reduction in starvation resistance was also associated with the genotype without a functional hsp27 gene in Drosophila (Hao et al. 2007). Typically, overexpressing fly mitochondrial hsp22 increased resistance to oxidative stress and extended lifespan by 32%, whereas flies that were not expressing this small hsp22 had a 40%decrease in lifespan and showed a reduction of 30% in locomotor activity (Morrow et al. 2004). These outcomes suggest there is a key role for shsps in antioxidant activity,locomotor activity and aging in insects. Furthermore,Zhang C et al. (2015, 2016) indicated that hsp25.4 and hsp21.4 mRNA levels were accumulated in the fat body to a maximum level at 1.5 and 3 h, respectively, when exposed to nucleopolyhedroviruses in Antheraea pernyi, which suggests shsps may be promising candidates for boosting the innate immunity of insects. Additionally, shsps are necessary for development and reproduction in Grapholita molesta and Drosophila, respectively (Michaud et al. 1997; Zhang B et al. 2015). Hsp21.3 exhibited the highest fold change in diapause individuals, which was highly significant, and along with five other shsps (hsp18.9, hsp19.8b, hsp21.7, hsp24.8 and hsp31.8), showed up-regulated expression during G. molesta diapauses (Zhang B et al. 2015). While hsp23 and hsp27 were significantly expressed in the absence of stress in the testes of Drosophila, the results implied that these two shsps play a vital role during spermatogenesis(Michaud et al. 1997). Thus, shsps not only are essential for different stress responses but have many physiological functions, such as innate immunity, development and reproduction in the insects.

Most recently, ten shsp genes have been identified in Tribolium castaneum, the red flour beetles, which is a notorious worldwide pest of stored grain and farinaceous materials for human consumption (Li Z et al. 2009; Contreras et al. 2013). Tchsp27, which is one of the T. castaneum species-specific shsps, was renamed as Tchsp18.3 based on its predicted molecular weight (Altincicek et al. 2008; Li Z et al. 2009). The amino acid sequences of Tchsp18.3 demonstrated the highest homology to hsp23.0 and hsp27.0 of Drosophila and hsp20.4 of B. mori (Xie et al. 2017).Hsp23.0 played an important role in hypoxia tolerance and muscle-specific over-expression of hsp23.0 protected muscle from heat shock stress in flies (Azad et al. 2009;Kawasaki et al. 2016). Meanwhile, Drosophila hsp27.0 was involved in the defense response against fungi, bacteria and prolonged exposure to dichlorvos (Pandey et al. 2016;Zhang L et al. 2016). In Bombyx mori non-diapause eggs,hsp20.4 displayed sudden rise in expression level from 42 to 48 h after oviposition (Ponnuvel et al. 2010). These results implied that Tchsp18.3, which was likely similar to its homologues in other insect species, possesses the functional diversity. Currently, it’s reported that Tchsp18.3 was significantly induced by the septic wounding and ultraviolet A (UV-A) exposure, indicated it was involved in immunity and stress responses of T. castaneum (Altincicek et al. 2008; Sang et al. 2012). This phenomenon was also observed after 1 h of exposure to UV-A light in the cell line BCIRL-TcA-CLG1 of T. castaneum. The outcome suggested that Tchsp18.3 plays a vital role in coping with deleterious stimuli to maintain cell normal growth and viability in vivo and in vitro (Garcia-Reina et al. 2017). Additionally, our previous study illustrated that Tchsp18.3 promoted positive effects in response to a high temperature as well as on resistance to starvation stresses in T. castaneum (Xie et al. 2017). Larval RNA interference (RNAi) of Tchsp18.3 caused significant lethality at the pupal stage, which seriously hindered the insects eclosion (Xie et al. 2017). The lacking of Tchsp18.3 at the pupal stage also decreased the eggs production of the red flour beetles (Xie et al. 2017). However, it is still unclear how Tchsp18.3 performs these functions in an insect model.To further explore the signaling systems of Tchsp18.3 involved in insect physiology, RNA sequencing technology was introduced in the present study, and these findings also shed new light into molecular regulatory mechanisms of other relative genes in T. castaneum.

2. Materials and methods

2.1. Beetle rearing

The T. castaneum strains Georgia-1 (GA-1) were reared at 30°C in whole wheat flour enriched with 5% yeast powder under standard conditions (Chen et al. 2016; Song et al.2017).

2.2. Double-strand RNA (dsRNA) synthesis and injection

For dsRNA synthesis, primers containing Tchsp18.3-specific sequences (sense primer: 5´-CGACCAGCT CATCGTTTCCT-3´, antisense primer: 5´-TTGTTGCACC GCTGGTGTA-3´) and the T7 polymerase promoter(5´-TAATACGACTCACTATAGGG-3´) at the 5´-end of both the sense primer and antisense primer was designed to amplify dsDNAs of Tchsp18.3. Tchsp18.3 dsRNAs synthesis and injection were also performed as described in a previous study (Sang et al. 2016). In this study, late larvae (20 days old) injected with Tchsp18.3 dsRNA were denoted as ds-Tchsp18.3 groups, whereas larvae injected with an equal volume of physiological buffer were set as control groups.The physiological buffer solution was prepared by combining 372.75 mg KCl and 38 mg Na3PO4·12H2O, then diluting it to a total volume of 1.0 L with ultrapure water. The pH of this solution will be 7.0.

2.3. RNA isolation and sequencing

Total RNAs were isolated from a pool of the six larvae of ds-Tchsp18.3 and control samples on the fifth day after injection with the RNAiso Plus reagent (TaKaRa, Dalian,China) according to the manufacturer’s instructions. The quality and quantities of total RNAs were determined with 1% agarose gels using an Agilent 2100 Bioanalyzer and by spectrophotometry on a NanoDrop Spectrophotometer(Thermo Scientific, Waltham, MA, USA), respectively. After extracting the total RNAs from the samples, poly(A) mRNA was isolated by using oligo(dT) beads and then treated with the fragmentation buffer. The cleaved RNA fragments(approximately 200 bp) were purified with a QIAquick PCR Extraction Kit (Qiagen, Germany) and then transcribed into first-strand cDNA using reverse transcriptase and random hexamer primers. This procedure was followed by second strand cDNA synthesis using buffer, dNTPs, DNA polymerase I and RNaseH following the manufacturer’s recommendation. The double-stranded cDNA was further purified with a QiaQuick PCR Extraction Kit and washed with EB buffer. Finally, sequencing adaptors were disrupted into the fragments and the required fragments were enriched by PCR amplification. Single-end RNA sequencing libraries of control and ds-Tchsp18.3 samples were prepared and sequenced on the Illumina HiSeq™ 2000 Platform. Raw sequence data are submitted to the Short Read Archive(SRA) of the NCBI database, and the accession numbers are SRR3087513 (control) and SRR5753969 (ds-Tchsp18.3),respectively.

2.4. Quality control and functional annotation

Firstly, ERANGE Software (http://woldlab.caltech.edu/gitweb/) was used to calculate the normalized gene locus expression level by assigning reads to their origin sites and counting them (Wang Z et al. 2009). The expression level of a gene from the RNA sequencing was normalized with the RPKM method (Reads Per Kilobase of exon model per Million mapped reads) (Li et al. 2014). The cutoff value for determining the gene transcriptional pro filing of each gene was determined based on a 95% confidence interval for all of the RPKM values (Mortazavi et al. 2008). Then,all of the detected beetle genes were assigned to different function categories using the software Blast2GO (http://www.blast2go.org/) (Conesa et al. 2005). Then, the Gene Ontology (GO) classi fication was carried out using WEGO(Web Gene Ontology Annotation Plot) (http://wego.genomics.org.cn/cgibin/wego/index.pl) to classify these GO terms by using Fisher’s exact test with a false discovery rate (FDR)correction to obtain an adjusted P-value between certain test gene groups and the whole annotation genome (Ye et al.2006). Finally, the KEGG biological pathway analysis was further performed by Cytoscape software (version 2.6.2)(http://www.cytoscape.org/) with the ClueGO plugin (http://www.ici.upmc.fr/cluego/cluegoDownload.shtml).

Sequencing produced raw image data that were transformed by base calling into sequence data, which were defined as raw reads. All the raw sequences generated from the sequencing machines were cleaned by discarding the adaptor sequences, containing N sequences and low-quality sequences before data analysis. The clean reads were then aligned to the reference genome of T. castaneum(http://www.beetlebase.org/) using SOAPaligner/soap2 and allowing up to two base mismatches (Li R et al. 2009;Liu et al. 2011). Failed mapped reads were progressively trimmed off, one base at a time, from the 3´-end and mapped to the genome again until a match was found (unless the read had been trimmed by ＜27 bases) (Xiong et al. 2017).The obtained data were further analyzed for functional annotation of all the detected genes.

2.5. ldentification and annotation of differentially expressed genes (DEGs)

All DEGs and transcript abundances were calculated with the Program R and the DEGseq packaged with a Random Sampling model based on the read count for each gene in the ds-Tchsp18.3 and control samples (Wang et al.2010). The combinations of FDR≤0.001 and the absolute value of log2Ratio≥1 were the thresholds to judge the significance of each gene expression difference (Storey and Tibshirani 2003). The GO functional annotation and KEGG enrichment analyses were performed to identify which DEGs were significantly enriched in GO terms and KEGG biological pathways compared with the wholetranscriptome background using the formula described in a previous study (Liu et al. 2011). The calculated P-value was further subjected to a Bonferroni correction using a corrected P-value≤0.05 as a threshold. GO terms fulfilling this condition were defined as significantly enriched GO terms in DEGs (Chen et al. 2016). We used the Blastall Program to annotate the pathways of DEGs against the KEGG database. Subsequently, we adopted the same formula as the formula in GO analysis to carry out DEG pathway enrichment analysis. In this study, KEGG pathways with a Q-value≤0.05 were significantly enriched pathways with DEGs (Li et al. 2014).

1.2.3 统计学方法 应用SPSS 13.0统计软件包，计量资料用均数、百分比进行描述，组间比较采用t检验。

2.6. Quantitative real-time PCR (qRT-PCR)

This study identified 569 DEGs from the ds-Tchsp18.3 and control samples. These DEGs suggested that Tchsp18.3 was not only vital for various stress responses but also involved in innate immunity indirectly through viperin,dorsal, Hdd11, PGRP2, defensin1 and defensin2 as well as other immune-related genes. Additionally, knockdown of Tchsp18.3 suppressed the antioxidant activity process,which most likely modulated the effects of Tchsp18.3 in development and reproduction of T. castaneum. The DEGs encoding Blimp-1, Gld, Drm, Kinesin-14, Pthr2, Delta(11)-like and EGF-like domain protein 2 were associated with the development and reproduction of ds-Tchsp18.3.Furthermore, the serine protease signaling pathway was activated by Tchsp18.3 to participate in the beetle’s stress responses and innate immunity as well as development and reproduction. These results provide clues to screen for the target genes that are involved in important biological processes of insects and could lead to the discovery of more suitable and sustainable management for pest control.

3. Results

3.1. lllumina sequencing and reads assembly

To uncover the differences in T. castaneum gene expression between the control group and ds-Tchsp18.3 group, a pooled cDNA sample from two groups was sequenced with the Illumina sequencing platform. We con firmed that Tchsp18.3 was knocked down in the ds-Tchsp18.3 samples compared to the control groups before sequencing (Appendix B). RNA sequencing generated a total of 14 154 435 sequence reads,including 7 315 113 clean reads and 6 839 322 clean reads from the control and ds-Tchsp18.3 of the red flour beetles,respectively (Table 1; Appendix C). In the total control reads, 65.97% mapped either to a unique (63.03%) or to multiple (2.94%) genomic locations, whereas 64.75% of the ds-Tchsp18.3 reads showed either a unique match (62.31%)or a multiple-position match (2.44%) (Table 1). A total of 12 761 and 12 586 genes were then detected in the control and ds-Tchsp18.3 samples, respectively, which ranged from 81 to 63 354 bp. These genes could be divided into seven categories based on their lengths (Fig. 1-A). Approximately half of the genes were distributed between 600 to 1 800 bp,and an average of 31.99% of the genes ranged from 600 to 1 200 bp, whereas 21.68% ranged from 1 200 and 1 800 bp(Appendix D). There was a similar distribution pattern for the detected genes in each category between the two samples (Fig. 1-A). Additionally, RPKM values were further applied to normalize and quantify the gene expression.The majority gene RPKM values ranged from 10 to 50 for both the control (5 263 genes, 41.24%) and ds-Tchsp18.3(5 202 genes, 41.23%) groups, whereas fewer RPKM values were over 1 000, which suggested that there was only small proportion of genes with a high expression level (Fig. 1-B; Appendix E). No significant differences were observed from the analysis of the RPKM values between the two groups (Fig. 1-B). To analyze whether the number of the detected genes increased proportionally with the sequencing amount, sequence saturation analyses were performed. The results showed that the number of detected genes stopped increasing when the number of sequence reads reached 1.5 million (Appendix F). For mRNA expression, heterogeneity and redundancy are two key characteristics. A small proportion of mRNA is highly expressed, while the majority of mRNA is expressed at low levels. Thus, the distribution of tag expression was used to check the normality of the RNA sequencing data. The results showed that there were nearly equal proportions with more than half of detected genes having a coverage ＞80% between the two libraries (Fig. 2-A and B;Appendix G). Overall, these two libraries yielded 13 299 genes (Appendix D) and provided abundant data for analyzing Tchsp18.3 regulatory systems.

要想充分保证校企合作的积极效果，就必须在展开人才培养的过程中，切实做好教学与实践的有机结合，不断强化学校和行业企业之间的沟通、合作与交流。相对于普通高等院校的教育教学活动而言，职业院校的体育教学也表现出了一定的独特性，具体来说就是属性的特殊性，是具备双重属性的，既有普通教育的属性，也有职业专门性教育的属性。

Table 1 Summary of reads from the RNA sequencing data1)

1) Of the total (100%) control reads, 65.97% matched either a unique (63.03%) or multiple (2.94%) genomic locations, whereas 64.75%of the total (100%) ds-Tchsp18.3 reads showed either a unique match (62.31%) or a multiple-position match (2.44%). The remaining 34.03 and 35.25% of reads were unmapped in the control and ds-Tchsp18.3, respectively.

Total reads Total mapped reads Unique matches Multi-position matches Total unmapped reads Control 7 315 113 (100%) 4 825 890 (65.97%) 4 610 640 (63.03%) 215 250 (2.94%) 2 489 223 (34.03%)ds-Tchsp18.3 6 839 322 (100%) 4 428 282 (64.75%) 4 261 293 (62.31%) 166 989 (2.44%) 2 411 040 (35.25%)

Fig. 1 The range in gene lengths (bp) and expression levels(RPKM) of control and ds-Tchsp18.3. A, length distribution of the genes detected by RNA sequencing. B, RPKM range of genes of the control and ds-Tchsp18.3.

3.2. GO classification and KEGG analysis of all detected genes

Additionally, 18 DEGs and 9 DEGs were involved in the development process (Table 4) and reproduction (Appendix J) under biological processes (Fig. 3), respectively. The DEGs encoding Blimp-1, Gld, Drm, Kinesin-14, Pthr2,Delta(11)-like and EGF-like domain protein 2 were related to the effects of Tchsp18.3 on development and reproduction in T. castaneum. specifically, Blimp-1 was up-regulated and the other 6 genes were down-regulated by RNAi for Tchsp18.3 (Table 4; Appendix J). The Blimp-1 expression was significantly induced by 20-hydroxyecdysone (20E)and its product existed during high-ecdysteroid periods and turned over rapidly in Drosophila (Agawa et al. 2007).Additionally, its prolonged expression resulted in a reduction of ftz transcription factor 1 (ftz-f1) expression, which is essential for insect development and reproduction (Agawa et al. 2007). A temporally specific ftz-f1 knockdown using RNAi during the pupal stage exhibited defects in eclosion and in the eye, the antennal segment, the wing and the leg in Drosophila (Sultan et al. 2014). Additionally, the loss of ftz-f1 gene expression led to a complete block in egg production through reducing vitellogenin transcripts and affecting oocyte maturation in T. castaneum (Xu et al. 2010). The results indicated that over-expression of Blimp-1 inhibited the expression of ftz-f1 and produced developmental abnormalities as well as reduced fertility.Apart from its function in immune responses, Gld is also required for insect eclosion as well as normal sperm storage and utilization (Cavener and MacIntyre 1983; Cavener et al. 1986; Iida and Cavener 2004). The Drosophila Gld mutants displayed a recessive lethal trait at the late pupal stage, and these mutants could be rescued by excising the anterior end of the pupal case 0–2 days prior to the normal adult emergence time. This outcome suggested that the function of Gld in fly pupae is to aid in the degradation of the puparium cuticle in preparation for eclosion of the adult(Cavener and MacIntyre 1983). Moreover, a knockdown of female Drosophila Gld expression decreased the amount of sperm stored in the spermathecae and caused a highly asymmetrical sperm distribution in the two spermathecal capsules, and this storage defect was especially severe when the mutant females were crossed with a Gld mutant male that had previously mated a few hours before the experimental cross, which showed that Gld facilitates sperm uptake and release through the spermathecal ducts (Iida and Cavener 2004). Together with the RNAi phenotypes from the Tchsp18.3 analyses, we showed that knockdown of Tchsp18.3 promoted and repressed the levels of Blimp-1 and Gld expression, respectively, and led to developmental and reproductive defects of T. castaneum. Furthermore, an odd family gene Drm encodes a zinc finger protein, and it is necessary for the development of the Drosophila embryonic small intestine, an anterior domain of the ectodermal hindgut(Uddin et al. 2011). In fly Drm mutants, the elongation of the embryonic hindgut was interrupted and it could not form the small intestine, which resulted in embryogenesis defects(Iwaki et al. 2001). In addition, most (91%) of the Kinesin family member mRNA levels changed abundantly throughout gametophyte development with 52% of Kinesin mRNAs enriched during the mitotic phase and 39% enriched during differentiation in Marsilea vestita. The lack of Kinesin-14 in Marsilea vestita led the nuclei of sperm cells to become elliptic and abnormal, which demonstrated that Kinesin-14 is essential for spermatogenesis (Tomei and Wolniak 2016). Indeed, Drm and Kinesin-14 were both decreased after RNAi of Tchsp18.3, which may be the reason why ds-Tchsp18.3 females were completely sterile and the males had low reproduction rates in T. castaneum. More importantly, RNAi against beetle Gld and Drm as well as Pthr2, Delta(11)-like and EGF-like domain protein 2 caused development arrest, although the extent of the lethality varied significantly according to the affected differential genes (Table 4; Appendix J). The loss of Drm, Kinesin-14 and Pthr2 at the T. castaneum pupal stage also induced embryonic defects (Table 4; Appendix J) (Donitz et al.2015). In general, these results suggested that Tchsp18.3 participated in development and reproduction in cooperation with Blimp-1, Gld, Drm, Kinesin-14, Pthr2, Delta(11)-like and EGF-like domain protein 2 in T. castaneum.

Fig. 2 The distribution of genes coverage of control (A) and ds-Tchsp18.3 (B). The percentages in pie chart and the numbers in parentheses indicate the percent of genes and the numbers of genes falling in this coverage range, respectively. The percentage interval on the right indicates the range of genes coverage. The gene coverage is the percentage of a gene covered by the reads,and this value is equal to the ratio of the base number in a gene covered by unique mapping reads to the total base number for that gene.

Fig. 3 Histogram of Gene Ontology (GO) classifications. The results are summarized in three main categories: biological process,cellular component and molecular function. The green- and red-colored numbers of the Y-axis on the right represent the number of differentially expressed genes (DEGs) and all of the detected genes, respectively.

3.3. ldenti fication and annotation of differentially expressed genes

Our previous study showed that Tchsp18.3 was involved in stress resistance, development and reproduction of the red flour beetles (Xie et al. 2017). Therefore, we used RNAsequencing technology to further identify the regulatory systems of Tchsp18.3 in T. castaneum. Knockdown of Tchsp18.3 in late larvae (20 days old) led to 312 DEGs being significantly up-regulated and 257 DEGs being significantly down-regulated (Fig. 4-A). Based on GO functional annotation and KEGG analysis, 246 DEGs were classified into 47 functional GO terms (Appendix H) and 282 DEGs were clustered in 147 biological pathways (Appendix I).The results provided fundamental clues for determining the functions and signaling regulatory mechanisms of Tchsp18.3 in T. castaneum.

Fig. 4 The change and validation of differentially expressed genes (DEGs). A, the number of genes up-regulated and down-regulated in ds-Tchsp18.3 relative to the control are indicated. B, a comparative analysis of the expression trends of the selected DEGs as determined by RNA sequencing and qRT-PCR. Ap, antimicrobial peptide; Duox, dual oxidase; Hr46,hormone receptor-like in 46; Z9desB, Z9 acyl coenzyme-A desaturase B; CuZnSOD, copper- and zinc-containing superoxide dismutase; Cpr62Bc, cuticular protein 62Bc; TwdlT,TweedleT; NiPp1, nuclear inhibitor of protein phosphatase 1;Cyp301a1, cytochrome P450 301a1; KCP, kielin/chordin-like protein; MCORP, multicopper oxidase relative protein; H115,serine protease H115.

GO function and KEGG pathway enrichment analyses were performed for all DEGs. According to GO functional analysis, 246 DEGs were assigned to three main GO functional categories and were then divided into 47 subcategories (P-value≤0.05) (Fig. 3). There were 24 subcategories for biological process, 14 for cellular components and 9 for molecular function (Appendix H). The metabolic process category (GO: 0008152) with 111 DEGs (45.87%)was the most prominent category for biological process. The cell part (GO: 0044464) and cell (GO: 0005623) categories both included 41 DEGs (16.94%) that were dominant in the main categories of cellular components. Catalytic activity(GO: 0003824) with 118 DEGs (45.87%) was dominant in the main category of molecular function (Fig. 3; Appendix H). A high percentage of DEGs were also classified as a response to stimulus (GO: 0050896) with 31 DEGs,a development process (GO: 0032502) with 18 DEGs,reproduction (GO: 0000003) with 9 DEGs and antioxidant activity (GO: 0016209) with 5 DEGs. Additionally, KEGG analyses showed that 282 DEGs were annotated in 147 biological pathways (Appendix I). Among them, 16 pathways were significantly enriched in DEGs compared to the whole-transcriptome background (Q-value≤0.05),including pancreatic secretion (ko04972), protein digestion and absorption (ko04974) as well as glycine, serine and threonine metabolism.

3.4. DEGs related to stress resistance, development and reproduction

In Drosophila, over-expressing hsp22 enhanced resistance to oxidative stress and extended lifespan by 32%, while absence of the expression of this small hsp22 led to a 40%decrease in lifespan and displayed a reduction of 30% in locomotor activity (Morrow et al. 2004). This result implied that Tchsp18.3, which may be similar to hsp22, is also involved in antioxidant activity and the aging of T. castaneum.The RNAi of Tchsp18.3 promoted the expression levels of Duox and CuZnSOD but repressed the Cat, Pex and Pexlike gene expression. These differentially expressed genes were simultaneously assigned to the response to stimulus(Table 2) and antioxidant activity (Table 3) GO terms.Reactive oxygen species (ROS), including superoxide anions and hydrogen peroxide produced in cell metabolism,have been shown to have a role in insect innate immunity as a potent pathogen-killing agent (Bahia et al. 2013; Deng and Zhao 2014). Under normal growth conditions, the Duox contributed to the generation of ROS, which maintained the homeostasis of microorganisms in the Drosophila gut (Kim and Lee 2014). However, higher ROS levels caused by over-expressed Duox also resulted in disruption of cellular function and structure as well as led to serious oxidative damage of DNA, RNA, proteins and lipids (Zhu et al. 2016).The residual/excessive ROS can be eliminated concurrently by antioxidant systems such as Cat, Pex and Pex-like, which aim to protect hosts against oxidative stress (Ahn et al. 2012;Zhao et al. 2013). Silencing of the Cat gene in the fourthinstar stage of Spodoptera litura larvae increased mortality because of induced ROS production, which suggested that ROS scavenging by Cat is important for S. litura survival(Zhao et al. 2013). Additionally, the over-expression and low-expression of Pex in Tetraodon nigroviridis differentially regulated the intracellular levels of ROS, and the latter could not reduce H2O2 levels (Sun et al. 2015). Meanwhile,qRT-PCR analysis showed that Pex-like was also induced in response to H2O2 stimulation in A. cerana cerana, and this result was con firmed by Western blot assay (Zhu et al.2016). The results indicated that Pex and Pex-like may play a crucial role in the antioxidant defense against ROS damage. Generally, the increases of CuZnSOD in vivo could also promote resistance to oxidative stress and the results reported that mRNA levels of the CuZnSOD gene were upregulated instantaneously after exposure to H2O2 treatment in Brachionus calyci florus (Van Raamsdonk and Hekimi 2009; Yang et al. 2013). These data exhibited that the loss of expression of Cat, Pex and Pex-like caused by reduction of Tchsp18.3 that was compensated for by up-regulation of CuZnSOD in response to excessive ROS, which further suggests that Tchsp18.3 participated in antioxidant stress responses through cross-talk with these antioxidant activityrelated genes.

4. Discussion

To further determine which of the 13 299 genes were significantly and differentially expressed between the ds-Tchsp18.3 and control samples, a likelihood ratio test was applied to compare the RPKM-derived read count (Marioni et al. 2008). We required a fold-change≥2 and FDR≤0.001,which screened out 569 DEGs. Compared to the control group, 312 DEGs were up-regulated and 257 DEGs were down-regulated in the ds-Tchsp18.3 sample (Fig. 4-A).Additionally, 16 DEGs were randomly selected for qRTPCR to validate the quality of RNA sequencing results. The results showed that expression patterns of these selected DEGs were similar between qRT-PCR and RNA sequencing(Fig. 4-B), which con firmed that RNA sequencing is a robust method to identity and quantify gene expressions.

Table 2 Differentially expressed genes (DEGs) involved in the stress response (significant at P-value≤0.05)

1) T, ds-Tchsp18.3; C, control. 2) FDR, false discovery rate. 3) Hsp, heat shock protein; PGRP2, peptidoglycan-recognition protein 2; Beta-1,3-gbp 2, beta-1,3-glucan-binding protein 2; Rttpk,receptor type tyrosine phosphatase kappa; Duox, dual oxidase; Z9desB, Z9 acyl-CoA desaturase B; CuZnSOD, copper- and zinccontaining superoxide dismutase; Pte-1, protein trapped in endoderm-1; FadD1, fatty acid desaturase D1; Delta(11)-like, acyl-CoA delta(11) desaturase-like; Pex-like, peroxidase-like; Drm, drumstick protein; Cat, catalase; Pthr2, parathyroid hormone-related peptide receptor. –, there is no protein being annotated. 4) –, there is no evidence of phenotype at present. 5) 11 days after pupal injection (includes death as pupa and adult). 6) 14 days after female pupal injection. 7) 11 days after larval injection (includes death as larva, prepupa and pupa).

Gene ID Log2Ratio (T/C)1)Regulation (T/C)1) FDR2) Protein3) Phenotype after RNAi of DEGs in Tribolium castaneum4)TC002669 8.09 Up 2.40E-77 Dusky –TC009126 4.41 Up 1.28E-09 Hsp68b –TC010517 4.36 Up 5.84E-153 Defensin2 –TC006250 3.40 Up 2.41E-97 Defensin1 –TC003541 3.34 Up 6.76E-05 Hsp22.2 –TC010172 3.14 Up 2.10E-66 Hsp70A1 –TC009706 2.94 Up 2.47E-06 Hsp68a –TC002386 2.73 Up 3.49E-04 Wnt-8a –TC010611 2.58 Up 3.08E-29 PGRP2 –TC007620 2.40 Up 2.32E-25 Beta-1,3-gbp 2 –TC006793 1.95 Up 1.69E-04 Hsp20.8 –TC009605 1.83 Up 2.04E-55 Rttpk –TC000521 1.74 Up 2.36E-14 – –TC004593 1.60 Up 1.72E-22 Duox –TC016415 1.24 Up 2.10E-20 Z9desB –TC011675 1.13 Up 6.11E-09 CuZnSOD –TC004716 1.05 Up 4.75E-05 Pte-1 –TC001409 1.02 Up 3.58E-06 Grapes –TC015383 –7.12 Down 1.06E-39 FadD1 –TC015349 –5.43 Down 2.00E-45 Delta(11)-like 20.0% lethality5)TC006593 –3.33 Down 7.81E-69 Hsp20.7b –TC011222 –3.25 Down 3.07E-106 Pex-like –TC011614 –2.97 Down 9.28E-36 Viperin 50.0% lethality5)TC006347 –2.53 Down 3.52E-05 Drm 10.0% lethality5) and more than 80.0% of embryonic and larval musculature patterns became irregular6)TC010105 –1.83 Down 2.22E-125 Hsp21.8 –TC005527 –1.54 Down 1.33E-18 Hsp20.7a –TC008096 –1.38 Down 8.62E-10 Dorsal –TC011090 –1.31 Down 3.33E-06 Cat –TC010267 –1.22 Down 2.27E-04 Pthr2 20.0% lethality5) and 30.0% of embryonic and larval musculature patterns became irregular6)TC004661 –1.21 Down 7.88E-49 Pex 10.0% lethality7) and 30.0% lethality5)TC013059 –1.06 Down 8.39E-50 Hdd11 –

4.1. The cross-talk exists between the immunity and stress response in T. castaneum

Typically, antimicrobial peptides (AMPs), such as cecropin and defensin, were over-expressed at 36°C and underexpressed at 15 and 32°C in Aedes aegypti larvae, which indicates there is an intimate relationship between the immune and stress responses (Muturi et al. 2012). This interaction between the innate immune and stress responses was also reported for the Drosophila Malpighian tubule(Davies et al. 2012). Septic wounding resulted in a stronginduction of both potential immune and stress genes in the red flour beetles, which further confirmed that there might be cross-talk between the immunity and stress response(Altincicek et al. 2008). More importantly, Tchsp18.3 increased in response to septic injury in T. castaneum,which suggests it plays a vital role in innate immunity(Altincicek et al. 2008). According to the GO functional annotations and references, 31 DEGs involved in a response to stimulus were identified (Table 2). Among these DEGs were viperin, dorsal, Hdd11, PGRP2, defensin1 and defensin2, which were also found to be associated with an immune response. Viperin has been shown to be active against a wide range of DNA and RNA viruses (Hee and Cresswell 2017), and when following the challenge with the viral mimic, its expression level was significantly induced in Branchiostoma japonicum (Lei et al. 2015). The amphioxus cells transfected with viperin can kill lymphocystis disease virus (LCDV) or inhibit its propagation, and co-incubation of viperin with white spot syndrome virus (WSSV) markedly attenuates its infectivity (Lei et al. 2015). The results indicated that viperin can promote resistance against viral infection by participating in immune system cell signaling transduction in vitro and in vivo. Meanwhile, the NF-κB/rel-related gene dorsal in Drosophila was expressed in larval and adult fat bodies, where its mRNA expression was noticeably enhanced after bacterial challenge (Reichhart et al. 1993; Lemaitre et al. 1995). The genes encoding inducible antibacterial peptides in fruit flies contain NF-κB/rel-related nucleotide sequences and showed that the dorsal protein can bind to such motifs and transactivate a reporter gene in co-transfection experiments with a fly cell line in a sequence-specific manner (Reichhart et al. 1993). The methuselah gene was also regulated by dorsal to participate in resistance to different kinds of stresses, including high temperatures, starvation and dietary paraquat in Drosophila(Kim et al. 2006). These findings implied that dorsal was related to the insect immune response as well as stress resistance. Additionally, the expression of the Hdd11 gene was undetectable in naive larvae when measured with Northern blot hybridization, but it was strongly induced in fat body after injecting bacterial cells or peptidoglycan into the larvae of Samia cynthia ricini (Bao et al. 2003), which indicated that constitutively increased induction of Hdd11 could be responsible for increased levels of immunity. As expected, our RNA-sequencing data showed that viperin,dorsal and Hdd11 were down-regulated in T. castaneum after RNAi of Tchsp18.3 (Table 2), which suggests that ds-Tchsp18.3 resulted in weakened immunity. Furthermore,PGRP2, defensin1 and defensin2 were up-regulated after RNAi of Tchsp18.3 in the red flour beetles (Table 2), while silencing of Drosophila PGRP-LB caused activation of the immune deficiency (IMD) signaling pathway and knocking down fly PGRP-SC led to significantly increased expression of diptericin and attacin (Wang J et al. 2009; Gao et al.2015), which suggests that PGRP2 probably inhibited the IMD signaling pathway and decreased the synthesis of AMPs as well as reduced the immunity of T. castaneum.Additionally, defensin1 was increased in response to bacterial challenges in Bombus terrestris and defensin2 was also induced following Cry3Ba spore-crystal treatment in thered flour beetles (Erler et al. 2011; Contreras et al. 2015).The over-expression of defensin1 and defensin2 are likely compensated for the loss of Tchsp18.3 and even viperin,dorsal and Hdd11 expression in insects, which maintained the homeostasis of microorganisms in vivo. Furthermore,PGRP2, defensin1 and defensin2 were also up-regulated in the fourth-instar larvae of A. aegypti when exposed to high temperature stress (Muturi et al. 2011), which suggests that innate immunity was tightly coupled with stress response. In this study, we offer the first evidence to show that cross-talk between the immunity and stress response was regulated by Tchsp18.3 and the decline in the immune response

might reduce the stress resistance of ds-Tchsp18.3 insects.However, knockdown of Tchsp18.3 can facilitate Tchsp20.8,Tchsp22.2, Tchsp68a, Tchsp68b and Tchsp70A1 gene expression as well as suppress Tchsp20.7a, Tchsp20.7b and Tchsp21.8 expression (Table 2), which shows that Tchsps do not act independently of each other but further study is needed to determine what kind of interactions do they have.

Table 3 Differentially expressed genes (DEGs) associated with antioxidant activity (significant at P-value≤0.05)

1) T, ds-Tchsp18.3; C, control. 2) FDR, false discovery rate. 3) Duox, dual oxidase; CuZnSOD, copper- and zinc-containing superoxide dismutase; Pex-like, peroxidase-like; Cat, catalase. 4) –, there is no evidence of phenotype at present. 5) 11 days after larval injection (includes death as larva, prepupa and pupa). 6) 11 days after pupal injection (includes death as pupa and adult).

Gene ID Log2Ratio (T/C)1) Regulation (T/C)1) FDR2) Protein3) Phenotype after RNAi of DEGs in Tribolium castaneum4)TC004593 1.60 Up 1.72E-22 Duox –TC011675 1.13 Up 6.11E-09 CuZnSOD –TC011222 –3.25 Down 3.07E-106Pex-like –TC011090 –1.31 Down 3.33E-06 Cat –TC004661 –1.21 Down 7.88E-49 Pex 10.0% lethality5) and 30.0% lethality6)

Table 4 Differentially expressed genes (DEGs) related to developmental processes (significant at P-value≤0.05)

1) T, ds-Tchsp18.3; C, control. 2) FDR, false discovery rate. 3) Blimp-1, B lymphocyte-induced maturation protein-1; Absent MD naos, absent MD neurons and olfactory sensilla; Gld-like, glucose dehydrogenase-like; Rttpk, receptor type tyrosine phosphatase kappa; Hr46, hormone receptor in 46-like; CuZnSOD, copper- and zinc-containing superoxide dismutase; Drm, drumstick protein; Pthr2, parathyroid hormone-related peptide receptor. 4) –, there is no evidence of phenotype at present. 5) 11 days after pupal injection (includes death as pupa and adult). 6) 14 days after female pupal injection. 7) 11 days after larval injection (includes death as larva, prepupa and pupa).

Gene ID Log2Ratio (T/C)1) Regulation (T/C)1) FDR2) Protein3) Phenotype after RNAi of DEGs in Tribolium castaneum4)TC002669 8.09 Up 2.40E-77 Dusky –TC000163 5.42 Up 2.40E-250 Trynity –TC014741 4.21 Up 7.15E-49 Blimp-1 –TC003170 3.73 Up 2.01E-05 Absent MD naos –TC001440 3.29 Up 9.08E-152 Dusky-like –TC002386 2.73 Up 3.49E-04 Wnt-8a –TC016355 2.38 Up 5.49E-67 Gld-like –TC006564 2.19 Up 5.61E-06 Multiple-wing-hairs –TC005267 1.97 Up 1.48E-16 Forked –TC009605 1.83 Up 2.04E-55 Rttpk –TC008909 1.53 Up 9.60E-42 Hr46 –TC011675 1.13 Up 6.11E-09 CuZnSOD –TC001409 1.02 Up 3.58E-06 Grapes –TC011323 1.01 Up 3.20E-09 Stumps –TC009212 –4.40 Down 1.39E-124 EGF-like domain protein 2 10.0% lethality5)TC006347 –2.53 Down 3.52E-05 Drm 10.0% lethality5) and more than 80.0% embryonic and larval musculature patterns became irregular6)TC008096 –1.38 Down 8.62E-10 Dorsal –TC010267 –1.22 Down 2.27E-04 Pthr2 20.0% lethality7) and 30.0% of embryonic and larval musculature patterns became irregular6)

Fig. 5 Analyses of KEGG biological pathway enrichment with Q-values≤0.05 according to RNA sequencing. The numbers of the Y-axis on the right (e.g., 0, 0 and 282, 9 708) represent the number of differentially expressed genes (DEGs) and all of the detected genes, respectively.

Table 5 Differentially expressed genes (DEGs) associated with pancreatic secretion pathways (significant at Q-value≤0.05)

1) T, ds-Tchsp18.3; C, control. 2) FDR, false discovery rate. 3) SP, serine protease; CipA2-gamma, calcium-independent phospholipase A2-gamma; Plrp2-like, pancreatic lipase-related protein 2-like. –, there is no protein being annotated. 4) –, there is no evidence of phenotype at present. 5) 11 days after pupal injection (includes death as pupa and adult). 6) 14 days after female pupal injection. 7) 11 days after larval injection (includes death as larva, prepupa and pupa).

Gene ID Log2Ratio (T/C)1) Regulation (T/C)1) FDR2) Protein3) Phenotype after RNAi of DEGs in Tribolium castaneum4)TC005904 9.03 Up 1.61E-150 Carboxypeptidase A –TC008659 5.20 Up 2.33E-274 SP P87 –TC005905 4.96 Up 0 Carboxypeptidase B –TC030071 4.49 Up 3.08E-82 Notopleural –TC001417 4.12 Up 2.80E-11 Lipase member H-B-like –TC030609 4.09 Up 9.49E-04 SP stubble protein –TC008657 4.00 Up 2.06E-30 SP P84 –TC0152453.59Up0– –TC004622 2.75 Up 2.36E-32 SP H51 –TC013419 2.09 Up 8.72E-06 Phospholipase A2C –TC000548 1.73 Up 4.98E-05 SP-like protein –TC000249 1.41 Up 9.14E-47 SP H4 –TC000495 1.35 Up 2.14E-23 SP P8 –TC009386 1.01 Up 5.69E-13 CipA2-gamma –TC015181 –3.46 Down 9.91E-09 Phospholipase A2E 10.0% lethality5) and 30.0% of embryos without visible larval cuticle6)TC016121 –2.99 Down 3.45E-94 SP P167 20.0% lethality7)TC008554 –2.42 Down 3.58E-22 SP H82 –TC011935 –2.20 Down 3.10E-05 Plrp2-like 10.0% lethality5)TC010934 –1.90 Down 4.63E-20 SP H115 60.0% lethality7)TC015237 –1.42 Down 2.32E-33 SP P151 60.0% lethality5)TC014083 –1.39 Down 5.57E-17 SP P147 30.0% lethality5) and 30.0–50.0%embryonic and larval musculature patterns became irregular6)TC002767 –1.37 Down 8.84E-08 SP P38 10.0% lethality5)TC030074 –1.36 Down 1.32E-12 SP 9 20.0% lethality5)TC004027 –1.31 Down 5.21E-16 Carboxypeptidase A –TC005130 –1.30 Down 1.45E-05 SP P142 10.0% lethality5) and 30.0% of embryonic and larval musculature patterns became irregular6)TC004657 –1.09 Down 4.41E-16 Phospholipase A1-like 10.0% lethality7)TC005133 –1.03 Down 3.14E-06 Esterase FE4-like 20.0% lethality7) and more than 80.0% embryos lacked visible larval cuticle6)TC004026 –1.03 Down 2.32E-19 Carboxypeptidase A 30.0% lethality7)TC008504 –1.01 Down 1.91E-98 SP P80 30.0% lethality5) and more than 80.0% of embryonic and larval musculature patterns became irregular6)

4.2. Tchsp18.3 participated in development and reproduction by regulating antioxidant activity and other physiological processes

Multiple GO functional groups (Fig. 3), such as a response to stimulus (GO: 0050896), antioxidant activity (GO: 0016209),developmental process (GO: 0032502) and reproduction(GO: 0000003), were inextricably associated with Tchsp18.3 functions. The response to stimulus (GO: 0050896) category contains 31 DEGs (Table 2), which include viperin, dorsal,Hdd11, Peptidoglycan recognition protein 2 (PGRP2),defensin1 and defensin2 genes that are not only crucial for the stress response but also innate immunity. Dual oxidase(Duox), Copper and Zinc-containing Superoxide Dismutase(CuZnSOD), Catalase (Cat), Peroxidase-like (Pex-like)and Pex were simultaneously assigned to the response to stimulus (GO: 0050896) (Table 2) and antioxidant activity(GO: 0016209) categories (Table 3), which were related to the antioxidant stress response and probably contributed to the development and reproduction of T. castaneum.We also found that there were 18 DEGs and 9 DEGs that participated in the developmental process (GO: 0032502)and reproduction (GO: 0000003), respectively (Table 4;Appendix J). Those DEGs included B lymphocyte-induced maturation protein-1 (Blimp-1), Glucose dehydrogenase(Gld), Drumstick protein (Drm), Kinesin-14, Parathyroid hormone-related peptide receptor 2 (Pthr2), acyl-CoA Delta(11) desaturase-like (Delta(11)-like) and EGF-like domain protein 2 genes that play significant roles in insects eclosion and fertility. Additionally, in the KEGG pathways(Fig. 5), there were 29 DEGs and 24 DEGs annotated as pancreatic secretion (ko04972) (Table 5) as well as protein digestion and absorption (ko04974) (Appendix K)processes, respectively. Twenty-one of these DEGs are the same in these two processes with 16 DEGs that are members of the serine protease family, which are likely responsible for the stress response and innate immunity as well as development and reproduction of the T. castaneum.The results showed that there are considerable differences in the physiological processes between the ds-Tchsp18.3 and control groups.

In addition to the role of lifespan regulation, antioxidant activity has been found to be involved in development and reproduction (Taracena et al. 2015; Deng et al. 2016).Silencing of Cat gene expression in Rhodnius prolixus resulted in an increase in mortality rates of 24–30% in the first-instars, molting rates reduced by 80% in third-instars,oviposition rates decreased by 84% and embryogenesis arrested in 20% of laid eggs (Taracena et al. 2015). The results also showed that loss of the Cat gene in female Lutzomyia longipalpis led to a significant increase in mortality and a reduction in developing oocytes (Diaz-Albiter et al. 2011). The Sitobion avenae nymphs were continuously exposed to dietary dsRNA of Pex through an artificial diet that reduced the survival rate and the ecdysis index revealed that Pex is important to maintain the growth and development of S. avenae (Deng et al.2016). In Caenorhabditis elegans, extracellular Pex acted in embryonic development and in curbing the regeneration of adult axons (Gotenstein et al. 2010). Moreover, depletion of Pex at the T. castaneum larval and pupal stages gave rise to 10 and 30% lethality, respectively (Tables 2 and 3)(Donitz et al. 2015). Consistent with our previous study,RNAi of Tchsp18.3 caused development abnormities before adulthood and reduced egg production in T. castaneum (Xie et al. 2017). Therefore, these results suggest that the effects of ds-Tchsp18.3 on the development and reproduction of the red flour beetles might be consequences of decreased antioxidant activity.

To facilitate the transcription analysis of gene expression, all the detected beetle genes were assigned to three function categories (i.e., biological process, molecular function and cellular component). Of the 13 299 detected genes, 6 831 were classified into 25, 14 and 10 functional groups of the three main categories, respectively, based on sequence homology (Fig. 3). It was clear that the dominant distributions are from the cellular process category (GO: 0009987) under biological process with 3 678 genes (53.84%), followed by the cell part (GO: 0044464) and cell (GO: 0005623)categories under cellular component with 3 183 genes(46.60%), and the binding (GO: 0005488) category under molecular function with 3 388 genes (49.60%) (Appendix H). In addition, we noticed that some genes were annotated in response to stimulus (GO: 0050896) and antioxidant activity (GO: 0016209) as well as development processes (GO: 0032502) and reproduction (GO: 0000003)(Fig. 3). To investigate biological behavior further, 13 299 genes were mapped to the reference canonical pathways in the KEGG, and KEGG annotation was obtained for 9 708 of them and enriched with 147 pathways. The pathways with the greatest representation included metabolic pathways(1 537 genes, ko01100), pathways in cancer (293 genes,ko05200), ubiquitin-mediated proteolysis (289 genes,ko04120), purine metabolism (284 genes, ko00230) and regulation of the actin cytoskeleton (270 genes, ko04810)(Appendix I).

4.3. Knockdown of Tchsp18.3 amplified the serine protease signaling pathway

Furthermore, combining the KEGG biological pathways(Fig. 5) for pancreatic secretion (Table 5) as well as protein digestion and absorption (Appendix K), the 9 members of the serine protease (SP) family were down-regulated, but 7 members were up-regulated after RNAi of Tchsp18.3.In addition to digestive functions, SPs have significant roles in multiple physiological processes, including stress responses and innate immunity as well as development and reproduction (Gorman et al. 2000; Zou et al. 2006; Wang et al. 2016). SP3, SPH42 and SP49 transcript levels were increased after the saline or Escherichia coli injection, and the pathogen challenge also led to a much stronger induction of SP1, SP2, SP6 and SP41 gene transcription in A. mellifera(Zou et al. 2006). Furthermore, ablation of larvae SPP167 and SPH155 resulted in development arrest with 20%and 60% lethality (includes deaths as larvae and pupae),respectively, in the red flour beetles (Table 5; Appendix K) (Donitz et al. 2015). While loss of SPP151, SPP147,SPP38, SP9, SPP142 and SPP80 at the T. castaneum pupal stage gave rise to morphological deficiency with 60,30, 10, 20, 10 and 30% lethality (includes deaths as pupae and adult), respectively (Table 5; Appendix K) (Donitz et al.2015). Additionally, we found that 30–50%, 30% and over 80% of the beetle embryonic musculature patterns became irregular after female pupae injection of ds-TcSPP147, ds-TcSPP142 and ds-TcSPP80, respectively (Table 5; Appendix K) (Donitz et al. 2015), which was consistent with silencing of SPP2c in silkworm embryos that significantly reduced the degradation rate of residual yolk proteins and further suppressed embryogenesis (Wang et al. 2016). The results indicated that the SP signaling pathways were activated and amplified by Tchsp18.3, which subsequently affected the stress and immune responses as well as development and reproduction of T. castaneum. Interestingly, the lack of these SPs in beetles is probably compensated by increasing the expression of the other major SPs (Table 5; Appendix K), at least partially. The result was similar to outcomes with SP inhibitors in which the expression of some SPs was reduced,whereas other SPs expression levels were significantly increased as part of an integrated compensation response in T. castaneum (Oppert et al. 2010; Perkin et al. 2017).Overall, these results implied that knockdown of Tchsp18.3 changed large-scale SP expression patterns and amplified the SP signaling pathways to further regulate the beetle’s diverse physiological functions.

5. Conclusion

To validate RNA sequencing data, total RNAs were extracted from six larvae of T. castaneum on the fifth day after injection of buffer and ds-Tchsp18.3 as described above. Additionally,1 μg of total RNA was converted to cDNAs by HiScript II Reverse Transcriptase (Vazyme, Nanjing, China) (Li et al.2017). Sixteen DEGs were randomly selected to check the expression patterns between the control and ds-Tchsp18.3 groups with qRT-PCR, which was performed with Power SYBR Green Master Mix (Roche) by the Step One Plus Real-Time PCR System (Applied Biosystems, Life Technologies,USA) and programmed at 95°C for 10 min, followed by 40 cycles of 95°C for 15 s and 60°C for 1 min, and finally run at 95°C for 15 s, 60°C for 1 min and 95°C for 15 s. The data are expressed as relative mRNA levels that were standardized to a control gene, T. castaneum ribosomal protein S3 (Tcrps3) through using the ΔΔCT method (Livak and Schmittgen 2001). Primers for these DEGs are shown in Appendix A. All the data were analyzed using the one way analysis of variance (ANOVA) program in SPSS version 13.0 (SPSS Inc., Chicago, IL, USA) and are presented as the means±standard error.

“医疗风险无处不在”已经是世界医疗界的共识。新技术和新方法的出现，医生的决断能力，环境变得复杂多样都会给医疗带来误诊，引起事故[19-21]。目前教研室使用医学仿真模拟人，建立了完善的科学的培训体系，设置了相关内容，培养了学生的疾病诊断能力，以患者为中心，为患者着想，防止患者受到精神或身体的损伤，杜绝了教学的不必要纠纷。学生进入病房前，都接受过了医学模拟人的技能培训，考核通过。所以，医学仿真模拟人对医学生的自信心建立，医疗技能提升都有非常大的帮助。

Acknowledgements

This work was supported by the National Natural Science Foundation of China (31572326 and 31172146) and the Priority Academic Program Development (PAPD) of Jiangsu Higher Education Institutions, China.

根据大量研究显示，虽然村级财务管理工作受到了广大村民的热切关注，各级政府要花费大量的时间与人力和物力推动村务财务管理制度透明化建设。但是一些自然村仍旧受到计划思想的困扰，在村级财务管理方面暴露了大量的问题。不能够有效实现财务管理透明化。在村级管理中记账审核制度不够健全，由于农村财务管理人员普遍专业素养较低，对于村级财务工作不够熟悉，缺乏相应的原则性以及责任心，很容易导致村一级的开支失衡，在一些会计凭证处理上违反国家相关规定。同时又由于农村村级财务管理工作公开的不够具体和详细，特别是牵扯到一些实质性内容并未公开，只是公开一些细枝末节的不是很重要的财务报表。

Appendices associated with this paper can be available on http://www.ChinaAgriSci.com/V2/En/appendix.htm

本文所关注的解释变量是董事会结构，涉及四个方面的评估。董事会规模通常采用董事会人数进行度量，董事会独立性采用独立董事比例度量，两职分离和审计委员会的设立均为虚拟变量，若等于1则表示两职分离和设立审计委员会，两职兼任和未设立审计委员会时则等于0。

建构主义是学习理论中行为主义发展到认知主义以后的进一步发展，最早的提出者可追溯到瑞士心理学家让·皮亚杰（J.Piaget）和美国教育心理学家杰罗姆·布鲁纳(Jerome Seymour Bruner)。70年代末，以布鲁纳为首的美国教育心理学家将苏联教育心理学家维果斯基的思想介绍到美国，极大地推动了建构主义思想的发展。建构主义目前还未形成稳定的、统一的体系,其分类存在不同的方法。国内学者陈琦等归纳了六种不同倾向的建构主义：激进建构主义、社会性建构主义、社会文化认知的观点、信息加工的建构主义、社会建构论和控制论系统。

3.2.2 完善天然气管网。加快推进上海东侧天然气主干管网建设，尽早打通沪苏北通道。将宝钢专线并入全市天然气主干网，进一步提高上海东北区域的供气能力。

References

Aevermann B D, Waters E R. 2008. A comparative genomic analysis of the small heat shock proteins in Caenorhabditis elegans and briggsae. Genetica, 133, 307–319.

Agawa Y, Sarhan M, Kageyama Y, Akagi K, Takai M, Hashiyama K, Wada T, Handa H, Iwamatsu A, Hirose S, Ueda H. 2007.Drosophila Blimp-1 is a transient transcriptional repressor that controls timing of the ecdysone-induced developmental pathway. Molecular and Cellular Biology, 27, 8739–8747.

Ahn H M, Lee K S, Lee D S, Yu K. 2012. JNK/FOXO mediated PeroxiredoxinV expression regulates redox homeostasis during Drosophila melanogaster gut infection.Developmental and Comparative Immunology, 38, 466–473.

Altincicek B, Knorr E, Vilcinskas A. 2008. Beetle immunity:Identification of immune-inducible genes from the model insect Tribolium castaneum. Developmental and Comparative Immunology, 32, 585–595.

Arrigo A P. 2000. sHsp as novel regulators of programmed cell death and tumorigenicity. Pathologie Biologie (Paris),48, 280–288.

Azad P, Zhou D, Russo E, Haddad G G. 2009. Distinct mechanisms underlying tolerance to intermittent and constant hypoxia in Drosophila melanogaster. PLoS ONE,4, e5371.

Bahia A C, Oliveira J H, Kubota M S, Araujo H R, Lima J B,Rios-Velasquez C M, Lacerda M V, Oliveira P L, Traub-Cseko Y M, Pimenta P F. 2013. The role of reactive oxygen species in Anopheles aquasalis response to Plasmodium vivax infection. PLoS ONE, 8, e57014.

Bakthisaran R, Tangirala R, Rao C M. 2015. Small heat shock proteins: Role in cellular functions and pathology.Biochimica Biophysica Acta, 1854, 291–319.

Bao Y, Mega K, Yamano Y, Morishima I. 2003. cDNA cloning and expression of bacteria-induced Hdd11 gene from erisilkworm, Samia cynthia ricini. Comparative Biochemistry Physiology (Part C: Toxicology Pharmacolgy), 136,337–342.

Cavener D R, MacIntyre R J. 1983. Biphasic expression and function of glucose dehydrogenase in Drosophila melanogaster. Proceedings of the National Academy of Sciences of the United States of America, 80, 6286–6288.

Cavener D R, Otteson D C, Kaufman T C. 1986. A rehabilitation of the genetic map of the 84B-D region in Drosophila melanogaster. Genetics, 114, 111–123.

Chen X, Xiong W, Li C, Gao S, Song X, Wu W, Li B. 2016.Comparative RNA-sequencing pro filing reveals novel Deltaclass glutathione S-transferases relative genes expression patterns in Tribolium castaneum. Gene, 593, 13–20.

Chen X, Zhang Y. 2015. Identification of multiple small heatshock protein genes in Plutella xylostella (L.) and their expression pro files in response to abiotic stresses. Cell Stress Chaperones, 20, 23–35.

Conesa A, Gotz S, Garcia-Gomez J M, Terol J, Talon M,Robles M. 2005. Blast2GO: A universal tool for annotation,visualization and analysis in functional genomics research.Bioinformatics, 21, 3674–3676.

Contreras E, Benito-Jardon M, Lopez-Galiano M J, Real M D,Rausell C. 2015. Tribolium castaneum immune defense genes are differentially expressed in response to Bacillus thuringiensis toxins sharing common receptor molecules and exhibiting disparate toxicity. Developmental and Comparative Immunology, 50, 139–145.

Contreras E, Rausell C, Real M D. 2013. Proteome response of Tribolium castaneum larvae to Bacillus thuringiensis toxin producing strains. PLoS ONE, 8, e55330.

Csermely P, Schnaider T, Soti C, Prohaszka Z, Nardai G.1998. The 90-kDa molecular chaperone family: structure,function, and clinical applications. A comprehensive review.Pharmacology Therapeutics, 79, 129–168.

Davies S A, Overend G, Sebastian S, Cundall M, Cabrero P,Dow J A, Terhzaz S. 2012. Immune and stress response‘cross-talk’ in the Drosophila Malpighian tubule. Journal of Insect Physiology, 58, 488–497.

Deng F, He Q, Zhao Z. 2016. Suppressing a peroxidase gene reduces survival in the wheat aphid Sitobion avenae.Archives of Insect Biochemistry Physiology, 93, 86–95.

Deng F, Zhao Z. 2014. Influence of catalase gene silencing on the survivability of Sitobion avenae. Archives of Insect Biochemistry Physiology, 86, 46–57.

Diaz-Albiter H, Mitford R, Genta F A, Sant’Anna M R, Dillon R J. 2011. Reactive oxygen species scavenging by catalase is important for female Lutzomyia longipalpis fecundity and mortality. PLoS ONE, 6, e17486.

Donitz J, Schmitt-Engel C, Grossmann D, Gerischer L, Tech M, Schoppmeier M, Klingler M, Bucher G. 2015. iBeetle-Base: A database for RNAi phenotypes in the red flour beetle Tribolium castaneum. Nucleic Acids Research, 43,D720–D725.

Erler S, Popp M, Lattorff H M. 2011. Dynamics of immune system gene expression upon bacterial challenge and wounding in a social insect (Bombus terrestris). PLoS ONE, 6, e18126.

Fu X, Jiao W, Chang Z. 2006. Phylogenetic and biochemical studies reveal a potential evolutionary origin of small heat shock proteins of animals from bacterial class A. Journal of Molecular Evolution, 62, 257–266.

Gao Y, Tang T, Gu J, Sun L, Gao X, Ma X, Wang X, Liu F,Wang J. 2015. Downregulation of the Musca domestica peptidoglycan recognition protein SC (PGRP-SC) leads to overexpression of antimicrobial peptides and tardy pupation.Molecular Immunology, 67, 465–474.

Garcia-Reina A, Rodriguez-Garcia M J, Ramis G, Galian J. 2017. Real-time cell analysis and heat shock protein gene expression in the TcA Tribolium castaneum cell line in response to environmental stress conditions. Insect Science, 24, 358–370.

Gorman M J, Andreeva O V, Paskewitz S M. 2000. Sp22D: A multidomain serine protease with a putative role in insect immunity. Gene, 251, 9–17.

Gotenstein J R, Swale R E, Fukuda T, Wu Z, Giurumescu C A,Goncharov A, Jin Y, Chisholm A D. 2010. The C. elegans peroxidasin PXN-2 is essential for embryonic morphogenesis and inhibits adult axon regeneration. Development, 137,3603–3613.

Hao X, Zhang S, Timakov B, Zhang P. 2007. The Hsp27 gene is not required for Drosophila development but its activity is associated with starvation resistance. Cell Stress Chaperones, 12, 364–372.

Haslbeck M, Franzmann T, Weinfurtner D, Buchner J. 2005.Some like it hot: The structure and function of small heatshock proteins. Nature Structural Molecular Biology, 12,842–846.

Hee J S, Cresswell P. 2017. Viperin interaction with mitochondrial antiviral signaling protein (MAVS) limits viperin-mediated inhibition of the interferon response in macrophages. PLoS ONE, 12, e0172236.

Iida K, Cavener D R. 2004. Glucose dehydrogenase is required for normal sperm storage and utilization in female Drosophila melanogaster. Journal of Experimental Biology,207, 675–681.

Iwaki D D, Johansen K A, Singer J B, Lengyel J A. 2001.Drumstick, bowl, and lines are required for patterning and cell rearrangement in the Drosophila embryonic hindgut.Developmental Biology, 240, 611–626.

Kawasaki F, Koonce N L, Guo L, Fatima S, Qiu C, Moon M T,Zheng Y, Ordway R W. 2016. Small heat shock proteins mediate cell-autonomous and -nonautonomous protection in a Drosophila model for environmental-stress-induced degeneration. Disease Models & Mechanisms, 9, 953–964.

Kegel K B, Iwaki A, Iwaki T, Goldman J E. 1996. AlphaB-crystallin protects glial cells from hypertonic stress.American Journal of Physiologly, 270, C903-C909.

Kim H, Kim J, Lee Y, Yang J, Han K. 2006. Transcriptional regulation of the methuselah gene by dorsal protein in Drosophila melanogaster. Molecules and Cells, 21,261–268.

Kim K K, Kim R, Kim S H. 1998. Crystal structure of a small heat-shock protein. Nature, 394, 595–599.

Kim S H, Lee W J. 2014. Role of DUOX in gut inflammation:Lessons from Drosophila model of gut-microbiota interactions. Frontiers in Cellular and Infection Microbiology,3, 116.

Lee G J, Pokala N, Vierling E. 1995. Structure and in vitro molecular chaperone activity of cytosolic small heat shock proteins from pea. Journal of Biological Chemistry, 270,10432–10438.

Lei M, Liu H, Liu S, Zhang Y, Zhang S. 2015. Identification and functional characterization of viperin of amphioxus Branchiostoma japonicum: Implications for ancient origin of viperin-mediated antiviral response. Developmental and Comparative Immunology, 53, 293–302.

Lemaitre B, Meister M, Govind S, Georgel P, Steward R,Reichhart J M, Hoffmann J A. 1995. Functional analysis and regulation of nuclear import of dorsal during the immune response in Drosophila. The EMBO Journal, 14, 536–545.

Li C, Wu W, Sang M, Liu X, Hu X, Yun X, Li B. 2014.Comparative RNA-sequencing analysis of mthl1 functions and signal transductions in Tribolium castaneum. Gene,547, 310–318.

Li C J, Yun X P, Yu X J, Li B. 2017. Functional analysis of the circadian clock gene timeless in Tribolium castaneum.Insect Science, doi: 10.1111/1744-7917.12441

Li R, Yu C, Li Y, Lam T W, Yiu S M, Kristiansen K, Wang J.2009. SOAP2: An improved ultrafast tool for short read alignment. Bioinformatics, 25, 1966–1967.

Li Z W, Li X, Yu Q Y, 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 Evolution Biology, 9, 215.

Liang D, Benko Z, Agbottah E, Bukrinsky M, Zhao R Y. 2007.Anti-Vpr activities of heat shock protein 27. Molecular Medicine, 13, 229–239.

Lindquist S, Craig E A. 1988. The heat-shock proteins. Annual Review of Genetics, 22, 631–677.

Liu B, Jiang G, Zhang Y, Li J, Li X, Yue J, Chen F, Liu H, Li H, Zhu S, Wang J, Ran C. 2011. Analysis of transcriptome differences between resistant and susceptible strains of the citrus red mite Panonychus citri (Acari: Tetranychidae).PLoS ONE, 6, e28516.

Livak K J, Schmittgen T D. 2001. Analysis of relative gene expression data using real-time quantitative PCR and the 2-ΔΔCT method. Methods, 25, 402–408.

Marioni J C, Mason C E, Mane S M, Stephens M, Gilad Y.2008. RNA-seq: An assessment of technical reproducibility and comparison with gene expression arrays. Genome Research, 18, 1509–1517.

Martinez-Paz P, Morales M, Martin R, Martinez-Guitarte J L,Morcillo G. 2014. Characterization of the small heat shock protein Hsp27 gene in Chironomus riparius (Diptera) and its expression pro file in response to temperature changes and xenobiotic exposures. Cell Stress Chaperones, 19,529–540.

Michaud S, Marin R, Westwood J T, Tanguay R M. 1997. Cellspecific expression and heat-shock induction of Hsps during spermatogenesis in Drosophila melanogaster. Journal of Cell Science, 110, 1989–1997.

Van Montfort R, Slingsby C, Vierling E. 2001. Structure and function of the small heat shock protein/alpha-crystallin family of molecular chaperones. Advances in Protein Chemistry, 59, 105–156.

van Montfort R L M, Basha E, Friedrich K L, Slingsby C, Vierling E. 2001. Crystal structure and assembly of a eukaryotic small heat shock protein. Nature Structural Biology, 8,1025–1030.

Morrow G, Battistini S, Zhang P, Tanguay R M. 2004. Decreased lifespan in the absence of expression of the mitochondrial small heat shock protein Hsp22 in Drosophila. Journal of Biological Chemistry, 279, 43382–43385.

Mortazavi A, Williams B A, McCue K, Schaeffer L, Wold B.2008. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nature Methods, 5, 621–628.

Muturi E J, Kim C H, Alto B W, Berenbaum M R, Schuler M A. 2011. Larval environmental stress alters Aedes aegypti competence for Sindbis virus. Tropical Medicine &International Health, 16, 955–964.

Muturi E J, Nyakeriga A, Blackshear M. 2012. Temperaturemediated differential expression of immune and stressrelated genes in Aedes aegypti larvae. Journal of the American Mosquito Control Association, 28, 79–83.

Nakamoto H, Vigh L. 2007. The small heat shock proteins and their clients. Cellular and Molecular Life Sciences, 64,294–306.

Oppert B, Elpidina E N, Toutges M, Mazumdar-Leighton S.2010. Microarray analysis reveals strategies of Tribolium castaneum larvae to compensate for cysteine and serine protease inhibitors. Comparative Biochemistry Physiology(Part D: Genomics Proteomics), 5, 280–287.

Pandey A, Saini S, Khatoon R, Sharma D, Narayan G, Kar Chowdhuri D. 2016. Overexpression of hsp27 rescued neuronal cell death and reduction in life- and health-span in Drosophila melanogaster against prolonged exposure to dichlorvos. Molecular Neurobiology, 53, 3179–3193.

Perkin L C, Elpidina E N, Oppert B. 2017. RNA interference and dietary inhibitors induce a similar compensation response in Tribolium castaneum larvae. Insect Molecular Biology,26, 35–45.

Ponnuvel K M, Murthy G N, Awasthi A K, Rao G, Vijayaprakash N B. 2010. Differential gene expression during early embryonic development in diapause and non-diapause eggs of multivoltine silkworm Bombyx mori. Indian Journal of Experimintal Biology, 48, 1143–1151.

Van Raamsdonk J M, Hekimi S. 2009. Deletion of the mitochondrial superoxide dismutase sod-2 extends lifespan in Caenorhabditis elegans. PLoS Genetics, 5, e1000361.

Reichhart J M, Georgel P, Meister M, Lemaitre B, Kappler C,Hoffmann J A. 1993. Expression and nuclear translocation of the rel/NF-kappa B-related morphogen dorsal during the immune response of Drosophila. Comptes Rendus Academie Sciences III, 316, 1218–1224.

Ritossa F M, Vonborstel R C. 1964. Chromosome puffs in Drosophila induced by ribonuclease. Science, 145,513–514.

Sang M, Li C, Wu W, Li B. 2016. Identification and evolution of two insulin receptor genes involved in Tribolium castaneum development and reproduction. Gene, 585, 196–204.

Sang W, Ma W H, Qiu L, Zhu Z H, Lei C L. 2012. The involvement of heat shock protein and cytochrome P450 genes in response to UV-A exposure in the beetle Tribolium castaneum. Journal of Insect Physiology, 58, 830–836.

Singh M K, Tiwari P K. 2016. Cloning & sequence identification of Hsp27 gene and expression analysis of the protein on thermal stress in Lucilia cuprina. Insect Science, 23,555–568.

Song X, Huang F, Liu J, Li C, Gao S, Wu W, Zhai M, Yu X, Xiong W, Xie J, Li B. 2017. Genome-wide DNA methylomes from discrete developmental stages reveal the predominance of non-CpG methylation in Tribolium castaneum. DNA Research, 24, 445–457.

Storey J D, Tibshirani R. 2003. Statistical significance for genomewide studies. Proceedings of the National Academy of Sciences of the United States of America,100, 9440–9445.

Sultan A R, Oish Y, Ueda H. 2014. Function of the nuclear receptor FTZ-F1 during the pupal stage in Drosophila melanogaster. Development Growth & Differentiation, 56,245–253.

Sun C C, Dong W R, Zhao J, Nie L, Xiang L X, Zhu G, Shao J Z. 2015. Cysteine-independent catalase-like activity of vertebrate Peroxiredoxin 1 (Prx1). Journal of Biological Chemistry, 290, 19942–19954.

Taracena M L, Oliveira P L, Almendares O, Umana C,Lowenberger C, Dotson E M, Paiva-Silva G O, Pennington P M. 2015. Genetically modifying the insect gut microbiota to control Chagas disease vectors through systemic RNAi.PLoS Neglected Tropical Diseases, 9, e0003358.

Tomei E J, Wolniak S M. 2016. Transcriptome analysis reveals a diverse family of kinesins essential for spermatogenesis in the fern Marsilea. Cytoskeleton (Hoboken), 73, 145–159.

Tsvetkova N M, Horvath I, Torok Z, Wolkers W F, Balogi Z,Shigapova N, Crowe L M, Tablin F, Vierling E, Crowe J H,Vigh L. 2002. Small heat-shock proteins regulate membrane lipid polymorphism. Proceedings of the National Academy of Sciences of the United States of America, 99, 13504–13509.

Uddin S N, Yano M, Murakami R. 2011. The drumstick gene acts cell-non-autonomously and triggers specification of the small intestine in the Drosophila hindgut. International Journal of Developmental Biology, 55, 945–952.

Wang D, Zhang Y, Dong Z, Guo P, Ma S, Guo K, Xia Q, Zhao P.2016. Serine protease P-IIc is responsible for the digestion of yolk proteins at the late stage of silkworm embryogenesis.Insect Biochemistry Molecular Biology, 74, 42–49.

Wang J, Wu Y, Yang G, Aksoy S. 2009. Interactions between mutualist Wigglesworthia and tsetse peptidoglycan recognition protein (PGRP-LB) in fluence trypanosome transmission. Proceedings of the National Academy of Sciences of the United States of America, 106, 12133–12138.

Wang L, Feng Z, Wang X, Zhang X. 2010. DEGseq: An R package for identifying differentially expressed genes from RNA-seq data. Bioinformatics, 26, 136–138.

Wang Z, Gerstein M, Snyder M. 2009. RNA-Seq: A revolutionary tool for transcriptomics. Nature Reviews Genetics, 10,57–63.

Waters E R, Aevermann B D, Sanders-Reed Z. 2008.Comparative analysis of the small heat shock proteins in three angiosperm genomes identi fies new subfamilies and reveals diverse evolutionary patterns. Cell Stress Chaperones, 13, 127–142.

Waters E R, Rio florido I. 2007. Evolutionary analysis of the small heat shock proteins in five complete algal genomes.Journal of Molecular Evolution, 65, 162–174.

Xie J, Hu X X, Zhai M F, Yu X J, Song X W, Gao S S, Wu W, Li B. 2017. Characterization and functional analysis of hsp18.3 gene in the red flour beetle, Tribolium castaneum. Insect Science, doi: 10.1111/1744-7917.12543

Xiong W, Zhai M, Yu X, Wei L, Mao J, Liu J, Xie J, Li B.2017. Comparative RNA-sequencing analysis of ER-based HSP90 functions and signal pathways in Tribolium castaneum. Cell Stress Chaperones, 23, 29–43.

Xu J, Tan A, Palli S R. 2010. The function of nuclear receptors in regulation of female reproduction and embryogenesis in the red flour beetle, Tribolium castaneum. Journal of Insect Physiology, 56, 1471–1480.

Yang J, Dong S, Jiang Q, Si Q, Liu X. 2013. Characterization and expression of cytoplasmic copper/zinc superoxide dismutase (CuZnSOD) gene under temperature and hydrogen peroxide (H2O2) in rotifer Brachionus calyci florus.Gene, 518, 388–396.

Ye J, Fang L, Zheng H, Zhang Y, Chen J, Zhang Z, Wang J,Li S, Li R, Bolund L. 2006. WEGO: A web tool for plotting GO annotations. Nucleic Acids Research, 34, W293–W297.

Zhang B, Zheng J, Peng Y, Liu X, Hoffmann A A, Ma C S. 2015.Stress responses of small heat shock protein genes in Lepidoptera point to limited conservation of function across phylogeny. PLoS ONE, 10, e0132700.

Zhang C, Dai L, Wang L, Qian C, Wei G, Li J, Zhu B, Liu C.2015. Inhibitors of eicosanoid biosynthesis influencing the transcripts level of sHSP21.4 gene induced by pathogen infections, in Antheraea pernyi. PLoS ONE, 10, e0121296.

Zhang C, Zhu B, Dai L S, Liu C, Luo X. 2016.Nucleopolyhedroviruses (NPV) induce the expression of small heat shock protein 25.4 in Antheraea pernyi. Gene,591, 327–332.

Zhang L, Kang H, Jin S, Zeng Q T, Yang Y. 2016. Hsp27 gene in Drosophila ananassae subgroup was split by a recently acquired intron. Journal of Genetics, 95, 257–262.

Zhang Y, Liu Y, Guo X, Li Y, Gao H, Xu B. 2014. sHsp22.6,an intronless small heat shock protein gene, is involved in stress defence and development in Apis cerana cerana.Insect Biochemistry and Molecular Biology, 53, 1–12.

Zhao H, Yi X, Hu Z, Hu M, Chen S, Muhammad R U, Dong X, Gong L. 2013. RNAi-mediated knockdown of catalase causes cell cycle arrest in SL-1 cells and results in low survival rate of Spodoptera litura (Fabricius). PLoS ONE,8, e59527.

Zhu M, Zhang W, Liu F, Chen X, Li H, Xu B. 2016.Characterization of an Apis cerana cerana cytochrome P450 gene (AccCYP336A1) and its roles in oxidative stresses responses. Gene, 584, 120–128.

Zou Z, Lopez D L, Kanost M R, Evans J D, Jiang H. 2006.Comparative analysis of serine protease-related genes in the honey bee genome: Possible involvement in embryonic development and innate immunity. Insect Molecular Biology,15, 603–614.