Population‐specific transcriptional differences associated with freeze tolerance in a terrestrial worm

Abstract Enchytraeus albidus is a terrestrial earthworm widespread along the coasts of northern Europe and the Arctic. This species tolerates freezing of body fluids and survives winters in a frozen state. Their acclimatory physiological mechanisms behind freeze tolerance involve increased fluidity of membrane lipids during cold exposure and accumulation of cryoprotectants (glucose) during the freezing process. Gene regulatory processes of these physiological responses have not been studied, partly because no gene expression tools were developed. The main aim of this study was to understand whether the freeze tolerance mechanisms have a transcriptomic basis in E. albidus. For that purpose, first the transcriptome of E. albidus was assembled with RNAseq data. Second, two strains from contrasting thermal environments (Germany and Greenland) were compared by mapping barcoded RNAseq data onto the assembled transcriptome. Both of these strains are freeze tolerant, but Greenland is extremely freeze tolerant. Results showed more plastic responses in the Greenland strain as well as higher constitutive expression of particular stress response genes. These altered transcriptional networks are associated with an adapted homeostasis coping with prolonged freezing conditions in Greenland animals. Previously identified physiological alterations in freeze‐tolerant strains of E. albidus are underpinned at the transcriptome level. These processes involve anion transport in the hemolymph, fatty acid metabolism, metabolism, and transport of cryoprotective sugars as well as protection against oxidative stress. Pathway analysis supported most of these processes, and identified additional differentially expressed pathways such as peroxisome and Toll‐like receptor signaling. We propose that the freeze‐tolerant phenotype is the consequence of genetic adaptation to cold stress and may have driven evolutionary divergence of the two strains.


| INTRODUCTION
The majority of Enchytraeids are soil and litter dwelling oligochaetes widespread in Arctic and temperate regions. In particular, in the colder regions, the soil layers in which enchytraeids are living freeze every winter and temperatures may drop far below the melting point of their body fluids increasing the risk of internal ice formation. It is, therefore, of interest to understand how animals like enchytraeids can survive winters in these cold areas. Cold hardiness in enchytraeids has been reported in several studies, and they survive sub-zero temperatures using either of the strategies freeze tolerance or cryoprotective dehydration (Kähler, 1970;Pedersen & Holmstrup, 2003;Slotsbo, Maraldo, Malmendal, Nielsen, & Holmstrup, 2008;Sømme & Birkemoe, 1997).
Freeze tolerance can be defined as the process where the animal tolerates internal, but only extracellular, freezing of body fluids, which induces cellular dehydration (Zachariassen, 1985). In contrast, cryoprotective dehydration is a strategy where the entire animal becomes dehydrated when exposed to sub-zero temperatures, which enables it to avoid freezing of body fluids (Holmstrup, Bayley, & Ramløv, 2002).
Enchytraeus albidus (Clitellata; Enchytraeidae) can be found in organic matter-rich soils (e.g., in compost) as well as in decaying seaweed at the supralittoral zone along the coasts of Western Europe, Scandinavia, and Greenland (Fisker, Overgaard, Sørensen, Slotsbo & Holmstrup, 2014). This species is freeze tolerant and has a wide physiological tolerance range to freezing temperatures in its habitat (Slotsbo et al., 2008;. In general, freeze tolerance of ectothermic animals is promoted by adaptive modifications of membrane lipid composition that reduce the risk of deleterious phase transitions of cellular membranes during cold exposure, that is, homeoviscous adaptation (Kostal, 2010), and the accumulation of low molecular weight cryoprotectants, that can reduce the ice fraction at sub-zero temperatures (Zachariassen, 1979). Clearly, a wide range of other molecular mechanisms is important for the ability to tolerate freezing of extracellular fluids in ectothermic organisms, but relatively little is known of this (Storey & Storey, 2013). Freeze-tolerant animals promote freezing at high sub-zero temperatures using ice nucleating proteins or inoculative freezing through the skin for controlled triggering of ice growth, and further may employ ice active (antifreeze) proteins that inhibit ice recrystallization. As cellular dehydration is an important consequence of freeze tolerance, aquaporins and other trans-membrane channels that allow movement of cryoprotectants are important as well as cytoskeleton remodeling preserving cellular integrity. Further, anti-oxidant defenses, heat shock proteins, and other chaperones are important for survival of freezing (for a review, see Storey & Storey, 2013).
In E. albidus, freeze tolerance involves active mobilization of glucose (from glycogen reserves) starting when the freezing process initiates and is therefore more akin to freeze tolerance of vertebrate ectotherms like frogs than of insects where cryoprotectant loading takes place months before the actual freezing of hemolymph starts Holmstrup, Costanzo, & Lee, 1999;Slotsbo et al., 2008). Strains of this species inhabiting different geographic locations have different freeze tolerance capacities. Arctic strains are-not surprisingly-clearly more freeze tolerant than strains of temperate regions and this is likely related to the capacity to adjust membrane phospholipid composition to low temperature, and the ability to accumulate large glycogen reserves during autumn.
Studies at the transcriptomic level using E. albidus were first available in a lower throughput fashion in 2011, for example, the first in house microarray was developed by Amorim et al. (2011) to study the effects of chemical and natural stressors. This was further developed and used to explore the ecotoxicological mechanisms of several contaminants (Gomes, Soares, Scott-Fordsmand, & Amorim, 2013;Novais, De Coen, & Amorim, 2012). Here, we attempt to obtain a more mechanistic understanding of freeze tolerance in E. albidus by studying differential transcriptomic responses between a freeze-tolerant strain and a freeze-sensitive strain. After de novo assembly and annotation of a consensus E. albidus transcriptome, two strains were compared using RNA sequencing; one originating from Greenland and one from Germany. The strain from Greenland tolerates much lower temperatures and endures freezing for much longer periods than worms from Germany Silva, Holmstrup, Kostal, & Amorim, 2013). Henceforth, we term the Greenland strain "freeze tolerant," while the Germany strain is termed "freeze sensitive." Our main aim was to identify key transcriptomic pathways linked to physiological alterations associated with freeze tolerance.

| Species
Enchytraeus albidus (Clitellata; Enchytraeidae) from Germany were obtained from a commercial supplier who originally collected the worms from garden compost (Büchner Zierfischfutter, Jena, Germany; coordinates: 51°51′N, 9°50′E). In our laboratory, the worms were cultured for 1 year in the laboratory in agricultural (loamy) soil at 5°C before being used in experiments. The worms were fed weekly with rolled oats mixed with dried and crushed macroalgae (predominantly Fucus spp., collected near Aarhus, Denmark). Worms from Greenland were collected in 2010 from decaying seaweed near the seashore at Kobbefjord about 20 km south-east of Nuuk (coordinates: 64°8′N, 51°23′W). These worms were also kept in the laboratory as described for Germany worms for about 1 year prior to experiments. Two experts in enchytraeid taxonomy verified (based on morphology and anatomy characters) that both strains indeed were E. albidus (B. Christensen, personal communication; R. Schmelz, personal communication).

| Exposures
To ensure that we had broad representation of stress-related genes in our reference transcriptome, we sequenced a pooled sample including mRNA from worms exposed to heat, cold, desiccation, pesticides, heavy metals, photoperiods, as well as worms at several developmental stages (adults, juveniles). Table S1 summarizes the exposures with relevant detailed exposure levels (Supporting Information). Worms were exposed following the standard guideline (OECD 2004) adapted to specific test treatment where needed. In short, ten adult worms with well-developed clitellum were introduced into glass vessels (50 ml) containing 25 g of test soil (moistened to 50% of the water holding capacity). Four replicates per condition were used. The tests were run at 20 ± 1°C with a 16-hr:8-hr light:dark photoperiod (except when testing the different temperatures and photoperiod treatments, see Table S1).

| Freezing experiment for RNA sequencing
Worms were either exposed to −5°C (hereafter termed "frozen") or kept at 2°C as control (hereafter termed "control"). Groups of ten adult enchytraeids were placed in 9-ml vials containing approximately 5-g fresh soil moistened with 0.22 ml tap water g −1 dry soil. The vials were closed with perforated lids (to allow ventilation) and kept at 2°C for 1 week prior to the freezing exposure. Freezing of the enchytraeids was accomplished by placing vials in a walk-in freezer set to −1.5 ± 0.2°C. After 1 hr, a small ice crystal was added to the soil surface to ensure a controlled freezing of soil and worms at a high sub-zero temperature (Slotsbo et al., 2008). After 16 hr, the soil and worms were frozen, and the vials were moved to a programmable cooling cabinet (Binder, Tütlingen, Germany) accurate to ± 0.2°C. The temperature was then gradually lowered to −5°C during 24 hr and kept at this temperature for 48 hr. With this freezing protocol, we simulated a harsh, but not unrealistic, freezing event in the natural habitat of Arctic E. albidus (Coulson et al., 1995). Frozen enchytraeids were sampled by quickly thawing the soil in cold tap water and blotting them dry with filter paper before they were placed in Eppendorf tubes. The tubes were snap frozen in liquid N 2 and kept at −80°C until RNA extraction. Unfrozen control worms received the same treatment in identical vials, but were kept at 2°C, before snap freezing in liquid nitrogen for RNA extraction. Treatments and controls were replicated three times for freeze-tolerant and freeze-sensitive strains.
It should be noted that both strains survive this freezing exposure . Each replicate contained ten pooled animals to obtain sufficient RNA for sequencing.

| RNA isolation, library prep, and sequencing
Total RNA was extracted using the TRIzol extraction method (Invitrogen, Belgium), followed by DNase treatment (Fermentas, Germany). The quantity and purity of the isolated RNA were measured spectrophotometrically on a NanoDrop ND-1000 Spectrophotometer (Thermo Scientific), and its quality was checked on a denaturing formaldehyde agarose gel. In the case of Illumina Miseq sequencing, aiming to provide the basis for a complete transcriptome assembly, a pool was prepared including two replicates of each condition (depicted in Table S1), each replicate contributed 1 μg of Total RNA to the pool. Subsequently, Evrogen (Moscow, Russia) applied their SMART approach (Zhu, Machleder, Chenchik, Li, & Siebert, 2001), followed by their DSN normalization method (Zhulidov et al., 2004) to obtain a normalized cDNA library. The TruSeq mRNA library prep kit was used to prepare a normalized library for sequencing on a MiSeq. The Miseq version 3 kit was used according to manufacturer's instructions to sequence the normalized library for 300-bp cycles paired end.
For the actual freezing experiment, aiming to obtain the freezerelated transcriptomics, three replicate samples per treatment were selected and prepared for sequencing using the Illumina TruSeq ® RNA Library Preparation Kit v2 according to the protocol supplied by the manufacturer. Briefly, in this procedure, total RNA is purified for mRNA using oligo-dT attached magnetic beads and mRNA is

| Assembly and data analysis
Sequence quality control on raw reads from both the MiSeq and the HiSeq runs was performed using FastQC. Raw read quality was high in general but the last 10 basepairs of the MiSeq reverse run were low in quality based on the calculated Phred score (inverse log of probability that the corresponding base call is incorrect). Consequently, low quality bases were trimmed using Trimmomatic with a quality cutoff set at 24 using Phred33 encoding according to Bolger, Lohse, & Usadel (2014). The transcriptome was assembled using Trinity RNA-Seq (Haas et al., 2013) using the default k-mer length settings. The two strains (Greenland and Germany) were assembled separately with each assembly containing the six strain-specific samples combined with the MiSeq data. Each set of assembled contigs was further clustered using Cap3 software (Huang & Madan, 1999) in order to merge and extend the original contigs generated by Trinity. We compared contigs from each strain by running best reciprocal BLASTn software locally to define a core set of transcripts, which was used as mapping reference for each individual strain (Camacho et al., 2009). This core set of transcripts was tested for completeness and assembly redundancy using Benchmarking Universal Single-Copy Orthologs (BUSCO) analysis against the metazoan database using default settings (Simão, Waterhouse, Ioannidis, Kriventseva, & Zdobnov, 2015). Subsequently, the core set was annotated using Blast2GO pro version 4.1 as follows. First, translated BLASTx was applied to identify homologs in the nonredundant NCBI sequence database, using expect-value cutoff of 10 −5 and at least 50% sequence similarity. Then, Gene ontology terms were mapped and annotated, including enzyme codes if available.
The mapping references for RNAseq expression analysis were generated using Bowtie2 (Langmead & Salzberg, 2012), while the actual mapping and read quantification was preformed with a combination of Bowtie2 and eXpress (Langmead & Salzberg, 2012;Roberts & Pachter, 2013). Statistical analysis of differential gene expression was performed in R using the EdgeR package (McCarthy, Chen, & Smyth, 2012), employing a GLM on a simple contrast in which the freezing treatment for each strain was directly compared to the control conditions. Moreover, the treatment x strain interaction was statistically tested. Subsequent p-values were corrected for multiple testing using the False Discovery Rate (FDR) method (Benjamini & Hochberg, 1995). This yielded two sets of significant genes that were subjected to Gene Ontology (GO) analysis using the TopGO package in R (Alexa, Rahnenfuhrer, & Lengauer, 2006). For GO analysis, we focused on the Biological Process and Molecular Function components as these two are the most informative. GO terms which only contained one transcript were removed from the significant GO term list. Additional KEGG ortholog annotation was performed using the KEGG Automated Annotation Server (http://www.genome.jp/tools/kaas/), and gene expression pathway analysis was performed using the R-packages GAGE in conjunction with Pathview using standard settings (Luo, Friedman, Shedden, Hankenson, & Woolf, 2009). We chose the four expression fold change ratios as determined by the gene expression interaction analysis from EdgeR as input for the KEGG pathway analysis in the GAGE package.

| RESULTS AND DISCUSSION
The assembly and annotation of a reference transcriptome based on both MiSeq data and HiSeq data are presented first. The datasets were combined in order to generate a comprehensive transcriptome covering as many transcripts as possible. Subsequently, raw reads from the RNA sequencing experiment (HiSeq data) were mapped against the reference transcriptome to study differential gene expression. Finally, differential expression networks were associated with physiological alterations, using annotation information as well as pathway analyses.

| Transcriptome assembly and annotation
For the assembly of the reference transcriptome, we made use of two Illumina datasets: The HiSeq run on the twelve barcoded treatment samples (i.e., two strains, frozen and control, three replicates) and a MiSeq run that contained a normalized pool of stress-response samples (Table S1). Initially, we pooled all the data for a joined assembly, resulting in ~270 million quality-controlled reads to be used in Trinity-RNAseq, yielding in total 292,933 contigs with a median length of 413 basepairs. However, mapping of the two individual strain reads showed that the assembly was highly strain specific. As this strainspecific assembly confounded gene expression interaction effects between strain and treatment, we decided to assemble the strains separately and use best reciprocal BLAST-hit (Moreno-Hagelsieb & Latimer, 2007) to define a core transcriptome. In this core transcriptome, genes have a specific sequence for each strain but are linked between strains as orthologs. The strain-specific assembly yielded 56,050 contigs (median length of 846 bp) for the Germany strain and 51,918 contigs (median length of 845 bp) for the Greenland strain.
Subsequent best reciprocal BLASTn analysis between the two strains yielded a final core transcriptome set of 40,002 genes with an average median length of 950 bp, which was used as mapping reference for the differential gene expression analysis (Table 1).
The two gene sets were annotated with a combination of BLASTx and Blast2GO and yielded 13,231 genes with a blast result (34%) for the Germany strain and 12,699 genes with a blast result (32%) for the Greenland strain (Table 1). About 32% of the genes which contained a BLAST-hit in the Germany strain were fully annotated and usable for Gene Ontology (GO) analysis. The transcriptome of the Greenland strain contained 31% of fully annotated genes, and KEGG pathway information could be retrieved for 25% of all assembled contigs in both strains ( Table 1). As the two sets were considered ortholog sets with highly similar annotation information, we merged the annotation results for GO analysis and pathway analysis to obtain an unbiased set of annotation terms per gene to be used for both strains. The combined annotation also yielded a total of 1,136 unique enzyme codes which could be mapped to metabolic pathways using IPATH to determine the completeness of metabolic pathways represented in the transcriptome (Letunic, Yamada, Kanehisa, & Bork, 2008). This is visualized as red edges in Fig. S1, suggesting that most of the metabolic pathways are represented except for metabolism of terpenoids and polyketides, which are usually more associated to microbes and plants. The carbohydrate and fatty acid metabolic pathways, which are involved in mechanisms of freeze tolerance, are well represented. Finally, transcriptome completeness was determined by applying BUSCO analysis, T A B L E 1 Summarized properties of assembly, annotation, and quality control of the core transcriptome used for mapping of the RNA sequencing data for Germany strain (G) and Greenland strain (GR) using a curated reference set of single-copy core genes present in all metazoans. This assessment shows that about 84% of this core set is identified in the current E. albidus transcriptome, indicating a high level of completeness (Table 1). We have to note, however, that about 9% of these core set genes are present as duplicates. This suggests that the assembly process generated some redundancy within the final set of contigs (Simão et al., 2015).

| Differential gene expression analysis
Each RNAseq library generated on average 20.3 million of paired-end raw sequence reads per sample, of which 97% passed quality control.
After mapping of the RNAseq reads of the Germany strain to the assembled contigs, we identified 3,569 transcripts to be differentially expressed between unfrozen/control conditions (+2°C) and the frozen treatment (−5°C). From these, 2,426 transcripts were significantly upregulated, while 1,143 were significantly down-regulated (Figure 1). In case of the Greenland strain, we identified only 940 transcripts to be significantly differentially expressed between unfrozen (controls) and frozen worms. Only 215 transcripts were significantly up-regulated, while 725 were significantly down-regulated ( Figure 1). Figure 2 shows the smear plots, visually illustrating the differences between the two strains, which already indicates that the transcriptome profile of the freeze-tolerant strain represents a phenotype that is much less perturbed due to freezing, when compared to the freeze-sensitive Germany strain. This lack of perturbed transcriptional profile has been observed regularly both in animals and plants (Roelofs, Aarts, Schat, & Van Straalen, 2008), and thus seems to represent a general pattern of tolerance. It reflects the ability of a tolerant phenotype to maintain homeostasis under circumstances that would be considered as stressful for nonadapted organisms. Several examples show that this is a consequence of the alteration to increased constitutive expression of otherwise stress-inducible genes coding for enzymes/peptides that protect the organism against environmental stress (Halimaa et al., 2014;Roelofs et al., 2009).
Functional annotation and GO term enrichment analysis of significantly different transcripts show that the up-regulated genes in the Germany strain can be linked to changes in lipid metabolism, while significantly down-regulated genes are mainly related to mRNA processing and splicing (Table S2). In contrast, significantly up-regulated genes in the Greenland strain are mainly associated with transport of sodium and drugs as well as more general processes related to transcriptional regulation. In several arthropods, sodium/calcium transporters are up-regulated in response to cold (Des Marteaux, Mckinnon, Udaka, Toxopeus, & Sinclair, 2017;Gao et al., 2009). In situ hybridization showed that increased sodium transport was localized at the nephridial canal, a region previously associated with ion reabsorption (Gao et al., 2009). Cold acclimation in trout is also associated with increased sodium pump activity in erythrocytes. The authors suggested that this compensates for the reduced trans-membrane gradients of sodium and potassium in erythrocytes during cold temperature regimes (Raynard & Cossins, 1991). In insects, however, cold exposure can result in chill coma and chilling injury, which is related to reduced active ion transport. This causes leaking of cellular K + and depolarization of the resting potential, which ultimately can result in influx of Ca ++ , apoptosis and necrosis (Overgaard & MacMillan, 2017).
In freeze-tolerant E. albidus, these physiological constraints are also a threat for survival, but there might be additional problems with pumping ions in a frozen state when cells are dehydrated and much of the extracellular fluids are changed to ice. We cannot say whether the differential expression of genes coding for the mentioned ion transporters we observed is directly related to cold injuries, but the results do suggest that differences in freeze tolerance, in addition to many other processes, may also somehow be associated with the expression of ion pumps. More targeted studies are needed to look into this.

| Population × freeze treatment interaction
By comparing the differential expression profiles of the freezetolerant Greenland strain with the freeze-sensitive Germany strain, we were able to identify genes and molecular processes that changed during freezing and were different in the two strains. These gene sets may have undergone transcriptional regulatory evolution in order to genetically adapt to more frequent and more severe freezing conditions in the Greenland strain. To that end, we statistically analyzed transcripts that show an interaction between genotype (Greenland versus Germany) and response to freezing (control versus frozen). In total, 1,374 transcripts showed such significant interaction (Table S3). Gene ontology enrichment analysis was performed on the annotated transcripts belonging to each of the three clusters. Transcripts of cluster 1 could be grouped into blood coagulation, chromatin condensation, and positive regulation of transcription. Transcripts belonging to cluster 2 were linked to lipid metabolism (long-chain fatty acid biosynthesis), oxidation reduction, and metabolism and transport of anions, lipids, and drugs. Unfortunately, the level of annotation in cluster 3 was too low to provide significantly enriched processes. However, a member of the HSP40/DnaJ family (DnaJB2 homolog) was identified, which is down-regulated in the Germany strain, slightly up-regulated in the Greenland strain, but shows a high constitutive expression level in the Greenland strain (transcript EAID 31697, Table S3). This family of molecular chaperones engages in a complex with Heat shock protein 70, thereby stimulating its adenosin triphosphate domain (ATPase) providing substrate specificity of Hsp70 (Ohtsuka & Hata, 2000). Interestingly, we identified a transcript in cluster 2, homologous to hsp70 (transcript EAID_14676,

| Linking differential expression to physiological alterations
When comparing gene expression of frozen and unfrozen control animals recall that the "frozen" treatment will potentially include mRNA transcribed both during the freezing process itself (at −1.5°C) and/or during the subsequent 72 hr where animals were kept frozen at −5°C.
The mRNA abundance therefore represents a response to the freezing process where ice formation and cell shrinking takes place, but apparently inducing transcriptional changes even though cells were dehydrated. Previous studies showed that ice fraction at −5°C is relatively low due to accumulation of high glucose concentrations allowing for a significant metabolic rate at this temperature, and therefore probably also allowing for a significant rate of differential transcription Patrício Silva, Enggrob, Slotsbo, Amorim, & Holmstrup, 2014). Still, we are careful with the interpretation of the observed differential transcription. For example, we can only assume that this mRNA has been translated into protein. Alternatively, the mRNA was regulated in anticipation of the processes associated with thawing, but our study does not allow for teasing apart these two different possibilities. oxidative stress defense mechanisms were previously suggested to be involved in freeze tolerance in this species. We will show results and discuss these processes in the context of molecular pathways identified by the interaction analysis between treatment and strain within the RNAseq data (Figure 3).

| Fatty acid metabolism
The GO enrichment analysis already indicated that fatty acid metabolism was significantly affected (  in several geographically different strains of E. albidus. They showed that the preponderance of shorter unsaturated PLFAs (e.g., linoleic acid) is associated with increased freeze tolerance of strains from arctic regions, a pattern that has also been observed in Drosphila species (Slotsbo et al., 2016). Here, we provide evidence that several fatty acid elongases are being suppressed in the freeze-tolerant Greenland strain upon freezing treatment, while these genes remain unchanged in the German strain. This suggests that the freeze-tolerant Greenland strain potentially is more flexible in controlling PFLA length over a wide temperature gradient. However, as we have no evidence that transcription of these genes actually results in protein synthesis responsible for PLFA modifications in the frozen state, our results cannot be directly used as molecular support to observations showing that freeze-tolerant strains of E. albidus have higher proportions of short chain PLFAs than freeze-sensitive strains (Fisker et al., 2015).

| (Poly)sugar metabolism and transport
We annotated four genes involved in glucose metabolism and transport in cluster 2 (Figure 3), among which are two sodium-dependent glucose transporters, glucose dehydrogenase, and glucose isomerase (involved in glycolysis). Their expression patterns are highly comparable to the expression patterns observed among genes involved in fatty acid metabolism. During freeze treatment these genes are slightly down-regulated in the Greenland strain, but the absolute levels of gene expression are mostly significantly higher when compared to the Germany strain.
In contrast, glycogenin-1, a protein involved in glycogen biosynthesis, is less expressed in the Greenland strain when compared to the Germany strain. Upon freezing treatment, this gene is strongly down-regulated in the Greenland strain, while expression levels in the Germany strain remain high.
Glucose is a well-known cryoprotectant in oligochaetes, because it is the primary blood sugar of these animals (Holmstrup et al., 1999).

Gene expression
Gene expression seems to be decreased as well. This is in line with previous observations from Slotsbo et al. (2008), that glucose levels are accumulating in the Greenland strain at the expense of glycogen levels to serve as a cryoprotectant. This conclusion is further underpinned by a severe down-regulation of glycogenin-1 in the tolerant strain (Figure 4b). This protein is essential in glycogen biosynthesis, indicating that glycogen biosynthesis is down-regulated to keep glucose levels high for optimal cryoprotection. It also shows again that Greenland worms respond much more flexible to temperature changes when compared to the Germany strain with regard to physiological processes (PFLA composition, glucose content) essential in coping with below zero temperatures. The level of plasticity in these responses may be an important genetic adaptation to withstand fluctuating freezing temperatures, which maybe a realistic scenario under natural conditions.

| Oxidative stress
Genes involved in oxidative stress were observed both within cluster 1 and cluster 2. Cluster 1 (Figure 3) contained gluthatione-Stransferases, thioredoxins, and glucuronosyltransferase. Overall, gene expression levels were always significantly lower in the Greenland strain when compared to the Germany strain, and the freezing treatment caused an increase in this differential expression pattern. Cluster 2 contained several ABC transporters and Cyp p450s (Table S3, Figure   5). These genes were all slightly down-regulated due to freezing in the Greenland strain, while being slightly up-regulated in the Germany strain. They suggest that the Greenland strain encounters less oxidative stress than the Germany strain, although ABC transporters and Cyp p450s may not necessarily always be associated with oxidative damage. Nevertheless, this effect becomes even more pronounced under freezing conditions.
Response to oxidative stress is a very general process associated with almost all types of stress (Tomanek, 2015).  showed a higher activity of superoxide dismutase as well as high levels of glutathione in freeze-tolerant Greenland worms, when compared to Germany strain. This indicates a higher capacity of tolerant phenotypes to withstand oxidative stress during freezing.
Here, we observed genes coding for thioredoxins and glutathione-S-transferases to be down-regulated in the Greenland strain, when compared to the Germany strain ( Figure 5). This may support the observation of  that these worms suffer much less from oxidative stress when frozen, but we did not investigate the abundance and activity of the enzymes and the actual balance of ROS production. Moreover, the amplitude of transcriptional variance in response to freezing is much higher in freeze-tolerant phenotypes when compared to freeze-sensitive animals ( Figure 5), suggesting that plasticity of transcriptional regulation of freeze tolerance-associated gene networks is an adaptive feature of freeze-tolerant phenotypes.
Interestingly, two metallothionein (mt) transcripts were significantly up-regulated only in the freeze-tolerant Greenland strain.
Metallothionein is mainly involved in metal detoxification, but can also be activated by oxidative stress (Haq, Mahoney, & Koropatnick, 2003;Ling et al., 2016). This was recently observed in the earthworm, Dendrobaena octaedra, in response to freezing (Fisker, Holmstrup, & Sorensen, 2016), and the authors proposed two scenarios to explain this result. Metal ion concentrations are expected to increase rapidly upon freeze-induced dehydration of the cell, which could trigger mt activation. Alternatively, mt could be induced by hypoxia as a result of reduced diffusion of oxygen into the worm body. Anoxic conditions lead to increased ROS, while two previous studies showed that mt can act as anti-oxidant agent (English & Storey, 2003;Viarengo et al., 1999). Thus, mt induction in the Greenland strain may be an indirect response to oxidative stress and/or increased metal ion concentrations at the cellular level.

| Pathway analysis
In total, 56 pathways showed a significant interaction between population and freeze treatment (Table S4). The majority of these (51)  Toll-like receptor pathway are significantly up-regulated in Greenland under freezing conditions. Several previous studies suggest that insects increase immunity after cold exposure (Sinclair, Ferguson, Salehipour-shirazi, & MacMillan, 2013). Our data suggest that innate immunity may play an important role in freeze tolerance in E. albidus.
Whether this is a nonadaptive by-products of a generalized stress response, or the result of a coadaptation linking pathogens and cold remains to be elucidated in the future. In summary, the pathway analysis suggests that the freeze-tolerant strain is able to maintain important metabolic functions, such as energy (ATP) production as well as innate immune responses even under freezing conditions.

| FUTURE DIRECTIONS
Our study is one of the first employing a global transcriptomics approach to the elucidation of freeze tolerance physiology in F I G U R E 6 Differential expression of genes representing components of the "Peroxisome" KEGG pathway, as an example of pathways that are significantly differentially expressed between freeze-tolerant Greenland strain and freeze-sensitive German strain, in this case under frozen conditions (−5°C invertebrates. We identified several major responses to freezing that deserve further investigation in more targeted studies. Recent studies suggest that tolerance of freeze-thaw events in enchytraeids has a significant genetic component ).
The present study suggests these genetic differences have a basis in the expression of molecular chaperones, and pathways involved in fatty acid metabolism, cell signaling, energy metabolism, innate immunity, and reduction of oxidative stress caused by freezing and possibly thawing. Such issues would be important when assessing the potential threats of climate change driving increased frequency of freezing and thawing events which may threaten soil biota in the Arctic and subarctic areas.

| CONCLUSION
Overall, we observed two patterns of gene expression regulation in freeze-tolerant animals. One pattern resembles constitutive overexpression of otherwise highly inducible stress response genes like Hsp40/DnaJ and Hsp70. Such patterns have been described in several previous studies dealing with genetic adaptation to stress tolerance. Secondly, the freeze-tolerant phenotype shows plastic responses in transcriptional regulation of genes involved in pathways that were earlier described to tolerate freezing conditions in physiological studies (PFLA composition, sugar metabolism, and oxidative stress). This flexibility is proposed to maintain homeostasis under freezing conditions, which is observed at the physiological level as well as the transcriptome level. This apparently results in much less overall transcriptional perturbation in the freeze-tolerant phenotype, while strong perturbation of the complete transcriptional machinery is observed in reference to E. albidus animals under realistic freezing conditions. We speculate that this may be the consequence of the type of stressor. Temperature fluctuations are of a seasonal and diurnal nature, requiring highly plastic phenotypic response. As such, increased plasticity in protective cascades may be much more effective to counteract the physiological challenges incurred by freezing.

DATA ACCESSIBILITY
All raw RNAseq data as well as all assemblies are submitted to NCBI under BioProject ID code PRJNA386732.