Comparative Analysis of Transcriptome Responses to Cold Stress in Galeruca daurica (Coleoptera: Chrysomelidae).

Galeruca daurica (Joannis) has become a new insect pest in the Inner Mongolia grasslands since 2009, and its larvae and eggs have strong cold tolerance. To get a deeper insight into its molecular mechanisms of cold stress responses, we performed de novo transcriptome assembly for G. daurica by RNA-Seq and compared the transcriptomes of its larvae exposed to five different temperature treatments (-10, -5, 0, 5, and 25°C for 1 h and then recovered at 25°C for 1 h), respectively. Compared with the control (25°C), the numbers of differentially expressed genes (DEGs) decreased from 1,821 to 882, with the temperature declining from 5 to -10°C. Moreover, we obtained 323 coregulated DEGs under different low temperatures. Under four low temperatures (-10, -5, 0, and 5°C), a large number of genes were commonly upregulated during recovery from cold stresses, including those related to cuticle protein, followed by cytochrome P450, clock protein, fatty acid synthase, and fatty acyl-CoA reductase; meanwhile, lots of genes encoding cuticle protein, RNA replication protein, RNA-directed DNA polymerase, and glucose dehydrogenase were commonly downregulated. Our findings provide important clues for further investigations of key genes and molecular mechanisms involved in the adaptation of G. daurica to harsh environments.

Galeruca daurica (Joannis) is found mainly in Mongolia, Russia (Siberia), Korea, and China. In china, this beetle is distributed in Inner Mongolia, Xinjiang, and Gansu province (Yang et al. 2010), mainly feeding on onion plants such as Allium mongolicum, A. polyrhizum, and A. ramosum (Hao et al. 2014). Extensive outbreaks of this pest since 2009 have caused great losses to pasture in the Inner Mongolian grasslands and levels of damage continue to increase (Zhou et al. 2016). Galeruca daurica has one generation each year and overwinters as eggs, and its larvae begin to appear as early as in early April when temperature fluctuates acutely and even the lowest temperature decreases below −10°C in its natural habitats (Hao et al. 2015). The overwintering eggs are strongly tolerant to cold with the lower medial lethal temperature (LT50) of −33.08°C (Gao et al. 2015), and the larvae in spring are also highly tolerant to cold with the LT50 of −8.5 to −10.1°C (Li et al. 2014). Rapid cold hardening could significantly enhance the supercooling capacity of the G. daurica larvae (Li et al. 2015). However, to date, a few investigations have been conducted on the underlying molecular mechanisms of cold resistance in G. daurica.

Temperature is a key factor affecting the distribution and performance of insects (Chown and Nicolson 2004). Low temperature in winter can be a threat to the persistence of insect populations. An insect’s cold hardiness is critical in determining its survival and its ability to develop, reproduce, disperse and cause crop losses in the following season (McDonald et al. 2000). Therefore, there has been longstanding interest in the mechanisms by which insects are able to survive in the most severe of winter climates (Bale and Hayward 2010). Most recently, omic technologies have revolutionized how environmental physiologists investigate stress response pathways in insects. Transcriptomics has been used to characterize gene expression changes in response to a great variety of environmental stressors, and this helps not only to identify candidate genes underpinning stress responses but also to predict putative function through the association of coexpressed or differently expressed genes (Hayward 2014). Especially, the fast development of novel high-throughput sequencing technologies, such as RNA-Seq, has provided an opportunity to investigate cold hardiness and signaling-associated genes in different species by de novo assembly and also facilitated quick identification and analysis of the vast majority of transcriptomes (Zhang et al. 2015). Dunning et al. (2013) first used RNA-Seq to identify cold-responsive genes in a New Zealand alpine stick insect, and then further compared the transcriptional responses to low temperature among different populations of alpine and lowland species of New Zealand alpine stick insects, and suggested that cuticle modification may have accompanied colonization of low-temperature alpine environments and the development of a more cold-hardy phenotype (Dunning et al. 2014). Since then, RNA-Seq has been used to investigate the gene expression changes in response to thermal stress in several insect species, such as Drosophila virilis (Parker et al. 2015), Cryptolaemus montrouzieri (Zhang et al. 2015), Microdera punctipennis (Tusong et al. 2016), Nilaparvata lugens, Sogatella furcifera, and Laodelphax striatellus (Huang et al. 2017). These studies show that cold stress can change the expression levels of more than hundreds of genes associated with transcription, metabolism, and cuticular organization, especially enzyme-related genes. Nevertheless, the above authors only used one temperature, 5, 0, or −5°C, as stress or acclimation temperature in their study, but insects may response differentially to different low temperatures. The LT50 of G. daurica larvae is lower down to −10.1°C (Li et al. 2014), so we chose −10°C as the lowest stress temperature in this study.

In this study, we used RNA-Seq to obtain and characterize the global transcriptomes of G. daurica larvae exposed to five different temperatures and further analyzed the differentially expressed genes (DEGs) and identified commonly regulated DEGs during cold stress. We also performed quantitative real-time PCR to explore the expression profiles of several candidate thermal-responsive genes in G. daurica.

Materials and Methods

Experimental Insect and Sample Preparation

G. daurica was originally collected in April 2014 from Xilingol, Inner Mongolia, China (44.215835 N, 114.000067 E, elevation 1,002 m) and reared in the laboratory at 25 ± 1°C with 14:10 (L:D) h cycle. The experiment included five temperature treatments, i.e., −10°C (T1), −5°C (T2), 0°C (T3), 5°C (T4), and 25°C (T5). The treatment at 25°C (T5) was used as a control. For each treatment, 40 second-instar larvae were collected in a 5-ml centrifuge tube and placed in a low-temperature chamber (LRH-100CB, Shanghai Yiheng Scientific Instrument Company, China) under treatment temperature for 1 h, then allowed to recover at 25°C for 1 h, and frozen quickly in liquid nitrogen and stored at −80°C until RNA extraction.

RNA Extraction

Total RNA was extracted from 40 second-instar larvae of G. daurica using the MiniBEST Universal RNA Extraction Kit (Takara, Japan) according to the manufacturer’s instructions. The purity, concentration, and integrity of RNA samples were determined using Nanodrop (Thermo Scientific, Wilmington, DE), Agilent 2100, and Qubit 2.0. The mRNA samples were purified and fragmented using the TruSeq RNA Sample Preparation Kit v2-Set A (Illumina, San Diego, CA) to remove rRNA. After additional quality control, cDNA library construction, Illumina sequencing, and de novo transcriptome assembly were performed in the Beijing Biomarker Biotechnology Company (Beijing, China). Each temperature treatment included one biological replicate (40 second-instar larvae).

Illumina Sequencing and De Novo Assembly

The G. daurica transcriptomes were sequenced on an Illumina HiSeq 2500 platform and performed de novo assembly from the short paired-end reads. Before transcript assembly, all raw reads were processed to remove low-quality and adaptor sequences to obtain the high-quality clean reads. The clean reads were then assembled into contigs, transcripts, and unigenes using the Trinity platform (Grabherr et al. 2011). The raw data have been submitted to the NCBI Sequence Read Archive with accession number SRP150886.

Gene Functional Annotation

BLAST software (Altschul et al. 1997) was used to obtain homologous sequences (E-value ≤ 10−5) by searching several public protein databases, including the National Center for Biotechnology Information (NCBI) nonredundant protein (NR) database (Deng et al. 2006), the UniProt/Swiss-Prot (Apweiler et al. 2004), Gene Ontology (GO; Ashburner et al. 2000), Cluster of Orthologous Groups of proteins (COG; Tatusov et al. 2000), and the Kyoto Encyclopedia of Genes and Genomes (KEGG; Kanehisa et al. 2004).

Differential Gene Expression Analysis

Bowtie was used to align reads reference sequences (Langmead et al. 2009), and the expression levels were estimated by RSEM (Li and Dewey 2011). Gene expression level was measured as numbers of clean reads mapped on the reference sequence and normalized to fragments per kilobase per million (FPKM) fragments. Differential expression analysis was conducted in pair-wise comparisons using the method of Benjamini–Hochberg using the EBSeq software (Leng et al. 2013). The false discovery rate (FDR) was used to determine the threshold of the P-value in multiple tests and analyses. The FDR < 0.01 and the absolute value of log2 (FC) ≥ 1 were used as thresholds to define significantly different gene expressions (Benjamini and Yekutieli 2001). Functional classification of the DEGs was conducted with BLAST analysis as above. The software topGO was used to carry out GO enrichment analysis. For pathway enrichment analysis, all DEGs were aligned to terms in KEGG database and searched for significantly enriched KEGG terms. The gene sequences were also aligned to the COG database to predict and classify functions (Tatusov et al. 2000).

Quantitative RT–PCR Verification

To validate the reliability of RNA-Seq, 11 genes were selected for the qRT-PCR assay. The purified RNA samples were reverse-transcribed using the PrimeScript first-strand cDNA Synthesis Kit (Takara, Japan) following the manufacturer’s protocol. Gene-specific qRT–PCR primers were designed using Primer 3.0 software (http://primer3.ut.ee/). The primer pairs used in qRT-PCR reactions are listed in Supp Table 1 (online only). Succinate dehydrogenase (SDHA) gene was used as a reference gene for expression analyses (Tan et al. 2017a). Negative controls without either template or transcriptase were included in each experiment. qRT–PCR was carried out on the FTC-3000 Real-Time Quantitative Thermal Cycler (Funglyn Biotech Inc., Canada) using Go Taq q PCR Master Mix (Promega) according to the manufacturer’s instructions. Experiments were performed in a 10-µl reaction mixture (1 µl cDNA, 0.2 µl 10 µmol/L forward primer, 0.2 µl 10 µmol/L reverse primer, 5 µl GoTaq qPCR Master Mix, and 3.6 µl RNase free water). All reactions used the following conditions: denaturation at 95°C for 10 min followed by 40 cycles of 95°C for 15 s, 60°C for 1 min and a dissociation at the end of each reaction. Three biological replicates with three technological replicates each were performed for each reaction. The relative abundance of target genes was calculated relative to the reference genes according to the 2−ΔΔCT method (Livak and Schmittgen 2001). The Pearson correlation between the digital gene expression and qPCR results was analyzed using SPSS 13.0 software package (SPSS, Chicago, IL).

Results

Sequencing, RNA-Seq Assembly, and Functional Annotation

In total, 26.86 Gb of clean data passed the Illumina quality filter after transcriptome sequencing of five cDNA samples with Q30 > 85% (Supp Table 2 [online only]). We pooled all high-quality reads (115,022,842) from the five samples to perform the de novo transcriptome assembly (Supp Table 3 [online only]). Using paired-end joining and clustering according to the similarity of contigs, these contigs were further assembled into 142,994 transcripts with a mean length of 1,110 bp and an N50 of 2,109 bp, and 72,352 unigenes with a mean length of 793 bp and an N50 of 1,527 bp. 27,796 unigenes (38.42%) were longer than 500 bp and 15,414 unigenes (21.30%) were longer than 1 kb. This indicates a high-quality assembly, providing a sequence basis for future studies.

To validate and annotate the assembled unigenes, using BLASTx with a cut off E-value of 10−5, they were searched against five public databases, including NR, UniProt/Swiss-Prot, GO, COG, and KEGG. In total, 31,513 unigenes (43.56%) were successfully annotated in at least one of these databases (Supp Table 4 [online only]). Overall, the unigene sequences were most similar to gene sequences from Tribolium castaneum (36.54%), Dendroctonus ponderosae (9.95%), Acyrthosiphon pisum (9.14%), Bombyx mori (5.03%), Camponotus floridanus (3.34%), and Hydra magnipapillata (2.70%) via BLASTx matches (Supp Fig. 1 [online only]).

Differentially Expressed Genes Under Cold Stresses

To investigate the gene expression under cold stresses, we used the FPKM mapped reads method to calculate the expression levels of the unigenes from T1, T2, T3, and T4 samples compared with the control (T5). As shown in Table 1, 882, 1,053, and 1,441 differential expressed genes (DEGs) were identified between the control (25°C) and cold stresses at −10, −5, 0, and 5°C, respectively. The numbers of DEGs declined from 1,821 to 882 when the cold stress temperatures decreased from 5 to −10°C.

Table 1.

Number of differentially expressed genes

DEG set	All DEG	Upregulated	Downregulated
T1_vs_T5	882	539	343
T2_vs_T5	1,053	493	560
T3_vs_T5	1,441	416	1,025
T4_vs_T5	1,821	434	1,387

Differentially Expressed Genes Coregulated in Response to Cold Stresses

Common DEGs were those that were simultaneously upregulated or downregulated under various temperature treatments when compared with the control (25°C). Across four cold treatments (−10, −5, 0, and 5°C), 323 common DEGs (142 upregulated and 181 downregulated) by cold stresses were identified. Of 323 common DEGs, 268 genes (128 upregulated and 140 downregulated) were annotated in the NCBI database (Supp Table 5 [online only]). The top five most differentially upregulated genes were putative glycine-rich protein (c52471.graph_c0) and four hypothetical proteins (c43224.graph_c0, c64046.graph_c0, c43224.graph_c0, and c65962.graph_c0), whereas the top five most differentially downregulated included those encoding two cuticle proteins (c74868.graph_c0 and c68355.graph_c0) and three RNA replication proteins (c38194.graph_c0, c33509.graph_c0, and c75103.graph_c0). Moreover, a large number of genes were upregulated, including those related to cuticle protein, followed by cytochrome P450, clock protein, fatty acid synthase, and fatty acyl-CoA reductase; meanwhile, lots of genes encoding cuticle protein, RNA replication protein, RNA-directed DNA polymerase, and glucose dehydrogenase were highly downregulated.

To further investigate the changes in organism function across cold treatments, we performed functional analysis using the GO and KEGG enrichment. Of 323 common DEGs under cold stress, 95 (29.41%) DEGs (41 upregulated and 54 downregulated) and 36 (11.15%) DEGs (10 upregulated and 26 downregulated) were annotated in the GO and KEGG databases, respectively. According to the GO classification (Table 2), the most frequent assignments of the DEGs coregulated under cold stress were oxidation–reduction process, peptide cross-linking, and nucleoside metabolic process in the biological processes. In the molecular function category, the coregulated DEGs were mostly assigned to oxidoreductase activity, structural constituent of cuticle, binding, monooxygenase activity, and structural constituent of chitin-based cuticle. For cellular components, only extracellular region was significantly enriched.

Table 2.

Significantly enriched GO terms in the DEGs coregulated by cold stresses

GO term	DEGs	Unigenes	Corrected P-value
Biological process
Oxidation–reduction process	20 (13.61%)	278 (3.64%)	7.51E-5
Peptide cross-linking	3 (2.04%)	5 (0.07%)	0.0166
Nucleoside metabolic process	3 (2.04%)	6 (0.08%)	0.0328
Cellular component
Extracellular region	13 (18.84%)	106 (2.51%)	3.81E-7
Molecular function
Oxidoreductase activity	23 (10.18%)	249 (2.46%)	1.06E-6
Structural constituent of cuticle	19 (8.43%)	26 (0.26%)	1.69E-24
Heme binding	9 (3.98%)	48 (0.47%)	1.37E-4
Iron ion binding	8 (3.54%)	50 (0.49%)	0.0019
Chitin binding	5 (2.21%)	23 (0.23%)	0.0190
Monooxygenase activity	5 (2.21%)	25 (0.25%)	0.0290
Structural constituent of chitin-based cuticle	4 (1.77%)	13 (0.13%)	0.0219

DEGs (differentially expressed genes).

KEGG analysis revealed that a large number of coregulated DEGs under cold stresses were significantly enriched in the following pathways: drug metabolism-other enzymes, purine metabolism, starch and sucrose metabolism, metabolism of xenobiotics by cytochrome P450, drug metabolism-cytochrome P450, and fatty acid biosynthesis (Fig. 1).

Fig. 1.

Significantly enriched KEGG pathways in the differentially expressed genes (DEGs) coregulated by cold stresses. The y-axis indicates the pathways, and the x-axis indicates the percentage of DEGs in each pathway. The numbers on the right side of bars indicate the numbers of DEGs significantly enriched in the relevant pathways.

Validation of Gene Expression Profiles by qRT–PCR

To confirm the accuracy and reproducibility of the Illumina RNA-Seq results, eleven DEGs were chosen for qRT-PCR. These DEGs included glucose dehydrogenase (GLD), tyrosine 3-monooxygenase (PLE), lactase-like protein (LCTL), gamma-glutamyltranspeptidase 1 (GGT1), glycine N-methyltransferase (GNMT), aromatic-l-amino-acid decarboxylase (DDC), facilitated trehalose transporter Tret1 (TRET1), cytochrome P450 307a1 (SPO), probable chitinase 3 (CHT3), phosphosulfate synthase (PAPSS), and 40S ribosomal protein S6 (RPS6). The results showed that ten genes displayed the same expression patterns in the qRT-PCR assays as in the RNA-Seq data (R > 0.878, P < 0.05) except PLE (R = 0.3763, P > 0.05; Supp Fig. 2 [online only]). Thus, the RNA-Seq results were considerably reliable for the identification of DEGs under cold stresses in this study.

Discussion

Galeruca daurica is a new insect pest in the Inner Mongolia grasslands, northern China. Our previous studies show this pest has stronger cold tolerance (Li et al. 2014, 2015; Zhou et al. 2016). However, due to the lack of its molecular information, the molecular mechanisms underlying cold hardiness are unknown. Comprehensive investigation of gene expression regulation under cold stress is very important to understand the biochemical and physiological adaptation process in G. daurica. In this study, we de novo assembled five transcriptomes from the G. daurica larvae exposed under five different temperatures, −10, −5, 0, 5, and 25°C, respectively, by RNA-Seq. KEGG analysis revealed that most of the cold-regulated DEGs were enriched in the following pathways: drug metabolism-other enzymes, purine metabolism, starch and sucrose metabolism, metabolism of xenobiotics by cytochrome P450, drug metabolism-cytochrome P450, and fatty acid biosynthesis (Fig. 1). There are some similar results in the investigation on transcriptome responses to cold stress in the ladybird, C. montrouzieri (Zhang et al. 2015) and the desert beetle, M. punctipennis (Tusong et al. 2016). However, the most representative pathway was ribosome, followed by oxidative phosphorylation, spliceosome, and amino sugar and nucleotide sugar metabolism for the seabuckthorn carpenter moth, Eogystia hippophaecolus (Cui et al. 2017). This difference may be due to the different cold exposure durations endured by the tested insects. The E. hippophaecolus larvae were kept at 5 and −5°C for 10 h, whereas the insects were only exposed for 1–2 h in two former studies.

In this study, we identified 11 cuticular protein genes from 268 coregulated DEGs after different cold stresses (Supp Table 5 [online only]). Cold-responsive cuticular protein genes have been identified in many other insects such as flies (Qin et al. 2005), wasps (Colinet et al. 2007), beetles (Carrasco et al. 2011), locusts (Wang et al. 2012), stick insects (Dunning et al. 2013), rice planthoppers (Huang et al. 2017), and seabuckthorn carpenter moth (Cui et al. 2017), suggesting that the change in insect cuticle may play an important role in adaptation to low temperature. Nevertheless, the physiological role of cuticular proteins in insect cold hardiness has not yet been determined (Dunning et al. 2014) and still needs to be further characterized (Huang et al. 2017).

The cytochrome P450 family is a superfamily associated with oxidative metabolism. It has been previously reported that cytochrome P-450 is related to temperature regulation (Nakashima et al. 1996). Huang et al. (2017) found that some cytochrome P450 genes were upregulated under both low and high temperatures, and they hypothesized that cold and heat stress increase oxidative stress in the insect body, and the antioxidant defense system could be induced to resist oxidative stress damage to the organism. In the present study, several genes encoding cytochrome P450 were also found to be significantly upregulated under cold stresses. Zhang et al. (2015) and Cui et al. (2017) also found that many DEGs under cold stresses were significantly enriched in the KEGG pathways related to cytochrome P450. The above results support the hypothesis of Huang et al. (2017).

Trehalose is the major blood sugar in insects playing a crucial role as an instant source of energy and in dealing with abiotic stresses. Trehalose transporter in Anopheles gambiae was reported to discharge the trehalose from the fat body into the hemolymph, which played important roles in coping with stressful conditions (Clark et al. 2009). Cold stress could induce not only an increase of the trehalose phosphate synthase gene in the mRNA levels but also the trehalose concentration in Helicoverpa assulta (Cha and Lee 2016). In this study, two coregulated DEGs encoding trehalose transporter were upregulated and one coregulated DEG was downregulated under different cold stresses (Supp Table 5 [online only]). Our previous study showed that a trehalose-6-phosphate synthase gene GdTPS was induced to upregulate by low temperatures from 0 to −14°C in the G. daurica larvae (Lu et al. 2017). This suggests that the G. daurica larvae may be more adaptive to cold environment, which is consistent with the phenomenon that this pest only occurs in the cold northern regions in China (Hao et al. 2015).

As a highly conserved molecular chaperone, Hsps play important roles in cold hardiness and organism responses to other stressors that can affect the folding and functional conformation of proteins (Clark and Worland 2008, King and MacRae 2015). Many studies have proved that Hsp genes can be induced to express by cold and heat (King and MacRae 2015). Zhang et al. (2015) found that many HSP genes were differentially expressed under cold (4°C for 2 h) in ladybirds. Tusong et al. (2016) also found three Hsp90 and four Hsp70 upregulated by cold stress at 4°C for 3 h in M. punctipennis. Unexpectedly, compared with the control at 25°C, no Hsp gene was regulated during recovery (25°C for 1 h) from cold stress (−10, −5, 0, and 5°C) in this study (Supp Table 5 [online only]). However, Tan et al. (2017b, 2018) found by qPCR analysis that three Hsp genes from G. daurica were significantly upregulated under cold stress for 1 h, such as GdHsp70 at −10°C, GdHsp10 at −5 and 0°C, and GdHsp60 at −10, −5, 5, and 10°C. This difference might be due to the differences in Hsp, insect species, and intensity of temperature treatment. It has been reported that the degree of Hsp induction depends on sex, intensity of treatment, and strain in Bombyx mori (Li et al. 2012).

Circadian clock is related to temperature stress in several insect species (Yoshii et al. 2007, Guerra and Reppert 2013, Chu et al. 2016, Huang et al. 2017). In this study, five clock genes were upregulated and two were downregulated by cold stress (Supp Table 5 [online only]). Clock genes have a large influence on the regulation of metabolic processes in Drosophila (Xu et al. 2011, Sahar and Sassone-Corsi 2012). However, it is not clear if the changes in these genes have a direct influence on the increase in cold tolerance. Therefore, just as suggested by Parker et al. (2015), clock genes are good candidates for orchestrating the changes in metabolic profile during cold acclimation.

In conclusion, we present a global survey for transcriptome profiles in G. daurica under different cold stresses using RNA-Seq. Comparative transcriptome analysis identified many genes and pathways common under different cold stresses, implicating molecular basis of cold stress in this grassland pest. Many cold-regulated genes were first identified in this study, including some functionally unknown genes. These newly found genes may be important to G. daurica in withstanding extremely cold environment in the Inner Mongolia grasslands. Our data will facilitate further molecular investigations on cold tolerance of G. daurica and provide new insights into insect adaptation to harsh environments.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.