Metabolomic and transcriptomic responses of ticks during recovery from cold shock reveal mechanisms of survival

Ticks are hematophagous ectoparasites and are vectors for a myriad of disease-causing pathogens in vertebrate hosts throughout their range (Sonenshine and Roe, 2013). The American dog tick, Dermacentor variabilis, is one of the most widely distributed ticks in North America and is a prolific vector for Rocky Mountain spotted fever and tularemia (de la Fuente et al., 2008; Goethert and Telford, 2009). In recent years, the geographic distribution of D. variabilis has changed, including an expansion of this species northward (Dergousoff et al., 2013). This augmentation of geographic range is occurring in nearly all ixodid ticks (Daniel et al., 2003; Lindgren et al., 2000; Raghavan et al., 2019; Sonenshine, 2018), and is likely increasing the risk of tick-borne diseases as a result of altered human and animal interaction with ticks.

There are a multitude of factors that influence the distribution of ticks; however, low temperatures in winter are one factor that limits the latitudinal and altitudinal distribution of various tick species (Daniel et al., 2003; Dantas-Torres and Otranto, 2011; Dergousoff et al., 2013; Lindgren et al., 2000). Dermacentor variabilis overwinters in any of its life stages particularly at the northern limits of its range (Belozerov et al., 2002; Burg, 2001; Smith and Cole, 1941; Sonenshine, 1993). To enhance overwintering survival, ticks employ a variety of strategies including accumulation of cryoprotectants (Neelakanta et al., 2010; Yu et al., 2014), entering a diapause-like state (Yoder et al., 2016) and seeking sheltered hibernacula (Burks et al., 1996; Rosendale et al., 2016a). These approaches likely improve survival by reducing the direct and/or indirect effects of low temperature as D. variabilis, like most tick species, are chill susceptible and experience mortality above their super-cooling points (Burks et al., 1996; Dautel and Knülle, 1996; Needham et al., 1996; Rosendale et al., 2016a; Yu et al., 2014). Another potential strategy for survival of sub-zero temperatures is employing the rapid cold-hardening (RCH) response, which has been shown to occur in several tick species, including D. variabilis (Rosendale et al., 2016a; Wang et al., 2017; Yu et al., 2014).

During RCH, arthropods can dramatically improve their cold hardiness after a brief exposure to low temperatures (Lee et al., 1987). In contrast to long-term cold acclimation, which occurs over days to weeks, RCH occurs in minutes to hours and has been documented in a wide range of arthropods (Elnitsky and Lee, 2010). RCH is a critical component of ectotherm survival when daily fluctuating temperatures swing between above and below freezing. The mechanisms of RCH are incompletely understood; however, various physiological changes have been implicated in the process. RCH is regulated by several signaling events, including MAP kinase, apoptosis and calcium signaling (Teets and Denlinger, 2013). These signaling pathways, among others, contribute to the accumulation of various cryoprotectants, cause changes in the fatty acid composition of cell membranes, inhibit apoptotic pathways and cause changes in chaperone protein abundance (Teets and Denlinger, 2013). However, there is still much that is unknown about the underlying mechanisms of RCH, with even less known about this phenomenon in non-insect arthropods.

The role that transcriptional changes play in the RCH response is unclear because of discrepancies among species and study design (Teets et al., 2020). Although an upregulation of several genes has been observed during RCH (reviewed in Teets et al., 2020), most studies indicated that the RCH response occurs with little to no change in transcript expression (Sinclair et al., 2007; Teets et al., 2012; Vesala et al., 2011). However, the transcriptional profile of flies recovering from RCH and cold shock (CS) is distinct from that of untreated controls (Teets et al., 2012). Similarly, the accumulation of cryoprotectants during RCH exposure is minimal (Teets et al., 2020), whereas a recovery period following RCH and CS can dramatically impact the metabolic profile (Teets et al., 2012). These transcriptional and metabolomic adjustments during recovery from RCH and CS may contribute to cold injury repair and/or prepare the organism for subsequent cold exposure.

In ticks, there is a paucity of information on the physiological and molecular responses to low temperatures, including recovery from cold exposure. It is also unclear how important the RCH response is in tick overwintering. Dermacentor variabilis, Haemaphysalis longicornis and Dermacentor silvarum undergo RCH (Rosendale et al., 2016a; Wang et al., 2017; Yu et al., 2014), but Ixodes scapularis seems to lack this response (Vandyk et al., 1996). The bush tick, H. longicornis, responds to short-term cold acclimation through changes in water content and protein levels, but not glycerol; however, these alterations seem to be life-stage specific (Yu et al., 2014). The actual mechanisms of cold hardiness in ticks remain largely unknown; therefore, the purpose of this study was to examine the physiological and molecular responses of D. variabilis during recovery from cold and RCH exposure using a combined transcriptomic and metabolomic approach followed by targeted functional studies. Previous reports suggest that a recovery period is required to elicit transcriptional changes following RCH (reviewed in Teets et al., 2020); therefore, this study included a recovery period. With this experimental design, the hypothesis that gene expression during recovery from CS is modulated by RCH pretreatment can be tested.

Unfed Dermacentor variabilis (Say 1821) and Amblyomma maculatum Koch 1844 were acquired from the Oklahoma State University (OSU) Tick Rearing Facility (Stillwater, OK, USA). The laboratory colonies at OSU are kept under constant conditions [14 h:10 h, light:dark (L:D), 97% relative humidity (RH), and 25±1°C] year round. Upon arrival, groups of 10–15 ticks were transferred to 15 cm³ mesh-covered vials and placed in closed chambers containing a supersaturated solution of potassium nitrate, providing 93% RH (Winston and Bates, 1960). Ticks were kept in these chambers at 26±1°C and 15 h:9 h L:D until used in experiments, usually within 1 month. For treatments, ticks were randomly chosen from multiple rearing batches for each replicate.

To examine the effect of rapid cold hardening on cold tolerance, cold-shock survival was determined using a 2 h exposure to subzero temperatures. Groups of 10 D. variabilis or A. maculatum (N=7 groups) were placed in 1.5 cm³ mesh-covered tubes and these tubes were placed in 50 ml vials, which were plugged with insulating foam then suspended in a chilled water:ethylene glycol (40:60) solution. Temperatures were maintained (±0.1°C) using a programmable refrigerated bath (Arctic A25; Thermo Scientific, Pittsburgh, PA, USA). The temperature inside the tick-containing tubes equilibrated with the bath temperature within approximately 10 min, so an incubation time of 2 h 10 min was used for all trials. Cold-shocked ticks were transferred directly from 26°C to the subzero temperature, whereas RCH ticks were exposed to 4°C for 2 h prior to the subzero temperature. Following treatment, ticks were returned to rearing conditions and assessed for survival after 48 h recovery and were considered alive if they could self-right and/or move several body lengths spontaneously or in response to being handled.

For RNA sequencing (RNA-Seq) and metabolomics analyses, female D. variabilis were individually transferred to 1.5 cm³ mesh-covered tubes and exposed to the following temperature treatments using the same general method described above: control, 26°C for 6 h; CS, 26°C for 2 h followed by −5°C for 2 h, and a 2 h recovery at 26°C; RCH, 4°C for 2 h, followed by −5°C for 2 h, and a 2 h recovery at 26°C. At the end of the treatment, ticks were flash frozen and stored at −80°C. Following low-temperature exposure, some ticks show signs of life after 2 h of recovery, only to succumb to death several hours later. Therefore, a temperature that resulted in no mortality (−5°C) was chosen for RNA-Seq treatments to prevent the confounding effects on transcriptional changes in dying ticks.

To obtain adequate RNA concentrations and reduce individual variation, each RNA sample consisted of 8 D. variabilis that were removed from −80°C and immediately homogenized together. Ticks were manually cut into small pieces, placed in 1 ml of chilled TRIzol^® reagent (Invitrogen, Carlsbad, CA, USA) and homogenized using a BeadBlaster 24 microtube homogenizer (Benchmark Scientific, Edison, NJ, USA). Total RNA was extracted in TRIzol following the manufacturer's protocol and then treated with DNase I (Thermo Scientific) to eliminate potential genomic DNA contamination. RNA was concentrated using the GeneJET RNA Cleanup and Concentration Micro Kit (Thermo Scientific), and the final RNA concentration and purity were determined with a Nanodrop (Thermo Fisher Scientific, Waltham, MA, USA).

The DNA Sequencing and Genotyping Core at the Cincinnati Children's Hospital Medical Center (CCHMC) constructed the poly(A) library and performed the sequencing. RNA was quantified using the Qubit 3.0 Fluorometer (Life Technologies, Carlsbad, CA, USA), and 150–500 ng of total RNA was poly(A) selected and reverse transcribed using the TruSeq Stranded mRNA Library Preparation Kit (Illumina, San Diego, CA, USA). Samples were fitted with one of 96 adapters containing a different 8-base molecular barcode for high-level multiplexing. Following 15 cycles of PCR amplification, completed libraries were sequenced on a HiSeq 2500 sequencing system (Illumina) in Rapid Mode. Reads were single-end and 75 bases in length with approximately 30 million high-quality reads being generated per sample. Raw RNA-Seq data were uploaded to the National Center for Biotechnology Information's (NCBI) Sequence Read Archive: BioProject PRJNA657863.

Sequences generated by Illumina sequencing were trimmed for ambiguities (0 ambiguities were allowed) and quality (using a Phred quality score limit of 0.05). The 5′ and 3′ termini were trimmed to remove 5 and 8 nucleotides from the sequence, respectively, and resulting sequences less than 40 nucleotides in length were discarded. The quality of the resulting cleaned sequences was verified with the FastQC package (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Transcripts were de novo assembled using three separate assembly programs: Trinity (Grabherr et al., 2011), Velvet-Oases (Schulz et al., 2012) and CLC Genomics Workbench (Qiagen, Redwood City, CA, USA), with a minimum contig length of 200. Contigs generated from these programs were combined into a single assembly and redundant contigs were removed by analyzing the assembly with the CD-HIT-EST (Huang et al., 2010) clustering algorithm using a 95% similarity threshold. Quality of the assembly was assessed by the identification of orthologs based on BUSCO and CEGMA gene sets (Parra et al., 2007; Simão et al., 2015).

Reads were mapped to the de novo assembly using CLC Genomics with at least 70% of the read having at least 80% identity with the reference and a mismatch cost of 2. Reads were allowed to align to no more than 10 different places on the reference assembly. Expression values were measured as reads per kilobase of transcript, per million mapped reads (RPKM). Significant changes in expression among samples were analyzed using the ‘differential expression for RNA-Seq tool’ in CLC Genomics. This analysis utilizes a negative binomial Generalized Linear Model and differences between treatment groups are determined with a Wald test. Transcripts were considered differentially expressed when the false discovery rate (FDR) P-value was ≤0.05 and the fold change was ≥|2|.

To identify differentially expressed contigs, sequences were searched (BLASTx) with an expect value (E-value) threshold of 0.001 against NCBI's arthropod, non-redundant (nr) protein database, the SwissProt protein database, and the reference protein sequences (RefSeq) for fly (Drosophila melanogaster) and the black-legged tick Ixodes scapularis. The highest scoring BLAST hit in each search was used to assign a gene ID to each contig. Enrichment of functional pathways was identified by submitting hits from the I. scapularis search to the PANTHER Gene List Analysis (Mi et al., 2017), DAVID functional annotation database (Huang et al., 2009) and g:Profiler Gene Group Functional Profiling (Reimand et al., 2016). Gene ontology (GO) terms were analyzed for over-representation, and functional annotations were considered significantly enriched when the corrected P-value was ≤0.05. A full list of GO terms that were enriched in two or more of the analyses was generated for each treatment and the number of significant GO terms was reduced by using REVIGO (Supek et al., 2011) to eliminate redundant terms. Lastly, specific unidentified contigs were searched (BLASTx) with an E-value threshold of 0.001 against the recently completed tick genomes (Jia et al., 2020) to identify tick-specific and Dermacentor-specific contigs.

To validate the RNA-Seq results, qPCR validation was performed as described in Rosendale et al. (2016b) with a separate group of D. variabilis (N=5–6) exposed to CS or RCH treatments with recovery periods as described for RNA-Seq. Ticks were individually homogenized in 1 ml of chilled TRIzol^® reagent and RNA was extracted and concentrated as described for RNA-Seq. cDNA was made from approximately 300 ng RNA using the DyNAmo cDNA Synthesis Kit (Thermo Scientific). qPCR reactions consisted of KiCqStart SYBR Green qPCR ReadyMix (Sigma Aldrich, St Louis, MO, USA), 300 nmol l⁻¹ forward and reverse primers, cDNA diluted 1:50 and nuclease-free water. Primers were designed based on sequences obtained from the RNA-Seq analysis (Table S1). These contigs were chosen as their expression levels in RNA-Seq were neither particularly high nor low, and their fold-changes did not exceed the median fold-change for differentially expressed contigs. qPCR reactions were analyzed using an MX3005p Real-time PCR System (Agilent Technologies). PCR consisted of activation for 3 min at 95°C followed by 40 cycles of denaturation (10 s at 95°C), annealing/extension (30 s at 55°C) and denaturation (10 s at 95°C). Following amplification, a melt curve analysis was performed from 55 to 95°C with 0.5°C increments every 15 s. Samples were run in triplicate and the average C_q value was determined. The ΔΔC_q method (Schmittgen and Livak, 2008) was used to determine expression levels, using β-actin to normalize genes of interest. Actin was chosen as the reference gene as it has been used successfully as a reference gene in other tick stress studies (Rosendale et al., 2016b); additionally, there was no difference in the expression of this gene among our samples (control versus CS FDR P=0.99; control versus RCH FDR P=0.92). The fold-change in these genes was determined and the logarithmic fold-change was plotted against the corresponding value from the RNA-Seq analysis, and a Pearson correlation coefficient (r) was determined.

A nuclear magnetic resonance (NMR)-based metabolomics approach was used to identify metabolites present in ticks. NMR experiments were performed at the NMR-based Metabolomics Core at CCHMC. A total of 15 D. variabilis per sample were lyophilized overnight and the dried samples were weighted into 2 ml standard tubes containing 2.8 mm metal beads (Bertin Corp., Rockville, MD, USA). Ticks were homogenized 3 times for 30 s at 4000 rpm with a Minilys homogenizer (Bertin Corp.) and polar metabolites were extracted using a modified Bligh and Dyer extraction (Bligh and Dyer, 1959). Briefly, cold methanol and water were added to the samples in bead tubes and homogenized for 30 s. The samples were transferred into glass tubes containing cold chloroform and water (methanol:chloroform:water ratio of 2:2:1.8). The mixture was vortexed, incubated on ice for 10 min, and centrifuged at 2000 g for 5 min. The polar phase was transferred into a 1.5 ml tube and dried by vacuum centrifugation for 2–3 h at room temperature. The dried metabolites were re-suspended in 0.6 ml of NMR buffer containing 100 mmol l⁻¹ phosphate buffer (pH 7.3), 1 mmol l⁻¹ TMSP [3-(trimethylsilyl)2,2,3,3-d4 propionate], and 1 mg ml⁻¹ sodium azide prepared in deuterium oxide. A final volume of 550 μl of each sample was transferred into a 5 mm NMR tube (Norell Inc., Marion, NC, USA) for NMR data acquisition.

Data were acquired using one-dimensional ¹H Nuclear Overhauser Effect Spectroscopy (NOESY) NMR experiments on a 600 MHz INOVA NMR spectrometer (Agilent Technologies, Santa Clara, CA, USA). Representative samples were analyzed with two-dimensional ¹H-¹³C heteronuclear single quantum correlations (HSQC) and ¹H-¹H Total Correlation Spectroscopy (TOCSY) for metabolite annotation. Spectral data were phase corrected, baseline corrected, and aligned to the internal reference standard peak, TMSP, using TopSpin 3.1 (Bruker, Billerica, MA, USA).

Metabolite identity was assigned using Chenomx NMR Suite 7.8 (Chenomx Inc., Edmonton, AB, Canada), two-dimensional NMR experiments, and reference spectra found in databases such as the Human Metabolome Database (Wishart, 2007), the Madison metabolomics consortium database (Cui et al., 2008) and the biological magnetic resonance data bank (Markley et al., 2008). Metabolite concentrations were determined using the Chenomx NMR Suite integration tool and normalized to the dry mass of the ticks. Normalized metabolite concentrations were compared using a one-way ANOVA with Tukey post hoc tests.

To evaluate their cryoprotective effects, several of the metabolites that were found to be elevated following recovery from RCH were injected into female D. variabilis and ticks were subsequently exposed to subzero temperatures. Dermacentor variabilis (N=10 ticks per group) were injected with either insect Ringer's solution (187 mmol l⁻¹ NaCl, 21 mmol l⁻¹ KCl, 7 mmol l⁻¹ CaCl₂ and 1 mmol l⁻¹ MgCl₂) or insect Ringer's solution with 1, 0.1 or 0.01 mol l⁻¹ alanine, betaine, glycine or valine. Additionally, male and female D. variabilis and female A. maculatum were injected with Ringer's or 1 mol l⁻¹ betaine solution to allow for a sex and species comparison, respectively. The cuticle was pierced with a 28.5-gauge needle on the ventral idiosoma just posterior to the fourth coxa, a tapered polycarbonate tip (0.19 mm outer diameter) was inserted into the body cavity, and the solution (∼0.5 μl) was dispensed by a manual microdispenser (Drummond, Broomall, PA, USA). After injection, ticks recovered for 2 h at 26°C and 100% RH and were then exposed to subzero temperatures (−13.5°C for D. variabilis and −12.0°C for A. maculatum). Following CS, ticks were returned to 26°C and 100% RH and survival was assessed 48 h post-treatment. Survival was compared among Ringer's solution- and metabolite-injected ticks using a Student's t-test or ANOVA.

To determine whether adult D. variabilis demonstrate a typical RCH response, survival of several subzero temperatures was examined. There was no difference (F_2,15=0.18, P=0.84) in the dry mass of ticks in the control, CS or RCH groups (3.11±0.09, 3.07±0.16, 3.17±0.11 mg, respectively, means±s.e.m.), nor did treatments affect (F_2,15=0.25, P=0.78) water content of the ticks (approximately 1.3 mg water mg⁻¹ dry mass, 57%, for all groups). Survival was high (100%) for both CS and RCH groups exposed to −5°C; however, survival of the CS group dropped below 50% at temperatures lower than −8°C (Fig. 1). Overall, survival of RCH ticks was greater (F_1,72=48.4, P<0.0001) than that of CS ticks, with a 2 h pre-exposure to 4°C improving survival of −8 to −12°C by 50–30%.

Nine cDNA libraries were generated from the total homogenates of control ticks and those that recovered from CS or RCH (3 per group). These libraries were sequenced using Illumina HiSeq technology, resulting in 391,281,556 cleaned reads. Sequences from all samples were used in the de novo transcriptome assembly using the Trinity, Velvet-Oases and CLC Genomics Workbench assembly programs. Contigs from each of these three assemblies were combined into a single assembly, and duplicate contigs were removed with CD-HIT-EST (Table 1). Contig sizes ranged from 200 to ∼14,000 bases with the number of contigs decreasing with increasing contig length. Contigs were annotated by searching (BLASTx) NCBI's arthropod nr and SwissProt databases as well as RefSeq sequences from D. melanogaster and the Ixodes genus. The majority of the matches were to the black-legged tick, I. scapularis (Fig. 2). BUSCO and CEGMA scores were comparable with previous D. variabilis results (Rosendale et al., 2016a,b, 2017), where over 96% of BUSCO genes from the I. scapularis predicted gene set were present in our assembly.

Differential expression of transcripts was analyzed among the control ticks and those that recovered from CS or RCH. Seventy-three percent (62–80%) of the reads from control ticks, 80% (78–84%) from CS and 78% (77–79%) of reads from RCH ticks mapped to the reference assembly. We found 99, 1860 and 1376 genes with expression levels significantly different between control and CS, control and RCH, and CS and RCH ticks, respectively (Table 2). When compared with control ticks, most (∼65%) of the genes differentially expressed in CS ticks also showed expression changes when the CS was preceded with RCH (Fig. 3). The full set of differentially regulated genes is presented in Table S2. To validate our RNA-Seq results, qPCR was used to measure the expression of several genes, normalized to β-actin. qPCR and RNA-Seq results were similar based on the Pearson correlation coefficient (r=0.952; Fig. S1).

In ticks exposed to CS, there was an upregulation of several genes related to glucose metabolism, chaperone and repair; with most genes being salivary gland proteins. When analyzed for GO terms, only carbohydrate biosynthetic process and monosaccharide biosynthetic process were enriched in the upregulated genes of CS ticks . There was no GO term enrichment in the downregulated genes of CS ticks. In contrast to the CS group, ticks exposed to the RCH conditions showed a large number of differentially expressed genes. For genes upregulated during recovery from RCH and a subsequent CS, there was an enrichment of a variety of GO terms (Fig. 4; Fig. S2). GO terms that were enriched in RCH ticks included biological processes such as protein metabolic process, transport and oxidation–reduction process; molecular functions such as catalytic activity and binding; and cellular components such as ribosomes and cytoskeleton. For genes that were downregulated in response to RCH, there was an enrichment of various GO terms including metabolic processes, binding and organelles. The full list of GO terms that were differentially expressed in two or more analyses (see Materials and Methods) for all comparisons is shown in Table S3.

NMR-based metabolomics was used to identify differences in metabolite concentrations among control ticks and those that recovered from CS or RCH. Using spectra from ¹H NMR, the concentrations of 29 metabolites were identified (Table S4), of which six showed a significant difference between the control and RCH groups (Table 3). There were no differences in metabolite concentration between the CS group and the other two groups.

Based on our combined metabolomic and transcriptomic results, several metabolites were identified that accumulated during the RCH recovery process. To determine whether these metabolites have cryoprotective effects, ticks were injected with different concentrations (1, 0.1 or 0.01 mol l⁻¹) of alanine, betaine, glycine or valine solutions and subjected to subzero temperatures. Injection of 0.5 μl of these solutions would theoretically increase the concentration of these metabolites by approximately 150, 15 or 1.5 nmol mg⁻¹ dry mass; these lower concentrations are within the range of observed increases of these metabolites (Table 3). Survival of exposure to −13.5°C was significantly (P<0.05) affected by injection of betaine or valine as compared with that of ticks injected with insect Ringer's solution (control; Fig. 5A–C). However, there was no effect of concentration of these metabolites as D. variabilis injected with either betaine or valine at 1, 0.1 and 0.01 mol l⁻¹ showed no statistical difference in survival (Fig. 5A–C). In D. variabilis, female survival was ∼4-fold higher (t₂₀=5.07, P<0.0001) in the betaine-injected groups, while in males survival was ∼2-fold higher (t₂₂=2.26, P=0.034) (Fig. 5D). Amblyomma maculatum exposed to −12°C showed a classic RCH response, and cold hardiness was also improved through injection of betaine (Fig. 5E). Amblyomma maculatum receiving RCH treatment showed an 8-fold higher (P<0.01) survival than the CS group, whereas betaine injection improved survival 2.7-fold (P<0.05) over the Ringer's solution (control) group.

The American dog tick can transmit various diseases, including Rocky Mountain spotted fever, and has a very wide geographic distribution, making it an important disease vector in North America (de la Fuente et al., 2008). Augmentation of the geographic range of D. variabilis suggests an increase in the impact of this species as a vector (Dergousoff et al., 2013), and changes in distribution are occurring in a variety of tick species (Sonenshine, 2018). The geographic range of ticks is often limited by climate, including the ability to survive low winter temperatures (Dantas-Torres and Otranto, 2011; Dautel et al., 2016; Randolph, 2004); therefore, understanding winter survival may be important in predicting population dynamics. However, little is known about the mechanisms underlying low temperature survival in ticks. In the present study, we examined the response of ticks to CS and RCH and found changes at the transcriptome and metabolome levels that revealed multiple biological pathways that are likely to be important in the recovery from these conditions. We also demonstrated the ability of betaine to improve cold hardiness in ticks.

By exposing ticks to a brief (2 h) period at 4°C prior to exposure to subzero temperatures, cold hardiness was dramatically improved. Adult D. variabilis displayed the classic RCH response, which is common in arthropods and has been documented in multiple classes (Teets et al., 2020). Although RCH has been most thoroughly examined in insects, it has been documented in Acari in both mites (Broufas and Koveos, 2001; Ghazy and Amano, 2014) and ticks (Rosendale et al., 2016a; Wang et al., 2017; Yu et al., 2014). However, in ticks, the RCH response does not seem to be ubiquitous. Ixodes scapularis seems to lack a RCH response (Vandyk et al., 1996), whereas H. longicornis (Yu et al., 2014) and A. maculatum (present study) undergo RCH. Even within a genus (Dermacentor), there is variation in the RCH response; D. variabilis (present study; Rosendale et al., 2016a,b) and D. silvarum (Wang et al., 2017) have both been shown to undergo RCH, whereas RCH is not found in D. albipictus (Holmes et al., 2018). This discrepancy may be an artifact of methodologies used in these various studies as certain acclimation temperatures and/or durations are more effective at eliciting a RCH response in ticks than others, and the physiological state of ticks can greatly impact stress tolerance (Rosendale et al., 2016a, 2017; Wang et al., 2017; Yu et al., 2014). In H. longicornis, RCH results in increased water content (Yu et al., 2014); however, in adult D. variabilis (present study) and D. silvarum (Wang et al., 2017), water content does not change, suggesting that alteration in water content is not a critical component of the RCH response in ticks.

Following recovery from CS and RCH, a suite of genes that likely contribute to cold hardiness was differentially regulated. For the CS ticks, a relatively small number of genes (99) were differentially regulated as compared with the much larger number that changed in the RCH group (1860), and there were many differences between CS and RCH ticks (1376 genes). This is different from what was observed in flesh flies exposed to similar CS and RCH treatments (Teets et al., 2012). In Teets et al. (2012), gene expression profiles after 2 h recovery were nearly identical between flies directly exposed to a CS and those given RCH prior to CS. Because of our experimental design, it is impossible to determine whether the transcriptomic changes in the RCH group occurred during the brief acclimation period or during the recovery phase; however, studies on flies suggest that it is the recovery period that is critical for the changes in gene expression (Sinclair et al., 2007; Teets et al., 2012). Although many of the changes may occur during recovery, our data indicate that in ticks, a brief acclimation period triggers a much more robust transcriptomic response than CS alone. The actual mechanism that accounts for this more robust response remains unclear. It is possible that the RCH acclimation is necessary to trigger a transcriptomic response following recovery from cold. Alternatively, the RCH treatment may allow ticks to more quickly restore physiological functioning to alter gene expression following cold exposure, whereas ticks directly exposed to CS may require a longer recovery period.

This transcriptomic response included an upregulation of 728 genes that contained a variety of genes and/or GO terms that likely contribute to enhanced cold hardiness (Fig. 6). Many of these genes were categorized as binding, including calcium binding and iron ion binding, both of which are important mechanisms in the RCH response (Gerken et al., 2015; Teets et al., 2020). Adjustment to membrane transport properties is important in both seasonal acclimation and RCH (Armstrong et al., 2012; Koštál et al., 2006; Teets and Denlinger, 2013), and long-term acclimation leads to an upregulation of transport-related genes (Enriquez and Colinet, 2019). In ticks that recovered from RCH, 20 genes related to transporter activity were upregulated, including several specifically related to ion transport, which could be important in restoring ion homeostasis during CS recovery (MacMillan et al., 2012; Teets and Denlinger, 2013). Structural components of the cell, including the plasma membranes and cytoskeleton, are vulnerable to damage from cold, and an important RCH mechanism is modification to these structures (Kim et al., 2006; Li and Denlinger, 2008; Teets and Denlinger, 2013). In our ticks that recovered from RCH, several GO terms related to the cytoskeleton were enriched in the upregulated genes and 44 genes related to the membrane were upregulated. In addition, several transcripts with similarity to heat-shock protein 70 (Hsp70) were upregulated. Hsps act as molecular chaperones to protect proteins, contributing to cold tolerance, and Hsp70 is specifically important in surviving low temperature (Rinehart et al., 2007).

Among the groups of ticks that recovered from RCH and CS, there were 56 common transcripts that were upregulated; for these genes, the GO terms carbohydrate biosynthetic process and gluconeogenesis were enriched. Genes of note include glucose 6-phosphatase (G6P) and phosphoenolpyruvate carboxykinase (PEPCK). Although accumulation of carbohydrate-based cryoprotectants is a hallmark response to long-term cold acclimation, RCH does not always result in increased cryoprotectant levels (Teets et al., 2012). Glycerol has been suggested as a cryoprotectant in ticks (Wang et al., 2017; Yu et al., 2014); however, neither glycerol nor glucose (another common cryoprotectant) changed levels in our experiments (discussed below). Alternatively, this pathway could be a common stress response in ticks as PEPCK is upregulated in dehydrated ticks (Rosendale et al., 2016b) and contributes to tolerance of oxidative stress in tick cells (Della Noce et al., 2019).

In overwintering D. silvarum, molecular functions such as catalytic activity and binding, cellular components including membrane, and biological processes such as oxidation–reduction process and metabolic process were highly enriched GO terms (as compared with summer ticks) (Yu et al., 2020). For most of the GO terms enriched in winter D. silvarum, there were genes that were upregulated in our RCH recovery ticks that matched those categories (Fig. 7). Enrichment of these terms in response to both long- and short-term cold exposure suggests the importance of these pathways in survival of low temperature.

Accumulation of various molecules during the cold-hardening process can serve a protective role against cold injury (Teets and Denlinger, 2013). Dermacentor variabilis showed several changes to the metabolome during recovery from RCH that were supported by the transcriptome data. Levels of several amino acids, including alanine, glycine and valine, were significantly increased following recovery from RCH. This augmentation of amino acid levels was potentially the result of increased protein catabolism, as there was an upregulation of proteolysis genes. Amino acids such as proline can act as potent cryoprotectants when they are accumulated prior to cold exposure (Koštál et al., 2011). Neither alanine nor glycine impacted CS survival of D. variabilis in our experiments; however, when valine was exogenously supplied prior to CS treatment, survival improved.

Valine has been suggested as one of multiple amino acids that improves seasonal cold hardiness in insects (Feng et al., 2016), and our data support the cryoprotective nature of this amino acid. Increases in valine concentration in arthropods have been noted both in response to long-term, low-temperature acclimation (Storey et al., 1986) and seasonally in overwintering insects (Feng et al., 2016; Qiang et al., 2012). Valine has been suggested to improve cold hardiness by contributing to the lower super-cooling point in winter-acclimated Lepidoptera (Qiang et al., 2012). However, while proline has been extensively studied in cold and freeze tolerance (Koštál et al., 2011), less is known about the specific cryoprotective mechanisms linked to other amino acids, such as valine. Teets et al. (2012) noted that RCH treatment alone was not enough to elicit an accumulation of valine in flies. If the valine in ticks is accumulated during the recovery period, it may serve a role in the restoration of homeostasis and/or may contribute to survival of subsequent cold exposure. Alternatively, it is possible that the protein breakdown is due to a disruption in homeostasis caused by CS damage, or proteolysis may be a common stress response in ticks as dehydrated D. variabilis (Rosendale et al., 2016b) and overwintering D. silvarum (Yu et al., 2020) also show an upregulation of this pathway.

The accumulation of betaine in ticks that recovered from our RCH treatment prompted investigation of this molecule as a potential cryoprotectant in ticks. Betaine was investigated more thoroughly than valine for a cryoprotective role as it has not been described as a cryoprotectant in arthropods, whereas some amino acids are known to have a cryoprotective effect (Teets and Denlinger, 2013). The increase in concentration of betaine was coordinated with an upregulation of several transcripts with high similarity to betaine-related genes, including a betaine-aldehyde dehydrogenase and a sodium- and chloride-dependent betaine transporter. Betaine is a common compatible osmolyte that has a protective role during osmotic stress by stabilizing protein structure and membrane integrity (Yancey, 2005). In addition to being an osmoprotectant, betaine also has cryoprotective properties in plants (Xing and Rajashekar, 2001). In the freeze-tolerant Alaskan beetle Upis ceramboides, a betaine-like substance accumulated under winter conditions, although its direct role as a cryoprotectant was not established (Walters et al., 2009). The positive effect of exogenously supplied betaine over a 100-fold range of concentration on the cold hardiness of D. variabilis suggests that this molecule can have a cryoprotective effect if accumulated prior to cold exposure. It has been suggested that betaine acts as a cryoprotectant through both colligative and non-colligative mechanisms in plants (Xing and Rajashekar, 2001). Further study is needed to confirm that the levels of betaine accumulation observed in this study are high enough to have colligative effects and what, if any, non-colligative benefits betaine may provide to cold-hardiness in ticks.

Administration of exogenous valine and betaine improved the cold hardiness of D. variabilis 2- to 3-fold as compared with controls. Neither of these metabolites was as efficacious in improving low-temperature survival as RCH, even at relatively high doses, indicating that a multifaceted response is necessary during RCH. There are several possible explanations for the lack of a concentration effect in our treatments. The cryoprotective effect of metabolites, including amino acids, can have non-colligative mechanisms (Toxopeus et al., 2019), and it is possible that this is the case for valine and betaine in these ticks. Alternatively, valine and betaine may show colligative effects on cold hardiness, but the lowest concentration administered (similar to physiological levels) in our experiments represents the maximal colligatively induced improvement in tick cold hardiness. Further study is needed to elucidate the cryoprotective mechanisms of valine and betaine.

Survival of low temperatures by ticks seems to be related to the accumulation of molecules such as glycerol (Wang et al., 2017; Yu et al., 2014); however, we did not observe a significant change in glycerol levels. There was also no change in glucose, another common cryoprotectant (Overgaard et al., 2007; Teets and Denlinger, 2013). It is possible that glycerol plays a role in the cold hardiness of D. variabilis and no changes in the levels of this molecule were detected as a result of the treatment methods. In H. longicornis, glycerol levels were not affected by 2 h at 0°C but did increase after 10 days at low temperature (Yu et al., 2014). Similarly, D. silvarum showed no change in glycerol levels after 2 h at 0°C, and only females increased glycerol in response to 2 h at −3°C (Wang et al., 2017).

Ticks spend the vast majority of their lives off-host, and establishment of a population is dependent on overwintering survival (Dergousoff et al., 2013; Gray et al., 2009). RCH substantially improves the ability of D. variabilis to survive CS and likely contributes to survival of ticks under ecologically relevant conditions. We found a suite of molecular and biochemical changes in D. variabilis that received RCH treatment prior to CS and a recovery period, and these changes were distinct from those in ticks that experienced CS and recovery only. These changes include an upregulation of genes related to cell structure, transport, signaling and metabolism that are common across a wide range of arthropods (Teets et al., 2020). Because of our experimental design, it is unclear whether these changes occurred during the brief acclimation or during the recovery period. If these changes occurred prior to CS, then they likely contributed to cold hardiness; however, if the recovery period is required, then these changes potentially contribute to subsequent cold exposures and/or contribute to repair and re-establishment of homeostasis. Enhanced cold hardiness in D. variabilis is facilitated by the supplementation of valine and betaine, which improved cold survival in both D. variabilis and A. maculatum. Further research is warranted to determine the cryoprotective mechanisms of physiological levels of betaine. Overall, ticks respond to RCH conditions by activating various pathways to reduce cold-induced damage, some of which directly contribute to survival of cold exposure while others may prepare ticks for subsequent chilling events.

We are grateful to two anonymous reviewers who contributed constructive comments on the paper.