Advertisement
Browse Subject Areas
?

Click through the PLOS taxonomy to find articles in your field.

For more information about PLOS Subject Areas, click here.

  • Loading metrics

Open Access

Peer-reviewed

Research Article

Transcriptomic analysis to uncover genes affecting cold resistance in the Chinese honey bee (Apis cerana cerana)

Correction

31 Aug 2017: Xu K, Niu Q, Zhao H, Du Y, Jiang Y (2017) Correction: Transcriptomic analysis to uncover genes affecting cold resistance in the Chinese honey bee (Apis cerana cerana). PLOS ONE 12(8): e0184468. https://doi.org/10.1371/journal.pone.0184468 View correction

Abstract

The biological activity and geographical distribution of honey bees is strongly temperature-dependent, due to their ectothermic physiology. In China, the endemic Apis cerana cerana exhibits stronger cold hardiness than Western honey bees, making the former species important pollinators of winter-flowering plants. Although studies have examined behavioral and physiological mechanisms underlying cold resistance in bees, data are scarce regarding the exact molecular mechanisms. Here, we investigated gene expression in A. c. cerana under two temperature treatments, using transcriptomic analysis to identify differentially expressed genes (DEGs) and relevant biological processes, respectively. Across the temperature treatments, 501 DEGs were identified. A gene ontology analysis showed that DEGs were enriched in pathways related to sugar and amino acid biosynthesis and metabolism, as well as calcium ion channel activity. Additionally, heat shock proteins, zinc finger proteins, and serine/threonine-protein kinases were differentially expressed between the two treatments. The results of this study provide a general digital expression profile of thermoregulation genes responding to cold hardiness in A. c. cerana. Our data should prove valuable for future research on cold tolerance mechanisms in insects, and may be beneficial in breeding efforts to improve bee hardiness.

Citation: Xu K, Niu Q, Zhao H, Du Y, Jiang Y (2017) Transcriptomic analysis to uncover genes affecting cold resistance in the Chinese honey bee (Apis cerana cerana). PLoS ONE 12(6): e0179922. https://doi.org/10.1371/journal.pone.0179922

Editor: Olav Rueppell, University of North Carolina at Greensboro, UNITED STATES

Received: March 1, 2017; Accepted: June 6, 2017; Published: June 26, 2017

Copyright: © 2017 Xu et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Data Availability: All relevant data are within the paper and its Supporting Information files.

Funding: This work is supported by the National Nature Science Foundation of China. The award number of this foundation is 31372382. The recipient of this foundation is Qingsheng Niu who is the second author of this paper. The funders were undertaken to explore cold tolerance mechanism of honey bee and had an important role in study design, decision to publish, and preparation of the manuscript.

Competing interests: The authers have declared that no competing interests exist.

Introduction

Temperature is one of the most important abiotic factors influencing insect activity and geographic distribution [1]. Insects may face daily or seasonal temperature fluctuations in temperate climates and subzero temperatures in frigid zones, severely compromising metabolic and physiological activity [2]. Although differences in population distribution patterns [3–5], developmental stages [6] and diapause status [7] all contribute to the considerable variation in cold tolerance among insects, they have evolved a series of consistent physiological and behavioral strategies to mitigate cold-related damage. Behaviorally, some insects (e.g., monarch butterflies Danaus plexippus) perform long-distance migration to warm habitats [8], while others take shelter in thermally buffered microclimates, such as tree bark crevices or under snow and frozen soils [1, 9]. Physiologically, insects are broadly divided into freeze-tolerant and freeze-susceptible species, based on their adaptation mechanisms; the former group can survive after a percentage of their extracellular fluids freeze, whereas the latter group cannot [10]. Winter survival rates and cold resistance is improved in freeze-susceptible insects through the regulation of physiological and metabolic pathways that lower their supercooling point (SCP), the temperature below which insects will freeze [11, 12]. The SCP is thus an indicator of cold resistance capacity. Insect SCP can be altered through extreme dehydration, purging the gut of ice nucleators, and producing multimolar concentrations of cryoprotectants [13, 14].

Cryoprotectants are divided into high molecular weight antifreeze-associated proteins and small-molecule cryoprotectants. The former group includes heat shock proteins (HSPs), antifreeze proteins, antifreeze glycolipids, and other cold-induced stress proteins, which all exhibit distinctive structures and properties that are critical to overcoming low temperatures [15–17]. In insects, the synthesis of small cryoprotectants is substantially common; these compounds include polyhydric alcohols, sugars, and glycerol [18]. Studies on cold-shock and chilling tolerance in Drosophila have shown that energetic reserves (e.g., glycogen, triacylglycerols, and proline) may be important in cold-coping mechanisms [19]. In the cold-hardy gall fly (Eurosta solidaginis), for example, the larvae exhibit high levels of glycerin, sorbitol, glucose, and trehalose biosynthesis [20]. Other overwintering insects also synthesize glycoproteins and various amino acids; these molecules likely play a key role in lowering the SCP as part of a cold-stress response [21, 22].

The vast majority of research on insect cold resistance has focused on solitary species, with social insects being relatively ignored. The honey bee (Apis spp.) is an eusocial insect and major pollinator, with high economic value worldwide; an understanding of cold-tolerance mechanisms in the honey bee thus has important practical implications. Honey bees have been anthropogenically introduced to a wide range of tropical and temperate geographical regions. Physiological and behavioral thermoregulation is essential for the ectothermic honey bees at both individual and colony levels [23, 24]. Bee colonies tend to have high stable temperatures in the central core area (where the brood is located) under a variety of environmental conditions [25], achieved via the clustering and layering of worker bees to form an insulating shell [23, 26]. Because of this ability to maintain a constant temperature, studies on thermoregulation in honey bees tend to focus on the colony. However, individual forager bees, for example, must still survive colder environmental conditions in early spring and late autumn to successfully gather food [27]. Currently, we know that honey bees maintain their body temperature above ambient conditions through burning fuel during flight [27]. However, when temperatures drop below 12°C, individual honeybees are unable to fly, and they become comatose at temperatures below 10°C [28]. More work is required to understand the molecular regulatory mechanisms of cold resistance in honey bees.

The Chinese honey bee (Apis cerana cerana) is an important endemic species that is critical to environmental improvement and ecosystem balance in China [29]. A. c. cerana has a lower SCP [30] and is more cold-hardy than the Western honey bee (A. mellifera), displaying high forager activity below 14°C and pollinating winter-flowering plants in southern China [31]. Thus, A. c. cerana should be an effective model organism for studying transcriptional regulation of cold resistance in honey bees. In this study, we characterized gene expression in cold-treated A. c cerana adults, using Illumina sequencing technology to examine their transcriptomes. Through identifying key temperature-regulated genes and gene networks, this research should offer important insight into the mechanisms of cold tolerance.

Materials and methods

Insect treatments

Chinese honey bees (A. c. cerana) were maintained in an experimental apiary from the College of Animal Science and Veterinary Medicine of Shanxi Agriculture University (Taigu, China). Winter bees were bred on September 15, 2015 and emerged on October 5, 2015. Newly emerged adults were marked with paint. Adults 20-day post-emergence were collected from three hives on October 25, 2015, when the environmental temperature in this region ranged from 6°C to 15°C. Three hundred adults (n = 100/hive) were collected and each individual bee was placed into 15 mL centrifuge tubes with small pores. Bees were then divided into two groups (n = 150/group) and treated for 2 h at either 0°C (cold treatment; ZOT) or 25°C (normal-temperature treatment; ZRT) in constant-temperature incubators. After 2 h, 60 randomly selected bees (n = 30/group) were allowed to recover for 30 min before mortality was calculated as an indication of phenotypic (cold-tolerance) variation. The remaining individuals (120/group) were immediately frozen in liquid nitrogen for total RNA extraction and RNA-seq analysis.

Library construction

Total RNA was extracted using TRIzol reagent (Ambion, Foster City, CA, USA) following the manufacturer’s protocol. Total RNA concentration was determined using a Qubit® RNA Assay Kit and the Qubit® 2.0 Fluorometer (Life Technologies, Carlsbad, CA, USA), and RNA sample quality was assessed with agarose gel electrophoresis. Novogene Bioinformatics Technology Co., Ltd, Beijing, China (http://www.novogene.cn/) performed RNA library construction. Sequencing libraries were generated using the NEBNext® Ultra™ Directional RNA Library Prep Kit for Illumina® (NEB) in accordance with the manufacturer’s protocol. After clustering index-coded samples, each library preparation was sequenced on an Illumina HiSeq 2000 platform and 50-bp single-end (SE) reads were generated.

Read mapping to the reference genome

Raw reads in the fastq format were first processed using in-house Perl scripts. they were then cleaned by removing of low-quality reads and those containing adapters or poly-N (quality limit 0.05). High-quality clean reads were used for downstream analyses. The Q20, Q30, and GC contents of the clean data were calculated. Next, clean tags were retained and mapped to the reference genome of Apis cerana (https://www.ncbi.nlm.nih.gov/genome/?term=Apis+cerana). TopHat was the selected as the mapping tool because it can generate a splice junction database based on the gene-model annotation file, leading to better mapping results than other non-splice-site mapping tools.

Bioinformatic analysis of RNA-seq data

Reads mapped to each gene were counted in HTSeq (version 0.5.4p3) [32]. Fragments per kilobase of transcript per million fragments sequenced (FPKM) of each gene were calculated from gene length and mapped read count [33]. Differential expression analysis of the two temperature treatments was performed using the DEGSeq R package v1.12.0 [34]. The Benjamini and Hochberg method was used to adjust P-values [35]. The significance threshold for differential expression was set at a corrected P-value < 0.05 and log2 (fold-change) of 0.5.

GO and KEGG pathway enrichment analyses

Gene ontology (GO) enrichment analysis of the differentially expressed genes (DEGs) was performed using the GOseq R package (version 1.10.0), with a correction for gene length bias. GO terms with a corrected P <0.05 were considered significantly enriched by DEGs. Significant enrichment (corrected P < 0.05) in a Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis was tested with KOBAS (http://kobas.cbi.pku.edu.cn/).

qRT-PCR analyses

Differential gene expression was verified with qRT-PCR (quantitative reverse transcription polymerase chain reaction). First-strand cDNA was synthesized from 1 μg of total RNA, using the ReverTra Ace-a First-strand cDNA Synthesis Kit (TaKaRa, Shiga Prefecture, Japan). Primers designed for each gene are given in S1 Table. The qRT-PCR thermocyling program was as follows: 95°C for 30 s; followed by 45 cycles of 95°C for 5 s and 60°C for 30 s. The β-actin gene was used as an internal standard for normalization. All samples were performed in triplicate. The results were normalized to the expression level of the constitutive Lbr. Relative quantification of PCR results was determined using the 22DDCt method [36].

Results

Differential phenotypes of ZOT and ZRT

In present study, we subjected two groups of bees to different temperatures (0°C and 25°C) for 2 h. Comparisons of activities over time revealed that physiological function of the 0°C group declined gradually, whereas it did not in the 25°C group. After 20 min, the honeybees under 0°C went into cold coma (abdominal stop peristalsis). At the end of different temperature treatments, we found no mortality in either group (n = 30/group).

Mapping of RNA-seq reads to the A. cerana genome

The three 0°C libraries (ZOT-1, ZOT-2, and ZOT-3) and three 25°C libraries (ZRT-1, ZRT-2, and ZRT-3) generated 41.56–59.70 million (M) raw reads per sample (Table 1). After quality control, the two groups averaged 44 667 109 and 49 380 011 clean reads, respectively. All short reads in both groups mapped perfectly to the A. cerana genome, and 70.91–73.92% of the tags in all six libraries could be uniquely mapped to the reference genome. Approximately 45% of reads were mapped to known exons and 27% were located in predicted intergenic or intronic regions.

thumbnail
Download:
Table 1. Statistics for filtering and mapping reads.

https://doi.org/10.1371/journal.pone.0179922.t001

Analysis of DEGs

Before DEG analysis, an RNA-seq correlation test (Pearson’s) and FPKM distribution were used to check sample variation across the two treatments. Strong positive correlations (R2 > 0.8; Fig 1A) were found between samples, indicating reasonable sample selection. The FPKM and density distributions showed that most genes were in the same groups (Fig 1B and 1C).

thumbnail
Download:
Fig 1. Bioinformatic analyses of RNA-seq data.

(A) Pearson correlation between samples. (B) FPKM density distribution of two groups. (C) FPKM distribution of two groups. ZOT: 0°C treatment; ZRT: 25°C treatment.

https://doi.org/10.1371/journal.pone.0179922.g001

Based on FPKM values, we detected 10 651 genes in the A. cerana genome. To better decipher the mechanisms driving the response to low temperatures, the DESeq tool was used to calculate differential expression. A loose transcript detected with p-adjusted < 0.05 and log2FoldChange was used to compare the ZOT library to the ZRT library (S2 Table). The ZOT group contained 269 significantly upregulated DEGs and 232 significantly downregulated DEGs (Fig 2A). The 501 DEGs and annotations found across the two groups are shown in S3 Table. Venn diagram analysis revealed that 8614 genes overlapped between the two groups, while 308 genes were expressed only in ZOT and 165 were expressed only in ZRT. Cluster analysis showed that DEGs between ZOT and ZRT were divided into three groups: upregulated in the first, downregulated in the second, and with varied expression in the third (Fig 2C).

thumbnail
Download:
Fig 2. Bioinformatic analyses of DEGs (differentially expressed genes).

(A) DEGs distribution of the two treatment groups (ZOT: 0°C, ZRT: 25°C). (B) Venn diagrams showing the number of genes expressed in the two groups. (C) Cluster analysis of DEGs.

https://doi.org/10.1371/journal.pone.0179922.g002

Gene annotation identified differential expression between the two treatments included nine genes encoding heat shock proteins (HSPs), seven encoding serine/threonine-protein kinase genes (STKs), and twelve encoding zinc finger protein (ZFPs) (S4 Table).

GO enrichment analysis of DEGs

The GOseq tool was used for GO analysis of DEGs, based on the Wallenius non-central hyper-geometric distribution [37]. Upregulated DEGs mostly enriched DNA binding, gated channel activity, as well as biosynthesis and metabolism regulation (Fig 3A). Down-regulated DEGs were mostly involved in anion binding, oxidoreductase activity, mitochondrial matrix, lipid metabolism, polysaccharide biosynthesis, and steroid metabolic processes (Fig 3B). More cold-tolerance-related GO terms were associated with upregulated than down-regulated DEGs. These cold-tolerance-related terms were associated with biological processes (energy biosynthesis and metabolism [of glycogen, nucleic acids, amino acids, polyol]; stimulus response, movement, ion transport, signal transduction), molecular function (binding; activities of gate channels, calcium channels, kinases, transmembrane transporters), and cellular components (microtubule cytoskeleton, cell, cell part, protein complex).

thumbnail
Download:
Fig 3. Bar graph showing gene ontology (GO) of DEGs between ZOT (0°C) and ZRT (25°C).

(A) The most enriched GO terms for upregulated DEGs between the two treatment groups. (B) The most enriched GO terms for downregulated DEGs between the two treatment groups.

https://doi.org/10.1371/journal.pone.0179922.g003

KEGG analysis of DEGs

KEGG pathway enrichment analysis [38] of DEGs revealed that these transcripts were assigned to 63 reference pathways (S5 Table). Among them, metabolic pathways were the most enriched. These DEGs enriched glycine, serine, and threonine metabolism; carbon metabolism; purine metabolism; glycerophospholipid metabolism; inositol phosphate metabolism; fatty acid metabolism; alanine, aspartate, and glutamate metabolism; amino sugar and nucleotide sugar metabolism; as well as cysteine and methionine metabolism. Additionally, many enriched pathways were involved in signal transduction pathways related to cell differentiation and proliferation, repair of injured organs, as well as regulation of cytosolic calcium concentration (e.g., Wnt signaling pathway, Hippo signaling pathway–fly, FOXO signaling pathway, Hedgehog signaling pathway, and phosphatidylinositol signaling system) [39, 40].

Validation of RNA-seq by qRT-PCR

To validate the expression profiles identified from RNA-seq, we performed qRT-PCR on several cold-stress-related DEGs (nine HSPs, seven STKs, and twelve ZFPs). The comparative analysis showed that qRT-PCR-detected expression patterns supported RNA-seq results, with only one DEG (in the STK protein family) being inconsistent (Table 2). However, the fold change of expression patterns from RNA-seq and qRT-PCR still deviated slightly, probably due to methodological differences [41].

thumbnail
Download:
Table 2. Results of qRT-PCR and RNA-seq on genes encoding cold-resistance proteins.

https://doi.org/10.1371/journal.pone.0179922.t002

Discussion

Insects vary considerably in their ability to survive low temperatures [14]. Insect susceptibility to physiological experimentation increases the difficulty of exploring cold-resistance mechanisms in the honey bee. However, with the availability of the Eastern honey bee (A. cerana) whole genome [42] and the development of RNA-seq, we now have effective alternatives for studying the molecular mechanisms of cold resistance.

This study identified DEGs between honey bees treated under two temperature conditions. Our data are in line with existing research on other insects. Previously, a DNA microarray study in Drosophila melanogaster identified 36 transcripts associated with low-temperature treatments (0°C for 2 h) [18]. Furthermore, an RNA-seq study in New Zealand stick insects showed that genes expressed under cold exposure (-5°C for 1 h) included HSPs, cuticular genes, as well as others involved in carbohydrate metabolism and physiological variation [43]. Transcriptomic and proteomic analysis in Chinese white wax scale insects (Ericerus pela) from different climatic regions found 2386 DEGs, involved in energy metabolism, catalytic activity, and response to stimuli [13].

To further investigate the molecular mechanism of cold resistance in A. c. cerana, we used RNA-seq to compare DEG transcript levels. We demonstrated that the number of upregulated DEGs (269) was slightly higher than the number of downregulated DEGs (232) between the two treatment groups. Next, KEGG analysis revealed that metabolic pathways (including carbon, amino acid, and nucleic acid metabolism) were the most enriched. Thus, energy metabolism was critical to cold resistance in honey bees, in line with the observation that honey bees consume honey and pollen to produce heat during low temperatures [23, 26]. Data from physiological studies of various insects [20, 44] also support these results. Finally, signal transduction pathways are important in cold stress responses [45], and in this study, four such pathways were found to be enriched.

Based on GO analysis, the DEGs were most enriched in biological processes, and specifically in metabolic processes. As the main energy sources for honey bees are glucose, fructose, maltose, and sucrose, it was unsurprising that the detected DEGs were involved in glucose metabolism. These included mitogen-activated protein kinase-binding protein 1 (MAPKBP1), poly(ADP-ribose) glycohydrolase (PARG), UDP-glucuronosyltransferase 1–8 (UTG 1–8), and maltase A3 (Mal-A3). MAPKBP1 is a key transcription factor in the NF-κB signaling pathway [46]. The role of insects MAPKBP1remains unclear, but a previous study found that endogenous MAPKBP1 and scaffold protein genes co-localize to stress granules; thus, the expression of this gene might be related to cold tolerance in honey bees [47]. RAPG is a hydrolyzing enzyme and regulator of poly(ADP-ribose) activity, producing free ADP-ribose residues to affect several important biological functions, such as DNA repair, chromosome stability, and maintenance of normal cell function [48]. Like other UTG1 members, UTG 1–8 is involved in glucuronidation, a process that enhances the conversion of lipophilic xenobiotics and endobiotics to more water-soluble compounds. For example, UTG 1–8 can split UDP-glucuronic acid and catalyze glucuronic acid transfer to a small hydrophobic molecule [49]. Finally, Mal-A3 is a maltase that catalyzes maltose hydrolysis to glucose [50].

Free amino acids have a strong affinity to small molecular compounds, leading to increased cell fluid concentration, as well as improved cell water retention and stability [1]. Indeed, high accumulation of free amino acids in the hemolymph appears to be a major biochemical characteristic of insect cold-stress response [20–22]. Corroborating those findings, we identified DEGs associated with serine, lysine, glycine, threonine, and methionine metabolism. Two upregulated and zero downregulated genes were linked to methionine and threonine, while four downregulated and zero upregulated genes were linked to serine, lysine, and glycine. These results suggest that methionine and threonine may act as cryoprotectants in honey bees. Interestingly, although lipids are another major energy source, lipid metabolism was most strongly enriched in downregulated DEGs in this study. Lipids are an essential energy source in cell membrane synthesis, and energy from lipid metabolism is primarily allocated to long-distance travel in insects [51]. Therefore, lipid metabolism may be downregulated in the cold because honey bee flight ability is hampered. However, lipids may power honey bees to flight at the appropriate temperature, while maintaining membrane structure in the cold.

A major biological signal to cold response is variation in intracellular calcium ion concentration [52]. Upregulated GO terms related to calcium channel activity included voltage-dependent calcium channel subunit alpha-2/delta-3 (CACNA2D3) and small conductance calcium-activated potassium channel protein (SK). CACNA2D3 is an auxiliary member of the alpha-2/delta subunit family in the voltage-dependent calcium channel complex; it regulates calcium ion influx upon membrane polarization [53]. SKs are gated solely through an increase in internal calcium concentration and are therefore extremely sensitive to calcium. Our observation of SK and CACNA2D3 expression indicated that greater calcium ion activity and increased intracellular calcium concentration improves cold tolerance in honey bees. Our results are similar to previous findings in D. melanogaster [52].

The DEGs related to metabolic processes and the calcium channel activity terms identified indicate that some small cryoprotectants (including glucose, methionine, threonine and calcium) may play an important role in the cold-coping mechanisms of the honey bee. We can apply the information gained from our study to optimize honey bee diet in the winter. Currently, the most important supplement for honey bees in the winter in most China northern states is sucrose (processed white table sugar), while the addition of free amino acids and minerals in diets has often been overlooked. The results of this study highlight the importance and necessity of dietary supplementation with small cryoprotectants. The partial substitution of invert sugar (its main component sugars are glucose and fructose) for sucrose and extra methionine, threonine and calcium added in the production of stored, it could improve overwintering success and decrease the loss of honey bees.

Antifreeze-associated proteins, unlike small-molecule cryoprotectants, can induce ideal antifreeze effects through catalyzing specific biochemical reaction and affecting the conformation of the resulting molecules [54–56]. In this research, some DEGs were linked to three such proteins (HSPs, STKs, and ZFPs) and can be considered candidate genes for subsequent studies on honey bee cold-tolerance mechanisms. Additionally, based on these candidate genes, we can generate the corresponding proteins in vitro as a protectant against further cold stress.

HSPs are a super family of chaperone proteins rapidly biosynthesized in response to various environmental stress factors [57, 58]. We found nine DEGs associated with HSPs: two upregulated (heat shock factor 2-binding protein [HSF2BP] and heat shock factor protein 5 [HSF5]) and seven downregulated (heat shock protein 10 kDa [HSP10], heat shock protein sHSP24.2a [sHSP24.2a], heat shock protein sHSP22.6 [sHSP22.6], heat shock protein 60 kDa [HSP60], heat shock protein 90 kDa [HSP90], heat shock 70 kDa protein cognate 3 [HSC 70–3], heat shock 70 kDa protein cognate 4 [HSC 70–4]). As in other insects, the expression of HSPs can be induced quickly with low temperatures. Recently, a study in Sitodiplosis mosellana suggested that HSP90 and HSP70 expression at low temperatures (-5°C) was higher than under normal conditions [59]. Further, rapid cold hardening significantly upregulated three HSP genes (hsp40, hsp23, and hsp10) in Folsomia candida [60]. However, after prolonged cold exposure, HSPs expression levels would drop [61], indicating that the resultant state of the honey bees in the 0°C treatment was not induced by any new stimuli, and that a self-defense mechanism may have been established.

Unlike typical HSPs, HSC70 function varies across multiple species and treatments. Some studies suggest that HSC70 is expressed constitutively and thus, not stress-inducible [62]. Others report that HSC70 is involved in thermoregulation, as well as protein synthesis, folding, transport, and degradation [63]. Our results for HSC70 indicated that its expression is also sensitive to low temperature. Additionally, HSFs mediate HSPs induction through binding to the heat shock element of HSP genes, explaining why HSFs mRNA levels clearly increased in adult fruit flies exposed to 4°C [64]. HSF5 is part of a protein family that modulates HSPs expression; therefore, it is likely related to thermoregulation. Indeed, HSF5 expression was significantly higher in winter than in summer [65]. Further, HSF2BP primarily modulates HSF2 activation; an increase in its transcripts can be induced following a heat-shock treatment, but it also appears to be the most insensitive to temperature change among the HSPs [66]. Overall, our results show that lower temperatures downregulated HSPs and upregulated HSFs in honey bees, similar to findings in D. melanogaster [18]. These outcomes suggest that HSPs expression and activity varies in response to cold stress.

The zinc finger protein family (so-called due to the presence of the zinc finger domain) is large and widely distributed across taxa (plants, animals, microorganisms) [54]. These proteins are involved in growth and development as well as in abiotic stress responses (e.g., to salt, water, heat, and cold) [67, 68]. Recently, three differentially upregulated ZFPs were identified in heat-treated silkworm eggs [41]. This study also identified twelve ZFPs, including seven upregulated genes and five downregulated genes. Our results together with those of previous studies demonstrate that ZFPs are important in the cold-stress response of insects and that high expression of some specific ZFPs improves cold tolerance. Additionally, we found seven DEGs related to STKs, a major protein kinase group [69] that is implicated in modulating abiotic stress responses (including the cold) [70, 71]. However, research on insect STKs and ZFPs has been limited. Thus, our current results and their implications for honey bees should be considered with caution. The roles of differentially expressed STK and ZFP genes in the molecular mechanisms underlying honey bee cold-stress response warrant further study.

Conclusions

In conclusion, 501 DEGs were identified across two temperature treatments. The results of GO analysis revealed that DEGs significantly enriched terms associated with sugar and amino acid biosynthesis and metabolism, as well as with calcium ion channel activity. Various temperature-sensitive proteins (e.g., HSPs, STKs, and ZFPs) were also differentially expressed across temperature treatments. Our study provides insight into the molecular mechanisms of cold resistance in A. c. cerana and offers a scientific foundation for exploring candidate genes linked to cold stress in social insects.

Supporting information

S1 Table. Primers used for qRT-PCR analysis on differentially expressed genes (DEGs).

https://doi.org/10.1371/journal.pone.0179922.s001

(DOC)

S2 Table. DEGs between ZOT and ZRT, based on corrected P values and log2FoldChange.

https://doi.org/10.1371/journal.pone.0179922.s002

(XLS)

S3 Table. Annotation of DEGs in NR, NT, and SwissProt databases.

https://doi.org/10.1371/journal.pone.0179922.s003

(XLS)

S4 Table. Annotation of HSPs, STKs, and ZFPs.

https://doi.org/10.1371/journal.pone.0179922.s004

(XLS)

S5 Table. The most enriched KEGG pathway for DEGs between ZOT and ZRT.

https://doi.org/10.1371/journal.pone.0179922.s005

(XLS)

Acknowledgments

We would like to acknowledge the other members of the lab, Shuang Yang and Lina Guo, for their help with the experiments. We greatly appreciate the assistance of the scientists of Novogene Bioinformatics Technology Co., Ltd in sequencing the DEG libraries. We also would like to thank Editage for English language editing.

Author Contributions

  1. Conceptualization: KX YJ.
  2. Data curation: KX YD.
  3. Formal analysis: KX HZ.
  4. Funding acquisition: QN YJ.
  5. Investigation: KX YJ.
  6. Methodology: KX YJ.
  7. Project administration: QN YJ.
  8. Resources: QN YJ.
  9. Software: KX HZ.
  10. Supervision: QN YJ.
  11. Validation: QN KX.
  12. Visualization: KX.
  13. Writing – original draft: KX.
  14. Writing – review & editing: KX.

References

  1. 1. Doucet D, Walker VK, Qin W. The bugs that came in from the cold: molecular adaptations to low temperatures in insects. Cellular and molecular life sciences: CMLS. 2009;66(8):1404–18. Epub 2009/01/09. pmid:19129970.
  2. 2. Colinet H, Renault D, Roussel D. Cold acclimation allows Drosophila flies to maintain mitochondrial functioning under cold stress. Insect Biochem Mol Biol. 2016;80:52–60. pmid:27903433.
  3. 3. Hazell SP, Groutides C, Neve BP, Blackburn TM, Bale JS. A comparison of low temperature tolerance traits between closely related aphids from the tropics, temperate zone, and Arctic. J Insect Physiol. 2010;56(2):115–22. pmid:19723528.
  4. 4. Sorensen JG, Norry FM, Scannapieco AC, Loeschcke V. Altitudinal variation for stress resistance traits and thermal adaptation in adult Drosophila buzzatii from the New World. J Evol Biol. 2005;18(4):829–37. pmid:16033554.
  5. 5. Wang XH, Qi XL, Kang L. Geographic differences on accumulation of sugars and polyols in locust eggs in response to cold acclimation. J Insect Physiol. 2010;56(8):966–70. pmid:20416314.
  6. 6. Carrillo MA, Heimpel GE, Moon RD, Cannon CA, Hutchison WD. Cold hardiness of Habrobracon hebetor (Say) (Hymenoptera: Braconidae), a parasitoid of pyralid moths. J Insect Physiol. 2005;51(7):759–68. pmid:15885701.
  7. 7. Goto M, Sekine Y, Outa H, Hujikura M, Koichi S. Relationships between cold hardiness and diapause, and between glycerol and free amino acid contents in overwintering larvae of the oriental corn borer, Ostrinia furnacalis. J Insect Physiol. 2001;47(2):157–65. Epub 2000/11/07. pmid:11064022.
  8. 8. Larsen KJ, Lee RE. Cold tolerance including rapid cold-hardening and inoculative freezing of fall migrant monarch butterflies in Ohio. Journal of Insect Physiology. 1994;40(10):859–64. http://dx.doi.org/10.1016/0022-1910(94)90019-1.
  9. 9. Berkvens N, Bale JS, Berkvens D, Tirry L, De Clercq P. Cold tolerance of the harlequin ladybird Harmonia axyridis in Europe. J Insect Physiol. 2010;56(4):438–44. pmid:19951708.
  10. 10. Kostal V, Simek P, Zahradnickova H, Cimlova J, Stetina T. Conversion of the chill susceptible fruit fly larva (Drosophila melanogaster) to a freeze tolerant organism. Proc Natl Acad Sci U S A. 2012;109(9):3270–4. Epub 2012/02/15. pmid:22331891; PubMed Central PMCID: PMCPmc3295325.
  11. 11. Zhao YX, Kang L. Cold tolerance of the leafminer Liriomyza sativae (Dipt., Agromyzidae). Journal of Applied Entomology. 2000;124(3–4):185–9.
  12. 12. Kang L, Chen B, Wei JN, Liu TX. Roles of thermal adaptation and chemical ecology in Liriomyza distribution and control. Annual review of entomology. 2009;54:127–45. Epub 2008/08/20. pmid:18710304.
  13. 13. Yu SH, Yang P, Sun T, Qi Q, Wang XQ, Chen XM, et al. Transcriptomic and proteomic analyses on the supercooling ability and mining of antifreeze proteins of the Chinese white wax scale insect. Insect Sci. 2016;23(3):430–7. pmid:26799455.
  14. 14. Clark MS, Worland MR. How insects survive the cold: molecular mechanisms-a review. Journal of comparative physiology B, Biochemical, systemic, and environmental physiology. 2008;178(8):917–33. Epub 2008/06/28. pmid:18584182.
  15. 15. Dahlhoff EP, Rank NE. The role of stress proteins in responses of a montane willow leaf beetle to environmental temperature variation. J Biosci. 2007;32(3):477–88. Epub 2007/05/31. pmid:17536167.
  16. 16. King AM, MacRae TH. Insect heat shock proteins during stress and diapause. Annual review of entomology. 2015;60:59–75. Epub 2014/10/24. pmid:25341107.
  17. 17. Stetina T, Kostal V, Korbelova J. The Role of Inducible Hsp70, and Other Heat Shock Proteins, in Adaptive Complex of Cold Tolerance of the Fruit Fly (Drosophila melanogaster). PLoS One. 2015;10(6):e0128976. Epub 2015/06/04. pmid:26034990; PubMed Central PMCID: PMCPmc4452724.
  18. 18. Qin W, Neal SJ, Robertson RM, Westwood JT, Walker VK. Cold hardening and transcriptional change in Drosophila melanogaster. Insect Mol Biol. 2005;14(6):607–13. pmid:16313561.
  19. 19. Chen C-P, Walker VK. Cold-shock and chilling tolerance in Drosophila. Journal of Insect Physiology. 1994;40(8):661–9. http://dx.doi.org/10.1016/0022-1910(94)90093-0.
  20. 20. Holden CP, Storey KB. 6-phosphogluconate dehydrogenase from a freeze tolerant insect: Control of the hexose monophosphate shunt and NADPH production during cryprotectant synthesis. Insect Biochemistry & Molecular Biology. 1994;24(2):167–73.
  21. 21. Chen YJ, Sun XG, Zhang WG, Mu ZG, Guo GZ. Relation Between Variation of Water, Fat, Glycerol in vivo of Over-wintering Diaphania pyloalis Walker Larvae and Cold-hardiness. Canye Kexue. 2005;31(2):111–6.
  22. 22. Huang GY, Wang YC, You ZP. Study on Cold Tolerance Functions of Agrotis segetum. Journal of Zhejiang Forestry College. 1990;(2):140–6.
  23. 23. Heinrich B. Mechanisms of body-temperature regulation in honeybees, Apis mellifera. II. Regulation of thoracic temperature at high air temperatures. Physics Letters A. 1980;29(9):551–.
  24. 24. Stabentheiner A, Kovac H, Hagmüller K. Thermal behavior of round and wagtail dancing honeybees. Journal of Comparative Physiology B. 1995;165(6):433–44.
  25. 25. Owens CD. The thermology of wintering honey bee colonies. Technical Bulletins. 1971.
    • 26. Ritter W. Experimenteller Beitrag zur Thermoregulation des Bienenvolks (Apis mellifera L.). Apidologie. 1982;13(2):169–95. https://doi.org/10.1051/apido:19820205.
    • 27. Coelho JR. Heat Transfer and Body Temperature in Honey Bee (Hymenoptera: Apidae) Drones and Workers. Environmental Entomology. 1991;20(6):1627–35. https://doi.org/10.1093/ee/20.6.1627.
    • 28. Blumberg H, Silbersweig D, Martinez D, Ricketts S, Asis JD, Mcbride A, et al. Thermoregulation of African and European honeybees during foraging, attack, and hive exits and returns. Breast Cancer Research & Treatment. 1979;139(1):163–76.
    • 29. Liu Z, Yao P, Guo X, Xu B. Two small heat shock protein genes in Apis cerana cerana: characterization, regulation, and developmental expression. Gene. 2014;545(2):205–14. Epub 2014/05/20. pmid:24835315.
    • 30. LI ZY, Xue YB, Wang Z, Chang ZG. Study on the Super cooling Point of Honeybees and the Relationship between the Super cooling Point and the Cold Resistance. Journal of Bee. 2006;26(11):12–3.
    • 31. Yang GH. Harm of introducing the western honeybee Apis mellifera L. to the Chinese honeybee Apis cerana F. and its ecological impact. Acta Entomologica Sinica. 2005;48(3):401–6.
    • 32. Feng Y, van der Veeken J, Shugay M, Putintseva EV, Osmanbeyoglu HU, Dikiy S, et al. A mechanism for expansion of regulatory T-cell repertoire and its role in self-tolerance. Nature. 2015;528(7580):132–6. Epub 2015/11/26. pmid:26605529; PubMed Central PMCID: PMCPmc4862833.
    • 33. Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, et al. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010;28(5):511–5. Epub 2010/05/04. pmid:20436464; PubMed Central PMCID: PMCPmc3146043.
    • 34. Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010;11(10):R106. Epub 2010/10/29. pmid:20979621; PubMed Central PMCID: PMCPmc3218662.
    • 35. Benjamini Y, Hochberg Y. Controlling The False Discovery Rate—A Practical And Powerful Approach To Multiple Testing. Journal of the Royal Statistical Society. 1995;57(57):289–300.
    • 36. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method. Methods. 2001;25(4):402–8. Epub 2002/02/16. pmid:11846609.
    • 37. Young MD, Wakefield MJ, Smyth GK, Oshlack A. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 2010;11(2):R14. Epub 2010/02/06. pmid:20132535; PubMed Central PMCID: PMCPmc2872874.
    • 38. Kanehisa M, Araki M, Goto S, Hattori M, Hirakawa M, Itoh M, et al. KEGG for linking genomes to life and the environment. Nucleic Acids Res. 2008;36(Database issue):D480–4. Epub 2007/12/14. pmid:18077471; PubMed Central PMCID: PMCPmc2238879.
    • 39. Patel J, Keith RA, Salama AI, Moore WC. Role of calcium in regulation of phosphoinositide signaling pathway. Journal of Molecular Neuroscience. 1991;3(1):19. pmid:1653000
    • 40. Logan CY, Nusse R. The Wnt signaling pathway in development and disease. Annual Review of Cell & Developmental Biology. 2004;20(20):781.
    • 41. Liu P, Wang Y, Du X, Yao L, Li F, Meng Z. Transcriptome Analysis of Thermal Parthenogenesis of the Domesticated Silkworm. PLoS One. 2015;10(8):e0135215. pmid:26274803; PubMed Central PMCID: PMCPMC4537240.
    • 42. Park D, Jung JW, Choi BS, Jayakodi M, Lee J, Lim J, et al. Uncovering the novel characteristics of Asian honey bee, Apis cerana, by whole genome sequencing. BMC genomics. 2015;16:1. Epub 2015/01/03. pmid:25553907; PubMed Central PMCID: PMCPmc4326529.
    • 43. Dennis AB, Dunning LT, Sinclair BJ, Buckley TR. Parallel molecular routes to cold adaptation in eight genera of New Zealand stick insects. Scientific reports. 2015;5:13965. Epub 2015/09/12. pmid:26355841; PubMed Central PMCID: PMCPmc4564816.
    • 44. Chen CP, Walker VK. Increase in cold-shock tolerance by selection of cold resistant lines in Drosophila melanogaster. Ecological Entomology. 1993;18(3):184–90.
    • 45. Janska A, Marsik P, Zelenkova S, Ovesna J. Cold stress and acclimation—what is important for metabolic adjustment? Plant biology (Stuttgart, Germany). 2010;12(3):395–405. Epub 2010/06/05. pmid:20522175.
    • 46. Fu L, Shi J, Hu K, Wang J, Wang W, Ke X. Mitogen-activated protein kinase binding protein 1 (MAPKBP1) is an unfavorable prognostic biomarker in cytogenetically normal acute myeloid leukemia. Oncotarget. 2015;6(10):8144–54. Epub 2015/04/30. pmid:25924238; PubMed Central PMCID: PMCPmc4480741.
    • 47. Cohen-Katsenelson K, Wasserman T, Darlyuk-Saadon I, Rabner A, Glaser F, Aronheim A. Identification and analysis of a novel dimerization domain shared by various members of c-Jun N-terminal kinase (JNK) scaffold proteins. The Journal of biological chemistry. 2013;288(10):7294–304. Epub 2013/01/24. pmid:23341463; PubMed Central PMCID: PMCPmc3591637.
    • 48. Hanai S, Kanai M, Ohashi S, Okamoto K, Yamada M, Takahashi H, et al. Loss of poly(ADP-ribose) glycohydrolase causes progressive neurodegeneration in Drosophila melanogaster. Proc Natl Acad Sci U S A. 2004;101(1):82–6. Epub 2003/12/17. pmid:14676324; PubMed Central PMCID: PMCPmc314142.
    • 49. King CD, Rios GR, Green MD, Tephly TR. UDP-glucuronosyltransferases. Current drug metabolism. 2000;1(2):143–61. Epub 2001/07/24. pmid:11465080.
    • 50. Wanke V, Vavassori M, Thevelein JM, Tortora P, Vanoni M. Regulation of maltose utilization in Saccharomyces cerevisiae by genes of the RAS/protein kinase A pathway. FEBS Lett. 1997;402(2–3):251–5. Epub 1997/02/03. pmid:9037205.
    • 51. Nation JL. Insect physiology and biochemistry. Insect Physiology & Biochemistry. 2002.
      • 52. Takeuchi K, Nakano Y, Kato U, Kaneda M, Aizu M, Awano W, et al. Changes in temperature preferences and energy homeostasis in dystroglycan mutants. Science. 2009;323(5922):1740–3. Epub 2009/03/28. pmid:19325118.
      • 53. Gong HC, Hang J, Kohler W, Li L, Su TZ. Tissue-specific expression and gabapentin-binding properties of calcium channel alpha2delta subunit subtypes. The Journal of membrane biology. 2001;184(1):35–43. Epub 2001/11/01. pmid:11687876.
      • 54. Miller J, McLachlan AD, Klug A. Repetitive zinc-binding domains in the protein transcription factor IIIA from Xenopus oocytes. The EMBO journal. 1985;4(6):1609–14. Epub 1985/06/01. pmid:4040853; PubMed Central PMCID: PMCPmc554390.
      • 55. Molle V, Kremer L. Division and cell envelope regulation by Ser/Thr phosphorylation: Mycobacterium shows the way. Mol Microbiol. 2010;75(5):1064–77. Epub 2010/05/22. pmid:20487298.
      • 56. Bahar MH, Hegedus D, Soroka J, Coutu C, Bekkaoui D, Dosdall L. Survival and hsp70 gene expression in Plutella xylostella and its larval parasitoid Diadegma insulare varied between slowly ramping and abrupt extreme temperature regimes. PLoS One. 2013;8(9):e73901. Epub 2013/09/17. pmid:24040110; PubMed Central PMCID: PMCPmc3765401.
      • 57. Bahar MH, Hegedus D, Soroka J, Coutu C, Bekkaoui D, Dosdall L. Survival and hsp70 gene expression in Plutella xylostella and its larval parasitoid Diadegma insulare varied between slowly ramping and abrupt extreme temperature regimes. PLoS ONE. 2013;8(9):-.
      • 58. Liberek K, Lewandowska A, Zietkiewicz S. Chaperones in control of protein disaggregation. Embo Journal. 2008;27(2):328–35. pmid:18216875
      • 59. Cheng W, Li D, Wang Y, Liu Y, Zhu-Salzman K. Cloning of heat shock protein genes (hsp70, hsc70 and hsp90) and their expression in response to larval diapause and thermal stress in the wheat blossom midge, Sitodiplosis mosellana. J Insect Physiol. 2016;95:66–77. Epub 2016/10/21. pmid:27639943.
      • 60. Waagner D, Holmstrup M, Bayley M, Sorensen JG. Induced cold-tolerance mechanisms depend on duration of acclimation in the chill-sensitive Folsomia candida (Collembola). The Journal of experimental biology. 2013;216(Pt 11):1991–2000. Epub 2013/02/09. pmid:23393277.
      • 61. Sivan A, Shriram AN, Muruganandam N, Thamizhmani R. Expression of heat shock proteins (HSPs) in Aedes aegypti (L) and Aedes albopictus (Skuse) (Diptera: Culicidae) larvae in response to thermal stress. Acta tropica. 2016;167:121–7. Epub 2016/12/28. pmid:28024869.
      • 62. Mahroof R, Yan ZK, Neven L, Subramanyam B, Bai J. Expression patterns of three heat shock protein 70 genes among developmental stages of the red flour beetle, Tribolium castaneum (Coleoptera: Tenebrionidae). Comparative Biochemistry & Physiology Part A Molecular & Integrative Physiology. 2005;141(2):247–56.
      • 63. Yang Jing, Zhao Yang, Sheng Ying-Fang, et al. Identification of heat shock cognate protein 70 gene (Alhsc70) of Apolygus lucorum and its expression in response to different temperature and pesticide stresses. Insect Science. 2016;23(1):37–49. pmid:25448821
      • 64. Fujikake N, Nagai Y, Popiel HA, Kano H, Yamaguchi M, Toda T. Alternative splicing regulates the transcriptional activity of Drosophila heat shock transcription factor in response to heat/cold stress. Febs Letters. 2005;579(17):3842–8. pmid:15978579
      • 65. Lal SV, Brahma B, Gohain M, Mohanta D, De BC, Chopra M, et al. Splice variants and seasonal expression of buffalo HSF genes. Cell Stress and Chaperones. 2015;20(3):545–54. pmid:25655489
      • 66. Fu S, Rogowsky P, Nover L, Scanlon MJ. The maize heat shock factor-binding protein paralogs EMP2 and HSBP2 interact non-redundantly with specific heat shock factors. Planta. 2006;224(1):42–52. pmid:16331466
      • 67. Feurtado JA, Huang D, Wicki-Stordeur L, Hemstock LE, Potentier MS, Tsang EW, et al. The Arabidopsis C2H2 zinc finger INDETERMINATE DOMAIN1/ENHYDROUS promotes the transition to germination by regulating light and hormonal signaling during seed maturation. Plant Cell. 2011;23(5):1772–94. pmid:21571950; PubMed Central PMCID: PMCPMC3123948.
      • 68. Ben Saad R, Zouari N, Ben Ramdhan W, Azaza J, Meynard D, Guiderdoni E, et al. Improved drought and salt stress tolerance in transgenic tobacco overexpressing a novel A20/AN1 zinc-finger "AlSAP" gene isolated from the halophyte grass Aeluropus littoralis. Plant Mol Biol. 2010;72(1–2):171–90. Epub 2009/10/20. pmid:19838809.
      • 69. Faucher SP, Viau C, Gros PP, Daigle F, Le Moual H. The prpZ gene cluster encoding eukaryotic-type Ser/Thr protein kinases and phosphatases is repressed by oxidative stress and involved in Salmonella enterica serovar Typhi survival in human macrophages. FEMS Microbiol Lett. 2008;281(2):160–6. Epub 2008/03/04. pmid:18312574.
      • 70. Rudrabhatla P, Reddy MM, Rajasekharan R. Genome-wide analysis and experimentation of plant serine/ threonine/tyrosine-specific protein kinases. Plant Mol Biol. 2006;60(2):293–319. Epub 2006/01/24. pmid:16429265.
      • 71. Rudrabhatla P, Rajasekharan R. Developmentally regulated dual-specificity kinase from peanut that is induced by abiotic stresses. Plant Physiol. 2002;130(1):380–90. Epub 2002/09/13. pmid:12226517; PubMed Central PMCID: PMCPmc166570.