De novo transcriptome sequencing of Isaria cateniannulata and comparative analysis of gene expression in response to heat and cold stresses

Isaria cateniannulata is a very important and virulent entomopathogenic fungus that infects many insect pest species. Although I. cateniannulata is commonly exposed to extreme environmental temperature conditions, little is known about its molecular response mechanism to temperature stress. Here, we sequenced and de novo assembled the transcriptome of I. cateniannulata in response to high and low temperature stresses using Illumina RNA-Seq technology. Our assembly encompassed 17,514 unigenes (mean length = 1,197 bp), in which 11,445 unigenes (65.34%) showed significant similarities to known sequences in NCBI non-redundant protein sequences (Nr) database. Using digital gene expression analysis, 4,483 differentially expressed genes (DEGs) were identified after heat treatment, including 2,905 up-regulated genes and 1,578 down-regulated genes. Under cold stress, 1,927 DEGs were identified, including 1,245 up-regulated genes and 682 down-regulated genes. The expression patterns of 18 randomly selected candidate DEGs resulting from quantitative real-time PCR (qRT-PCR) were consistent with their transcriptome analysis results. Although DEGs were involved in many pathways, we focused on the genes that were involved in endocytosis: In heat stress, the pathway of clathrin-dependent endocytosis (CDE) was active; however at low temperature stresses, the pathway of clathrin-independent endocytosis (CIE) was active. Besides, four categories of DEGs acting as temperature sensors were observed, including cell-wall-major-components-metabolism-related (CWMCMR) genes, heat shock protein (Hsp) genes, intracellular-compatible-solutes-metabolism-related (ICSMR) genes and glutathione S-transferase (GST). These results enhance our understanding of the molecular mechanisms of I. cateniannulata in response to temperature stresses and provide a valuable resource for the future investigations.


Introduction
Entomopathogenic hyphomycete fungi, such as Beauveria bassiana and Metarhizium spp., are very common and abundant in ecosystems and many are well known to control insects, nematodes, and plant pathogens [1,2]. However, they are very susceptible to extreme temperatures (heat and cold) [3][4][5]. High temperatures affect the efficiency of entomopathogenic fungi as pest control agent by hindering fungal germination [4], growth [6,7], and virulence [8]. The shelf life of hyphomycetes-based formulations can also be reduced by those factors [9]. Low temperature does not kill conidia of entomopathogenic fungi, but it stops or delays conidial germination [5], and facilitates storage of fungal formulations [3].
Entomopathogenic hyphomycetes have evolved complex molecular mechanisms to protect themselves from extreme environmental temperatures. Many stress-related genes and metabolic pathways have been identified in response to environmental stimuli. Since the fungal cell wall is the first line of defense against harsh environments, cell-wall-major-componentsmetabolism-related (CWMCMR) genes have been shown to play a critical role [10,11]. The gene encoding β-1,3-glucan synthase (MaFKS) in M. acridum plays a role in the maintenance of cell wall integrity, hyperosmotic pressure tolerance, and conidiation [12]. Two GPI (glycosylphosphatidylinositol)-anchored protein Ecm33 orthologous genes, Bbecm33 in B. bassiana and Mrecm33 in M. robertsii, have been reported to be essential for cell wall integrity and multi-stress tolerance [13]. Both heat shock proteins (Hsps) and intracellular compatible solutes (e.g., trehalose, mannitol, glycerol) are important anti-stress agents [14,15]. Under heat stress, in M. anisopliae, 10 Hsp genes and three orthologous trehalose-6-phosphate synthase genes are up-regulated [16]. The expression of Hsp70 gene in B. bassiana is up-regulated by high (38˚C) and low (4˚C) temperatures, or UV stress [17]. Overexpression of the endogenous Hsp25 gene in M. robertsii improves its tolerance to several stresses, including heat and cold [18]. Both Mas5 and Mdj1, the members of Hsp40 in B. bassiana, play a very important role in environmental adaptation and host infection [19,20]. The deletion of mannitol-1-phosphate dehydrogenase (MPD) reduces the content of mannitol in B. bassiana, resulting in a decline of conidial tolerance to diverse stresses, including heat [21]. The ribosome pathway, endocytosis pathway and proteasome pathway are active in M. anisopliae under conditions of heat shock stress [16].
Fortunately, de novo transcriptome assembly and digital gene expression (DGE) sequencing using the Illumina offer the most cost-effective choice for characterizing non-model organisms without a reference genome [31]. In this study, we present the first comprehensive transcriptome characterization of I. cateniannulata and explore the influence of heat and cold temperatures on gene expressions. Using the Illumina HiSeq2000, we generated over 29 million clean reads and these reads were assembled into 17,514 unigenes without a reference genome. DGE analysis was then employed to further validate the genes and metabolic pathways in response to high or low temperatures in I. cateniannulata. This first transcriptomic report of I. cateniannulata functional genes provides an essential database for developing strategies to enhance resilience of entomopathogenic hyphomycetes to temperature variation.

Culture of I. cateniannulata
The wild-type strain ICBS918 of I. cateniannulata was used for this study. This strain was isolated from cadaver of Homona coffearia with a high virulence to H. coffearia and Adoxophyes honmai larvae [24]. This strain was cultured on potato dextrose agar (Difco) for 7 d at 26˚C under the light/dark cycle of 14h/10h. After this period, three replicates of the culture were exposed separately to three temperatures: 4˚C (cold or LT), 26˚C (control or NT), and 39˚C (heat or HT) for 4 h (for a total of nine replicates). Mycelia were then harvested for total RNA extraction.

RNA isolation, library preparation, and sequencing
Total RNA was extracted using the SV total RNA isolation kit (Promega) according to the manufacturer's protocol. The contamination and degradation of the extracted total RNA were preliminarily assessed on 1% agarose gels. The RNA integrity was checked with an Agilent Bioanalyzer 2100 system, and the purity and amount determined using an IMPLEN Nano-Photometer1 spectrophotometer. Total RNA was finally quantified using the Qubit RNA Assay Kit in a Qubit 2.0 Flurometer according to the manufacture protocol (Life Technologies, CA, USA).
Equimolar quantities of the total RNA from the nine samples (1 μg from each sample) were combined into one pool as the input material for the transcriptome library, which acts as a reference library. To build DGE libraries for the nine samples (three replicates in each of the three treatments), 3 μg of extracted RNA were used. The cDNA library was constructed using the NEBNext Ultra RNA Library Prep Kit for Illumina (New England Biolabs, Ipswich, MA, USA) and the final cDNA library was selectively enriched by PCR and purified with the AMPure XP system (Beckman Coulter, Beverly, USA). The cDNA library was then sequenced on the Illumina HiSeq 2000 platform by the Novogene Bioinformatics Technology Co. Ltd (Beijing, China). Finally, 100-bp paired-end raw reads for transcriptome and 100-bp singleend raw reads for DGE were generated. All raw read sequences were deposited in the NCBI Short Read Archive (SRA) database under the accession number of SRP073968.

De novo assembly and annotation
Prior to transcriptome assembly, raw sequencing reads were filtered through in-house Perl scripts to discard dirty reads, which included adaptors, unknown 'N' nucleotides (where the 'N' ratio was more than 10%), and low quality sequences (where the quality score was less than 5). Based on the clean reads, de novo transcriptome assembly into transcripts without a reference genome was carried out using short reads assembling program Trinity with min_kmer_cov set to 2 and all other parameters set by default [32]. For removing redundancy, if a component contained more than one assembled transcript, only the longest one was regarded as a unigene. To annotate the unigenes, we used NCBI BLAST 2.2.28+ with an E-value threshold of 10 −3 in the database of eukaryotic orthologous groups (KOG), and with an E-value threshold of 10 −5 in the NCBI non-redundant protein sequences (Nr) database, NCBI nucleotide collection (Nt) database and the Swiss-Prot protein database. Protein family (Pfam) was assigned using the HMMER 3.0 package. The database categories of Kyoto Encyclopedia of Genes and Genomes (KEGG) were assigned to the unigenes according to the KEGG Automatic Annotation Server (KAAS) online [33]. Based on BLASTX hits against the Nr and Swiss-Prot protein database (E value 10 −5 ), the gene ontology (GO) annotation of unigenes was obtained using Blast2GO v2.5 [34] and GO functional classification was performed using WEGO software [35]. The aligning results were compared against the Nr and Swiss-Prot databases with a priority order of Nr > Swiss-Prot to decide the unigenes' direction and coding sequences (CDSs). When unigenes were unaligned to any of the above databases, ESTScan software [36] was used to predict the coding regions and sequence orientation.

Analysis of DGE tags
Raw reads of DGE sequences were deposited in the NCBI Short Read Archive (SRA) database (http://www.ncbi.nlm.nih.gov/sra). Raw sequence data of the libraries for DGE profiling analyses were filtered to remove reads containing adaptors and reads with more than 10% unknown nucleotides, and reads with more than 50% of low-quality base (value 5). Clean reads were mapped into the assembled transcriptome reference database and the read counts for each gene derived from the mapping results were achieved by RSEM [37]. All read counts were normalized using expected number of fragments per kilobase of transcript sequence per millions base pairs (FPKM), as described by Trapnell et al. [38]. The DESeq R package (1.10.1) was used for differential expression analysis of unigenes between the control (26˚C) and the treatments (4˚C or 39˚C). To evaluate the reproducibility of DGE library sequencing, a Pearson correlation analysis was performed using the replicates between the various treatments. A corrected P-value of <0.05 (P values adjusted using the Benjamini & Hochberg method, P adj ) was used as the threshold to judge the significant differences in gene expression. The heatmap of the differentially expressed genes (DEGs) was constructed using the R packages of ggplot2 and pheatmap.
In order to analyze the biological functions and metabolic pathways of the DEGs, DEGs were subject to GO functional analysis using Blast2GO v2.5 [34] and KEGG enrichment analysis using the KOBAS software [39]. Pathways with P value cutoff of 0.05 were considered significantly enriched by DEGs.

Quantitative real-time PCR analysis
To confirm the DGE results, 18 genes were randomly selected for qRT-PCR verification. Specific primers for qRT-PCR were designed using the Primer Premier 5.0 software (Premier, Canada) and listed in the S1 Table. RNA for qRT-PCR analysis was extracted from the mycelia under the extreme temperature, as described above, using the SV total RNA isolation kit (Promega, USA). Reverse transcription was carried out according to the manufacturer's protocol of GoScript™ Reverse Transcription System (Promega, USA). The qRT-PCR was run in the CFX96 Touch TM Real-Time PCR Detection Systems (Bio-Rad, USA) using BRYT-Greenbased PCR reaction. PCR amplification was performed in a total reaction mixture of 20 μL, containing 20 ng cDNA, 10 μL 2 × GoTaq1 qPCR Master Mix (Promega, USA) and 0.2 μM of each primer. PCR was run with the standard thermal cycle conditions using the two-step qRT-PCR method: an initial denaturation at 95˚C for 30s, followed by 40 cycles of 3 s at 95˚C and 30 s at 60˚C. Each reaction was run in triplicate, and the average threshold cycle (C t ) was calculated for each replicate. The assembled glyceraldehyde-3-phosphate dehydrogenase unigene (GenBank accession no. KT944290) was employed as the internal gene and the relative expression of gene was calculated using the 2 -ΔΔCt method. The correlation of the fold change of the gene expression ratios between qRT-PCR and RNA-seq was checked by linear regression analysis in SPSS 17.0 software.

Assembled transcriptome
As shown in Table 1, a total of 2.92 Gb nucleotides were obtained with a Q20 percentage of 95.5%. The percentage of unassigned base "N" was 5.95% and the average GC content, 55.17% (Table 1). By overlapping information in high-quality reads, 24,707 transcripts were generated, with an average length of 1,581 bp and an N50 of 2,664 bp. From these transcripts, 17,514 unigenes were obtained with a mean size of 1,197 bp ( Table 2).

DGE library sequencing and mapping to the reference transcriptome
Nucleotides of 0.74-1.41 Gb were obtained in each of the nine DGE libraries with a Q20 percentage of over 95%. The percentage of unassigned base "N" was below 0.01% and the average GC content was 55.01%-56.30% (Table 1). The clean reads of the nine DGE libraries were mapped to the above-constructed transcriptome reference database for each sample and the total mapped reads were above 96% for each replicate (S5 Table), which suggested the transcriptome was a reliable reference. The square of the Pearson correlation coefficient (R 2 ) varied during 0.765-0.918 in HT treatment, 0.827-0.923 in LT treatment, and 0.882-0.933 in NT, indicating good operational stability and reliability of DEG library sequencing (Fig 4, S2 Fig).

Differentially expressed genes and qRT-PCR verification
There were a total of 5,686 DEGs in response to extreme temperature stresses, and the expression profile in the heatmap (Fig 5) showed significant differences among the three temperature treatments. In HT treatment, the expression of 4,483 genes was significantly changed when compared with the control (NT), with 2,905 genes being up-regulated and 1,578 genes downregulated. In LT treatment, the expression of 1,245 genes was up-regulated and 682 genes down-regulated when compared with the control (Fig 6, S6 Table). Overlapped transcripts between heat and cold stresses showed that there were 724 DEGs responding to both heat and cold stresses, 3,759 DEGs only responding to heat stress, and 1,203 DEGs only responding to LT treatment (Fig 7).
To verify the authenticity and reproducibility of the RNA-Seq results, 18 unigenes were randomly selected for qRT-PCR analysis and the expression profiles of the candidate unigenes by qRT-PCR were similar to the results of transcriptome analysis (Fig 8). The linear regression analysis showed significantly positive correlation of the relationship between gene expression ratios of qRT-PCR and RNA-seq was significantly positive (r = 0.75, P = 0.003) (Fig 9), confirming our transcriptomic data validity.

DEGs in response to heat treatment
A total of 4,483 DEGs including 2,905 up-regulated and 1,578 genes down-regulated were identified in HT. GO functional classification of 4,483 DEGs showed that they could be categorized into 50 functional groups, belonging to three main GO domains: biological processes (20), cellular components (16), and molecular functions (14). Among these groups, we found that cellular process, metabolic process, and single-organism process in the biological process KEGG enrichment analysis showed that 8 pathways were significantly enriched (S7 Table). Among them, the ribosome biogenesis in eukaryotes was the most significantly enriched pathway, which contained 29 DEGs (28 down-regulated DEGs and 1 up-regulated DEG); the aminoacyl-tRNA biosynthesis was the second most significantly enriched pathway, which contained 18 DEGs (10 down-regulated DEGs and 8 up-regulated DEGs). However, we focused on the genes that were involved in endocytosis pathway, because of its ability to destroy misfolded and damaged proteins. After HT treatment, 12 DEGs related to clathrin-dependent endocytosis (CDE) were observed (Fig 10A). Among these 12 genes, 6 genes were up-regulated including PLD (comp7126_c0), dynamin (comp7305_c0), Hsc70 (comp7450_c0), VPS22 (comp2053_c0), CHMP4 (comp8162_c0) and CHMP3 (comp7520_c0); 6 genes were downregulated including PLD (comp6558_c0), AP-2 (comp8686_c0 and comp6289_c0), Hsc70 (comp5575_c1), ArfGEF (comp4739_c0) and VPS45 (comp4862_c0).
In this study, four categories of DEGs functioning as temperature sensors were observed, including CWMCMR genes, Hsps, ICSMR genes and glutathione S-transferase (GST) genes. DEG analysis showed that there were 14 CWMCMR unigenes in the up-regulated unigenes, but no CWMCMR unigene was observed in the down-regulated unigenes (Table 4). In the Hsp superfamily, 25 unigenes from 6 Hsp families (Hsp100, Hsp90, Hsp70, Hsp60, Hsp40 and sHsp) were differentially regulated under heat treatments, including 24 up-regulated unigenes and 1 down-regulated unigene. Hsp40 family was the largest group with 8 family members and Hsp70 family was the second most abundant group with 6 family members (Table 5).  Eight ICSMR genes including trehalose metabolism-related genes, mannitol metabolismrelated genes, glycerol metabolism-related genes and arabinitol metabolism-related genes were differentially regulated (Table 6). Six GST genes were found to be up-regulated, and no GST gene was observed to be down-regulated (Table 7).

DEGs in response to cold treatment
In comparison with NT, 1,927 DEGs were found in LT including 1,245 up-regulated genes and 682 down-regulated genes. All DEGs were categorized into 49 functional groups by GO functional classification, including biological processes (20), cellular components (14), and molecular functions (15). Among these groups, we found cellular process, metabolic process, and single-organism process in the biological process ontology, cell, cell part, and membrane in the cellular component ontology, and binding and catalytic activity in the molecular function ontology were the dominant clusters (S4 Fig).

DEGs in response to both heat and cold treatments
A total of 724 co-regulated DEGs (Fig 7) including 496 up-regulated and 228 down-regulated genes were identified in both heat and cold treatment. All DEGs were categorized into 45 functional groups by GO functional classification, including biological processes (21), cellular components (15), and molecular functions (9). Among these groups, we found cellular process, metabolic process, and single-organism process in the biological process ontology, cell, cell part, and organelle in the cellular component ontology, and binding and catalytic activity in the molecular function ontology were the dominant clusters (S5 Fig). The co-regulated DEGs were highly enriched in glycolysis / gluconeogenesis, regulation of mitophagy-yeast, sulfur metabolism, glycerolipid metabolism and ABC transporters by KEGG analysis (S7 Table). Besides, there were two co-regulated DEGs, PLD (comp7126_c0) and CHMP4 (comp8162_c0), related to endocytosis pathway. DEG analysis revealed that 5 CWMCMR genes (Table 4), 4 Hsp genes (Table 5), 3 ICSMR genes (Table 6) and 1GST gene ( Table 7) were both up-regulated under heat and cold treatments. Table 4. Variation in the expression of CWMCMR genes when exposed to heat and cold treatments in I. cateniannulata. comp2210_c0 GPI anchored cell wall protein ☆ " GPI-anchored protein """ indicates up-regulated-genes and "☆" indicates unigenes no differentially expressed.

Unigene ID Annotation HT/NT LT/NT Category
https://doi.org/10.1371/journal.pone.0186040.t005 Table 6. Variation in the expression of ICSMR genes when exposed to heat and cold treatments in I. cateniannulata.

Discussion
As the genomic information of I. cateniannulata was unavailable, and the molecular response mechanisms of the species to environmental variation not well understood, we used next-generation sequencing technology and transcriptome analysis as an alternative for in-depth analysis of molecular responses of I. cateniannulata to high and cold temperatures. In the present study, the transcriptome characterization of I. cateniannulata generated 29.2 million clean reads and assembled into 17,514 unigenes with a mean size of 1,197 bp. While 12,095 (69%) unigenes were successfully annotated using the public databases, 5,419 (31%) unigenes had no homologous sequences to public databases. This indicates that I. cateniannulata may contain several species-specific unigenes. The variation in number of unigenes and the expression profiles of DEGs under heat and cold temperature stresses suggests that the molecular responses of I. cateniannulata to high and low temperatures may greatly vary as reported in Pyropia yezoensis [40]. Though the cold temperature slows down the metabolic rates and results in the dormancy in organism, we still identified many DEGs in metabolic pathways under cold temperature treatment in our study. It maybe the cold temperature treatment (4˚C for 4h) was not devastating to I. cateniannulata, and there were many DEGs have been activated. Our result is well in accordance with the study of Pyropia yezoensis in response to chilling and freezing stress [40], Chrysanthemum nankingense under low temperature [41], and Anthurium andraeanum under cold stress [42]. Endocytosis involves cells take up molecules (such as proteins) via vesicles, which is essential for cell-to-cell communication and cellular responses to external stimuli in all eukaryotic cells [43,44]. Endocytosis pathway is significantly enriched in lotus (Nelumbo Adans) [45] and filamentous fungus (M. anisopliae) [16] under conditions of heat shock stress. And this pathway is also induced by cold stress in amur carp (Cyprinus carpio haematopterus) [46]. Endocytosis pathway contains two major categories: CDE and CIE [47]. CDE is well characterized and highly conserved cellular process from humans to fungi [48]. After heat treatment, the pathway of CDE was activated in I. cateniannulata as reported in M. anisopliae [16]. Dynamin is one of critical factors for different stages of CDE [49]. It was significantly up-regulated in I. cateniannulata under heat stress. Here, one Hsc70 gene (comp7450_c0) up-regulation and the down-regulation of comp5575_c1 were observed after heat stress. Hsc70 play a very important role in the ATP-dependent dissociation of clathrin from clathrin coated vesicles (CCVs) and other key processes in CDE [50]. Late endosomes (also known as multivesicular bodies, MVBs) is one of the major trafficking hub of the endocytic pathway, which formation is essential for cells to destroy misfolded and damaged proteins [47]. The endosomal sorting complex required for transport (ESCRT) plays a very important role in the MVB sorting pathway [47]. In this study, three genes involved in the ESCRT pathway, including ESCRT-Ⅱ (VPS22) and ESCRT-Ⅲ (CHMP3 and CHMP4) were up-regulated under heat stress. It was speculated that Table 7. Variation in the expression of GSTs when exposed to heat and cold treatments in I. cateniannulata.

Unigene ID Annotation HT/NT LT/NT Category
glutathione S-transferase domain-containing protein " " """ indicates up-regulated-genes and "☆"indicates unigenes no differentially expressed. https://doi.org/10.1371/journal.pone.0186040.t007 the activation of CDE would facilitate the elimination of misfolded and damaged proteins, which caused by heat shock in I. cateniannulata. Apart from CDE, CIE pathway have been found in eukaryotic cells, which mediates many biological processes [51]. In this study, CIE was activated when I. cateniannulata was exposed to cold treatment. The gene encoding AMSH (associated molecule with a Src homology 3 domain), a deubiquitinating enzyme, was up-regulated as well as, three genes involved in the ESCRT pathway, ESCRT-0 (STAM (signal transducing adaptor molecule)) and ESCRT-Ⅲ (CHMP (charged multivesicular body proteins) 2 and 4). These results suggest that the activation of CIE pathway would accelerate the removal of denatured proteins caused by low temperature stress in I. cateniannulata. A total of 20 CWMCMR genes were up-regulated in response to high and/or low temperatures. Cell wall plays a very important role in protecting fungi against a variety of harsh environments such as heat, cold, desiccation and osmotic stress [52]. In this study, 1 β-1,3-glucan synthase orthologous gene (comp34522_c0) was up-regulated under cold temperature. We suggest that the cell may synthesize more β-1,3-glucan to enhance cell wall tensile strength against cold temperature. β-1,3-glucan has a coiled spring-like structure that confers a degree of elasticity and tensile strength to the cell wall [11]. β-1,3-glucan synthase (MaFKS)-RNAi transformants of entomopathogenic fungus M. acridum are more sensitive to agents that disturb the cell wall or cell membrane and to hyperosmotic stress in comparison with the wild type [12].
GPI-anchored protein orthologous gene (comp2210_c0) was up-regulated in response to cold temperature. GPI-anchored protein may not only play a role in maintaining the fungal cell wall integrity [13,53], but also contribute to their multi-stress tolerance [13]. So the upregulation of orthologous genes (comp2210_c0) of GPI-anchored protein may enhance the cell wall integrity and offer more resistance to cold temperature. Other up-regulated genes such as glucan-metabolism-related genes, chitin/chitosan-metabolism-related genes, and mannoprotein-metabolism-related genes, might have also been involved in maintaining cell wall integrity and increasing tensile strength against heat and cold. However, the exact functions of these genes are elusive.
Twenty-seven Hsp genes from 6 major families (Hsp100, Hsp90, Hsp70, Hsp60, Hsp40 and sHsp) were differentially regulated when exposed to heat and/or cold in I. cateniannulata. Hsps are ubiquitous in all prokaryotes and eukaryotes and can be induced by several kinds of stresses, including extremes temperatures, desiccation, and toxic substances [54][55][56]. In the filamentous fungus M. anisopliae, 10 Hsp genes are up-regulated in heat-treated conidia [16]. Here, one hsp98 homolog (comp7457_c0) was up-regulated in response to both heat and cold and one hsp78 homolog (comp7486_c0) was up-regulated in response to cold treatment. Hsp100 family is a major heat-regulated protein family in several species [57]. The Neurospora crassa Hsp98 protein, a member of hsp100 family, is highly expressed in response to heat shock [58]. Hsp78, a member of the Clp/Hsp100 family in S. cerevisiae, is required for the maintenance of mitochondrial function under heat stress [59]. Hsp90 homolog (comp5388_c0) was up-regulated upon exposure to heat. All eukaryotic cells produce prominent heat-shock protein with the molecular weight ranging from 80 to 90 kD, which are classified as the Hsp90 family [14,60]. Hsp90 transcription is significantly induced when exposed to various abiotic stresses such as heat, cold, and oxygen deprivation [61,62].
Interestingly, 5 members of Hsp70 family in I. cateniannulata were up-regulated in response to heat and might be invovled with the removal of denatured proteins. This is in agreement with other studies. For example, Hsp70 can prevent the aggregation of unfolding proteins and even refold aggregated proteins under heat [63]. The fungus B. emersonii has 10 putative Hsp70 homologs, and all Hsp70 genes (except for Hsp70-4 and Hsp70-6) are induced to different degrees upon exposure to heat [64]. The yeast S. cerevisiae has 14 Hsp70 genes, in which cytosolic Hsp70-Ssa (Ssa1, Ssa2, Ssa3 and Ssa4) is heat-inducible [65]. The transcriptional level of Hsp70 is up-regulated in B. bassiana when exposed to heat (38˚C), cold (4˚C), or UV stress [17].
Three Hsp60 homologs were up-regulated under heat in I. cateniannulata. The Hsp60 gene expression levels are also up-regulated in the Aspergillus fumigatus, A. terreus and Scedosporium apiospermum under heat [66]. In addition, 10 Hsp40/DnaJ homologs were up-regulated under heat and/or cold temperatures. Hsp40 family, also known as J-protein family, is the largest class of Hsp70 cofactors, which bind nonnative proteins and deliver them to Hsp70 [67,68]. Both Mas5 and Mdj1, the members of Hsp40 in B. bassiana, have been shown to be indispensable for the environmental adaptation and virulence [19,20].
In our study, 2 Hsp30 homologs (comp2221_c0 and comp6572_c0) were up-regulated under both heat and cold temperatures, and 1 Hsp20 (comp1380_c0) and 1 Hsp10 homologs (comp1717_c1) up-regulated when exposed to heat. The expression of Hsp23, a small heatshock protein gene in Trichoderma virens is increased when the fungus is grown at extreme temperatures (4, 10 or 41˚C) [69]. A small heat shock protein gene hsp25 is up-regulated in M. robertsii in response to extreme temperatures (4, 35, and 42˚C), and overexpression of hsp25 improves the growth of M. robertsii when exposed to heat [18]. The sHsps of A. nidulans have been shown to take part in resisting adverse conditions, including heat and cold as well as oxidative/osmotic stresses [70].
Three trehalose-metabolism-related genes were up-regulated after heat and/or cold treatments in I. cateniannulata. Trehalose is present in a wide range of organisms, and could serve as a stabilizer and protectant of proteins and cellular membranes against a variety of stresses such as heat, cold, oxidation, and desiccation [71]. Under heat or chemical stress, the increasing of trehalose in the cell, which associated with up-regulation of the trehalose-6-P phosphatase transcript in arbuscular mycorrhizal (AM) fungi Glomus intraradices was observed [72]. Several trehalose accumulation-related genes are up-regulated in M. anisopliae in response to heat [16].
Here, two mannitol 1-phosphate dehydrogenase (MPD) orthologous genes (comp7547_c0 and comp13549_c0) showed up-regulation in response to heat. It was speculated that the up-regulation of MPD would increase the content of mannitol in I. cateniannulata. One D-arabinitol 2-dehydrogenase orthologous gene (comp5307_c0) was up-regulated under heat treatment. Darabitol is one of the polyols found most frequently in fungi [73], which may act as adversity protectant [74]. In our study, 2 glycerol-metabolism-related genes were up-regulated when exposed to cold treatment, and 2 glycerol-metabolism-related genes were up-regulated in response to both heat and cold treatment in I. cateniannulata. The up-regulation of glycerol-metabolism-related genes may contribute to the increasing glycerol content against heat and cold temperatures.
GSTs play a very important role in response to oxidative stress by removing reactive oxygen species and regenerate S-thiolated proteins [75]. The GSTs involved in protecting Schizosaccharomyces pombe cells from damage causing by oxidative stress [76]. In our study, 6 GST genes were up-regulated under heat treatment, and 1 GST gene was up-regulated under cold treatment. We suggest that GST genes up-regulation might protect I. cateniannulata cells against damage resulting from oxidative stress induced by heat and cold treatments.
When I. cateniannulata was exposed to cold stress, the homologue of CSP (comp1755_c0) and GRP (comp7081_c0) were up-regulated. This is well in accordance with M. anisopliae's CSP (CRP1) and GRP (CRP2) homologue, which play a key role in against cold stress [77].

Conclusions
The combination of RNA-seq and DGE analysis based on next generation sequencing technology provided comprehensive information on gene expression of I. cateniannulata, an entomopathogenic fungus for which little genomic information was available. Many DEGs of I. cateniannulata were identified under heat and cold temperatures with significant differences in molecular responses. In this study, we mainly focused on endocytosis pathway and identified several genes that were either up or down regulated when exposed to changing temperatures. Candidate stress-related genes may be useful tools for improvement of strain tolerance against extreme environmental temperatures in I. cateniannulata.