Article ID: pjab.100.023
Tardigrades are microscopic animals that are renowned for their capabilities of tolerating near-complete desiccation by entering an ametabolic state called anhydrobiosis. However, many species also show high tolerance against radiation in the active state as well, suggesting cross-tolerance via the anhydrobiosis mechanism. Previous studies utilized indirect DNA damaging agents to identify core components of the cross-tolerance machinery in species with high anhydrobiosis capacities. However, it was difficult to distinguish whether transcriptomic changes were specific to DNA damage or mutual with anhydrobiosis. To this end, we performed transcriptome analysis on bleomycin-exposed Hypsibius exemplaris. We observed induction of several tardigrade-specific gene families, including a previously identified novel anti-oxidative stress family, which may be a core component of the cross-tolerance mechanism. We also identified enrichment of the tryptophan metabolism pathway, for which metabolomic analysis suggested engagement of this pathway in stress tolerance. These results provide several candidates for the core component of cross-tolerance, as well as possible anhydrobiosis machinery.
Tardigrades are microorganisms capable of entering an ametabolic state known as anhydrobiosis when faced with environmental desiccation.1) In this state, tardigrades are able to tolerate extreme environments, ranging from low to high temperatures,2)-4) pressure,5) organic solvents,6) and high dosage of ultraviolet7),8) or radiation.9)-17) Current knowledge of the anhydrobiosis machinery consists of both mechanisms that are broadly conserved across eukaryotes and tardigrade-specific proteins.18) Protection from potential damage during desiccation has been the focus of anhydrobiosis studies, and anhydrobiosis is known to sometimes cause DNA double strand breaks in tardigrades19) and also in sleeping chironomids.20) One example is the nucleus-localizing protein Damage suppressor (Dsup), which has been identified in Ramazzottius varieornatus and H. exemplaris to have DNA protective capabilities,21)-25) and human cell lines expressing Dsup have increased X-ray tolerance. However, the persistence of DNA damage in recovered samples suggests that this protection alone is not sufficient.
Studies have observed high tolerance to gamma radiation in the active state in many species throughout the phylum Tardigrada, sometimes even higher than during anhydrobiosis, suggesting that these organisms have high DNA repair capabilities, possibly a cross-tolerance acquired through anhydrobiosis mechanisms.26) The core machinery of tardigrade DNA repair does not differ significantly from that of other eukaryotes,27) with the exception of the duplication of several genes in R. varieornatus, e.g. MRE11.28) Other DNA repair factors regulated during anhydrobiosis include the regulation of RAD51 in Milnesium cf. tardigradum.29) Transcriptome-based studies have identified both non-homologous end joining repair and homologous recombination repair to be up-regulated during anhydrobiosis.30),31) Other studies have also observed regulation of DNA damage repair pathways in specimens exposed to ionizing and non-ionizing radiation; gamma rays, and ultraviolet C.8),32),33) In addition, two novel gene families related to stress response have been identified: Anhydrobiosis-related Mn (manganese)-dependent Peroxidase (AMNP), a novel Mn \(^{2+}\) -dependent peroxidase widely conserved throughout the phylum Tardigrada, and the Rv.g241 gene family, a potentially horizontally transferred gene with similarity to the C-terminus of bacterial peroxidases.31),33) However, the mechanism of anhydrobiosis is a complex union of various pathways, and DNA damage caused in the conditions used above was mostly based on the indirect attack of reactive oxygen species. This makes it difficult to classify whether the differentially expressed genes were regulated as a response to DNA damage or oxidative stress. Therefore, a more focused approach was required to identify key factors in the tardigrade DNA repair process.
To this end, we performed a transcriptomic analysis in Hypsibius exemplaris exposed to bleomycin, a commonly used compound that induces DNA double strand breaks. H. exemplaris is a relatively weak anhydrobiote that requires 48 h of preconditioning at high relative humidity for successful anhydrobiosis.34) We have previously assembled the genome of this species to identify key factors for anhydrobiosis,31) and genetic toolkits (e.g., RNAi) have also been developed for this species.35) In contrast, the stronger anhydrobiosis species R. varieornatus does not require preconditioning, making it difficult to distinguish core components of cross-tolerance using comparative transcriptomics. H. exemplaris shows tolerance to gamma rays in both embryonic and adult active stages at extremely high doses (3,000–4,000 Gy of \(^{137}\) Cs gamma ray),12) suggesting the existence of a prominent DNA protection, response, or repair machinery. UV light functions via electron excitation and low-LET ionizing radiation such as gamma rays function by both excitation and ionization, which may have multifactorial indirect effects on DNA, proteins, membranes, and other small molecules. In this study, we focused on the cross-tolerance specifically on DNA damage, through the use of DNA damage inducing drug bleomycin. An activated bleomycin molecule can decompose, producing reactive oxygen species, which can attack cellular molecules or bind to DNA through intercalation (or interaction with the minor grove of the double stranded DNA molecule) and cause cleavage.37) By exposing H. exemplaris to 100 μM bleomycin, we mimicked exposure to a high dosage of radiation, and we tested for genes that were induced in desiccation response as well. We observed a three-phased response, of which the early response had the most commonalities with anhydrobiosis entry through extensive time-course transcriptome sequencing. We found that several early responders that were observed in our previous analysis33) (e.g., CAHS, AMNP, Rv.g241) were included in the bleomycin early response, highlighting the importance of these gene families. We observed enrichment of the tryptophan metabolism pathway in the early response, but the metabolomic analysis was not consistent with the transcriptomic response. Taken together, we propose that the gene set identified here may comprise the core gene set required for cross-tolerance in H. exemplaris.
In order to analyze the transcriptomic response against exposure to bleomycin-induced DNA damage, we exposed H. exemplaris to 100 μM bleomycin for 24 h and collected specimens at 0, 1, 3, 6, 9, 12, and 24 h after bleomycin removal (Fig. 1A). A control sample without bleomycin exposure was also collected at the 0–h time point. These samples were subjected to RNA-Seq, where we generated 6–11M RNA-Seq reads (average 8.7M reads) for each sample, including an unexposed control sample, all in triplicate. Approximately 66–85% (average 79%) mapped to coding sequences (Supplementary Table S1). Approximately 70% of the transcripts had an average TPM greater than 1, consistent with previous studies. Initial clustering based on PCA analysis indicated that the transcriptome changed dynamically between 0–1 h and 12–24 h after bleomycin treatment (Fig. 1B). We detected 11,210 transcripts that were differentially expressed between any of the 24 time points. This is approximately three times more than our previous transcriptome analysis on the anhydrobiosis entry of H. exemplaris (1,422 genes significantly upregulated),31) suggesting a dynamic shift in the cellular state during recovery from 24–h bleomycin exposure.
(Color online) Transcriptome analysis of bleomycin-exposed H. exemplaris. [A] Experimental setup of the transcriptome analysis. H. exemplaris specimens were exposed to 100 μM bleomycin (BLM) for 24 h and were sampled at the corresponding time points. A control sample without bleomycin was also prepared. [B] PCA of transcriptome profiles. Transcripts with average expression \(>\) 1 were subjected to PCA. [C] Clustering of expression profiles. Transcripts with average expression \(>\) 1 were clustered based on Spearman’s correlation and split into eight groups based on total within sum of squares. A dynamic change was observed between 0–1 h, 9–12 h, and 12–24 h. [D] An illustration of the broader clustering of the eight groups. E, M, and L indicate the Early, Middle, and Late responses.
To obtain an overview of the transcriptional response, we clustered the transcriptome profiles of differentially expressed genes (DEGs), resulting in eight groups (Fig. 1C). Gene Ontology enrichment analysis of each group (Supplementary Table S2) showed that terms related to DNA repair were enriched in Groups 4 and 7, of which only Group 4 contained terms related to double strand break repair, in particular non-homologous end joining (NHEJ). This group also included mismatch repair (MMR) as well. In contrast, group 7 included nucleotide-excision repair (NER).
These eight clusters could be grouped into three larger clusters depending on the time point when they showed the highest expression, designated Clusters Early (E; 0–3 h, Groups 2, 3, 8), Middle (M; 4–9 h, Groups 4, 5), and Late (L; 12–24 h, Groups 1, 6, 7) (Fig. 1D). The majority of Early group transcripts did not show a significant increase in expression compared with the non-exposed control samples (5316/5714). In comparison, the majority of the Middle 1013/1780) and Late (3322/4218) response transcripts had significant increases at 3–9 h and 12–24 h, respectively. Gene ontology enrichment analysis of these significantly upregulated groups (early, middle, and late induced transcripts) suggested a three-step regulation of cellular processes (Supplementary Fig. S1A). First, induction of metabolic processes focusing on the tryptophan pathway occurs. These processes are followed by mitochondrial regulation, indicating recovery from mitochondrial stress and double strand break repair by NHEJ (3–12 h). Finally, DNA replication and metabolism are regulated in the late stages to restore homeostasis.
We also validated the expression of previously identified tardigrade-specific anhydrobiosis-related genes, such as the CAHS, SAHS, MAHS, LEAM, Dsup, TDR1, AMNP, and Rv.g241 gene families (Supplementary Table S3).22),33),38)-40) In contrast to previous studies using stresses other than anhydrobiosis,31),33),41) most CAHS and SAHS copies showed a significant increase in expression. On the contrary, the expression of MAHS (BV898_03788) and Dsup (BV898_01301) orthologs was significantly decreased during the time course. The LEAM ortholog (BV898_03429) did not show any significant change. The tardigrade-specific peroxide AMNP family and the stress-responsive Rv.g241 gene family also showed an increase in expression (Supplementary Fig. S1BC). We observed extremely high expression of TDR1 at 0 h (TPM 328), of which approximately TPM 100 is observed until 12 h.
2.2.Comparison with anhydrobiosis entry and identification of novel stress-responsive transcripts.We then compared the bleomycin response with anhydrobiosis entry, namely, transcriptome sequencing data of H. exemplaris anhydrobiosis (GSE94295).31) Re-analysis of the anhydrobiosis transcriptome data revealed 2,738 DEGs, of which 2,089 transcripts were significantly upregulated. Anhydrobiosis entry shared the largest number of transcripts with the Early response cluster (1,152), followed by the Late cluster (692) and then the Middle cluster (449) (Fig. 2A). Fisher’s exact test for enrichment indicated the Early response cluster to be most significant (Early: 1.9E-70, Middle: 2.8E-47, Late: 6.1E-12, FDR corrected), suggesting that the Early response may partially resemble the anhydrobiosis entry trajectory.
(Color online) Comparison with anhydrobiosis entry. [A] UpSet plot showing the number of genes shared between each cluster and anhydrobiosis entry. [B] Number of GO BP terms significantly enriched between each response cluster and anhydrobiosis. [C] Features of the 102 NLS predicted transcripts without annotations, induced in anhydrobiosis and either one of clusters A, B, C. Transcripts are ordered based on the clustering preformed at Supplement Fig. S1C. Z-scaled expression profiles (Expression profile), DEG profiles (DEG), averaged IUPRED3 score for protein sequence (IUPRED, estimates disordered regions within protein sequences), conservation patterns across Tardigrada and relative lineages are shown as “Conservation” (species names omitted for visualization). Specific species names can be found in the Supplementary Table S6. The red line for IUPRED indicates the average of the group. A protein sequence IUPRED3 score \(>\) 0.5 typically infer strong disordered nature of the protein.
To obtain an overview of the transcripts common to anhydrobiosis entry and each response group, we subjected each subgroup to Gene Ontology enrichment analysis (Supplementary Table S4, Supplementary Fig. S2ABC). We found that most of the enriched terms were time specific (37, 16, and 31 terms for Early, Middle, and Late response, respectively, Fig. 2B), but the terms “proteolysis” and “oxidation-reduction process” were common to all subgroups, emphasizing the importance of antioxidative stress processes. We then focused on the 37 terms common to the Early response and anhydrobiosis and found many terms related to transmembrane transport and tryptophan metabolism to be contained in this list.
We next screened the DEGs to attempt to identify novel stress-responsive transcripts without SwissProt orthologs. First, we aimed to obtain Early responding gene families, induced in anhydrobiosis entry and bleomycin early-response that have orthologs induced in R. varieornatus UVC exposure early-response (Supplementary Table S5). Using these thresholds, we identified 9 tardigrade-specific gene families, including two CAHS families (OG0001789, OG0002138), the AMNP and Rv.g241 gene families (OG0000230, OG0000231). For the remaining gene families, protein structures were predicted with AlphaFold2 using ColabFold and similarity search with DALI against the AlphaFold database (Supplementary Table S5). We found protein families with similarity (Identity \(>\) 30 or Z-score \(>\) 10) to metalloendopeptidases (OG0000456), hormone receptors (OG0000592), abnormal cell migration proteins (OG0000940), and LYSM domain proteins (OG0028858). Only one protein family could not be annotated (OG0010176). Additional screening with fold change (Max TPM / Control) \(>\) 2 and high expression (Max TPM \(>\) 100) left only the CAHS, Rv.g12777, and Rv.g241 gene families, again, highlighting the importance of these gene families. The TDR1 gene did not have a significant increase in expression in our anhydrobiosis dataset (TPM 0.9 to 6.19, 6.87x increase, FDR \(=\) 0.36).
Hypothesizing that novel DNA protection- or repair-related proteins would localize to the nucleus, we also tested for nucleus-localizing stress-responsive proteins: tardigrade-specific proteins with predicted nuclear localization signals (NLSs) and induced by both anhydrobiosis and bleomycin (Fig. 2C, Supplementary Table S6). Approximately 100 transcripts fitted these thresholds. Several of these transcripts were highly duplicated in Eutardigrada and were also conserved in Heterotardigrada. By querying these AlphaFold2 modeled protein structures against the AlphaFold2 SwissProt database using FastSeek, several of the proteins were annotated. We found two proteins that we identified as related genes to be induced from the transcriptome analysis (Fig. 3AB); 5-hydroxytryptamine receptor 1E (BV898_15649.t01, AF-Q6VB83-F1-model_v4.pdb E-value = 6.22E-12, RMSD between 163 pruned atom pairs is 1.124 angstroms across all 353 pairs: 10.861) and Ecdysone-induced protein 78C (BV898_05254.t01, AF-P45447-F1-model_v4 E-value 4.48E-08, RMSD between 84 pruned atom pairs is 1.044 angstroms; across all 391 pairs: 32.283). In addition, we found 18 transcripts with high IUPRED3 scores above 0.9 (four with 1.0), indicating highly disordered proteins (Fig. 3CDEF). Only one transcript (BV898_07132.t01) showed homology to a Swiss-Prot gene (28 kDa heat- and acid-stable phosphoprotein). The Dsup gene (BV898_01301) also fit into this category, thus suggesting that these four may be novel nuclear protection/repair proteins. Additionally, we screened the 102 genes for those with known DNA binding motifs and found only two: BV898_13251 and BV898_10351. The BV898_13251 gene was found to harbor a “Helix-loop-helix DNA-binding domain”, as present in TFEB transcription factors, and BV898_10351 was expressed at low levels ( \(<\) TPM 10). We concluded that we were unable to identify potential DNA-binding protection related genes through our current dataset.
(Color online) AlphaFold structures of NLS-containing proteins. AlphaFold2 structures modeled from NLS-predicted proteins that had no SwissProt BLASTP hits (E-value \(<\) 1e-15). [AB] Structures for those that had a FoldSeek hit against the AlphaFold2 SwissProt database (E-value \(<\) 1E-05). The yellow and blue structures indicate the reference and the predicted H. exemplaris protein structures, respectively. [A] BV898_05254.p01 (5-hydroxytryptamine receptor 1E (AF-Q6VB83-F1-model_v4.pdb E-value \(=\) 6.22E-12, RMSD between 163 pruned atom pairs is 1.124 angstroms across all 353 pairs: 10.861). [B] Ecdysone-induced protein 78C (AF-P45447-F1-model_v4 E-value 4.48E-08, RMSD between 84 pruned atom pairs is 1.044 angstroms; across all 391 pairs: 32.283). [C–F] Transcripts that had an IUPRED3 score of 1.0. [C] BV898_03773.p01. [D] BV898_07132.p01. [E] BV898_12866.p01. [F] BV898_18010.p01.
To validate whether tryptophan metabolites are accumulated during anhydrobiosis entry and bleomycin response, mirroring the DEG analysis, we subjected active, anhydrobiotic, and bleomycin-exposed H. exemplaris specimens to metabolomic analysis using liquid chromatography-mass spectrometry (LC-MS) and capillary ion chromatography-mass spectrometry (IC-MS). We detected a total of 237 compounds (Supplementary Table S7), with relatively high Spearman correlations between replicates. PCA analysis of the metabolite profiles showed high consistency between control (active) and bleomycin-exposed (Bleomycin) samples, but with some variance in anhydrobiosis samples (Anhydrobiosis) (Supplementary Fig. S3A). From these metabolite profiles, we detected a total of 113 compounds with significant differences between conditions (pairwise T-test, FDR \(<\) 0.05); 92, 49, and 95 compounds for Control vs Anhydrobiosis (Fig. 4, Supplementary Fig. S3B), Control vs Bleomycin, and Anhydrobiosis vs Bleomycin, respectively. Among the compounds that showed higher concentration in Anhydrobiosis, or Bleomycin compared with Control, eight were common to both conditions, 17 and 32 were Bleomycin and Anhydrobiosis specific. These eight compounds were composed of 3-methylhistidine, 3-phenyllactate, ADP-ribose, CDP-choline, deamido-NAD \(+\) , lysine, riboflavin, and S7P. KEGG pathway enrichment analysis of each condition combination (e.g., Control vs Anhydrobiosis, Control vs Bleomycin, Anhydrobiosis vs Bleomycin), indicated enrichment of several KEGG pathways, of which only Control vs Bleomycin showed enrichment in tryptophan metabolism (FDR \(=\) 0.002264990). Mapping of transcriptome and metabolome data to the KEGG tryptophan metabolism map (ko00380) revealed accumulation of tryptophan, 5-hydroxy-L-tryptophan, L-kynurenine, and 3-hydroxy-anthranilate compounds and the regulation of enzymes contributing to this pathway (Fig. 5). Notably, we did not see accumulation of these compounds in anhydrobiosis samples, although the corresponding enzymes seem to be up-regulated.
(Color online) Metabolomic profiles of active, anhydrobiotic, and bleomycin-exposed H. exemplaris. Z-scaled concentration profiles for 113 metabolites with significant differences between Control-Anhydrobiosis, Control-Bleomycin, Anhydrobiosis-Bleomycin.
(Color online) KEGG pathway map of tryptophan metabolism. Data for two comparisons are indicated in a single circle (metabolite) or square (transcript); For metabolites, Left: Control vs Bleomycin; Right : Control vs Anhydrobiosis; For transcripts; Left: Control vs Bleomycin-0 h (24 h exposure, no recovery time), Right: Active vs Tun. Colors correspond to the following conditions: Gray: not detectable or not detected (0 μg/ μL for each condition), Blue: decrease in expression/concentration (lower in Bleomycin/Anhydrobiosis); White: detected, but no significant difference; Red: increase in expression/concentration (higher in Bleomycin/Anhydrobiosis).
Tardigrades in general have been observed to have high radiation tolerance, which has been hypothesized to be a byproduct of the acquisition of the anhydrobiosis mechanism. We previously performed transcriptome analysis of R. varieornatus exposed to gamma rays and UVC to identify core components contributing to this cross-tolerance.33) Building on this, we tested for such factors that could be induced using a chemical approach, in this study using bleomycin. Bleomycin causes both direct and indirect DNA damage, making it comparable to previous radiation-based methods that mainly cause indirect damage through reactive oxygen species.
In this study we used bleomycin at a concentration of 100 μM. In Chinese hamster DC-3F fibroblasts, this concentration is approximately equivalent to 1,500,000 Gy (1 nM to 15 Gy).42) However, considering that (1) bleomycin itself does not diffuse through the plasma membranes, (2) we used individual tardigrades, which is quite different from exposure to cell cultures, (3) tardigrades have a cuticular exoskeleton that may further impede cellular uptake, 100 μM would not directly correspond to this extremely high gamma ray dose. A concentration of 50 μM has been used in Caenorhabditis elegans,43) 10–50 μg/mL in Drosophila melanogaster (corresponding to 7–35 μM for 1415.55 g/mol bleomycin),44),45) so considering the higher tolerance to DNA damage observed in H. exemplaris, the concentration used in our study is neither extremely high nor low. A previous study used an exposure of 100 μM for four or five days against H. exemplaris and analyzed the transcriptomic response at the point of exposure end.40) Therefore, in order to induce DNA damage, but achieve near 100% survival, H. exemplaris specimens were exposed to 100 μM bleomycin for 24 h and were sampled over the time course of 24 h after removal of tardigrades from bleomycin solution.
We generated time-course transcriptome sequencing data on H. exemplaris specimens recovering from a 24–h exposure to 100 μM bleomycin. We observed reduction in movement during bleomycin exposure, similar to that observed in UVC-exposed R. varieornatus.33),46). We observed a dynamic shift in the transcriptome at 0–1 h and 9–24 h after bleomycin clearance, suggesting the existence of Early, Middle, and Late response phases. Although this three-step response was similar to that observed in R. varieornatus upon irradiation, the timing of the transcriptome shifts was different; 4–6 h and 24–36 h for UVC, 6–9 h in gamma-ray response.32),33) This may be due to several factors; the early response may be induced during 24–h bleomycin exposure, differences in the requirement of stress response factors between R. varieornatus and H. exemplaris, differences in the mechanism by which gamma-ray, UVC, and bleomycin cause DNA damage.
Clustering of the transcriptome profiles and detection of enriched gene ontologies revealed that MMR and NHEJ in the middle response (Group 4), and NER are induced in the Late response. This suggests that H. exemplaris may repair devastating double strand breaks in 1–3 h after drug clearance by NHEJ and MMR. Notably, the Ku70/Ku80 complex genes were already upregulated at the 0 h time point compared with the control samples, possibly indicating that DNA damage had already been detected, but the repair was not yet induced. The fact that cell replication is rare in active H. exemplaris,47) with the exception of the digestive tract, and that they are equipped with specialized protection mechanisms such as Dsup, suggests that double strand breaks may not be a priority in damage repair, but rather a focus on minor DNA damage that may cause disruption of gene expression. Further analysis in the Early, Middle, and Late responses showed that DNA replication and central carbon metabolism are regulated in the Late response, thus suggesting that DNA repair and cell replication may occur at least 12–24 h after exposure to bleomycin. The NER pathway induced during this period may be contributing to the repair of remaining single nucleotide level damage, possibly those that emerged during the error prone NHEJ or during DNA replication. In addition, we found regulation of mitochondria-related genes during the middle response, suggesting mitochondrial reorganization at this stage. An ultrastructural study in H. exemplaris anhydrobiosis has observed abnormal morphology in the mitochondria.48) In addition, regulation of the mitochondrial chaperone BCS1 has also been observed,31) suggesting the presence of mitochondrial stress.
We also validated the expression of previously identified anhydrobiosis genes and observed the induction of CAHS, AMNP, and Rv.g241 orthologs. These gene families were also found to be shared between anhydrobiosis, bleomycin recovery, and UVC response in R. varieornatus. In addition, a single AMNP ortholog (BV898_01090) and Rv.g241 ortholog (BV898_18007) were also found to be significantly induced in the proteomic data of D942 exposed H. exemplaris,49) suggesting regulation of this ortholog through AMPK50) and highlighting the importance of these two gene families in the stress response in H. exemplaris and R. varieornatus. AMNP and Rv.g241 gene families were identified as stress responding families through transcriptome analysis in UVC-exposed R. varieornatus.33) The AMNP gene family was identified as a Mn \(^{2+}\) -dependent peroxidase, while the function of Rv.g241 remains to be identified, but considering that it is conserved at the C-terminus of bacterial peroxidases, it may be some kind of cofactor that enhances the antioxidative functions of peroxidases. Surprisingly, we found that many CAHS orthologs were also induced. CAHS proteins have been found to form a high-order structures in a desiccated state, which has been hypothesized to protect cellular molecules.51) Such high-order structures may not be required during recovery from damage; therefore, CAHS proteins may have additional stress-responsive functionalities, presumably related to its disordered nature with highly abundant charged residues. A recent study employing transcriptome sequencing on bleomycin-exposed H. exemplaris identified a DNA damage response gene TDR1 (BV898_14257).40) Although the TDR1 gene was highly induced in our bleomycin dataset similar to the previous study, we did not observe induction in the anhydrobiosis course. This may imply that there is limited DNA damage during anhydrobiosis, or the induction of TDR1 requires a prolonged recovery period. We hypothesized that the R. varieornatus TDR1 ortholog may be induced in our previous gamma ray exposure transcriptome dataset32); however, we did not find a TDR1 ortholog in other accessible tardigrade genomes using BLASTp using the truncated BV898_14257 amino acid sequence (E-value \(<\) 1E-5).33) Synteny-based analysis suggested g12065.t1 as a candidate ortholog with no known functional domain annotation; however, the sequence similarity was extremely low, and AlphaFold2 prediction suggested opposite alpha-/beta- propensities. The expression of Rv.g12065 was moderate (approx. TPM 50) for the gamma ray exposure time course. These data may imply that this gene may be a Hypsibius-specific gene, or an ortholog pair with extremely low similarity, as observed between HeDsup and RvDsup. We have also identified approximately 100 novel genes that encode a nuclear localization signal suggesting a nucleolus-localizing protein. Many of these proteins are highly disordered, a feature similar to tardigrade-specific anhydrobiosis genes.22),38),39). These proteins may be candidates for future functional analyses.
We then compared the transcriptome response to anhydrobiosis entry. We found that transcripts regulated during anhydrobiosis entry shared the most transcripts with the Early response (42%), followed by the Late (25%) and Middle (16%) responses. Gene Ontology enrichment analysis identified 37 gene ontology terms to be enriched in the early response-anhydrobiosis gene set. This list included many terms related to tryptophan metabolism and the oxidation-reduction process. Tryptophan metabolism has recently emerged as an important pathway in human therapeutics, immunology, and neuronal function.36) A metabolome study in Polypedilum vanderplanki, an anhydrobiotic insect, has identified the tryptophan oxidation pathway as one of the most anhydrobiosis-responsive pathways, e.g., kynurenic acid and xanthurenic acids showed high accumulation during anhydrobiosis recovery.52) In our metabolomic analysis, we observed accumulation of tryptophan metabolites in only in bleomycin-exposed samples, not during anhydrobiosis entry. This was not consistent with the gene expression profiles; most enzymes in the tryptophan metabolism pathway are induced in both bleomycin exposure and anhydrobiosis entry. Comparison of the tryptophan metabolite profiles with P. vanderplanki indicated similar profiles; constant 3-hydroxy-kynurenine concentration, with no significant difference in expression of kynurenine aminotransferase (K00816) during bleomycin exposure (down-regulated during anhydrobiosis entry). We did not detect xanthurenic acid or kynurenic acid in any samples. These data suggest that while tryptophan metabolism may be an important factor in tardigrades and P. vanderplanki, there are differences in how these pathways contribute to anhydrobiosis itself.
In this study, we performed time-course transcriptome analysis of bleomycin-exposed H. exemplaris to identify factors that constitute the core element of the cross-tolerance machinery. Through comparative transcriptomics, we propose that previously identified stress responders (i.e., tardigrade-specific peroxidase family) as well as the tryptophan metabolism are such elements. We attempted to validate this through metabolomic analysis, suggesting different profiles between anhydrobiosis and chemical stress. Taken together, we propose that antioxidative functions through tryptophan metabolism may play a role in the tardigrade stress response, though to what extent remains to be clarified.
H. exemplaris Z151 strain was cultured in accordance with our previous studies.31),34) Briefly, specimens were placed on 2% agar gels using Volvic water, supplied with Volvic water containing Chlorella vulgaris. Individuals were collected and moved to a new agar plate every 7 days.
4.2. Transcriptome sequencing of bleomycin exposed individuals.Fifty individuals per sample were placed in 96-well flat bottom plates and exposed to 200 μL of 100 μM bleomycin (Funakoshi, Cat No. B4518, dissolved in Volvic water) for 24 hours (N \(=\) 5). A control series exposed to Volvic water was also set. After 24 hours of bleomycin exposure, the control specimens and T0 specimens were collected and placed in a microtube with the least amount of water possible and frozen at \(-\) 80℃. The remaining samples were washed twice with 300 μL Volvic water to remove bleomycin and were placed on 45 mm agar plates with 2 mL Volvic water supplemented with 10 μL Chlorella for feeding. These specimens were incubated at 18℃ until sampling. After 1,3, 6, 9, 12, and 24 hours, each sample was collected as described above and frozen at \(-\) 80℃.
Frozen samples were submitted for transcriptome sequencing (N \(=\) 3). Total RNA was extracted using a Direct-Zol Micro Kit and quantified with Qubit RNA HS. An Illumina sequencing library was constructed from 30 ng of total RNA using NEB Next Ultra RNA in accordance with the manufacturer’s protocol. The final library was sequenced on a NextSeq 500 instrument (Illumina). Multiplexed reads were demultiplexed with bcl2fastq (Illumina). The quality of RNA-Seq reads was validated using FastQC v0.11.3.53)
4.3. Informatics analysis.The genomes, coding, and amino acid sequences and annotations for H. exemplaris nHd3.1.5 and R. varieornatus v101 were downloaded from http://ensembl.tardigrades.org.31) In addition, we predicted nuclear localizing proteins using NLStradumus v1.8.54) Sequenced reads were mapped to CDS and quantified with RSEM v1.2.2655) and bowtie2 v2.2.956) using the Trinity align_and_estimate_abundance.pl utility.57) Differential expression was tested with DESeq2 v1.22.258) using Trinity run_DE_analysis.pl, and transcripts with adjusted \(p\) -values below 0.05 and fold change over 1.5 or below 2/3 were designated as differentially expressed genes. The transcriptome profile was clustered by the Ward method using Spearman correlation. Each DEG cluster was subjected to Gene Ontology enrichment analysis with GOstat v2.48.059) using InterProScan Gene ontology annotations.60)
To determine gene families that are commonly regulated in multiple stresses, we used the transcriptome sequencing data of H. exemplaris anhydrobiosis (GSE94295)31) as well as the time-course RNA-Seq data of UVC-exposed R. varieornatus (GSE152753).31) The transcriptome sequencing data of H. exemplaris were subjected to differential expression analysis using the same method described above. DEGs of R. varieornatus UVC exposure were obtained from our previous study.33) To identify ortholog families, the amino acid sequences of H. exemplaris, R. varieornatus, and other tardigrade genomes from our previous analysis33) were clustered using OrthoFinder v2.5.1.61) Ortholog clusters containing DEGs of H. exemplaris bleomycin exposure, anhydrobiosis entry, and R. varieornatus UVC exposure were detected, and the annotations were validated. For gene families that were not previously annotated, we submitted the amino acid sequences to ColabFold (Alphafold_batch https://colab.research.[1]google.com/github/sokrypton/ColabFold/blob/main/batch/AlphaFold2_batch.ipynb)62) and predicted the protein structure with AlphaFold.63) Only one model was predicted for amino acid sequences capable of prediction. The predicted protein structures were submitted to DALI v5 search64) or FoldSeek v6.29e255765) against the AlphaFold Protein Structure Database.66) Homology were determined by (1) identity \(>\) 30 or (2) z-scale \(>\) 10 for DALI and E-value \(<\) 1e-5 for FoldSeek searches. Tardigrade predicted proteomes were screened for TDR1 orthologs using BLASTp search (E-value \(<\) 1e-5). Additionally, we conducted synteny searches using JCVI McScan between R. varieornatus and H. exemplaris.67)
4.4. Metabolome analysis.Hydrophilic metabolites were analyzed using two different analytical methods. Cationic metabolites were measured using LC-MS, as described previously.68),69) Briefly, LC-MS analysis was performed using an Agilent 1290 Infinity LC (Agilent Technologies, Santa Clara, CA) equipped with a Q Exactive Orbitrap mass spectrometer (Thermo Fisher Scientific, San Jose, CA). Separations were performed on a HILIC-Z column (150 \(\times\) 2.1 mm, 2.7 μm; Agilent Technologies). The mobile phase was composed of 20 mM ammonium formate \(+\) 0.25% ( \(v/v\) ) formate (A) and 20 mM ammonium formate \(+\) 0.25% ( \(v/v\) ) formate in 90% ( \(v/v\) ) acetonitrile (B). The flow rate was 0.25 mL/min, and the following linear gradient was used: 0–15 min, 100% to 70% B; 15–20 min, 70% to 10% B; 20–23 min, 10% B, followed by equilibration with 100% B for 7 min. The injection volume was 1 μL, and the column temperature was maintained at 40℃. The Q Exactive mass spectrometer was operated in a heated electrospray ionization (HESI) positive-ion mode using the following source parameters: spray voltage \(=\) 3.5 kV, capillary temperature \(=\) 250℃, sheath gas flow rate \(=\) 40 (arbitrary units), auxiliary gas flow rate \(=\) 10 (arbitrary units), auxiliary gas temperature \(=\) 300℃, sweep gas flow rate \(=\) 0, S-lens \(=\) 35 (arbitrary units). Data were acquired using a combination of full MS scan mode and parallel reaction monitoring (PRM) mode. The parameters in full MS scan mode were as follows: resolution, 35,000; auto gain control target, 3 \(\times\) 10 \(^{6}\) ; maximum ion injection time, 200 ms; scan range, 50–700 \(m/z\) . The parameters in PRM mode were as follows: resolution, 17,500; auto gain control target, 2 \(\times\) 10 \(^{5}\) ; maximum ion injection time, 100 ms; inclusion \(m/z\) list: 104.0706 (GABA), 118.0863 (Val), 132.1019 (Ile, Leu), 166.0863 (Phe), 182.0482 (methionine sulfone).
Anionic metabolites were analyzed by capillary IC-MS as previously described.70) Briefly, capillary IC-MS was performed using a Dionex ICS-5000 \(+\) system coupled with a Q Exactive Orbitrap mass spectrometer. Separation of anionic metabolites was performed using a Dionex IonPac AS11-HC-4 μm column (250 \(\times\) 0.4 mm, 4 μm; Thermo Fisher Scientific). The injection volume was 0.4 μL, and the column temperature was maintained at 35℃. The eluent flow rate was 20 μL/min and the potassium hydroxide (KOH) gradient was as follows: 0–2 min, 1 mM; 2–16 min, 1 mM to 20 mM; 16–35 min, 20 mM to 100 mM; 35–40 min, 100 mM; followed by equilibration with 1 mM KOH for 5 min. Isopropanol containing 0.1% acetate was delivered as the make-up solution at 5 μL/min. The Q Exactive mass spectrometer was operated in ESI negative-ion mode. ESI parameters were as follows: sheath gas, 20 (arbitrary units); auxiliary gas, 10 (arbitrary units); sweep gas, 0; spray voltage, 4.0 kV; capillary temperature, 300℃; auxiliary gas temperature \(=\) 300℃; S-lens, 50 (arbitrary units). Data were acquired in full MS scan mode and the parameters were as follows: resolution, 70,000; auto gain control target, 3 \(\times\) 10 \(^{6}\) ; maximum ion injection time, 100 ms; scan range, 70–1,000 \(m/z\) .
Quantified metabolites were further normalized by the total amount of compounds detected to the average of all samples. Metabolite profiles were subjected to pairwise T-test in R, and compounds with FDR \(<\) 0.05 were considered significant differences. Compound names were converted to KEGG IDs using MetaboAnalyst 5.0,71),72) and those that did not have a correct match were further searched for in the PubChem, ChEBI, or KEGG databases. The corresponding compounds might not be found in any of the databases or were a peak containing two compounds that were not be separated with our setup, and were therefore omitted from enrichment analysis: “1-Amino-1-cyclopentanecarboxylate \(+\) Pipecolate”, “3PG \(+\) 2PG”, “4-Acetylbutyrate”, “4-Hydroxymethylimidazole”, “4-Oxohexanoate”, “4-Oxopentanoate”, “AMP \(+\) 3 \('\) -AMP”, “Citramalate \(+\) 2-Hydroxyglutarate”, “Isonicotinamide \(+\) Nicotinamide”, “SDMA”. Linking data between KEGG compound IDs and pathway IDs were obtained from the KEGG REST service for all analyzed metabolites. For each KEGG pathway, we preformed Fisher’s exact test to detect enrichment using FDR-corrected \(p\) -values \(<\) 0.05 as the threshold. Metabolite and transcriptome profiles were visualized on KEGG pathway maps with Pathview.73)
4.5. Data analysis.Statistics and data manipulation were performed with R and custom scripts using the G-language Genome Analysis Environment v.1.9.1.74),75) Biovenn was used to visualize area -proportional Venn diagrams.76)
4.6. Data availability.RNA-Seq reads have been uploaded to GEO under the accession ID SE168917. The metabolome data are included in Supplementary Table S7. The raw data for metabolomic analysis is available upon request from the authors.
We thank Naoko Ishii, Ayako Shirahata, and Yuki Takai for the tardigrade experiments and transcriptome sequencing. We also thank Satsuki Ikeda for the metabolome analysis. C. vulgaris used to feed the tardigrades was provided courtesy of Chlorella Industry. This work is supported by KAKENHI Grant-in-Aid for Transformative Research Areas (A) from the Japan Society for the Promotion of Science (JSPS, grant Number 21H05279), Joint Research by Exploratory Research Center on Life and Living Systems (ExCELLS program Nos. 19–501 and 22EXC601) and partly by research funds from the Yamagata Prefectural Government and Tsuruoka City, Japan.
Supplementary materials are available at https://doi.org/10.2183/pjab.pjab.100.023.
Edited by Takao SEKIYA, M.J.A
Correspondence should be addressed to: K. Arakawa, Institute for Advanced Biosciences, Keio University, 403-1, Nipponkoku, Daihouji, Tsuruoka, Yamagata 997-0017, JAPAN (e-mail: gaou@sfc.keio.ac.jp).
Anhydrobiosis-related Mn (manganese)-dependent Peroxidase
CAHSCytoplasmic Abundant Heat Soluble
DEGsDifferentially expressed genes
DsupDamage suppressor
HESIheated electrospray ionization
IC-MSion chromatography-mass spectrometry
KOHpotassium hydroxide
LC-MSliquid chromatography-mass spectrometry
MAHSMitochondrial Abundant Heat Soluble
MMRmismatch repair
NERnucleotide-excision repair
NHEJnon-homologous end joining
NLSnuclear localization signal
PRMparallel reaction monitoring
SAHSSecretory Abundant Heat Soluble
TDR1Tardigrade DNA damage Response protein 1
UVCUltraviolet C