XlONG Wen-feng, XlE Jia, WEl Lu-ting, ZHANG Si-si, SONG Xiao-wen, GAO Shan-shan, Ll Bin

Jiangsu Key Laboratory for Biodiversity and Biotechnology, College of Life Sciences, Nanjing Normal University, Nanjing 210023,P.R.China

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 inDrosophila 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; Csermelyet al.1998). All of these families are highly evolutionarily conserved, except for the sHSP family, which is more diverse than other four (Li Zet 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 (Kimet al.1998; Waterset 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 (Fuet al.2006). There is also a conserved β-sheet sandwich in a sHSP secondary structure that is responsible for dimer formation (van Montfort Ret al.2001; van Montfort R Let 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(Leeet 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 (Haslbecket 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 (Kegelet al.1996; Waterset 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;Tsvetkovaet al.2002; Haslbecket al.2005; Bakthisaranet al.2015). Furthermore, these sHSPs also participated in many other physiological processes, including cell cycle regulation and cell differentiation (Lianget al.2007;Bakthisaranet al.2015), interfering with apoptosis and defending against diseases, and enhancing the immune response (Arrigo 2000; Altinciceket al.2008).

Recent studies have showed that an increasing number of sHSPs can be found in different insect species and have various functions. Knockdown ofhsp22.6remarkably reduced temperature tolerance inApis cerana cerana, andAcchsp22.6was induced by abiotic and biotic stresses,including insecticides (cyhalothrin and pyridaben) as well as hormone (20-hydroxyecdysone) and bacteria(Ascosphaera apis) treatments (Zhanget al.2014).Similarly,hsp27transcript levels were dramatically upregulated in response to temperature changes and chemical exposure in bothChironomus ripariusandLucilia cuprina(Martinez-Pazet al.2014; Singh and Tiwari 2016). A significant reduction in starvation resistance was also associated with the genotype without a functionalhsp27gene inDrosophila(Haoet al.2007). Typically, overexpressing fly mitochondrialhsp22increased resistance to oxidative stress and extended lifespan by 32%, whereas fliesthat were not expressing this smallhsp22had a 40%decrease in lifespan and showed a reduction of 30% in locomotor activity (Morrowet al.2004). These outcomes suggest there is a key role forshspsin antioxidant activity,locomotor activity and aging in insects. Furthermore,Zhang Cet al.(2015, 2016) indicated thathsp25.4andhsp21.4mRNA levels were accumulated in the fat body to a maximum level at 1.5 and 3 h, respectively, when exposed to nucleopolyhedroviruses inAntheraea pernyi, which suggestsshspsmay be promising candidates for boosting the innate immunity of insects. Additionally,shspsare necessary for development and reproduction inGrapholita molestaandDrosophila, respectively (Michaudet al.1997; Zhang Bet al.2015).Hsp21.3exhibited the highest fold change in diapause individuals, which was highly significant, and along with five othershsps (hsp18.9,hsp19.8b,hsp21.7,hsp24.8andhsp31.8), showed up-regulated expression duringG.molestadiapauses (Zhang Bet al.2015). Whilehsp23andhsp27were significantly expressed in the absence of stress in the testes ofDrosophila, the results implied that these twoshspsplay a vital role during spermatogenesis(Michaudet al.1997). Thus,shspsnot only are essential for different stress responses but have many physiological functions, such as innate immunity, development and reproduction in the insects.

Most recently, tenshspgenes have been identified inTribolium castaneum, the red flour beetles, which is a notorious worldwide pest of stored grain and farinaceous materials for human consumption (Li Zet al.2009; Contreraset al.2013).Tchsp27, which is one of theT.castaneumspecies-specificshsps, was renamed asTchsp18.3based on its predicted molecular weight (Altinciceket al.2008; Li Zet al.2009). The amino acid sequences ofTchsp18.3demonstrated the highest homology tohsp23.0andhsp27.0ofDrosophilaandhsp20.4ofB. mori(Xieet al.2017).Hsp23.0played an important role in hypoxia tolerance and muscle-specific over-expression ofhsp23.0protected muscle from heat shock stress in flies (Azadet al.2009;Kawasakiet al.2016). Meanwhile,Drosophila hsp27.0was involved in the defense response against fungi, bacteria and prolonged exposure to dichlorvos (Pandeyet al.2016;Zhang Let al.2016). InBombyx morinon-diapause eggs,hsp20.4displayed sudden rise in expression level from 42 to 48 h after oviposition (Ponnuvelet al.2010).These results implied thatTchsp18.3, which was likely similar to its homologues in other insect species, possesses the functional diversity. Currently, it’s reported thatTchsp18.3was significantly induced by the septic wounding and ultraviolet A (UV-A) exposure, indicated it was involved in immunity and stress responses ofT.castaneum(Altinciceket al.2008; Sanget al.2012). This phenomenon was also observed after 1 h of exposure to UV-A light in the cell line BCIRL-TcA-CLG1 ofT.castaneum. The outcome suggested thatTchsp18.3plays a vital role in coping with deleterious stimuli to maintain cell normal growth and viabilityinvivoandinvitro(Garcia-Reinaet al.2017). Additionally, our previous study illustrated thatTchsp18.3promoted positive effects in response to a high temperature as well as on resistance to starvation stresses inT.castaneum(Xieet al.2017). Larval RNA interference (RNAi) ofTchsp18.3caused significant lethality at the pupal stage, which seriously hindered the insects eclosion (Xieet al.2017). The lacking ofTchsp18.3at the pupal stage also decreased the eggs production of the red flour beetles (Xieet al.2017). However, it is still unclear howTchsp18.3performs these functions in an insect model.To further explore thesignaling systems ofTchsp18.3involved 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 inT.castaneum.

2. Materials and methods

2.1. Beetle rearing

TheT.castaneumstrains Georgia-1 (GA-1) were reared at 30°C in whole wheat flour enriched with 5% yeast powder under standard conditions (Chenet al.2016; Songet al.2017).

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

For dsRNA synthesis, primers containingTchsp18.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 ofTchsp18.3.Tchsp18.3dsRNAs synthesis and injection were also performed as described in a previous study (Sanget al.2016). In this study, late larvae (20 days old) injected withTchsp18.3dsRNA were denoted as ds-Tchsp18.3groups, 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.3and 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.3samples 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

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 ofT.castaneum(http://www.beetlebase.org/) using SOAPaligner/soap2 and allowing up to two base mismatches (Li Ret al.2009;Liuet 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) (Xionget al.2017).The obtained data were further analyzed for functional annotation of all the detected genes.

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 Zet 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) (Liet 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 (Mortazaviet al.2008). Then,all of the detected beetle genes were assigned to different function categories using the software Blast2GO (http://www.blast2go.org/) (Conesaet 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 adjustedP-value between certain test gene groups and the whole annotation genome (Yeet 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).

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.3and control samples (Wanget 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 (Liuet al.2011). The calculatedP-value was further subjected to a Bonferroni correction using a correctedP-value≤0.05 as a threshold. GO terms fulfilling this condition were defined as significantly enriched GO terms in DEGs (Chenet 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 aQ-value≤0.05 were significantly enriched pathways with DEGs (Liet al.2014).

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

To validate RNA sequencing data, total RNAs were extracted from six larvae ofT.castaneumon the fifth day after injection of buffer and ds-Tchsp18.3as described above. Additionally,1 μg of total RNA was converted to cDNAs by HiScript II Reverse Transcriptase (Vazyme, Nanjing, China) (Liet al.2017). Sixteen DEGs were randomly selected to check the expression patterns between the control and ds-Tchsp18.3groups 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 ΔΔCTmethod (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.

3. Results

3.1. lllumina sequencing and reads assembly

To uncover the differences inT.castaneumgene expression between the control group and ds-Tchsp18.3group, a pooled cDNA sample from two groups was sequenced with the Illumina sequencing platform. We con firmed thatTchsp18.3was knocked down in the ds-Tchsp18.3samples 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.3of 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.3reads 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.3samples, 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 analyzingTchsp18.3regulatory systems.

Table 1 Summary of reads from the RNA sequencing data1)

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

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).

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

To further determine which of the 13 299 genes were significantly and differentially expressed between the ds-Tchsp18.3and control samples, a likelihood ratio test was applied to compare the RPKM-derived read count (Marioniet 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.3sample (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.

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

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 withTchsp18.3functions. The response to stimulus (GO: 0050896) category contains 31 DEGs (Table 2), which includeviperin,dorsal,Hdd11,Peptidoglycan recognition protein 2(PGRP2),defensin1anddefensin2genes 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)andPexwere 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 ofT.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 includedB 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) andEGF-like domain protein 2genes 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 theT.castaneum.The results showed that there are considerable differences in the physiological processes between the ds-Tchsp18.3and control groups.

4. Discussion

Our previous study showed thatTchsp18.3was involved in stress resistance, development and reproduction of the red flour beetles (Xieet al.2017). Therefore, we used RNAsequencing technology to further identify the regulatory systems ofTchsp18.3inT.castaneum. Knockdown ofTchsp18.3in 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 ofTchsp18.3inT.castaneum.

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

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 inAedes aegyptilarvae, which indicates there is an intimate relationship between the immune and stress responses (Muturiet al.2012). This interaction between the innate immune and stress responses was also reported for theDrosophilaMalpighian tubule(Davieset 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(Altinciceket al.2008). More importantly,Tchsp18.3increased in response to septic injury inT.castaneum,which suggests it plays a vital role in innate immunity(Altinciceket 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 wereviperin,dorsal,Hdd11,PGRP2,defensin1anddefensin2, which were also found to be associated with an immune response.Viperinhas 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 inBranchiostoma japonicum(Leiet al.2015). The amphioxus cells transfected withviperincan kill lymphocystis disease virus (LCDV) or inhibit its propagation, and co-incubation ofviperinwith white spot syndrome virus (WSSV) markedly attenuates its infectivity (Leiet al.2015). The results indicated thatviperincan promote resistance against viral infection by participating in immune system cell signaling transductionin vitroandin vivo. Meanwhile, the NF-κB/rel-related genedorsalinDrosophilawas expressed in larval and adult fat bodies, where its mRNA expression was noticeably enhanced after bacterial challenge (Reichhartet al.1993; Lemaitreet al.1995). The genes encoding inducible antibacterial peptides in fruit flies contain NF-κB/rel-related nucleotide sequences and showed that thedorsalprotein can bind to such motifs and transactivate a reporter gene in co-transfection experiments with a fly cell line in a sequence-specific manner (Reichhartet al.1993). Themethuselahgene was also regulated bydorsalto participate in resistance to different kinds of stresses, including high temperatures, starvation and dietary paraquat inDrosophila(Kimet al.2006). These findings implied thatdorsalwas related to the insect immune response as well as stress resistance. Additionally, the expression of theHdd11gene 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 ofSamia cynthia ricini(Baoet al.2003), which indicated that constitutively increased induction ofHdd11could be responsible for increased levels of immunity. As expected, our RNA-sequencing data showed thatviperin,dorsalandHdd11were down-regulated inT.castaneumafter RNAi ofTchsp18.3(Table 2), which suggests that ds-Tchsp18.3resulted in weakened immunity. Furthermore,PGRP2,defensin1anddefensin2were up-regulated after RNAi ofTchsp18.3in the red flour beetles (Table 2), while silencing ofDrosophila PGRP-LBcaused activation of the immune deficiency (IMD) signaling pathway and knocking down flyPGRP-SCled to significantly increased expression of diptericin and attacin (Wang Jet al.2009; Gaoet al.2015), which suggests thatPGRP2probably inhibited the IMD signaling pathway and decreased the synthesis of AMPs as well as reduced the immunity ofT.castaneum.Additionally,defensin1was increased in response to bacterial challenges inBombus terrestrisanddefensin2was also induced following Cry3Ba spore-crystal treatment in thered flour beetles (Erleret al.2011; Contreraset al.2015).The over-expression ofdefensin1anddefensin2arelikely compensated for the loss ofTchsp18.3and evenviperin,dorsalandHdd11expression in insects, which maintained the homeostasis of microorganismsin vivo. Furthermore,PGRP2,defensin1anddefensin2were also up-regulated in the fourth-instar larvae ofA.aegyptiwhen exposed to high temperature stress (Muturiet 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 byTchsp18.3and the decline in the immune response

might reduce the stress resistance of ds-Tchsp18.3insects.However, knockdown ofTchsp18.3can facilitateTchsp20.8,Tchsp22.2,Tchsp68a,Tchsp68bandTchsp70A1gene expression as well as suppressTchsp20.7a,Tchsp20.7bandTchsp21.8expression (Table 2), which shows thatTchspsdo 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)

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)

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

InDrosophila, over-expressinghsp22enhanced resistance to oxidative stress and extended lifespan by 32%, while absence of the expression of this smallhsp22led to a 40%decrease in lifespan and displayed a reduction of 30% in locomotor activity (Morrowet al.2004). This result implied thatTchsp18.3, which may be similar tohsp22, is also involved in antioxidant activity and the aging ofT.castaneum.The RNAi ofTchsp18.3promoted the expression levels ofDuoxandCuZnSODbut repressed theCat,PexandPexlikegene 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 (Bahiaet al.2013; Deng and Zhao 2014). Under normal growth conditions, theDuoxcontributed to the generation of ROS, which maintained the homeostasis of microorganisms in theDrosophilagut (Kim and Lee 2014). However, higher ROS levels caused by over-expressedDuoxalso resulted in disruption of cellular function and structure as well as led to serious oxidative damage of DNA, RNA, proteins and lipids (Zhuet al.2016).The residual/excessive ROS can be eliminated concurrently by antioxidant systems such asCat,PexandPex-like, which aim to protect hosts against oxidative stress (Ahnet al.2012;Zhaoet al.2013). Silencing of theCatgene in the fourthinstar stage ofSpodoptera lituralarvae increased mortality because of induced ROS production, which suggested that ROS scavenging byCatis important forS.liturasurvival(Zhaoet al.2013). Additionally, the over-expression and low-expression ofPexinTetraodon nigroviridisdifferentially regulated the intracellular levels of ROS, and the latter could not reduce H2O2levels (Sunet al.2015). Meanwhile,qRT-PCR analysis showed thatPex-likewas also induced in response to H2O2stimulation inA.cerana cerana, and this result was con firmed by Western blot assay (Zhuet al.2016). The results indicated thatPexandPex-likemay play a crucial role in the antioxidant defense against ROS damage. Generally, the increases ofCuZnSOD in vivocould also promote resistance to oxidative stress and the results reported that mRNA levels of theCuZnSODgene were upregulated instantaneously after exposure to H2O2treatment inBrachionus calyci florus(Van Raamsdonk and Hekimi 2009; Yanget al.2013). These data exhibited that the loss of expression ofCat,PexandPex-likecaused by reduction ofTchsp18.3that was compensated for by up-regulation ofCuZnSODin response to excessive ROS, which further suggests thatTchsp18.3participated in antioxidant stress responses through cross-talk with these antioxidant activityrelated genes.

In addition to the role of lifespan regulation, antioxidant activity has been found to be involved in development and reproduction (Taracenaet al.2015; Denget al.2016).Silencing ofCatgene expression inRhodnius prolixusresulted 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 (Taracenaet al.2015). The results also showed that loss of theCatgene in femaleLutzomyia longipalpisled to a significant increase in mortality and a reduction in developing oocytes (Diaz-Albiteret al.2011). TheSitobion avenaenymphs were continuously exposed to dietary dsRNA ofPexthrough an artificial diet that reduced the survival rate and the ecdysis index revealed thatPexis important to maintain the growth and development ofS.avenae(Denget al.2016). InCaenorhabditis elegans, extracellularPexacted in embryonic development and in curbing the regeneration of adult axons (Gotensteinet al.2010). Moreover, depletion ofPexat theT.castaneumlarval and pupal stages gave rise to 10 and 30% lethality, respectively (Tables 2 and 3)(Donitzet al.2015). Consistent with our previous study,RNAi ofTchsp18.3caused development abnormities before adulthood and reduced egg production inT.castaneum(Xieet al.2017). Therefore, these results suggest that the effects ofds-Tchsp18.3on the development and reproduction of the red flour beetles might be consequences of decreased antioxidant activity.

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 encodingBlimp-1,Gld,Drm,Kinesin-14,Pthr2,Delta(11)-likeandEGF-like domain protein 2were related to the effects ofTchsp18.3on development and reproduction inT.castaneum. specifically,Blimp-1was up-regulated and the other 6 genes were down-regulated by RNAi forTchsp18.3(Table 4; Appendix J). TheBlimp-1expression was significantly induced by 20-hydroxyecdysone (20E)and its product existed during high-ecdysteroid periods and turned over rapidly inDrosophila(Agawaet al.2007).Additionally, its prolonged expression resulted in a reduction offtz transcription factor 1(ftz-f1) expression, which is essential for insect development and reproduction (Agawaet al.2007). A temporally specificftz-f1knockdown using RNAi during the pupal stage exhibited defects in eclosion and in the eye, the antennal segment, the wing and the leg inDrosophila(Sultanet al.2014). Additionally, the loss offtz-f1gene expression led to a complete block in egg production through reducing vitellogenin transcripts and affecting oocyte maturation inT.castaneum(Xuet al.2010). The results indicated that over-expression ofBlimp-1inhibited the expression offtz-f1and produced developmental abnormalities as well as reduced fertility.Apart from its function in immune responses,Gldis also required for insect eclosion as well as normal sperm storage and utilization (Cavener and MacIntyre 1983; Caveneret al.1986; Iida and Cavener 2004). TheDrosophilaGldmutants 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 ofGldin flypupae 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 femaleDrosophila Gldexpression 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 aGldmutant male that had previously mated a few hours before the experimental cross, which showed thatGldfacilitates sperm uptake and release through the spermathecal ducts (Iida and Cavener 2004). Together with the RNAi phenotypes from theTchsp18.3analyses, we showed that knockdown ofTchsp18.3promoted and repressed the levels ofBlimp-1andGldexpression, respectively, and led to developmental and reproductive defects ofT.castaneum. Furthermore, an odd family geneDrmencodes a zinc finger protein, and it is necessary for the development of theDrosophilaembryonic small intestine, an anterior domain of the ectodermal hindgut(Uddinet al.2011). In flyDrmmutants, the elongation of the embryonic hindgut was interrupted and it could not form the small intestine, which resulted in embryogenesis defects(Iwakiet al.2001). In addition, most (91%) of theKinesinfamily member mRNA levels changed abundantly throughout gametophyte development with 52% ofKinesinmRNAs enriched during the mitotic phase and 39% enriched during differentiation inMarsilea vestita. The lack ofKinesin-14inMarsilea vestitaled the nuclei of sperm cells to become elliptic and abnormal, which demonstrated thatKinesin-14is essential for spermatogenesis (Tomei and Wolniak 2016). Indeed,DrmandKinesin-14were both decreased after RNAi ofTchsp18.3, which may be the reason why ds-Tchsp18.3females were completely sterile and the males had low reproduction rates inT.castaneum. More importantly, RNAi against beetleGldandDrmas well asPthr2,Delta(11)-likeandEGF-like domain protein 2caused development arrest, although the extent of the lethality varied significantly according to the affected differential genes (Table 4; Appendix J). The loss ofDrm,Kinesin-14andPthr2at theT.castaneumpupal stage also induced embryonic defects (Table 4; Appendix J) (Donitzet al.2015). In general, these results suggested thatTchsp18.3participated in development and reproduction in cooperation withBlimp-1,Gld,Drm,Kinesin-14,Pthr2,Delta(11)-likeandEGF-like domain protein 2inT.castaneum.

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 theserine protease(SP) family were down-regulated, but 7 members were up-regulated after RNAi ofTchsp18.3.In addition to digestive functions,SPshave significant roles in multiple physiological processes, including stress responses and innate immunity as well as development and reproduction (Gormanet al.2000; Zouet al.2006; Wanget al.2016).SP3,SPH42andSP49transcript levels were increased after the saline orEscherichia coliinjection, and the pathogen challenge also led to a much stronger induction ofSP1,SP2,SP6andSP41gene transcription inA.mellifera(Zouet al.2006). Furthermore, ablation of larvaeSPP167andSPH155resulted in development arrest with 20%and 60% lethality (includes deaths as larvae and pupae),respectively, in the red flour beetles (Table 5; Appendix K) (Donitzet al.2015). While loss ofSPP151,SPP147,SPP38,SP9,SPP142andSPP80at theT.castaneumpupal 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) (Donitzet 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-TcSPP142and ds-TcSPP80, respectively (Table 5; Appendix K) (Donitzet al.2015), which was consistent with silencing ofSPP2cin silkworm embryos that significantly reduced the degradation rate of residual yolk proteins and further suppressed embryogenesis (Wanget al.2016). The results indicated that the SP signaling pathways were activated and amplified byTchsp18.3, which subsequently affected the stress and immune responses as well as development and reproduction ofT.castaneum. Interestingly, the lack of theseSPsin beetles is probably compensated by increasing the expression of the other majorSPs(Table 5; Appendix K), at least partially. The result was similar to outcomes with SP inhibitors in which the expression of someSPs was reduced,whereas other SPs expression levels were significantly increased as part of an integrated compensation response inT.castaneum(Oppertet al.2010; Perkinet al.2017).Overall, these results implied that knockdown ofTchsp18.3changed large-scale SP expression patterns and amplified the SP signaling pathways to further regulate the beetle’s diverse physiological functions.

5. Conclusion

This study identified 569 DEGs from the ds-Tchsp18.3and control samples. These DEGs suggested thatTchsp18.3was not only vital for various stress responses but also involved in innate immunity indirectly throughviperin,dorsal,Hdd11,PGRP2,defensin1anddefensin2as well as other immune-related genes. Additionally, knockdown ofTchsp18.3suppressed the antioxidant activity process,which most likely modulated the effects ofTchsp18.3in development and reproduction ofT.castaneum. The DEGs encodingBlimp-1,Gld,Drm,Kinesin-14,Pthr2,Delta(11)-likeandEGF-like domain protein 2were associated with the development and reproduction of ds-Tchsp18.3.Furthermore, the serine protease signaling pathway was activated byTchsp18.3to 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.

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.

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

Aevermann B D, Waters E R. 2008. A comparative genomic analysis of the small heat shock proteins inCaenorhabditis elegansand 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-1is 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 mediatedPeroxiredoxinVexpression regulates redox homeostasis duringDrosophila melanogastergut 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 insectTribolium castaneum.Developmental and Comparative Immunology,32, 585–595.

Arrigo A P. 2000.sHspas 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 inDrosophila 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 inAnopheles aquasalisresponse toPlasmodium vivaxinfection.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-inducedHdd11gene 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 inDrosophila 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 inDrosophila 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 glutathioneS-transferases relative genes expression patterns inTribolium castaneum.Gene,593, 13–20.

Chen X, Zhang Y. 2015. Identification of multiple small heatshock protein genes inPlutella 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 castaneumimmune defense genes are differentially expressed in response toBacillus thuringiensistoxins 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 ofTribolium castaneumlarvae toBacillus thuringiensistoxin 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 theDrosophilaMalpighian tubule.Journal of Insect Physiology,58, 488–497.

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

Deng F, Zhao Z. 2014. Influence of catalase gene silencing on the survivability ofSitobion 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 femaleLutzomyia longipalpisfecundity 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 beetleTribolium 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 theMusca domesticapeptidoglycan 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 TcATribolium castaneumcell 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. TheC.elegansperoxidasin 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. TheHsp27gene is not required forDrosophiladevelopment 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.Viperininteraction 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 femaleDrosophila 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 theDrosophilaembryonic 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 aDrosophilamodel 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 inDrosophila 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 fromDrosophilamodel of gut-microbiota interactions.Frontiers in Cellular and Infection Microbiology,3, 116.

Lee G J, Pokala N, Vierling E. 1995. Structure andin vitromolecular 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 amphioxusBranchiostoma 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 inDrosophila.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 ofmthl1functions and signal transductions inTribolium castaneum.Gene,547, 310–318.

Li C J, Yun X P, Yu X J, Li B. 2017. Functional analysis of thecircadian clockgene timeless inTribolium 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 mitePanonychus 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-ΔΔCTmethod.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 proteinHsp27gene inChironomus 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 inDrosophila 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 inDrosophila.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 altersAedes aegypticompetence 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 inAedes aegyptilarvae.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 ofTribolium castaneumlarvae 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 ofhsp27rescued neuronal cell death and reduction in life- and health-span inDrosophila melanogasteragainst 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 inTribolium castaneumlarvae.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 silkwormBombyx 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 inCaenorhabditis 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 ofDrosophila.Comptes Rendus Academie Sciences III,316, 1218–1224.

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

Sang M, Li C, Wu W, Li B. 2016. Identification and evolution of two insulin receptor genes involved inTribolium castaneumdevelopment 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 beetleTribolium castaneum.Journal of Insect Physiology,58, 830–836.

Singh M K, Tiwari P K. 2016. Cloning & sequence identification ofHsp27gene and expression analysis of the protein on thermal stress inLucilia 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 inTribolium 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 inDrosophila 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 vertebratePeroxiredoxin 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 fernMarsilea.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 theDrosophilahindgut.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 ofhsp18.3gene 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-basedHSP90functions and signal pathways inTribolium 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 rotiferBrachionus 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 ofsHSP21.4gene induced by pathogen infections, inAntheraea 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 inAntheraea pernyi.Gene,591, 327–332.

Zhang L, Kang H, Jin S, Zeng Q T, Yang Y. 2016.Hsp27gene inDrosophilaananassae 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 inApis 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 ofSpodoptera litura(Fabricius).PLoS ONE,8, e59527.

Zhu M, Zhang W, Liu F, Chen X, Li H, Xu B. 2016.Characterization of anApis cerana ceranacytochrome 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.