Comparative analysis of the transcriptome of the overwintering desert beetle Microdera punctipennis

Accepted Manuscript Comparative analysis of the transcriptome of the overwintering desert beetle Microdera punctipennis Kuerban Tusong, Xiaoxing Guo, Shanshan Meng, Xiaoning Liu, Ji Ma PII:

S0011-2240(17)30069-X

DOI:

10.1016/j.cryobiol.2017.06.009

Reference:

YCRYO 3861

To appear in:

Cryobiology

Received Date: 22 February 2017 Revised Date:

27 June 2017

Accepted Date: 27 June 2017

Please cite this article as: K. Tusong, X. Guo, S. Meng, X. Liu, J. Ma, Comparative analysis of the transcriptome of the overwintering desert beetle Microdera punctipennis, Cryobiology (2017), doi: 10.1016/j.cryobiol.2017.06.009. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Comparative analysis of the transcriptome of the overwintering desert beetle Microdera punctipennis

RI PT

Kuerban Tusong, Xiaoxing Guo, Shanshan Meng, Xiaoning Liu, Ji Ma*

Xinjiang Key Laboratory of Biological Resources and Genetic Engineering, College

M AN U

830046, China

SC

of Life Science and Technology, Xinjiang University, 666 Shengli Road, Urumqi

*Corresponding author at: Xinjiang Key Laboratory of Biological Resources and

EP

TE

Urumqi 830046, China

D

Genetic Engineering, College of Life Science and Technology, Xinjiang University,

AC C

E-mail address: [email protected] (K. Tusong), [email protected] (J. Ma)

1

ACCEPTED MANUSCRIPT 1

Abstract The cold tolerance mechanisms of insect have been studied extensively on the model species Drosophila and a few other species at the transcriptional level. However studies on insects that inherit strong cold tolerance are limited. Cold hardy Tenebrionid beetle Microdera punctipennis is endemic to Gurbantonggut Desert, northwest of China. However，its genomic information is lacking. To investigate the overwintering mechanisms of M. punctipennis adult, RNA-seq was performed on the winter adults and the control adults that were kept in laboratory at 30 ℃. A total of 175,247 unigenes were acquired with an average length of 645 bp. By using DESeq package, we identified 3,367 unigenes that were up-regulated and 7,988 down-regulated in the winter adults compared with the controls. To further our understanding of these differentially expressed genes (DEGs), Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed. Pathway analysis showed that the “ECM-receptor interaction”, “PI3K-Akt signaling pathway”, “Estrogen signaling pathway”, “Tight junction”, and “Regulation of actin cytoskeleton”, etc. might play important roles in M. punctipennis overwintering. The DEGs results from the RNA-Seq were confirmed partially by qRT-PCR for 13 DEGs, which showed high consistence with a Pearson's correlation coefficient of 0.851. Overall, the sequence data will provide basic information for subsequent bioinformatical analysis and mining of the genes responsible for cold tolerance in M. punctipennis, as well as for understanding the molecular mechanisms of desert beetle overwintering.

23

Key words: Tenebrionidae, Transcriptome, Low temperature, Overwintering,

24

Beetle

25

1. Introduction

26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41

Temperature affects the seasonal activity and geographic distribution of insects [7], and hence can impact the dynamics of native insect populations, spread of invasive pests, and the success of species introduced as biological control [4]. Insects in desert are regularly exposed to adverse environmental conditions including cold temperatures [39-41]. Thus, cold hardiness is a primary component of insect fitness, and one of the determinants of insect distribution [3]. Desert insects can quickly adjust themselves to the change of environmental temperature by reducing metabolic activity, building-up metabolic wastes, as well as desiccating [49]. In subzero conditions some physiological activities are still active such as cryoprotectants synthesis [48, 56] and diapause development [20]. Among beetles, Tenebrionidae is one of the most successful animals in desert [37], they adopted seasonal behavioral changes to avoid inhospitable environment [25, 58]. The beetle Microdera punctipennis, one of the Tenebrionidae, is an endemic species to the Gurbantonggut Desert, Xinjiang, northwest of China [18]. The desert environment is characterized by dry and extreme temperatures with -40 ℃ of the minimum air temperature and -22 ℃ of the average low air temperature in winter.

AC C

EP

TE

D

M AN U

SC

RI PT

2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22

2

ACCEPTED MANUSCRIPT The mean annual air temperature varies between 6 ℃ to 10 ℃ [1]. M. punctipennis is flightless, and night active [58, 60]. It lives in the desert adverse environment by hiding under the litter and burrowing into the loose sand. This type of beetle is freeze avoidant, it could lower its super cooling point below -19.6 ℃ in January [17]. It passes winters in dry sand about 5 cm or deeper beneath the earth surface. The average temperatures of the soil surface and soil-in-5cm in January were -12℃ and -5℃, respectively [17]. We have studied several cold hardiness related genes in this insect [28, 32, 38, 46, 59], but the underlying molecular mechanisms that make the beetle overwintering are not fully understood. RNA-Seq technology has been successfully applied to discover new genes, expression profile, and to investigate comparative, functional and evolutionary genomics of non-model organisms that have limited genomic information. In recent years, RNA-Seq has been used to study low temperature transcriptome of a few insects, such as Drosophila virilis, Cryptolaemus montrouzieri and Micrarchus nov. sp. 2 [9, 35, 63], which were manually treated with low temperature in laboratories except for Ericerus pela [62]. Insects are generally most cold tolerant at their overwintering stages [7]. Thus we collected the overwintering M. punctipennis adults by digging into the sand about 5 cm beneath the earth surface in the field on January 11, 2016 and sent them for transcriptome sequencing. The control beetles were selected from those kept at 30 ℃ after collected from the field on October 3, 2015. In total, 11,355 differentially expressed genes (DEGs) were detected. Therefore， analysis of M. punctipennis transcriptome would help to study the mechanisms of the overwintering process of desert beetles.

65

2. Methods

66

2.1 Insect materials

67 68 69 70 71 72 73 74 75 76

M. punctipennis adults were collected from the southern edge of the Gurbantonggut Desert (N 44°24´, E 087°51´). The winter samples were collected on January 11, 2016, and were immediately frozen in liquid nitrogen. The soil at the collecting site was dry and not frozen at the time. The soil (5 cm beneath the earth surface) temperature was -8.8 ℃, the soil water content was 4.98%. Four of the winter individuals were sent for sequencing. The control samples were selected from those collected from the field on October 3, 2015, and then kept in incubator at 30 ℃ and 16:8 L: D photoperiod as described by Wang [58]. After about two weeks being kept in lab, the insects began to mate and the females laid eggs.

77

2.2 RNA extraction, library construction and RNA sequencing

78 79 80 81 82

Four insects (Mp_WF1, Mp_WF2, Mp_WM1 and Mp_WM2) were winter samples and three (Mp_30F1, Mp_30F2 and Mp_30M) were the control samples. Total RNA was extracted from each insect using TRIzol Reagent (Sangon Biotech, China). RNA concentration was measured by Qubit® RNA Assay Kit in Qubit® 2.0 Flurometer (Life Technologies, CA, USA), integrity was assessed by using RNA

AC C

EP

TE

D

M AN U

SC

RI PT

42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64

3

ACCEPTED MANUSCRIPT Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA). For each sample, 1.5 µg of RNA was used to generate sequencing library by using NEBNext® Ultra™ RNA Library Prep Kit for Illumina® (NEB, USA) following the manufacturer’s instructions. Library quality was assessed on AMPure XP system and Agilent Bioanalyzer 2100 system respectively. The index-coded sample clustering was performed with TruSeq PE Cluster Kit v3-cBot-HS (Illumia) on the cBot Cluster Generation System according to the manufacturer’s instructions. The library was sequenced on Illumina Hiseq platform and 150 bp paired-end reads were obtained. All samples were sent to Novogene Bioinformatics Technology (Beijing, China) for RAN sequencing.

94

2.3 De novo assembly and functional annotation

SC

RI PT

83 84 85 86 87 88 89 90 91 92 93

Raw reads were processed with in-house Perl script by Novogene Bioinformatics. Clean reads were acquired by removing the adapter and reads with ploy-N and low quality from raw reads. Phred quality score (Q20), GC content, and clean data sequence duplication level were calculated. The all left files (read1 files) from the sequenced libraries were pooled into one file (left.fq), and right files (read2 files) into another (right.fq). Transcriptome was assembled with left.fq and right.fq by Trinity (r20140413p1) [16] with all the other parameters were set to default and min_kmer_covfor to 2. Unigenes were searched against NCBI nr database (ftp://ncbi.nih.gov), Swiss-Prot database (http://web.expasy.org/docs/swiss-prot_guideline.html), COG and KOG (ftp://ncbi.nih.gov/pub/COG/COG), KEGG, (http://www.genome.jp/kegg/) with the E-value cut-off of 1e-5. Genes were tentatively identified based on the best alignment against known sequences. Blast2GO software was used to predict the sequence function and assign GO terms (http://www.geneontology.org/). KAAS (KEGG Automatic Annotation Server) was used to predict metabolic pathways on KEGG database. Unigenes were searched against Pfam database for domain/family prediction by HMMER 3.0 with E-value of 0.01.

112

2.4 Differential expression profiling and enrichment

113 114 115 116 117 118 119

Clean reads were mapped back to the assembly, and readcount was obtained. Differential expression analysis was performed with DESeq package in R environment. Genes with an adjusted P-value (Padj) < 0.05 were considered as differentially expressed. For enrichment analysis, all differentially expressed genes (DEGs) were mapped to GO terms and KEGG pathways, and then the significantly enriched categories were screened.

120

2.5 Validation of the RNA-Seq results by qRT-PCR

121 122 123

Six up-regulated and seven down-regulated DEGs in the overwintering transcriptome were randomly selected for verification of gene expression levels by using quantitative real-time PCR. The primer sequences for qRT-PCRs of these DEGs

AC C

EP

TE

D

M AN U

95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111

4

ACCEPTED MANUSCRIPT

136

3 Results

137

3.1 Sequencing and assembly

138 139 140 141 142 143 144 145

In the seven cDNA libraries, four came from the biological replicates of the winter samples (Mp_W) and three from the biological replicates of the control samples (Mp_30). The obtained clean paired-end reads of each sample ranged from 39,272,986 (Mp_WF1) to 56,907,376 (Mp_30M) (Table 1). The total clean nucleotides sequenced from each sample exceeded 5.89 Gb. In addition, three libraries of RNA-seq data (Mp_4) from our previous study (unpublished work) were also added to the transcriptomic assembling process.

146 147

Table 1. Reads of each sample generated from Illumina sequencing.

SC

M AN U

D

TE EP

Table 1

After Trinity assembly, total of 260,407 transcripts from all libraries were generated. During the de novo assembly, Trinity software produces potential alternative spliced isoforms. The isoforms originated from the same locus were considered to have the same “chrysalis component”, “butterfly sub-component” and share some paths of de Bruijn graph [57]. Finally, 175,247 unigenes in total were obtained with an N50 of 1,084 bp (i.e., 50% of the transcriptome is in unigenes of this size or larger), and an average length of 645 bp. Of the 175,247 unigenes, 50,527 (28.48%) were > 500 bp, 23,913 (13.65%) were > 1,000 bp and 10,511 (6.00%) were > 2,000 bp (Table 2).

AC C

148 149 150 151 152 153 154 155 156 157 158 159 160

RI PT

127 128 129 130 131 132 133 134 135

are listed in Supplementary file 1. The specificity of the primer sets was examined by running the PCR products on agarose gel to ensure a single band amplification. About 1μg of total RNA was used for the synthesis of cDNA with M-MLV Reverse Transcriptase (Progema, USA). qRT-PCR was performed in 25 l system, containing 1 l of forward and reverse primers of each, SYBR Green Real-Time PCR Master Mixes (Invitrogen, USA) and 1 l of cDNA template on STEP ONE PLUS™ Real-Time PCR System (Applied Biosystems, USA). The cycling conditions were as follows: 95 ℃ for 10 min followed by 40 cycles of 95 ℃ for 15 s, 60 ℃ for 30 s, and 72 ℃ for 15 s. Each sample had three individuals as replicates. The expression level of each DEG was normalized to that of the elongation factor 1-alpha (EF-1a), an internal reference gene used in our previous studies [32, 38]. Relative gene expression level was calculated with the method of 2-ΔΔCt [29].

124 125 126

Table 2. Statistics of the assembled transcriptome.

Table 2 5

ACCEPTED MANUSCRIPT 161

3.2 Functional annotation of the unigenes

163 164 165 166 167 168 169 170 171 172 173

To predict and analyze the function of the assembled sequences, similarity searches were performed by sequence- and domain-based alignments against several public databases with the aforementioned method 2.3. 64,634 unigenes (36.88% of the total unigenes) were significantly aligned a sequence at least in one of the public databases. Among of them, 48,495 (27.67%) were found in NR (Table 3). Only 6,157 (3.51%) had annotations in all databases. To predict the potential domains inside the assembled unigenes, open reading frame (ORF) was searched for each sequence, all transcripts with predicted ORFs were used to search against Pfam database using Profile Hidden Markov Models. Finally, 39,250 unigenes (22.40%) were annotated with Pfam domain information (Table 3).

174 175 176

Table 3. Statistics of the functional annotations of the unigenes in public databases.

M AN U

SC

RI PT

162

Table 3

177

3.3 Differential expression analysis

179 180 181 182 183 184 185 186 187 188 189 190 191

Unigenes Assembly were set as a reference, bowtie2 package was used for mapping clean reads from each sample to the reference respectively. The readcount of each unigene in each sample was calculated using RSEM package based on the bowtie2 results. Technical biases were normalized by FPKM (number of fragments per kilobase of transcript per million mapped fragments) transition. To identify genes with different expression levels, DESeq package was used to calculate the expression level for each gene based on the normalized readcount. Genes with FPKM > 0.3 are considered as expressed. In total, 11,355 genes (6.48%) were detected as differentially expressed with the absolute value of log2 fold-change (log2FC) ≥ 1 and Padj ＜ 0.05. Among these DEGs, 3,367 genes (29.65% of the total DEGs) were induced during winter, while 7,988 genes (70.35%) showed down-regulation (Supplementary file 2). The distribution of genes’ expression fold-change are shown in Figure 1.

AC C

EP

TE

D

178

Figure 1 192 193 194 195 196

Figure 1. Volcano plot of DEGs between group Mp_W and Mp_30. The x-axis corresponds to log2 fold-change value, and the y-axis displays the mean expression value of log 10 (Padj). Red dots represent the up-regulated unigenes (Padj ＜ 0.05), 6

ACCEPTED MANUSCRIPT green dots represent the down-regulated unigenes (Padj ＜ 0.05).

199

3.4 GO enrichment analysis of DEGs

200 201 202 203 204 205 206 207 208 209 210 211 212 213

In total, 40,167unigenes had GO annotations, which accounts for 22.92% of the whole unigenes (Table 3). Among of them 1,176 (34.93% of the total up-regulated genes) were up-regulated and 5,380 (67.35% of the total down-regulated genes) were down-regulated. After DEGs were mapped to Gene Ontology database, gene numbers were calculated from each GO term. Hypergeometric test was used to identify significantly enriched GO terms in DEGs against to the genomic background. When Padj was set to < 1e-5, there were 32 terms of biological process (BP) category, one term of molecular function (MF) category and seven terms of cellular component (CC) category that enriched significantly in the up-regulated genes; 17 terms of BP, 28 terms of MF and two terms of CC enriched significantly in the down-regulated genes (Supplementary file 3). These enriched GO terms were visualized at different levels by using topGO package in R environment. The deepest (or more specific) level of these enriched GO terms in the up- or down-regulated genes were shown in table 4.

214 215 216

Table 4. Significantly enriched GO terms at the deepest level in the up- or down-regulated genes.

M AN U

SC

RI PT

197 198

EP

TE

Among the up-regulated genes, 16 terms of BP, five terms of MF and three terms of CC were the most specific GO terms that were enriched in the overwintering M. punctipennis transcriptome. Among of them, the top five specific terms in BP were "animal organ development" (with Padj value of 1.13e-12), "negative regulation of multicellular organismal process" (1.04e-09), "nervous system development" (1.43e-09), "positive regulation of multicellular organismal process" (1.12e-08) and "negative regulation of response to stimulus" (1.12e-08). The top three MF were "structural constituent of cuticle" (2.41e-05), "insulin-like growth factor binding" (0.001053) and "ferric iron binding" (0.001508). The top three CC were "basement membrane" (2.41e-06), "integrin complex" (7.76e-06) and "extracellular membrane-bounded organelle" (3.41e-05). Among the down-regulated genes, eight terms of BP, 11 terms of MF and two terms of CC were the most specific terms. Among of them, the top three specific BP were "carbohydrate metabolic process" (6.03e-12), "peptide biosynthetic process" (6.59e-11) and "translation" (1.39e-10). The top three specific MF were "serine-type endopeptidase activity" (1.54e-12), "structural constituent of ribosome" (1.88e-09) and "hydrolase activity, hydrolyzing O-glycosyl compounds" (1.03e-07). The top two CC were "ribosome" (1.83e-10) and "Holliday junction helicase complex" (0.02688).

AC C

217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236

D

Table 4

7

ACCEPTED MANUSCRIPT

3.5 KEGG enrichment analysis of DEGs

238 239 240 241 242 243 244 245 246

To identify the biochemical pathways and the signal transduction pathways that were significantly changed in the overwintering adults, we performed metabolic pathway enrichment analysis on DEGs. 3,367 up-regulated and 7,988 down-regulated unigenes were identified to be affected by winter. These DEGs were significantly enriched into 11 metabolic pathways (Padj ＜ 0.05), including carbohydrate metabolism, translation, membrane transport, immune system, xenobiotics biodegradation and metabolism, transport and catabolism, lipid metabolism and metabolism of cofactors and vitamins (Figure 2).

Figure 2

Figure 2. Top 20 enriched pathways by DEGs. Rich factor is the ratio of

deferentially expressed genes to all genes annotated to the same pathway; q value (ranges from 0 to1) is the Padj value obtained by multiple hypothesis tests. As q value approaches to 0, the enrichment of the pathway becomes more significant.

M AN U

248 249 250 251 252 253 254

SC

247

RI PT

237

3.6 Expression level verification

256 257 258 259 260 261

To confirm the expression levels of the identified DEGs in the overwintering M. punctipennis. Six up-regulated and seven down-regulated DEGs were randomly selected for qRT-PCR. The results showed that the expression patterns of these genes were in concordance with the expression patterns obtained from the RNA-seq data (Table 5) with a high Pearson’s correlation coefficient of 0.851.

262 263

Table 5. Relative expression levels of 13 DEGs by qRT-PCR

TE

EP

AC C

264

D

255

Table 5

265

4 Discussion

266 267 268 269 270 271 272 273

In this study, changes of gene expression at the genomic level in the overwintering desert beetle M. punctipennis in contrast to that at 30 ℃ were examined by using RNA-seq technology. In total, 175,247 unigenes were finally obtained with an N50 of 1,084 bp in all samples. Only 64,634 (36.88% of the total unigenes) were annotated, suggesting that many genes are function-unknown in M. punctipennis. In total, 11,355 (6.48%) were detected as differentially expressed. Among these DEGs, 3,367 (29.65% of the total DEGs) were up-regulated in winter, and 7,988 (70.35%) showed down-regulation. Some DEGs were confirmed by qRT-PCR with a Pearson’s 8

ACCEPTED MANUSCRIPT correlation coefficient of 0.851. 70.35% of the DEGs were down-regulated when the insect was in an inactive overwintering state. This is similar to the transcriptomic analysis of the overwintering Chinese white wax scale insect, E. pela [62]. However, 29.65% of the DEGs were up-regulated during winter, indicating some physiological activities may remain active or get started [48, 56]. M. punctipennisadults can go through winter by decreasing its supercooling point (SCP) to lower than -19.6 ℃ via synthesizing antifreeze proteins and increasing haemolymph osmolality [17]. The genome-wide transcriptional dynamics of the overwintering M. punctipennis will help us to clarify the potential transcriptional regulation of the processes that may be relevant to overwintering.

284

4.1 Winter-associated genes

285 286 287 288 289 290 291 292 293 294 295 296

By blasting DEGs from M. punctipennis transcriptomic data with those in the overwintering E. pela, we found 40 up-regulated homologous transcripts (Table 6). Most of them belong to heat shock protein 70 (Hsp70), cuticle protein, myosin, tubulin, tropomyosin, actin, serine protease inhibitor, 40S ribosomal protein, ATP synthase, thioredoxin, elongation factor, enolase, annexin, glyceraldehyde-3phosphate dehydrogenase (GADPH), and antifreeze proteins, etc. GADPH and antifreeze proteins affect super cooling, the other genes may play roles in maintaining dormant state and the related physiological processes. These genes were up-regulated in M. punctipennis at the transcription level and abundantly translated into proteins in E. pela during winter. These similar gene expression pattern in both species may represent a ‘core set’ of genes which might contribute to the successful overwintering in these species.

297 298 299

Table 6. Winter-associated genes in M. punctipennis transcriptome.

SC

M AN U

D

TE

Table 6

EP

300

RI PT

274 275 276 277 278 279 280 281 282 283

4.2 Winter-associated GO terms

302 303 304 305 306 307 308 309 310 311 312 313

In molecular function category, structural constituent of cuticle was the most enriched term in the up-regulated genes. This is similar to the low temperature response of the New Zealand stick insect and the transcriptome analysis of the cold acclimated D. melanogaster [10, 30], indicating that fortification of cuticle structure integrity might be the common response of insects to low temperature. The up-regulated genes were also enriched in hormone activity, insulin-like growth factor binding and peptidase inhibitor activity. This is consistent with the enriched development-related biological processes, such as cell proliferation, cell growth, tissue development, anatomical structure morphogenesis, animal organ development, nervous system development, and cell migration, etc. In transcriptomes of the cold acclimated D. melanogaster and in the winter morph of Drosophila suzukii the up-regulated genes were also significantly enriched in terms related to cell and

AC C

301

9

ACCEPTED MANUSCRIPT

EP

TE

D

M AN U

SC

RI PT

tissue development [30, 45]. We observed in early spring (March 2, 2010) in the desert, when snow just began to melt, adult M. punctipennis couples had vigorously mated. So the up-regulation of the development related genes in the overwintering M. punctipennis suggested that the adults take advantage of winter time to prepare for reproduction in the following spring. The control adults were collected in late autumn, at that time they were ready for overwintering. After kept in laboratory at 30 ℃ for about two weeks, the females began to lay eggs. Before sequencing they had been kept in laboratory for three months. The analysis indicated that several spermatogenesis-associated genes were highly up-regulated (Supplementary file2) in them. In addition, a transcript of protamine was highly expressed only in control samples. These results were in consistent with the fact that the control insects were reproductive. In D. melanogaster, protamine expressed in late spermatogenesis and replaced the histones to condense chromatins [2]. In the up-regulated genes in winter samples, negative regulation of signaling, of response to stimulus and of cell communication was another group of enriched terms. This was consistent with the dormant state of the overwintering insects. "Response to organic substance" and "ferric iron binding" were also the significantly enriched terms in winter M. punctipennis. In the hypoxia-ischemia treated rat brain tissue, "response to organic substance" is the most enriched term [61], so this may hint that the overwintering dormant M. punctipennis is subject to hypoxic environment in sand during winter. Meanwhile the enrichment of "ferric iron binding" may emphasize this situation because iron has potential toxicity in dormant organisms which cannot actively repair molecules damage caused by reactive oxygen species (ROS) [51]. Cells deal with this problem by producing proteins with high affinity for ferric or ferrous ions [27], such as ferritin and transferrin. Ferritin is an anti-oxidant protein [36], while transferrin involves in insect innate immunity [27]. We found that transcripts of ferritin heavy chain and serotransferrin precursor were up-regulated in the overwintering samples (Supplementary file 2). Thus, ferric iron binding may involve in anti-oxidation and innate immunity in M. punctipennis during winter. Down regulation of genes in winter samples was prominent, which accounted for 70.35% of the DEGs. At cellular component category, genes were significantly enriched in "ribosome" and "holliday junction helicase complex", indicating that cellular components that participate in translation, genetic recombination and DNA repair were down-regulated during winter. Correspondingly, in biological process category, genes participating in translation, ribosome biogenesis and tRNA metabolic process were also down-regulated. In molecular function category, structural constituent of ribosome was also correspondingly down-regulated, which is similar to the result of Drosophila suzukii transcriptome analysis [45]. For overwintering insect, down regulating these energy consuming processes is economic, but may not directly contribute to low temperature protection. Another group of enriched GO terms in the down-regulated genes were related to processes associated with carbohydrate metabolism, peptide biosynthesis, proteolysis, etc. This may indicate that the winter insects alter their energy use-temperature

AC C

314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357

10

ACCEPTED MANUSCRIPT relationship by modifying thermal sensitivity and repressing metabolic rate [47, 50], since insects generally use lipids and carbohydrates, and some utilize storage proteins, to fuel metabolism [5]. A prominent group of GO terms in the down regulated genes were related to metabolisms, such as oxidoreductase activity, hydrolase activity, and carbohydrate derivative binding, etc. Down regulation of metabolism occurs in a wide range of organisms during being stressed [50]. Thus, down-regulation of oxidoreductase activity might be one of the results of the suppressed metabolism in the overwintering M. punctipennis. In D. drosophila that was cold-acclimated for five days at 6 ℃, the same GO terms were also enriched in the down-regulated genes [30]. These obvious overlaps between M. punctipennis and D. melanogaster suggests that some candidate mechanisms of cold response are not only conserved among the Drosophila phylogeny [30] but also extended to far related insect species. It is worth mentioning that the transcripts of some proteins may not be present in mid-winter, although the proteins themselves may be. Some proteins may have synthesized in autumn prior to the onset of winter. In the down-regulated genes, more than 4000 sequences have no or very little expression in winter samples. Some of them were significantly enriched in GO terms like ‘translation’, ‘ribosome biogenesis’ and ‘structural constituent of ribosome’, etc. Besides, among the unexpressed genes, there were also significant part of transcripts annotated as “unnamed protein product”, “uncharacterized protein”, and “hypothetical protein”, etc.

379

4.3 Winter-associated pathways

380 381 382 383 384 385 386 387 388 389 390 391 392 393

KEGG enrichment analysis revealed that 11 pathways were significantly enriched in the DEGs (Figure 2), and among of them six were overlapped to those in the cold-acclimated D. melanogaster transcriptome, such as amino sugar and nucleotide sugar metabolism, lysosome, and metabolism of xenobiotics by cytochrome P450, etc. [30]. Most of the DEGs enriched into these pathways in winter M. punctipennis were down-regulated. In some pathways all the DEGs were down-regulated (Supplementary file 4). To obtain information regarding the mechanisms that make M. punctipennis overwintering, we focused on the 3367 up-regulated genes in the transcriptome data, and found that nine KEGG pathways significantly enriched in overwintering M. punctipennis (Table 7).

AC C

EP

TE

D

M AN U

SC

RI PT

358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378

Table 7. KEGG pathways enriched in the up-regulated genes of the overwintering M. punctipennis transcriptome. Table 7

394 395 396 397 398

The extracellular matrix (ECM) receptor interaction (ko04512) was the most significantly enriched pathway in the up-regulated genes in this study. 18 up-regulated and two down-regulated genes were involved in this pathway. ECM consists of a complex structural and functional macromolecules, they provide structural support for 11

ACCEPTED MANUSCRIPT

EP

TE

D

M AN U

SC

RI PT

organs and tissues, for cell layers in the form of basement membranes [19]. In our study, basement membrane was also a significantly enriched GO term. Thus, this pathway may play an important role in overwintering M. punctipennis. Specific interactions between ECM and cells control cellular activities such as apoptosis, differentiation, migration, proliferation and adhesion [26, 34]. Along with ECM-receptor interaction pathway, focal adhesion (ko04510) pathway was significantly enriched. This pathway is essential for maintaining cellular physiology [55]. These pathways constitute a number of components that mediate cell signal transduction in cell survival, migration and proliferation [14, 33]. Some components of the focal adhesions are signaling molecules that function in reconstruction of the actin cytoskeleton [6, 13, 44]. Consistently, regulation of actin cytoskeleton (ko04810) pathway was also significantly enriched. In the proteomic data of the overwintering E. pela, this pathway was also detected [62]. Two actin genes were up-regulated in the midgut of Culex pipiens in response to low temperature and diapause [24]. In our study, 30 DEGs were related to this pathway. Among of them, half were up-regulated (Table 6). Hence the specific modulation of this pathway during winter might be complex. The mitogen-activated protein kinase (MAPK) signaling pathway (ko04010) and phosphoinositide 3-kinase/protein kinase(3K-Akt) signaling pathway (ko04151) were significantly enriched. Signal transduction pathways also play roles in responses to low temperature stress [21, 57]. Our previous study showed that Ras and JNK genes in MAPK pathway were up-regulated in M. punctipennis treated at 4℃ and -4℃ [38]. MAPK families play important roles in insect rapid cold hardening (RCH) [11, 22] and in the regulation of initiation and termination of temperature dependent insect diapause [12, 23]. From overwintering E. pela transcriptomic and proteomic data, many transcripts and proteins are also related to MAPK and 3K-Akt signaling pathways [62].Therefore, MAPK signaling path way may play roles in most low temperature related processes such as rapid cold hardening, cold acclimation, adult insect overwintering, as well as diapauses. Interestingly, estrogen signaling pathway (ko04915) was also enriched. In Drosophila estrogen-related receptor (ERR) defines carbohydrate metabolism [52] to promote the synthesis of amino acids, nucleotides and lipids, thereby to support the cellular proliferation [52]. However, in the winter M. punctipennis transcriptome carbohydrate metabolic process and peptide biosynthetic process as well as ERR transcript, were all down-regulated. These conflict results suggested that the role of estrogen signaling pathway in the overwintering M. punctipennis needs to be further studied. Tight junction (ko04530) pathway was significantly enriched. This pathway maintains the integrity of the epithelial cell layers by intercellular junctional complexes such as tight junction (TJ) [43]. TJ controls epithelial proliferation and differentiation [15, 31]. One member of the TJ protein families is claudin [43]. Claudin is one of the most up-regulated genes in the winter M. punctipennis transcriptomes (log2FC= 6.03). Claudin controls paracellular permeability of water and ions in knockout mice [53, 54]. Increase in modest temperature and osmotic stress

AC C

399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442

12

ACCEPTED MANUSCRIPT can increase paracellular permeability [8, 42]. Therefore, in the overwintering M. punctipennis tight junction pathway may mediate the proliferation of epithelial cells to form a barrier to protect the insect from low temperature stress.

446

5 Conclusion

447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467

In this study, the genomic level of gene expression changes in the overwintering desert beetle M. punctipennis was obtained by using RNA-seq technology. A total of 175,247 unigenes were assembled with an average length of 645 bp. 64,634 unigenes (about 36.88 %) matched to the known proteins. The remaining 66.12% were functionally unknown genes. 11,355 unigenes were differently expressed in the overwintering beetle with 3,367 (29.65%) up-regulated and 7,988 (70.35%) down-regulated. In view of the enrichment analysis, we primarily outlined the picture of gene expression changes in the overwintering dormant M. punctipennis adults. There were up-regulated genes involving in negative regulation of signaling, of response to stimulus and of cell communication; meanwhile there were down-regulated genes involving in energy consuming processes such as translation, genetic recombination and DNA repair. There were also down regulated genes related to metabolisms, such as oxidoreductase activity, hydrolase activity, etc. To maintain the basal energy supply during winter, M. punctipennis adults altered their energy use by reducing carbohydrate metabolism, peptide biosynthesis and proteolysis. In addition, they up-regulated genes related to fortification of cuticle structure integrity to protect their body from damage. They adopted anti-oxidation mechanism during winter to deal with oxidative stress caused by ferric iron. Importantly, during winter they were actively preparing for breeding in the coming spring by up-regulating the expression of many developmental related genes.

468

Data archiving

469 470 471 472

The transcriptomic sequencing data supporting the results of this paper have been submitted to the National Center for Biotechnology Information (NCBI). The data sets can be accessed in the Short Read Archive (SRA) with the accession number SRP090908.

473

Supplementary files

474 475 476 477 478 479 480 481 482 483

Supplementary file1: Primer sequences used for qRT-PCR. Sequence ID and their primers used in quantitative real-time PCR analysis were listed. Supplementary file2: Differentially expressed genes. List of genes that differentially expressed during winter. The gene ID, corresponding readcount, expression fold-change, as well as functional annotation of Nr databases were listed. Supplementary file3: GO enrichment analysis. GO terms enriched in up- and down-regulated DEGs with a corrected Padj ≤ 1e-5 were listed. Supplementary file4: KEGG pathway enrichment analysis. KEGG pathways enriched in DEGs with a corrected Padj ≤ 0.05 were listed. Supplementary file5: Extracellular matrix（ECM）receptor interaction. The most

AC C

EP

TE

D

M AN U

SC

RI PT

443 444 445

13

ACCEPTED MANUSCRIPT enriched KEGG pathway in up-regulated DEGs. Red boxes indicate the up-regulated genes and blue boxes indicate the down-regulated genes.

486

Conflict of interests

487 488

There is no conflict of interests in this study. All authors have agreed to submit the manuscript to Cryobiology Journal.

489

Acknowledgements

490 491 492 493 494 495 496

The research was funded by the National Natural Science Foundation of China (No. 31360527), and grant from Xinjiang Key Laboratory of Biological Resources and Genetic Engineering (No. XJDX0201-2014-03). We would like to acknowledge our fellow scholars, Fengjuan Zhang and Mengge Ruan for help on collection and identification of the insects used in this research. We also thank Qinglin Cai and Yaohua Li for their help during the verification of DEGs by real-time PCR.

497

References

498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524

[1] B. Abudureheman, H.L. Liu, D.Y. Zhang, K.Y. Guan, Y.K. Zhang, The responses of the quantitative characteristics of a ramet population of the ephemeroid rhizomatous sedge Carex physodes to the moisture content of the soil in various locations on sand dunes, Sci. World J. (2014) 120186. [2] Z.A. Alvi, T. Chu, V. Schawaroch, A.V. Klaus, Genomic and expression analysis of transition proteins in Drosophila, Spermatogenesis 5 (2015) e1178518. [3] J.L. Andersen, T. Manenti, J.G. Sørensen, H.A. MacMillan, V. Loeschcke, J. Overgaard, How to assess Drosophila cold tolerance: chill coma temperature and lower lethal temperature are the best predictors of cold distribution limits, Funct. Ecol. 29 (2015) 55–65. [4] J.S. Bale, S.A.L. Hayward, Insect overwintering in a changing climate, J. Exp. Biol. 213 (2010) 980–994. [5] T. Burmester, Evolution and function of the insect hexamerins, Eur. J. Entomol. 96 (1999) 213–225. [6] K. Burridge, M. Chrzanowska-Wodnicka, Focal adhesions, contractility, and signaling, Annu. Rev. Cell Dev. Biol. 12 (1996) 463–519. [7] D.L. Denlinger, R.E. Lee, Low temperature biology of insects, Cambridge University Press, Cambridge, 2010. [8] K. Dokladny, P.L. Moseley, T.Y. Ma, Physiologically relevant increase in temperature causes an increase in intestinal epithelial tight junction permeability, Am. J. Physiol-Gastr. L 290 (2006) G204–G212. [9] L.T. Dunning, A.B. Dennis, D. Park, B.J. Sinclair, R.D. Newcomb, T.R. Buckley, Identification of cold-responsive genes in a New Zealand alpine stick insect using RNA-Seq, Comp. Biochem. Phys. D 8 (2013) 24–31. [10] L.T. Dunning, A.B. Dennis, B.J. Sinclair, R.D. Newcomb, T.R. Buckley, Divergent transcriptional responses to low temperature among populations of alpine and lowland species of New Zealand stick insects (Micrarchus), Mol. Ecol. 23 (2014)

AC C

EP

TE

D

M AN U

SC

RI PT

484 485

14

ACCEPTED MANUSCRIPT

EP

TE

D

M AN U

SC

RI PT

2712–2726. [11] Y. Fujiwara, D.L. Denlinger, p38 MAPK is a likely component of the signal transduction pathway triggering rapid cold hardening in the flesh fly Sarcophaga crassipalpis, J. Exp. Biol. 210 (2007) 3295–3300. [12] Y. Fujiwara, K. Shiomi, Distinct effects of different temperatures on diapause termination, yolk morphology and MAPK phosphorylation in the silkworm, Bombyx mori, J. Insect Physiol. 52 (2006) 1194–1201. [13] F.G. Giancotti, E. Ruoslahti, Integrin signaling, Science 285 (1999) 1028–1032. [14] A.P. Gilmore, L.H. Romer, Inhibition of focal adhesion kinase (FAK) signaling in focal adhesions decreases cell motility and proliferation, Mol. Biol. Cell 7 (1996) 1209–1224. [15] L. Gonzalez-Mariscal, S. Lechuga, E. Garay, Role of tight junctions in cell proliferation and cancer, Prog. Histochem. Cytochem. 42 (2007) 1–57. [16] M.G. Grabherr, B.J. Haas, M. Yassour, J.Z. Levin, D.A. Thompson, I. Amit, X. Adiconis, L. Fan, R. Raychowdhury, Q.D. Zeng, Z.H. Chen, E. Mauceli, N. Hacohen, A. Gnirke, N. Rhind, F. di Palma, B.W. Birren, C. Nusbaum, K. Lindblad-Toh, N. Friedman, A. Regev, Full-length transcriptome assembly from RNA-Seq data without a reference genome, Nat. Biotechnol. 29 (2011) 644–U130. [17] F. Hou, J. Ma, X. Liu, Y. Wang, X.N. Liu, F.C. Zhang, Seasonal changes in antifreeze protein gene transcription and water content of beetle Microdera punctipennis (Coleoptera: Tenebrionidae) from gurbantonggut desert in central Asia, Cryoletters 31 (2010) 359–370. [18] R.X. Huang, W. Wu, X.F. Mao, H.Y. Hu, Z.T. Fan, Y.J. Hou, X.P. Li, C.H. Du, H.G. Shao, X. Huang, T. Ouyang, The fauna of the desert insects of Xiniang and its formation and evolution, Xinjiang science and technology publishing house, Urumqi, 2005（in Chinese）. [19] R.O. Hynes, The extracellular matrix: not just pretty fibrils, Science 326 (2009) 1216–1219. [20] J.T. Irwin, V.A. Bennett, R.E. Lee, Diapause development in frozen larvae of the goldenrod gall fly, Eurosta solidaginis Fitch (Diptera : Tephritidae), J. Comp. Physiol. B 171 (2001) 181–188. [21] A. Janská, P. Maršík, S. Zelenková, J. Ovesná, Cold stress and acclimation - what is important for metabolic adjustment?, Plant Biol. 12 (2010) 395–405. [22] J. Kelty, A. Ancevski, Rapid Cold Hardening in Drosophila requires MAPK signaling, FASEB J. 24 (2010). [23] K. Kidokoro, K. Iwata, M. Takeda, Y. Fujiwara, Involvement of ERK/MAPK in regulation of diapause intensity in the false melon beetle, Atrachya menetriesi, J. Insect Physiol. 52 (2006) 1189–1193. [24] M. Kim, R.M. Robich, J.P. Rinehart, D.L. Denlinger, Upregulation of two actin genes and redistribution of actin during diapause and cold stress in the northern house mosquito, Culex pipiens, J. Insect Physiol. 52 (2006) 1226–1233. [25] B. Krasnov, Y. Ayal, Seasonal changes in darkling beetle communities (Coleoptera: Tenebrionidae) in the Ramon erosion cirque Negev Highlands, Israel, J. Arid Environ. 31 (1995) 335–347.

AC C

525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568

15

ACCEPTED MANUSCRIPT

EP

TE

D

M AN U

SC

RI PT

[26] R.P. Lanza, R. Langer, W.L. Chick, Principles of tissue engineering, Academic press, Austin, 1997. [27] J.H. Law, Insects, oxygen, and iron, Biochem. Biophys. Res. Commun. 292 (2002) 1191–1195. [28] J.Q. Li, X.Y. Lu, J. Ma, Characterization and expression analysis of attacins, antimicrobial peptide-encoding genes, from the desert beetle Microdera punctipennis in response to low temperatures, Cryoletters 38 (2017) 65–74. [29] K.J. Livak, T.D. Schmittgen, Analysis of relative gene expression data using real-time quantitative PCR and the 2−ΔΔCt method, Methods 25 (2001) 402–408. [30] H.A. MacMillan, J.M. Knee, A.B. Dennis, H. Udaka, K.E. Marshall, T.J.S. Merritt, B.J. Sinclair, Cold acclimation wholly reorganizes the Drosophila melanogaster transcriptome and metabolome, Sci. Rep. 6 (2016) 28999. [31] K. Matter, S. Aijaz, A. Tsapara, M.S. Balda, Mammalian tight junctions in the regulation of epithelial differentiation and proliferation, Curr. Opin. Cell Biol. 17 (2005) 453–458. [32] S.S. Meng, J.Q. Li, M.G. Ruan, X.Y. Lu, J. Ma, Cloning and expression analysis of extracellular signal - regulated kinase ERK gene from the desert beetle Microdera punctipennis, J. Environ. Entomol. 37 ( 2015) 987–995（in Chinese）. [33] J.E. Meredith, B. Fazeli, M.A. Schwartz, The extracellular matrix as a cell survival factor, Mol. Biol. Cell 4 (1993) 953–961. [34] C.M. Nelson, M.J. Bissell, Of extracellular matrix, scaffolds, and signaling: Tissue architecture regulates development, homeostasis, and cancer, Annu. Rev. Cell Dev. Biol. 22 (2006) 287–309. [35] D.J. Parker, L. Vesala, M.G. Ritchie, A. Laiho, A. Hoikkala, M. Kankare, How consistent are the transcriptome changes associated with cold acclimation in two species of the Drosophila virilis group?, Heredity 115 (2015) 13–21. [36] D.Q.D. Pham, J.J. Winzerling, Insect ferritins: typical or atypical?, Bba-Gen. Subjects 1800 (2010) 824–833. [37] G.D. Ren, Y.Z. Yu, The desert and semi-desert Tenebrionids of China, Hebei university publishing house, Baoding, 1999（in Chinese）. [38] M.G. Ruan, J.Q. LI, S.S. Meng, J. Ma, Cloning and expression profiling in response to low temperature of Ras GTPase-activating protein gene MpRasGAP in the desert beetle Microdera punctipennis (Coleoptera: Tenebrionidae), Acta Entomol. Sin. 58 (2015) 367–374（in Chinese）. [39] L. Sømme, The effect of prolonged exposures at low temperatures in insects, Cryoletters 17 (1996) 341–346. [40] L. Sømme, Invertebrates in hot and cold arid environments, Springer-Verlag, Berlin, 1995. [41] L. Sømme, K.E. Zachariassen, Adaptations to low temperature in high altitude insects from Mount Kenya, Ecol. Entomol. 6 (1981) 199–204. [42] G. Samak, T. Suzuki, A. Bhargava, R.K. Rao, c-Jun NH2-terminal kinase-2 mediates osmotic stress-induced tight junction disruption in the intestinal epithelium, Am. J. Physiol-Gastr. L 299 (2010) G572–G584. [43] E.E. Schneeberger, R.D. Lynch, The tight junction: a multifunctional complex,

AC C

569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612

16

ACCEPTED MANUSCRIPT

EP

TE

D

M AN U

SC

RI PT

Am. J. Physiol-Cell Ph. 286 (2004) C1213–C1228. [44] M.A. Schwartz, V. Baron, Interactions between mitogenic stimuli, or, a thousand and one connections, Curr. Opin. Cell Biol. 11 (1999) 197–202. [45] P.W. Shearer, J.D. West, V.M. Walton, P.H. Brown, N. Svetec, J.C. Chiu, Seasonal cues induce phenotypic plasticity of Drosophila suzukii to enhance winter survival, BMC Ecol. 16 (2016) 11. [46] M. Shi, X.J. Liu, X.N. Liu, J. Ma, Effect of thermal stresses on hsp70 gene expression in larvae of the desert insect Microdera punctipennis ( Coleoptera: Tenebrionidae), J. Environ. Entomol. 36 (2014) 335–342（in Chinese）. [47] B.J. Sinclair, Linking energetics and overwintering in temperate insects, J. Therm. Biol. 54 (2015) 5–11. [48] K.B. Storey, J.G. Baust, J.M. Storey, Intermediary metabolism during low temperature acclimation in the overwintering gall fly larva, Eurosta solidaginis, J. Comp. Physiol. B 144 (1981) 183–190. [49] K.B. Storey, J.M. Storey, Freeze tolerance in animals, Physiol. Rev. 68 (1988) 27–84. [50] K.B. Storey, J.M. Storey, Metabolic rate depression and biochemical adaptation in anaerobiosis, hibernation and estivation, Q. Rev. Biol. 65 (1990) 145–174. [51] X.N. Tang, B. Zhou, Iron Homeostasis in Insects: Insights from Drosophila Studies, Iubmb Life 65 (2013) 863–872. [52] J.M. Tennessen, K.D. Baker, G. Lam, J. Evans, C.S. Thummel, The Drosophila estrogen-related receptor directs a metabolic switch that supports developmental growth, Cell Metab. 13 (2011) 139–148. [53] S. Tsukita, M. Furuse, M. Itoh, Multifunctional strands in tight junctions, Nat. Rev. Mol. Cell Bio. 2 (2001) 285–293. [54] C.M. Van Itallie, J.M. Anderson, Claudins and epithelial paracellular transport, Annu. Rev. Physiol. 68 (2006) 403–429. [55] S.J. Wala, J.R. Karamchandani, R. Saleeb, A. Evans, Q. Ding, R. Ibrahim, M. Jewett, M. Pasic, A. Finelli, K. Pace, E. Lianidou, G.M. Yousef, An integrated genomic analysis of papillary renal cell carcinoma type 1 uncovers the role of focal adhesion and extracellular matrix pathways, Mol. Oncol. 9 (2015) 1667–1677. [56] K.R. Walters, T. Sformo, B.M. Barnes, J.G. Duman, Freeze tolerance in an arctic Alaska stonefly, J. Exp. Biol. 212 (2009) 305–312. [57] X.C. Wang, Q.Y. Zhao, C.L. Ma, Z.H. Zhang, H.L. Cao, Y.M. Kong, C. Yue, X.Y. Hao, L. Chen, J.Q. Ma, J.Q. Jin, X. Li, Y.J. Yang, Global transcriptome profiles of Camellia sinensis during cold acclimation, BMC Genomics 14 (2013) 415. [58] Y. Wang, X.N. Liu, J. Zhao, K. Rexili, J. Ma, The rearing and biology of the desert beetle, Microdera punctipennis, under laboratory conditions, J. Insect Sci. 11 (2011) 39. [59] Y. Wang, L.M. Qiu, W.J. Xie, W. Huang, F. Ye, F.C. Zhang, J. Ma, Cold tolerance of transgenic tobacco carrying gene encoding insect antifreeze protein, Acta Agron. Sin. 34 (2008) 397–402（in Chinese）. [60] Y. Wang, M. Shi, X. Hou, S. Meng, F. Zhang, J. Ma, Adaptation of the egg of the desert beetle, Microdera punctipennis (Coleoptera: Tenebrionidae), to arid

AC C

613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656

17

ACCEPTED MANUSCRIPT environment, J. Insect Sci. 14 (2014). [61] L.J. Yang, D.Q. Ma, H. Cui, Proteomic analysis of immature rat pups brain in response to hypoxia and ischemia challenge, Int. J. Clin. Exp. Pathol. 7 (2014) 4645–4660. [62] S.H. Yu, P. Yang, T. Sun, Q. Qi, X.Q. Wang, X.M. Chen, Y. Feng, B.W. Liu, Transcriptomic and proteomic analyses on the supercooling ability and mining of antifreeze proteins of the Chinese white wax scale insect, Insect Sci. 23 (2016) 430–437. [63] Y.H. Zhang, H.S. Wu, J.Q. Xie, R.X. Jiang, C.S. Deng, H. Pang, Transcriptome responses to heat- and cold-stress in ladybirds (Cryptolaemus montrouzieri Mulasnt) analyzed by deep-sequencing, Biol. Res. 48 (2015) 66.

RI PT

657 658 659 660 661 662 663 664 665 666 667 668

SC

669

AC C

EP

TE

D

M AN U

670

18

ACCEPTED MANUSCRIPT

Table 1. Reads for each sample generated from Illumina sequencing. Sample

Raw reads

Clean reads

Percentage (%)

Clean bases

1

Mp_WF1

41,251,374

39,272,986

95.86

5.89G

2

Mp_WF2

43,339,334

41,272,108

95.89

6.19G

3

Mp_WM1

48,901,642

46,512,750

95.89

6.98G

4

Mp_WM2

42,718,822

40,694,974

95.93

6.1G

5

Mp_30F1

56,482,288

54,839,246

96.27

8.23G

6

Mp_30F2

55,370,960

53,833,008

95.99

8.07G

7

Mp_30M

58,467,426

56,907,376

96.37

8.54G

8

Mp_4F1

53,113,762

51,603,206

95.99

7.74G

9

Mp_4F2

56,499,854

54,966,268

96.03

8.24G

10

Mp_4M

56,901,292

55,241,870

95.72

8.29G

SC

RI PT

Sample No.

M AN U

Sequenced reads with quality score higher than 20 (Q20) were considered as clean reads.

Table 2. Statistics of the assembled transcriptome. Parameters

Transcripts

Unigenes

260,407

175,247

250,659,776

113,006,338

201

201

Maximum transcript length (bp)

30,181

30,181

Average transcript length (bp)

963

645

2,060

1,084

111,517

50,527

66,270

23,913

33,278

10,511

Number Assembly length (bp)

Length > 500bp

EP

Length > 1000bp

TE

N50 (bp)

D

Minimum transcript length (bp)

AC C

Length > 2000bp

Table 3. Statistics of the functional annotations ofthe unigenes in public databases. Public Database NR NT Swiss-Prot GO KEGG KOG PFAM Annotated in all databases

Number of Unigenes 48,495 22,825 33,128 40,167 17,557 23,113 39,250 6,157

Percentage (%) 27.67 13.02 18.90 22.92 10.02 13.19 22.40 3.51

ACCEPTED MANUSCRIPT All annotated unigenes Total unigenes

64,634 175,247

36.88 100.00

GO_accession

Term_

Description

number

type

Padj

RI PT

Table 4. Significantly enriched GO terms at the deepest level in DEGs. Up

All genes

BP

2.41e−06

23

216

GO:0042127

regulation of cell proliferation

BP

1.07e−07

16

83

GO:0001558

regulation of cell growth

BP

7.76e−06

11

36

GO:0023057

negative regulation of signaling

BP

7.04e−06

14

1.04e−09

14

45

BP

1.12e−08

BP

4.73e−07

18

114

BP

1.12e−08

21

111

BP

1.28e−07

22

147

SC

regulation of localization

70

17

78

GO:0009888

tissue development

GO:0009653

M AN U

GO:0032879

anatomical structure morphogenesis

BP

7.71e−08

43

423

GO:0010033

response to organic substance

BP

7.76e−06

26

244

GO:0048513

animal organ development

BP

1.13e−12

31

190

GO:0007399

nervous system development

BP

1.43e−09

26

163

BP

7.09e−06

14

71

BP

7.68e−07

16

97

BP

2.28e−06

10

35

MF

2.41e−05

19

137

GO:0048585

GO:0010648 GO:0016477 GO:0045597

regulation of multicellular organismal development

negative regulation of response to stimulus

negative regulation of cell communication cell migration

positive regulation of cell differentiation

structural constituent of cuticle

AC C

GO:0042302

organismal process

D

GO:2000026

positive regulation of multicellular

TE

GO:0051240

organismal process

EP

GO:0051241

negative regulation of multicellular

BP

GO:0030414

peptidase inhibitor activity

MF

0.004637

19

209

GO:0005179

hormone activity

MF

0.002302

24

273

GO:0005520

insulin-like growth factor binding

MF

0.001053

7

24

GO:0008199

ferric iron binding

MF

0.001508

8

36

GO:0005604

basement membrane

CC

2.41e−06

8

15

CC

3.41e−05

7

14

CC

7.76e−06

8

20

Padj

Down

All genes

GO:0065010 GO:0008305

extracellular membrane-bounded organelle integrin complex

GO_accession

Description

number

Term_ type

GO:0005975

carbohydrate metabolic process

BP

6.03e−12

389

1,839

GO:0006508

proteolysis

BP

1.13e−09

432

2,059

ACCEPTED MANUSCRIPT GO:0043043 GO:1901071

peptide biosynthetic process glucosamine-containing compound metabolic process

BP

6.59e−11

380

2,080

BP

0.000226

45

139

translation

BP

1.39e−10

400

2,241

GO:0043604

amide biosynthetic process

BP

2.49e−10

415

2,287

GO:0042254

ribosome biogenesis

BP

4.42e−08

295

1,653

GO:0006399

tRNA metabolic process

BP

0.000177

185

861

GO:0016491

oxidoreductase activity

MF

1.11e−06

537

2,981

GO:0003735

structural constituent of ribosome

MF

1.88e−09

247

1,352

GO:0097367

carbohydrate derivative binding

MF

9.27e−07

757

4,356

GO:0050662

coenzyme binding

MF

3.98e−06

178

767

GO:0043168

anion binding

MF

3.25e−07

836

4,804

MF

1.03e−07

137

532

3.07e−06

709

4,073 4,080

GO:0035639

hydrolase activity, hydrolyzing O-glycosyl compounds purine ribonucleoside triphosphate binding

MF

SC

GO:0004553

RI PT

GO:0006412

purine ribonucleoside binding

MF

3.66e−06

709

GO:0004252

serine-type endopeptidase activity

MF

1.54e−12

156

509

GO:0032555

purine ribonucleotide binding

MF

3.98e−06

709

4,084

GO:0030554

adenyl nucleotide binding

MF

5.01e−06

638

3,621

GO:0005840

ribosome

CC

1.83e−10

315

1,699

GO:0009379

Holliday junction helicase complex

CC

0.02688

40

141

M AN U

GO:0032550

TE

D

BP: biological process, MF: molecular function, CC: cellular component. DEGs: differentially expressed genes.

Table 5. Relative expression levels of 13 DEGs by qRT-PCR

c87092_g1

C-Jun N-terminal kinase [Microdera dzhungarica] PREDICTED: uncharacterized protein CG1785 isoform X1 [Tribolium castaneum]

AC C

c82561_g1

c84500_g1

c71258_g1

Expression level

Annotation (BLASTX)

EP

Unigene ID

Hypothetical protein YQE_05244, partial [Dendroctonus ponderosae]

PREDICTED: 40S ribosomal protein SA [Tribolium castaneum]

qRT-PCR

RNA-seq

1.09±0.16

1.2893

2.27±0.13

1.9163

3.10±2.00

1.8233

1.70±0.05

1.2685

3.25±0.38

1.5578

1.41±0.08

1.1807

0.32±0.17

-1.3211

0.34±0.02

-1.4705

PREDICTED: MAP kinase-interacting

c84564_g1

serine/threonine-protein kinase 1 isoform X1 [Tribolium castaneum]

c74055_g2 c76306_g1 c38941_g1

PREDICTED: 60S ribosomal protein L28 [Tribolium castaneum] -NADH dehydrogenase subunit 1 (mitochondrion) [Asbolus verrucosus]

ACCEPTED MANUSCRIPT

c88637_g1 c87663_g3 c89072_g1

c88870_g1

PREDICTED: transketolase-like protein 2 isoform X2 [Tribolium castaneum] PREDICTED: NAD(P) transhydrogenase, mitochondrial [Tribolium castaneum] PREDICTED: calreticulin [Tribolium castaneum] PREDICTED: uncharacterized protein LOC657194 [Tribolium castaneum] PREDICTED: probable helicase with zinc finger domain [Tribolium castaneum]

0.27±0.07

-1.4282

0.27±0.07

-1.8445

0.20±0.03

-1.6451

0.76±0.08

-2.0663

RI PT

c85197_g1

0.71±0.07

-1.9769

SC

Table 6. Winter-associated genes in M. punctipennis transcriptome. NR Description

c33749_g1

40S ribosomal protein S14, partial

c29988_g1

ATP synthase F0 subunit 6

c70387_g2

beta actin

c153327_g1

Chain A, crystal structure of the third thioredoxin domain of

E-value

log2FC

Padj

3.8683E-103

2.72

0.01439

7.87352E-41

3.79

0.020726

0

3.21

0.00014885

3.28584E-53

6.01

0.032975

5.80358E-17

2.26

0.0079741

0

2.92

0.00056024

4.13545E-98

1.92

0.019336

c69164_g1

cuticle protein precursor

c73450_g1

elongation factor 2

c77694_g1

heat shock protein 68a

c81827_g1

heat shock protein 70

M AN U

Gene_id

2.0961E-142

6.05

2.404E-25

c87054_g2

heat shock protein 70

9.02682E-66

6.12

7.0961E-23

c72456_g1

myosin regulatory light polypeptide 9

4.7786E-108

3.87

0.0013527

c76037_g1

myosin-9

0

2.70

0.00083348

c90102_g1

PREDICTED: actin, aortic smooth muscle isoform X3

0

6.45

0.000049

c31338_g1

PREDICTED: alpha-enolase isoform X1

3.3662E-117

4.51

0.0017192

c44888_g1

PREDICTED: annexin A5

0

3.36

0.0040545

c67667_g1

PREDICTED: cuticle protein 19-like

4.50552E-32

2.12

0.03756

c70059_g1

PREDICTED: cuticle protein 7

2.70013E-29

2.65

0.0085972

c71607_g1

PREDICTED: cuticle protein 7

3.82851E-48

2.28

0.02584

c73704_g1

PREDICTED:

0

2.73

0.0073608

AC C

EP

TE

D

human erp46

glyceraldehyde-3-phosphate

dehydrogenase

isoform X1

c83051_g1

PREDICTED: heat shock 70 kDa protein

0

2.09

0.0005392

c68678_g1

PREDICTED: heat shock 70 kDa protein cognate 2

5.0366E-149

1.67

0.010736

c87496_g1

PREDICTED: heat shock 70 kDa protein cognate 5

0

2.04

0.0012092

c77710_g3

PREDICTED: heat shock cognate 71 kDa protein

7.9821E-175

3.66

0.027664

c72170_g1

PREDICTED: heat shock protein 70 A1

1.8125E-96

5.85

9.045E-30

c69698_g1

PREDICTED: larval cuticle protein A2B

7.69993E-39

2.03

0.010626

c88445_g1

PREDICTED: myosin-11

4.0646E-158

1.33

0.0373

c82532_g1

PREDICTED: myosin-3

2.7958E-154

1.81

0.010782

c65412_g1

PREDICTED: plasminogen activator inhibitor 1 isoform X2

0

4.95

4.3101E-10

ACCEPTED MANUSCRIPT PREDICTED: protein takeout

9.16996E-42

1.75

0.0012581

c81170_g3

PREDICTED: protein takeout

2.5228E-110

2.97

0.000053518

c52429_g1

PREDICTED: tropomyosin alpha-1 chain isoform X7

3.4083E-123

3.08

0.0067125

c65396_g1

PREDICTED: tropomyosin alpha-4 chain isoform X1

1.3966E-112

3.43

0.0024483

c51699_g1

PREDICTED: tropomyosin beta chain isoform X1

3.7246E-100

5.88

6.1838E-06

c152169_g1

PREDICTED: tubulin beta-6 chain

2.2398E-142

6.15

0.026879

c63711_g1

PREDICTED: tubulin beta-7 chain-like, partial

0

3.97

0.0001029

c81753_g1

PREDICTED: unconventional myosin-Ixa

0

c72535_g1

Serine (or cysteine) peptidase inhibitor, clade A, member 1C

0

c51698_g1

Serine protease inhibitor A3K

0

c79869_g1

Tubulin alpha-1C chain

0

c50353_g1

Tubulin beta chain

1.36587E-93

c85405_g2

antifreeze protein isoform 1

4.34379E-56

RI PT

c79135_g1

0.0059969

4.88

0.015137

4.94

0.01702

2.42

0.016508

5.95

0.035048

3.30

8.7351E-06

M AN U

SC

1.85

Table 7. Enriched KEGG pathways in the up-regulated genes of the overwintering M. punctipennis transcriptome. Genes with pathway annotation

ID

Pathway

P-Value

Up

All

18

191

2.70E-05

25

426

0.00105

26

496

0.003662

ECM-receptor interaction

ko04510

Focal adhesion

ko04151

PI3K-Akt signaling pathway

ko04391

Hippo signaling pathway - fly

11

145

0.00458

ko04915

Estrogen signaling pathway

17

288

0.005862

ko04530

Tight junction

12

218

0.029205

ko04810

Regulation of actin cytoskeleton

15

299

0.032055

ko04390

Hippo signaling pathway

11

207

0.043746

16

341

0.044861

TE

EP

MAPK signaling pathway

AC C

ko04010

D

ko04512

CE ED

PT

SC

M AN U

RI

AC C EP TE D

SC

M AN U

RI PT