Differential expression profiling of heat stressed tardigrades reveals major shift in the transcriptome

Differential expression profiling of heat stressed tardigrades reveals major shift in the transcriptome


Introduction
Tardigrada constitutes a group of microscopic animals with a welldefined cephalic region and a trunk bearing four pairs of legs . Distributed worldwide, tardigrades, also known as water bears, populate highly diverse microhabitats in marine, freshwater and terrestrial ecosystems (McInnes and Pugh, 2018;Nelson et al., 2018). Tardigrades divide into two major evolutionary lineages represented by eutardigrades and the more diverse heterotardigrades Morek et al., 2020). They are renowned for their extreme stress tolerance (Møbjerg and Neves, 2021) and their ability to enter cryptobiosis-a reversible state of suspended animation-which provides them with an extreme tolerance towards severe environmental conditions, including complete desiccation, high external solute concentrations, extremely low temperatures, high toxicant concentrations and oxygen depletion (Clegg, 2001;Guidetti et al., 2011;Møbjerg et al., 2011;Hygum et al., 2017;Sørensen-Hygum et al., 2018). Importantly, tardigrades need to be surrounded by a film of water, in order to remain in their physiologically active state. During desiccation, when tardigrades may lose more than 95% of their water (Crowe, 1972), they contract longitudinally, retract head and legs, while forming a quiescent and barrel shaped form, the so-called "tun" (Crowe and Madin, 1974;Bertolani et al., 2004;Halberg et al., 2013).
Albeit their ability to endure a range of extreme conditions, tardigrades were recently reported to be sensitive to prolonged periods of high temperature exposure in active as well as desiccated tun states (Neves et al., 2020a(Neves et al., , 2020b. Specifically, thermotolerance experiments on desiccated tuns of Ramazzottius varieornatus, a eutardigrade frequently found in European roof-gutters, revealed a decrease in the 50% mortality temperature from 82.7 • C, over 63.1 • C to 56.1 • C following 1 h, 24 h and 1 week of high temperature exposures, respectively (Neves et al., 2020a(Neves et al., , 2020b. In fully hydrated and thus metabolically active specimens, the 50% mortality temperature was 37.1 • C following a 24 h exposure, with a small increase to 37.6 • C following a brief acclimation period (Neves et al., 2020a). Thus, although selected tardigrade species, such as Macrobiotus hufelandi, historically have been acknowledged for their high temperature tolerance of above 100 • C for up to 30 min (Doyère, 1842;Rahm, 1921), it is now also clear that extremotolerant tardigrades, such as R. varieornatus, may be vulnerable to elevated temperatures, and thus the potential effects of global warming (Neves et al., 2020a). Nevertheless, tardigrades occurring in limno-terrestrial and marine tidal environments are most likely sidestepping the damaging effects of heat exposure by undergoing desiccation induced cryptobiosis and, hence, being able to survive summer periods in a warmer Earth (Sørensen-Hygum et al., 2018;Vecchi et al., 2021). Accordingly, it remains to be seen to which extent the ongoing global warming will affect tardigrade populations in general.
The limited thermotolerance of R. varieornatus suggests an interference of high temperatures with tardigrade cell and macromolecule function (Møbjerg and Neves, 2021). Hence, bioprotectants and other macromolecules, which otherwise provide tardigrades with their extraordinary stress tolerance, might indeed be denatured following prolonged heat exposure. A variety of proteins with proposed cytoprotective function are found among tardigrades. For instance, heat shock proteins (HSPs) and late embryogenesis abundant (LEA) proteins may act as chaperones and molecular shields, preventing protein aggregation during desiccation, and HSPs might further participate in the repair of damaged proteins during rehydration (Schill et al., 2004;Schokraie et al., 2010). Noteworthy, the expression of HSP70 seems to be induced by heat stress in the eutardigrade Richtersius cf. coronifer (Jönsson and Schill, 2007) and in Milnesium inceptum the expression of a small α-crystallin heat-shock protein gene (Mt-sHsp17.2) increased 779-fold following heat stress (Reuner et al., 2010). More recently, it has been suggested that heat-soluble "tardigrade-unique proteins" (Yamaguchi et al., 2012;Tanaka et al., 2015;Hashimoto et al., 2016), also known as tardigrade-specific intrinsically disordered proteins, may provide structural stabilization to eutardigrade cells during desiccation (Boothby et al., 2017;Yagi-Utsumi et al., 2021). These intrinsically disordered proteins, including SAHS (secretory abundant heat-soluble proteins), CAHS (cytoplasmic abundant heat-soluble proteins), MAHS (mitochondrial abundant heat-soluble proteins) and RvLEAM (a group 3 LEA protein located in the mitochondria) are missing in the heterotardigrade lineage (Kamilari et al., 2019;Murai et al., 2021). Eutardigrades furthermore express an unstructured nuclear protein, Dsup, which supposedly binds to nucleosomes and protects DNA from reactive hydroxyl radicals (Chavez et al., 2019;Minguez-Toral et al., 2020). Importantly, transcriptome analyses have shown that extremotolerant tardigrade species seem to constitutively express their genes during cryptobiotic survival (Hashimoto et al., 2016;Murai et al., 2021).
Here, we investigate the transcriptomic response to heat stress in R. varieornatus with the aim of providing new insights into the molecular processes affected by high temperatures. Specifically, we compare the transcriptome of physiologically active, heat-exposed (35 • C) tardigrades to that of physiologically active controls kept at 5 • C. To the best of our knowledge, the current study presents the first analysis of RNAseq datasets obtained from heat stressed tardigrades. Our transcriptome profiling reveals a surprising shift in the transcriptome comprising 9634 differentially expressed genes, corresponding to ~35.6% of the transcriptome. Among the very large number of differentially expressed genes, several code for eutardigrade specific proteins such as SAHS, CAHS, MAHS, Dsup, and RvLEAM. In addition, two stressrelated heat-shock proteins were found to be upregulated, namely a HSP70 and a small HSP, while other heat-shock proteins were found to be downregulated. A gene ontology analysis suggests a change in metabolism, transcription, and possibly translation as a response to the heat-shock. Our findings support previous observations revealing that high temperatures (> 30 • C) are considerably more stressful for the eutardigrade R. varieornatus than extreme conditions, such as desiccation (Hashimoto et al., 2016) or freezing at very low sub-zero temperatures (-80 • C) (unpublished observations).

Methods
A graphical representation of the methods used to evaluate the effect of exposing fully hydrated, active tardigrades to 35 • C as compared to 5 • C is provided in Fig. 1. In the figure, grey text and arrows refer to a reference transcriptome available from the European Nucleotide Archive under the project PRJEB47629, which was used for mapping the reads and quantifying differential gene expression.

Sampling of specimens, experimental setup and total RNA extraction
Specimens of Ramazzottius varieornatus Bertolani and Kinchin, 1993 (Eutardigrada, Ramazzottiidae) were obtained from sediment collected from a roof gutter in Nivå, Denmark (55 • 56 ′ 36.53" N, 12 • 30 ′ 00.90 ′′ E), in February 2018, and identified using the original taxonomic description of the species (cf. Bertolani and Kinchin, 1993). The sediment sample was frozen under wet conditions, stored at − 20 • C for several months and subsequently transferred to − 80 • C. Two subsamples were thawed, diluted in ultrapure water (Millipore Milli-Q® Reference, Merck, Darmstadt, Germany) and acclimated to 5 • C in May and June 2020, respectively. The subsamples were recurrently examined for tardigrades between June and July 2020 with the aid of a stereomicroscope. Highly active, adult (large) specimens of R. varieornatus showing spontaneous movement of legs were collected from the subsamples and transferred to embryo dishes using hand pulled Pasteur pipettes, at room temperature (RT; 21-25 • C). Subsequently, three groups of ca. 100-180 tardigrades eachi.e. three replicateswere transferred to Eppendorf tubes with 1.5 ml of moderately hard reconstituted water (MHRW; Khanna et al., 1997) and then exposed to 35 • C using a AccuBlock Digital Dry Bath heating system (Labnet International, Edison, NJ) for approximately 24 h. After the heat shock, the tardigrades were briefly checked for locomotor activity and processed to extract total RNA. In addition, three groups of ca. 100-180 tardigrades each were kept refrigerated (5 • C) for 24 h as controls, and then processed to extract total RNA. Total RNA was extracted using the RNeasy Plus Universal Mini Kit (Qiagen, Hilden, Germany) according to the instructions of the manufacturer, with the following modifications: for disruption of the tardigrade cuticle, the specimens were processed in BeadBug™ tubes (Benchmark Scientific, Sayreville, NJ) containing 2.3-mm and 1-mm diameter zirconia/silica beads during 2 × 1 min homogenization at 50 Hz. Initial quantity and quality of the extracted RNA were evaluated with a NanoDrop® ND-1000 spectrophotometer (Peqlab Biotechnologie GmbH, Erlangen, Germany).

mRNA library construction and sequencing
The total RNA samples were sent to BGI Europe A/S for library preparation and sequencing. Concentration and purity of each RNA sample were assessed by BGI using an Agilent 4200 (Agilent Technologies, CA, USA) before mRNA was selected, reverse transcribed and sequenced using the DNA nanoball sequencing (DNBseq-G400) platform. A total of 376,485,664 (100 bp) paired end reads were generated (Table 1). We subsequently assessed the quality of the raw RNA-Seq data produced by BGI with the software FastQC (Andrews, 2010) and the produced reads were then subject to trimming of adapters and lowquality stretches using AdapterRemoval v2.0 (Schubert et al., 2016). The quality of the reads was once again evaluated after adapter trimming using FastQC to ensure that the adapter sequences were properly trimmed (Ewels et al., 2016).

Differential expression analyses
We used a reference transcriptome of the Danish population of R. varieornatus, recently made available (ENA PRJEB47629), for mapping sequence reads and a pseudo-aligner, Salmon version 1.1.0 (Patro et al., 2017), to quantify transcript abundance in each of the six RNA-Seq datasets (i.e. the three datasets from heat-exposed (35 • C) tardigrades and three datasets from non-exposed tardigrades, kept at 5 • C). The NumbRead metric estimated by Salmon was subsequently used for three differential expression analyses using the Bioconductor packages DESeq2 (Love et al., 2014), EdgeR (Robinson et al., 2010) and Limma (Ritchie et al., 2015). Only genes that were expressed in more than half of the RNA-Seq datasets were considered in the differential expression analyses (19,456 quantified transcripts were included, whereas 133 did not meet this criterion). The quasi-likelihood F-test was used for EdgeR, as this test is effective at controlling the false discovery rate (FDR) and Fig. 1|. Graphical representation of the methods used to evaluate the effect of exposing fully hydrated, active tardigrades to 35 • C. Active specimens were randomly pooled into groups (3 × ca. 100-180) at room temperature (RT), either exposed to 35 • C or kept at 5 • C (controls) for 24 h and subsequently used to extract total RNA for differential gene expression analyses. Sequenced reads were mapped against a reference transcriptome. Grey text and arrows concern the preparation of this reference transcriptome (ENA PRJEB47629), which was assembled and annotated for another study.

Table 1|
Summary data for RNA samples generated from active heat exposed (replicate #1-3) and control (control #1-3) specimens of Ramazzottius varieornatus. most appropriate for experiments with few replicates. An FDR adjusted p-value of ≤0.05 was used to determine whether a gene was up-or down regulated in the heat-exposed (35 • C) tardigrades as compared to the controls kept at 5 • C. The Venn diagram (Fig. 2) and MA-Plot (Fig. 3) were generated in R (R Core Team, 2016) and edited with Illustrator CS5 version (Adobe Inc.).

Enrichment analysis
We performed a Gene Ontology (GO) enrichment analysis on the transcripts identified in the differential expression analyses using TopGO version 2.42.0 (Alexa and Rahnenfuhrer, 2020), with GO terms derived from the annotated reference transcriptome (ENA PRJEB47629). Fisher's exact test was implemented to compute the significant GO enrichments at the α-level of p ≤ 0.05.

Results
Gene expression profiles for the heat-exposed (35 • C) and the nonexposed (5 • C) tardigrades revealed a substantial number of differentially expressed genes (DEGs), independently of the method used to perform the differential expression analysis. Indeed, a total of 11,015 different DEGs were identified across the three methods (DESeq2, EdgeR and Limma): 9634 of these transcripts (i.e., 35.6% of the 27,054 assembled unigenes in the reference transcriptome) were found by all three methods as shown in the Venn diagram represented in Fig. 2. Limma, EdgeR and DESeq2 exclusively identified 199, 178 and 30, respectively. 32 DEGs were shared only by EdgeR and DESeq2, 87 DEGs by Limma and DESeq2 and 855 DEGs were shared only by EdgeR and Limma. Out of the 9783 DEGs found by DESeq2, 4901 of the transcripts were upregulated and 4882 were downregulated in respect to the controls ( Fig. 3; Table 2). EdgeR identified 10,699 DEGs for the heating experiments and, out of these, 5382 were upregulated and 5317 were downregulated. Limma identified 10,775 DEGs of which 5489 were upregulated and 5286 were downregulated (Table 2). Seven GO molecular function categories were found to be overrepresented among the transcripts identified by the differential expression analyses (p ≤ 0.05), namely: ATPase activity, protein heterodimerization activity, sequencespecific DNA binding, structural constituent of ribosome, chitin binding, NADH dehydrogenase (ubiquinone) activity and translation initiation factor activity ( Table 3).
Proteins that have been shown previously to be involved in stress responses in tardigrades (e.g., heat shock and tardigrade disordered proteins), were sought after among the annotated DEGs. In our analysis, 12 Heat shock transcripts were identified as differentially expressed by all three methods. Two of these transcripts were upregulated, whereas the remaining 10 were downregulated. The two upregulated transcripts were a HSP70 and a small HSP (annotated as HSP67Bb), while the downregulated transcripts encompassed one HSP90, three HSP70s isoforms, one HSP68, one HSP60, and four small HSPs. Among the upregulated transcripts, two SAHS genes, i.e. SAHS1 and SAHS2, were identified by all three methods. Seven CAHS transcripts were annotated and identified as DEGs by all three methods (i.e. CAHS1, CAHS2, CAHS3, CAHS68135, CAHS77611, CAHS86272 and CAHS89226). All the identified CAHS transcripts were upregulated, with the exception of CAHS68135, which was downregulated. In addition, one MAHS transcript was identified, and all three analyses found it to be upregulated in the heat-exposed tardigrades when compared to the controls kept at 5 • C. Likewise, the group 3 LEA gene RvLEAM, a mitochondrial protectant like MAHS, was found to be differentially expressed and

Table 2|
Differential gene expression analyses of RNA-seq data generated from control and heat exposed tardigrades with the number of upregulated and downregulated genes found by each of the three methods (DESeq2, EdgeR, and Limma).  upregulated. Furthermore, a Dsup isoform was found to be upregulated in the differential expression analyses.

Differential gene expression
Our results revealed that major transcriptional pertubation occurs in active specimens of R. varieornatus as a response to temperature increment. The latter is in clear contrast to the constitutive expression observed during cryptobiotic survival in both eu-and heterotardigrades (Hashimoto et al., 2016;Murai et al., 2021). Importantly, high temperatures seem to transform tardigrade physiology, with 35.6% of the transcriptome being differentially expressed following 24 h heat stress. Among the 9634 differentially expressed transcripts identified by all three methods used to perform differential expression analyses (i.e., DESeq2, EdgeR and Limma), we found 12 DEGs related to heat-shock proteins, seven CAHS, two SAHS, one MAHS, one LEA and one Dsup gene. Of these transcripts, all DEGs coding for tardigrade intrinsically disordered proteins were upregulated, except for one CAHS. In contrast, only two of the heat-shock related transcripts were upregulated (a small heat-shock protein and a HSP70), whereas the remaining 10 were downregulated. We hypothesize that the small heat-shock protein could bind to denatured proteins preventing aggregation and that the HSP70 isoform could assist in renaturing dysfunctional proteins at the expense of ATP (Willsie and Clegg, 2002;Jagla et al., 2018;Hibshman et al., 2020), thus counteracting the effect of heating on tardigrade proteins. The upregulation of Dsup indicates that this damage suppressor could have an important role in protecting DNA during heat-stress. Specifically, Dsup has been shown to protect chromosomal DNA by binding to nucleosomes thereby shielding the DNA from reactive hydroxyl radicals (Hashimoto et al., 2016;Chavez et al., 2019;Minguez-Toral et al., 2020). It is thus reasonable to expect that Dsup also protects R. varieornatus' DNA from hydroxyl radicals that form as a result of increased temperature. In the category of intrinsically disordered proteins that were upregulated in response to heating, we furthermore found one LEA isoform (i.e. RvLEAM). LEA proteins have been shown to protect cells acting as antioxidants and membrane and protein stabilizers under diverse situations of water stress (Tunnacliffe and Wise, 2007). RvLEAM and MAHS were both found to be upregulated following heat stress, and thus may serve a protective role against heating related stress in the mitochondria of R. varieornatus. Their putative role might indeed be related with the protection against reactive oxygen species, which are a byproduct of electron transport in the mitochondria.
The intrinsically unstructured SAHS1 and SAHS2 are secreted by eutardigrade cells, and thus localize to the extracellular space (Yamaguchi et al., 2012). Recent studies have revealed that SAHS1 has a β-barrel motif with two potential carboxyl group binding sites, each of them allowing the binding of either two monocarboxylic acids or one dicarboxylic acid (Fukuda et al., 2017). Thus, it is not unlikely that SAHS1 serves a protective role involving the binding of carboxylic acids in the extracellular space of eutardigrades. This shielding effect would protect molecules from oxidative stress, thereby preventing the generation of reactive oxygen species in the extracellular environment. Additionally, these carboxylic binding sites might be able to bind to fatty acids and as such, SAHS1 might have the capability of storing or transporting fatty acids. Hence, the SAHS1 protein might have a similar function to that of albumin found in mammals, which is also characterized by conserved intrinsically disordered regions (Litus et al., 2018). SAHS proteins such as SAHS1 might thus have a dual function in transporting and storing energy rich molecules, such as fatty acids, ensuring an energy reserve and at the same time protecting them against oxidative stress in the extracellular environment.
In general, our findings reveal that the increment of temperature induced major transcriptional changes in the heat-exposed R. varieornatus, which is indicated by the 9634 DEGs found to be in common across the three differential gene expression analyses performed in this study. Many previously identified stress related tardigrade genes were found to be upregulated in the heat exposed, active state of R. varieornatus. Hence, we conclude that such genes contribute to stress tolerance in general, and that they are not necessarily exclusive to the cryptobiotic response.

Gene ontology analysis
ATPase activity is among the GO molecular function categories ( Table 3) that were found to be enriched among the differentially expressed transcripts in heat-exposed versus non-exposed R. varieornatus. ATPase activity represents a common molecular function utilized by a wide array of enzymes that hydrolyze ATP, providing the necessary energy to catalyze chemical reactions and e.g. transport ions and other molecules against their electrochemical gradients. NADH dehydrogenase (ubiquinone) activity is another overrepresented category in our analysis. This enrichment could reflect a change in mitochondrial inner membrane respiratory chain complex I activity. The mitochondrial respiratory complex I accepts electrons from NADH + H + and transfers these to the ubiquinone, while simultaneously moving protons out of the mitochondrial matrix, thereby fueling the electrochemical gradient that drives ATP production (Whitehouse et al., 2019). Thus, both the overrepresentation of ATPase activity and NADH dehydrogenase (ubiquinone) activity indicate a change in the metabolism of the heat stressed tardigrades. Given the higher temperature, it seems likely that chemical energy is required at a larger scale in the heatexposed, as compared to the non-exposed, tardigrades to catalyze cellular processes (e.g., DNA and RNA synthesis, intracellular signaling, transmembrane transport, etc.). Future studies could provide important further insight into the metabolic changes that take place in tardigrade cells following temperature increments.
Two other overrepresented GO molecular function categories, protein heterodimerization activity and sequence-specific DNA binding proteins, have a common denominator: transcription factors. Heterodimerization, a result of two non-identical proteins associating through non-covalent interaction to form a complex, is common among transcription factors, receptors and some enzymes (Jordan and Devi, 1999;Funnell and Crossley, 2012;Singh and Jois, 2018). Moreover, sequencespecific DNA binding proteins may also represent transcription factors that regulate gene expression (Suter, 2020). As such, these two overrepresented categories are in agreement with a change in transcription as a result of heating.
The overrepresented gene ontologies show a pattern revolving around a change, not only in transcription, but also translation. Indeed, overrepresented GO categories also include structural constituents of ribosome and translation initiation factor activity. Therefore, our analysis further suggests that the translation of proteins is being regulated as a response to the heating stimulus. This further supports the observation that heat-exposed tardigrades change their metabolism, transcription and translation to match the increase in temperature and counteract its effects.
A key feature in the body plan of tardigrades is the presence of a chitinous cuticle (Baccetti and Rosati, 1971;Bussers and Jeuniaux, 1973;Greven et al., 2016). As chitin binding is also an overrepresented GO term (Table 3), we hypothesize that heat-exposed R. varieornatus may be responding to the unfavorable change in their environment by preparing for encystment and diapause. Specifically, the observed regulation in transcription of chitin binding proteins could indicate that the heat-exposed tardigrades are in the process of shedding their cuticle. Encystation in tardigrades involves formation of new cuticle layers that are shed and subsequently provide protection against unfavorable environmental conditions (Guidetti and Møbjerg, 2018;Møbjerg and Neves, 2021). Encystment is not normally encountered in strong cryptobionts, such as R. varieornatus. Nevertheless, cysts of a Ramazzottius species were reported by Murray (1907). Interestingly, cyst formation was also reported in the strong cryptobiont Echiniscoides sigismundi, a marine tidal heterotardigrade, following experiments on osmotic stress tolerance (Clausen et al., 2014). Thus, although encystment is not the default strategy for strong cryptobionts, it may represent an alternative strategy, when they are presented with a stress factor (e.g. high temperatures or fluctuating salinities) that does not by itself induce cryptobiosis. Clearly more comprehensive studies on R. varieornatus and other strong cryptobionts as well as species well-known to form cysts are required to clarify whether high temperatures induce encystment.
Our results show that R. varieornatus reacts to heat-stress with a change in transcription and translation and with an adjustment of metabolism, and, possibly, by preparing for encystment and subsequently diapause. Thermotolerance in R. varieornatus is thus based on a complex response involving the concertation of a plethora of cellular processes. Importantly, distinct tardigrade taxa and lineages exhibit unique molecular adaptations (Kamilari et al., 2019) and it thus remains to be seen, whether closer and more distantly related tardigrade species have a similar response to elevated temperatures. In order to have a more comprehensive knowledge of the molecular processes associated with thermotolerance, future studies on different evolutionary lineages and taxa will be required, as well as investigations into the influence of habitat. In order to explore possible climate-correlated thermotolerance profiles, latitudinal comparisons of thermotolerance in species and populations inhabiting different climate regions is needed.

Concluding remarks
Our findings indicate that heat-shock proteins as well as eutardigrade specific proteins are involved in the heat stress response of R. varieornatus. Yet, these proteins only represent a minute fraction of the thousands of genes found to be differentially expressed following heat stress. A more detailed inspection of the large number of differentially expressed transcripts found in our analyses would likely uncover other genes with importance for tardigrade stress responses. Understanding the role of tardigrade stress-related genes, and the cellular processes they are involved in, can potentially be of great biomedical and biotechnological interest. Indeed, the potential translational application of tardigrade specific proteins was recently shown in studies revealing that the DNA-associating damage suppressor (Dsup) protein renders protection against UV-C and hydroxyl radicals to Dsuptransfected human HEK293 cells (Hashimoto et al., 2016;Ricci et al., 2021). Moreover, tobacco plants that express the codon-optimized Dsup gene seem to have an enhanced tolerance towards UV and X-rays (Kirke et al., 2020), revealing a potential for tardigrade specific proteins in the engineering of crops that better tolerate the more stressful environmental conditions associated with global warming.

Data availability
The RNA-Seq data generated and used for differential expression analyses in this study are available at the European Nucleotide Archive under project PRJEB49649.

Code availability
The bioinformatic pipeline used to analyse the RNA-Seq datasets is available in GitHub.

Declaration of Competing Interest
The authors report no conflict of interest.