Skip to main content

Sexual dimorphism in the tardigrade Paramacrobiotus metropolitanus transcriptome



In gonochoristic animals, the sex determination pathway induces different morphological and behavioral features that can be observed between sexes, a condition known as sexual dimorphism. While many components of this sex differentiation cascade show high levels of diversity, factors such as the Doublesex-Mab-3-Related Transcription factor (DMRT) are widely conserved across animal taxa. Species of the phylum Tardigrada exhibit remarkable diversity in morphology and behavior between sexes, suggesting a pathway regulating this dimorphism. Despite the wealth of genomic and zoological knowledge accumulated in recent studies, the sexual differences in tardigrades genomes have not been identified. In the present study, we focused on the gonochoristic species Paramacrobiotus metropolitanus and employed omics analyses to unravel the molecular basis of sexual dimorphism.


Transcriptome analysis between sex-identified specimens revealed numerous differentially expressed genes, of which approximately 2,000 male-biased genes were focused on 29 non-male-specific genomic loci. From these regions, we identified two Macrobiotidae family specific DMRT paralogs, which were significantly upregulated in males and lacked sex specific splicing variants. Furthermore, phylogenetic analysis indicated all tardigrade genomes lack the doublesex ortholog, suggesting doublesex emerged after the divergence of Tardigrada. In contrast to sex-specific expression, no evidence of genomic differences between the sexes was found. We also identified several anhydrobiosis genes that exhibit sex-biased expression, suggesting a possible mechanism for protection of sex-specific tissues against extreme stress.


This study provides a comprehensive analysis for analyzing the genetic differences between sexes in tardigrades. The existence of male-biased, but not male-specific, genomic loci and identification of the family specific male-biased DMRT subfamily provides the foundation for understanding the sex determination cascade. In addition, sex-biased expression of several tardigrade-specific genes which are involved their stress tolerance suggests a potential role in protecting sex-specific tissue and gametes.


Reproductive modes in animals are typically categorized into two major categories: asexual and sexual. Sexually reproducing animals produce sex-specific gametes, and genetic exchange between sexes leads to higher genetic diversity [1, 2]. Gonochoristic animals usually show sexual dimorphism, not only in gametes, but also in somatic tissues, physiology, and behavior within a single species, demonstrating dynamic intraspecies differentiation.

Although many aspects of the mechanisms for inducing sex differences remain to be elucidated, they are usually regulated by sex-specific differences in the genome and gene expression [3]. Gonochoristic animals must undergo sex determination through common well-studied mechanisms that underlie the development of sex-specific organs. The physiological systems of sex determination vary among species, but are generally categorized into two types: determination by sex-linked chromosomes and environmental cues [4, 5]. Species possess sex chromosomes that show different karyotypes depending on their sex. In contrast, environmental cues, including temperature, nutritional status, and population density, act as initial cues for sex determination [6]. Regardless of the mode of sex determination, several widely conserved genes play crucial roles in sex-specific organ development.

The transcription factor family Doublesex and Mab-3 Related Transcription factor (DMRT) comprises key regulators of somatic tissue development in various animals [7]. In animals utilizing sex chromosomal sex determination systems, not only DMRT orthologs on the sex chromosome, but also on the autosomes, are involved in sex determination cascades to regulate the growth of sex-specific tissues [8, 9]. In contrast to the chromosomal sex determination system, environmental cues induce the development of sex-specific tissues in normally parthenogenetic individuals through the expression of DMRT orthologs (e.g., Dsx1 in Daphnia magna), leading to genetic exchange through mating [10,11,12].The interplay between highly diverse and conserved components in generating different sexes to overcome environmental and genetic challenges presents a significant challenge to the understanding of the sex determination cascade.

The phylum Tardigrada, a member of the Ecdysozoa with over 1,400 species [13], is divided into three classes: Heterotardigrada, Eutardigrada, and nomen dubium Mesotardigrada [14, 15]. Tardigrades are renowned for their ability to tolerate extreme environments, and studies have identified tardigrade-specific proteins that mediate tolerance against nearly complete desiccation and anhydrobiosis using species from the family Echiniscidae (Heterotardigrada) [16], Hypsibiidae (Eutardigrada) [17], and Macrobiotidae (Eutardigrada) [18]. Asexual (parthenogenesis) and sexual reproduction have been observed within this phylum, with reported instances of both gonochorism and hermaphroditism in sexually reproducing species [19]. Sexual dimorphism in morphology and behavior during mating have been observed [20, 21]. In contrast, we lack knowledge on the molecular mechanisms that induce sexual dimorphism, as most molecular and genomic studies have focused on parthenogenetic species [22].

To this end, we conducted genomic and transcriptomic comparisons between males and females of the model gonochoristic tardigrade Paramacrobiotus metropolitanus to identify the molecular factors related to sexual dimorphism. This species, which is rich in ecological information, has a reported 170 Mbp genome, is relatively easy to culture, and shows a male-biased sex ratio (Male:Female = 7:3), but morphological sexual dimorphism, excluding testis/ovary, has not been described [18, 21, 23,24,25]. The results in this study lay the foundation for subsequent studies aimed at identifying a master regulator of the sex determination cascade and sex-dependent genetic differences in tardigrades.


Tardigrade culture condition and specimen preparation

The tardigrade P. metropolitanus TYO strain was cultured following methods described in a previous report [24]. The specimens were sexed using the method described by Sugiura et al. [23]. The eggs of P. metropolitanus were individually placed in an agar-coated dish, and hatched individuals were separated and reared separately to avoid sex contamination. These specimens were then grown until the development of sexual organs that were used for sexing.

RNA extraction and sequencing

Total RNA was extracted as described by Arakawa et al. [26]. Two hundred and fifty specimens of each sex were placed in a 1.5 mL tube with minimal water, and 100 µL of TRIzol reagent was added (Thermo Fisher Scientific). Total mRNA was extracted using the Direct-zol RNA kit (Zymo Research), and the samples were transported to Chemical Dojin for sequencing. The transcriptome sequencing libraries were prepared with poly A selection using the NEBNext Ultra II RNA Library Prep Kit for Illumina (New England Biolabs) and were sequenced using the NovaSeq 6000 instrument (Illumina, 150 bp PE). Four and three replicates were prepared for males and females, respectively.

External data and annotation

Genome data for P. metropolitanus were obtained from our previous study [18]. Raw gDNA-Seq reads used to assemble the genome and RNA-Seq data for the hydrated and desiccated samples (2d) were downloaded from SRA with prefetch and fasterq-dump from the sra-toolkit suite v2.10.1 (, Accession ID: DRR144969, DRR146886). We have added additional annotations to the protein sequences using NCBI Conserved Domain Search [27], DeepLoc2 [28], or InterproScan v5.62–94.0 [29]. Tardigrade-specific anhydrobiosis genes were annotated based on previous studies [18, 30,31,32,33]. Nucleotide sequences for the coding regions were extracted using gffread v0.12.7 [34]. Protein structures were predicted by ColabFold v1.5.3 ( [35] with default settings and visualized ChimeraX v.1.7.0 [36]. The chromosome-level genome of Hypsibius exemplaris was downloaded from DNAzoo [37], and the positions of the gene predictions from our previous study [32] were converted to the new genome with LiftOff v1.6.3. Genome and gene predictions for Ramazzottius varieornatus were obtained from our previous report [32].

Gene expression analysis

Raw RNA-Seq reads were mapped to the coding sequences and quantified using RSEM v1.3.3 [38]. The raw counts were then subjected to statistical testing using DESeq2 v1.38.0, within the from the Trinity pipeline v2.15.1 [39, 40]. Transcripts with FDR values < 0.05 were identified as differentially expressed genes (DEGs). Gene Ontology Enrichment Analysis (GOEA) was performed based on InterProScan GO annotations using GOstats v2.68.0 and GSEABase v1.64.0 [41, 42]. Gene ontologies with p-values < 0.05 were considered significant. Singleton terms were removed from the final list of enriched terms.

To extract genomic regions enriched in transcripts biased to either sex, we performed enrichment analysis based on the number of DEGs with more than 10 × fold change within 200 kbp windows (100 kbp steps). Genomic bins were created with BEDtools v2.31.1 [43] and the number of genes fitting the criteria was calculated using BEDtools intersect. An in-house Rscript was used to perform Fisher’s exact test for each bin, and p-values were corrected by BH method. Regions with Q-value < 0.01 were considered enriched.

We also performed transcriptome assembly through Trinity v2.15.1 and StringTie v2.2.1 [40, 44]. The RNA-Seq data were mapped to the genome using Hisat2 v2.1.0 [45] and assembled with genome-guided Trinity or StringTie. A non-genome-dependent assembly was also produced with Trinity. The assembled information was passed into PASA v2.5.3 for variant detection and merged with the original gene prediction using EvidenceModeler v2.0.0 for a comprehensive gene prediction set [46, 47]. This gene set was also subjected to PASA expansion to identify additional splice variants. SAM file conversion was performed using SAMtools v1.16.1 [48].

Phylogenetic analysis

To identify and analyze the expression patterns of DMRT genes, we first performed an exhaustive search for genes harboring Doublesex-Mab-3 related domains (DM domains). Initial candidates were extracted based on the InterProScan searches performed above, and the corresponding amino acid sequences were submitted to a BLASTP v2.2.22 [49] search against P. metropolitanus. In addition, the amino acid sequences of P. metropolitanus DMRT orthologs were subjected to a BLASTP search against amino acid sequences predicted from various tardigrade genomes [33]. The amino acid sequences for the tardigrade DMRT orthologs, metazoan orthologs provided in a previous study [50], and velvet worm DMRT ortholog were pooled (Table S1) and then aligned with MAFFT v7.450 [51] and subjected to phylogenetic tree construction using IQTREE2 v2.2.2.6 [52]. The phylogenetic tree was visualized in FigTree v.1.4.3 ( The expression patterns of H. exemplaris and R. varieornatus DMRT orthologs during developmental stages were obtained from a previous report by our group [53]. Additional alignments for the DMRT proteins were performed by MAFFT v7.450 and visualized using MView (

Similarly, we conducted a phylogenetic analysis of CAHS genes. We obtained annotated CAHS sequences from our previous report [54] and pooled the amino acid sequences of P. metropolitanus CAHS candidate orthologs [18]. A phylogenetic tree was constructed using the same procedure.

Genome extraction

Virgin P. metropolitanus was prepared by the method described above, and a single tardigrade was placed in a 0.2 ml tube after 1% penicillin/streptomycin treatment for 2 h to remove contamination. Genomic DNA was extracted and prepared using the method described by Arakawa et al. [26]. An individual was crushed by pressing it against a tube wall using a pipette tip. Genomic DNA was extracted with Quick-gDNA MicroPrep kit (Zymo Research) with three freeze–thaw cycles and then following the manufacturer’s protocol. The extracted DNA was sheared to 550 bp target fragments with Covaris M220 and an Illumina library was constructed with a Thruplex DNA-Seq kit (Takara BioRubicon Genomics). Quantification, quality, and library size selection were performed with Qubit Fluorometer (Life Technologies) and TapeStation D1000 ScreenTape (Agilent Technologies), respectively. Sequencing library fragments in the range of 400–1,000 bp were cut and purified with a NucleoSpin Gel and PCR Clean-up kit (Clontech) and sequenced using a NextSeq500 sequencer with HighOutputMode 75 cycles kit (Illumina). The reads were de-multiplexed, and adaptor sequences were removed using the bcl2fastq v2 software (Illumina).

Genome reassembly

Previously published ONT raw reads were submitted for reassembly using Canu v2.2 [55], NextDenovo v2.5.2 [56], Shasta v0.11.1 [57], Flye v2.9.2-b1786 [58], redbean v2.5 (wtgbt2) [59], GoldRush v1.1.0 [60], SPADES v3.15.5 [61], Pecat v0.0.3 [62], and Raven v1.8.3 [63]. Polishing was performed using NextPolish v1.4.1 [64] for the NextDenovo assembly. Each assembly was evaluated using compleasm v0.2.2 (metazoa and eukaryota lineage) or BUSCO v5.5.0 [65, 66]. Completeness was also evaluated for H. exemplaris and R. varieornatus published genomes as well [32]. The coverage for 10 kbp bins was calculated as previously stated, where we used BWA-MEM2 v2.2.1 [67] instead of BWA-MEM. DMRT orthologs were searched with TBLASTN v2.2.22, using P. metropolinatus DMRT protein sequences as the query (E-value < 1e-50) [18]. Additionally, we co-assembled male and female short reads produced in this study, along with the ONT and Illumina datasets using SPADES v3.15.5.

Sex specific region analysis

To identify candidate sex chromosome regions, we employed the Y chromosome genome scan (YGS) method [68], which was previously used to identify Drosophila melanogaster sex chromosome contigs. Briefly, reads from the same sex were pooled, and 15-mers were extracted with jellyfish count v2.2.4 or v2.2.10 [69]. Scripts from the YGS method v.11b (8 Oct 2012 10AM) were then used to calculate the percentage of validated single-copy unique k-mers (P_VSC_UK) for each contig. This was performed for the previously published genome, as well as for the SPADES assembly performed above. We also tested the coverage for both sexes calculated from the gDNA-Seq data. The raw gDNA-Seq reads were mapped to the genome using BWA-mem2 v2.2.1 and converted into BAM files using SAMtools v1.17. The genome was split into 10 kbp bins and the average coverage for each bin was calculated using BEDtools v2.31.0. The values were then normalized by the median of all bins for that sample, and the average for males and females was computed.

For gene-level synteny analysis, we employed the Python version of McScan in the JCVI suite v1.2.7 [70]. Gene prediction and coding sequences were prepared for H. exemplaris and P. metropolitanus, and syntenic regions were identified and visualized using default settings [18, 22]. To identify the Dsup ortholog, candidates were identified using gene-level synteny. Disorderness was analyzed using DISOPRED ( and IUPRED3 ( and the protein structure predicted ColabFold v.1.5.3 [35, 71, 72].

Genotyping for male specific regions

Virgin specimens were replaced with single-worm lysis buffer (50 mM KCl, 10 mM Tris pH 8.2, 2.5 mM MgCl2, 0.45% NP-40, 0.45% Tween20, 0.01% gelatin, 2 μg of Proteinase K) [73]. The specimen was then dissolved by freeze–thaw cycles (three times for liquid N2 and RT) and incubated at 60 °C for 1.5 h and 95 °C for 25 min. Genotyping PCR was performed using the following conditions: 94 °C for 3 min; 40 cycles of 94 °C for 30 s, 50 °C for 30 s, and 68 °C for 1 min; and a final extension at 68 °C for 5 min. Primer sequences were designed using Primer3 [74] from the nucleotide sequences of scaffold Parri_scaffold0000295 (Table S2). Quick-Taq (TOYOBO) was used for the polymerase with the concentrations of each reagent, following the manufacturer’s instructions. Electrophoresis was performed at 100 V for 20 min with 1.2% agarose gel/TAE (NacalaiTesque), and then the gel was strained with ethidium bromide for 20 min. The DNA bands were visualized using ChemiDoc (BioRad).


Transcriptomics of P. metropolitanus sexes

To identify sex-specific gene expression and genomic loci, we produced 10–20 M reads of RNA-Seq data for male and female specimens (Table S3) that mapped approximately 80–90% of the genome. Based on these data, we quantified and conducted differential gene expression analysis. Principal Component Analysis (PCA) of the expression profiles indicated a clear distinction between the male and female samples (Fig. 1A). A total of 9,015 transcripts were differentially expressed, with 4,685 and 4,329 transcripts showing higher expression in females and males, respectively (Fig. 1B). Gene ontology enrichment analysis of each gene set indicated enrichment of various pathways (Figure S1). For females, we observed enrichment of RNA processing, cellular component biogenesis, and negative regulation of biological processes. In contrast, terms related to cyclic nucleotide biosynthetic processes, aminoglycan metabolic processes, and monatomic ion transport were enriched in males.

Fig. 1
figure 1

Transcriptomic analysis of both sex. A PCA of expression profiles. B Scatterplot of the expression profiles. Red dots indicate differentially expressed transcripts (FDR < 0.05)

The sex determination cascade comprises multiple genes, forming a signaling cascade that causes differentiation between the sexes. We first focused on DMRT, a well-conserved gene family that regulates sex-specific tissue development and behavior. Initial BLAST analysis identified five DMRT orthologs (PARRI_0009851, PARRI_0001169, PARRI_0005877, PARRI_0003090, and PARRI_0003093) [18]. We observed that three genes, PARRI_0003090, PARRI_0003090, and PARRI_0005877, were upregulated in males, whereas PARRI_0003090 was moderately expressed (TPM > 30) in males.

A diverse array of lineage-specific upstream signaling factors (e.g., tra2, nix, fem) induce sex-specific splicing variants of the doublesex gene, transmitting signals to the downstream sex development cascade [6]. Although the master regulator of sex determination is highly variable, several components of the cascade are highly conserved, such as the DMRT orthologs. Likewise, the transformer-2 (tra2) gene is a DNA-binding protein coupled with the Tra protein, causing sex-specific splicing of dsx in insects [6]. The P. metropolitanus Tra2 (PmTra2; PARRI_0000692) exhibited sex-biased variants (Figure S2). The female-biased variant of PmTra2 (PmTra2F, evm.model.Parri_scaffold0000002.194, Female: 77.45, Male:38.95) has an additional intron in the 5’ UTR compared to the male-biased variant (PmTra2M, evm.model.Parri_scaffold0000002.194.3.65434fff, Female: 1.07, Male: 3.36).

We also sought to identify genes participating in the sex cascade in Drosophila, e.g., fruitless (fru) and sexlethal (sxl) genes. BLAST searches identified two candidates for sxl orthologs (PARRI_0002227 and PARRI_0007430) [18], which showed contrasting expression profiles. However, a phylogenetic analysis of these genes could not determine whether these orthologs were sxl or a gene family with relatively high similarity; therefore, we cannot conclude whether these genes are sxl orthologs (data not shown). No hits were found for fru in P. metropolitanus nor any tardigrade genomes. Thus, we concluded that fru is missing and sxl remains questionable in P. metropolitanus.

Genomic loci of the male-biased genes

We detected a peculiar population of genes that were expressed approximately > 25 higher in males (Fig. 1B). Hypothesizing that these male-biased expressed transcripts may be sex-specific genes located on the sex chromosome, we conducted a genomic enrichment analysis to determine genomic loci enriched in these highly biased genes. Using a genomic bin of 200 kbp (corresponding to roughly 30 genes per bin) against differentially expressed transcripts that had over 10 × fold change than the other sex (Female: 674, Male: 1724), we detected 325 (29 scaffolds) and 12 (three scaffolds) bins for males and females, respectively (Fig. 2A). We noted that approximately 2% of the male-biased genes had more than TPM 10 in females (11% for male-expression of female-biased genes), thus implying the specificity of male-biased genes. Gene ontology enrichment analysis of genes located in these bins indicated a high enrichment of transcripts related to sperm function (Table S4, S5). Interestingly, two of the three male-induced DMRT paralogs (PARRI_0003090 and PARRI_0003093) were located within a bin enriched for male-biased genes on the scaffold Parri_scaffold0000005 (Fig. 2B) [18]. We also observed that the genes within and in the surrounding regions of these bins were also expressed in females, suggesting that these genomic loci may not be male-specific (Fig. 2C).

Fig. 2
figure 2

Multiple male-biased regions within the P. metropolitanus genome and their synteny. A Genome-wide enrichment analysis of male- or female-biased transcripts. Scaffolds were ordered by size and colored green and purple to visualize the scaffolds. The threshold of FDR < 0.01 was used. B, C Characteristics of scaffold Parri_scaffold0000005 harboring the DMRT paralogs. B FDR values from the genomic loci enrichment analysis plotted against the bin’s position. Blue and red indicate the values for males and females, respectively. C Expression fold-change log2(male + 0.1 / female + 0.1) of the genes plotted against their locations along the scaffold. Colors indicate whether the gene was differentially expressed. D Macro-scale synteny analysis to identify orthologous genomic loci in male-biased scaffolds. Gray lines indicate syntenic blocks between H. exemplaris and P. metropolitanus scaffolds. The synteny block highlighted in green indicates the location of DMRT paralog loci. The numbers on the bar indicate the chromosome number or scaffold ID for each genome assembly. E Synteny region of the DMRT paralog loci on P. metropolitanus scaffold Parri_scaffold0000005 and H. exemplaris Chromosome 1. Orthologs between H. exemplaris and P. metropolitanus determined by synteny analysis are connected by a gray ribbon. The corresponding region in H. exemplaris Chromosome 1 also aligned to the P. metropolitanus scaffold scaffold0000002 3.79–3.84 Mb, as indicated in between He Chr1 10.02–10.23 Mb and Pm scaffold000005 1.96–2.37 Mb. The Dmrt3090/3093 complex orthologs (PmDmrt3090 and PmDmrt3093) are also indicated

To evaluate whether these genomic loci enriched for male-biased genes were on the same chromosome, we preformed synteny analysis with the recently reported chromosome-level genome assembly of H. exemplaris. Paramacrobiotus metropolitanus has been observed to have a 2n = 10 karyotype, similar to that of H. exemplaris [23, 75]. While the queried 29 male-biased scaffolds did not focus on a particular chromosome, we observed a slight bias toward chromosomes 1, 2, and 5 (Fig. 2D). Furthermore, we observed that the region harboring the paralogous DMRT loci and the surrounding region on Parri_scaffold0000005 were missing in H. exemplaris, where a genomic region on a different scaffold (Parri_scaffold0000002) was inserted into in H. exemplaris (Fig. 2E). These data suggest that this genomic region may have emerged in the P. metropolitanus lineage.

Emergence of a novel dmrt93B-like subfamily specific to Macrobiotidae

Given the importance of paralogous DMRT genes located on Parri_scaffold0000005, we focused on the characterization of the orthologs to determine the characteristics of these paralogs.

First, we submitted the amino acid sequences of the P. metropolitanus DMRT family for phylogenetic analysis, incorporating various tardigrade DMRT orthologs from genome and transcriptome assemblies. Careful examination of the PARRI_0003093 gene structure revealed a misassembly of a single nucleotide insertion, identified through gDNA- and RNA-Seq read mapping. This caused a frameshift in the 3′ terminus, leading to a truncated coding sequence. Therefore, manual curation for this gene was performed, resulting in 463 amino acid sequences. Phylogenetic analysis identified PmDmrt99B (PARRI_0009851), PmDmrt93B (PARRI_0005877), and PmDmrt11E (PARRI_0001169) orthologs, as well as two Dmrt93B paralogs (PmDmrt3090 PARRI_0003090; PmDmrt3093 PARRI_0003093; the 3090/3093 complex). The 3090/3093 complex contained DMRT genes only from Macrobiotidae species, suggesting the acquisition of this subfamily in this lineage (Fig. 3A, Table S1). We also observed a phylum-wide loss of the doublesex subfamily. Furthermore, we observed an Echiniscidae-specific Dmrt93B subfamily that was not included in the 3090/3093 complex, while the relatively lower bootstrap support of this branch (88) complicated the phylogenetic position of this clade (Fig. 3A, Table S1). We only found the DM domain, but not the CUE-DMA domain in these subfamily members through InterProScan analysis. Interestingly, phylogenetic analysis indicated that the dsx-like gene of the velvet worm branched into a Doublesex clade with arthropods, suggesting that dsx emerged after the divergence of Tardigrada (Fig. 3A, Table S1). We did not detect two copies of the 3090/3093 complex in several other Macrobiotidae species. A direct comparison between PmDmrt3090 and PmDmrt3093 amino acid sequences indicated that the first 30–180 aa sequences were extremely similar, but the intron nucleotide sequences were completely different (Figure S4AB). Furthermore, multiple nanopore reads spanned the entire length of each gene. Together, we suggest that the two copies were not the result of misassembly of these loci. We also noted that no ONT reads spanned both PmDmrt3090 and PmDmrt3093. However, the 3090/3093 complex region spanned more than 30 kbp and the N50 length of the ONT data was approximately 17 kbp. It is possible that there are no ONT reads spanning the entire region. While reassembly of the ONT reads using more recent assembly methods produced a more contiguous assembly (NextDenovo + NextPolish; Table S6), these two genes were predicted to be two separate genes.

Fig. 3
figure 3

Phylogenetic analysis and the expression of DMRT orthologs. A Phylogenetic analysis of DMRT orthologs detected in the tardigrade genomes. The DMRT families were classified based on the orthologs of the model species. Bootstrap values of less than 90% are shown on the branch. The blue colored region indicates the Macrobiotidae-specific orthologs. The P. metropolitanus and Echiniscidae-specific Dmrt93B orthologs (Cornechiniscus lobatus, Echiniscus spp. Pseudoechiniscus sp.), and the partial dsx of velvet worm (Euperipatoides kanangrensis) are indicated in bold. B Expression of PmDmrt orthologs. Triangle points indicate differentially expressed genes and circles indicate non-significant changes. The gray line indicates x = y. C Multiple alignments of DM and CUE-DMA domains. Dm indicates D. melanogaster

Based on these annotations, we identified PmDmrt3090, PmDmrt3093, and PmDmrt93B to be significantly expressed in males; thus, all three induced copies belong to the Dmrt93B clade (Fig. 3B). We then utilized our previously reported single specimen RNA-Seq data of the embryonic and juvenile life stages of the parthenogenetic tardigrades H. exemplaris and R. varieornatus to evaluate whether DMRT orthologs in others may be functional. Only females have been observed in both species, suggesting the lack of masculinization in these species. All three Dmrt11E, Dmrt93B, and Dmrt99B orthologs in H. exemplaris and R. varieornatus (RvDmrt11E: g5527, RvDmrt93B: g9000, RvDmrt99B: g7078; HeDmrt11E: BV898_08851, HeDmrt93B: BV898_13063, HeDmrt99B: BV898_01934.) were expressed during embryonic stages (Figure S3), where Dmrt11E preceded Dmrt99B in both species, and the three DMRT genes were expressed at lower levels in juvenile and adult stages.

We further investigated the functionality of the DMRT orthologs by functional domain detection (Fig. 3C). While all five DMRT copies harbored the DM domain at the N-terminus, they did not contain the dimerization domain known to exist in dsx proteins required for DNA binding and sex-specific splicing variants. Furthermore, we did not find the ubiquitin binding-related CUE-DMA domain in PmDmrt3090 by domain search analysis. However, multiple alignment of the five orthologs and the D. melanogaster DMRT sequence suggested the conservation of several residues within the region corresponding to the CUE-DMA domain, implying the conservation of this domain. By modeling the protein structure with AlphaFold2 and aligning the D. melanogaster Dmrt93B structure, we observed that the C-terminal region showed structural homology with CUE-DMA domain-like helices (RMSD: 0.276–0.574, Figure S4CDE), suggesting that PmDmrt3090 may also harbor the CUE-DMA domain.

Contradictory data from whole genome sequencing and PCR-based genotyping

Based on our observations of several male-biased but not male-specific genomic regions, we hypothesized that these regions were not sex-specific chromosome structures. To evaluate this, we sequenced the genomes of both sexes at low coverage. We produced approximately 50–60 M reads, corresponding to roughly 20–25 × coverage (Table S2). Approximately 80–90% of the reads were mapped to the genome, resulting in roughly 15–20 × coverage.

We first calculated the coverage of the 10 kbp bins genome wide. Initial PCA of the coverage profiles did not show a clear difference between males and females (Fig. 4A). We identified several bins with half of the average genome-wide coverage that were not found in females (Fig. 4B). These characteristics are similar to those of heterozygotic chromosomes, particularly the X chromosome of males in the XY sex determination system. All of the bins that were identified as male-biased by the transcriptome analysis had genome-wide average coverage, suggesting that all regions exist in females (Fig. 4C). We also evaluated whether we could detect male or female specific regions through k-mer based analysis using the YGS method. Scaffolds that have a high number of “percent validated single-copy unmatched k-mers” (P_VSC_UK) indicate sex-specificity. Although no scaffolds had a P_VSC_UK ratio of 100, we detected five scaffolds fitting the XY sex chromosome structure rather than the ZW scheme with an arbitrary threshold of P_VSC_UK > 80 (Fig. 4D, E, F). Similar profiles were observed by SPADES reassembly using all gDNA-Seq data from our and previous studies (Fig. 4G, H, I). We also noticed that many scaffolds from both assemblies had P_VSC_UK values of approximately 50% in both female-to-male (XY) and male-to-female (ZW) analyses (Fig. 4F, I), which indicates that the corresponding region is both male- and female-specific.

Fig. 4
figure 4

Paramacrobiotus metropolitanus lacks sex specific regions for both sexes. A PCA of genomic coverage profiles for male and female gDNA-Seq data. B Average coverage for 10 kbp bins, normalized by the median of all bins. The brown and red lines indicate the whole-genome average and half-genome average for males and females, respectively. The range between the 0 × and 10 × coverage ratio to the median is shown. C Genome coverage of male-biased bins. D, E, G, H YGS analysis for (D, G) female-to-male (XY system) and (E, H) male-to-female (ZW system) based on the (D, E) published P. metropolitanus assembly and (G, H) SPADES reassembly. The red line indicates a P-VSC-UK threshold of 80. F, I Scatter plot for P-VSC-UK values for male-to-female and female-to-male analysis for the (F) published genome and the (I) SPADES reassembly. Contigs shorter than 1,000 bp were removed from the SPADES plot. J gDNA-Seq coverage for all samples and the location of genotyping primers within contig Parri_scaffold0000295. Blue and red correspond to male and female samples, respectively. K Electrophoresis of genotyping primers designed for male specificity. Male specificity was not observed for any primer set. L gDNA-Seq coverage of PmDmrt3090/PmDmrt3093 harboring scaffold Parri_scaffold0000005 for all samples. Colors indicate samples for each sex

To evaluate the male specificity observed in the in-silico analysis, we designed several primers to amplify regions in the scaffold Parri_scaffold0000295 that were identified as male-specific (Fig. 4J, Table S2). Evaluating individual genomic coverage indicated that in this scaffold, a single male sample had near-zero coverage, in contradiction to the other two male samples (Fig. 4J). However, PCR genotyping indicated the existence of this region in females as well, which contradicted the results obtained from the in-silico analysis (Fig. 4K). Thus, we concluded that we could not derive sex chromosomes or male-specific regions, and the male-specific regions detected above may have been an artifact of differences between individuals. We also evaluated the male-specificity of the paralogous DMRT loci on Parri_scaffold0000005, where coverage analysis suggested that this region was not male-specific (Fig. 4L).

We noticed a relatively low level of RNA-Seq mappability to the reported genome (~ 90%), which led us to re-evaluate the current genome. Completeness analysis indicated 72.9% completeness with the most recent version of BUSCO. Furthermore, we observed several scaffolds with inconsistent coverage distribution in our sex-separated data, but not in Hara et al. Illumina data [18]. Therefore, we tested if recent assemblers would result in a more contiguous and complete assembly, compared to the Canu assembled current genome. However, we were not able to obtain a more complete genome, with the maximum being a 0.4% increase for the assembly derived with NextDenovo + NextPolish (Table S6). Other statistics had a large increase; N50 from 1.0 M to 1.3 M, longest scaffold length 4.48 M to 9.23 M. For comparison, we evaluated the completeness of other high-quality tardigrades genomes, namely R. varieornatus and H. exemplaris. Both BUSCO and compleasm resulted in completeness values similar to P. metropolitanus, R. varieornatus (C:74.6%) and H. exemplaris (C:73.3%). These data suggests either tardigrade genomes may lack some BUSCO genes, or the gene detection algorithm of the current BUSCO software may not fit the genome of tardigrades, resulting in lower BUSCO scores. Therefore, we used the current genome for P. metropolitanus for later analysis.

Sex-bias in anhydrobiosis-related genes

A major feature of tardigrades is their ability to survive environmental extremes, a phenomenon known as cryptobiosis [76]. Tolerance to near-complete desiccation is known as anhydrobiosis [77]. Several tardigrade-specific gene families, i.e. cytosolic-abundant heat soluble (CAHS) and Secretory Abundant Heat Soluble (SAHS), have been implicated in anhydrobiosis protection [22]. A recent study observed tissue-specific expression of anhydrobiosis genes [78]. Both males and females are capable of anhydrobiosis, in which protective genes are expressed in sex-specific organs, such as the testes or ovaries. Therefore, we hypothesized the presence of sex-biased anhydrobiosis genes.

We used our previously reported RNA-Seq data for the hydrated active state and the tun state, desiccated for two days, to identify genes induced during anhydrobiosis. We detected approximately 4,500 differentially expressed transcripts, slightly fewer than in our previous report, possibly due to the different methods used for differential expression analysis. We then compared the expression profiles of anhydrobiosis and between sexes and observed approximately 1,800 transcripts that were differentially expressed under both conditions (Fig. 5A). As hypothesized, we observed that three CAHS and one SAHS ortholog were sex-biased, possibly indicating tissue specificity (Fig. 5A). Interestingly, all three CAHS orthologs induced in males were the only three among the 13 CAHS orthologs that were not differentially expressed during anhydrobiosis (Table S7). Phylogenetic analysis indicated that these CAHS orthologs were CAHS1 (PARRI_0016931), putative CAHS5 (PARRI_0006576), and CAHS5 (PARRI_0002229) orthologs, following the proposed naming scheme of Fleming et al. [54]. In contrast, the SAHS ortholog, detected as differentially expressed, was induced in the females. We also found six orthologs of tardigrade-specific manganese-dependent peroxidase [33] to be highly expressed in males but not in females. Only four genes were found to be induced during anhydrobiosis.

Fig. 5
figure 5

Sexual bias in anhydrobiosis genes and identification of PmDsup ortholog. A Comparison of gene expression profiles between the sexes during anhydrobiosis. Log2 (Tun + 0.1) / (Active + 0.1) were plotted for the x-axis and for the y-axis log2 (Male + 0.1) / (Female + 0.1). Red dots indicate transcripts detected as differentially expressed in both comparisons. Sex-biased, but not anhydrobiosis induced, anhydrobiosis related genes have been noted by colored annotations: Blue: CAHS, Green: SAHS, Red: AMNP, and Black: Dsup. B Synteny analysis to identify orthologous genomic loci in H. exemplaris and P. metropolitanus. HeDsup and PmDsup have been annotated in the plot. C DISOPRED and IUPRED3 scores (D) Protein structure predicted by ColabFold. The N-terminus to the C-terminus shows gradient colors from blue to red. E Expression of PmDsup in both sexes

Based on the identification of the H. exemplaris ortholog of the Damage suppressor (Dsup, BV898_01301) gene, we also searched for a P. metropolitanus Dsup ortholog through gene synteny with H. exemplaris [31, 79]. We identified the gene PARRI_0005796 as a Dsup ortholog candidate (Fig. 5B). This protein was annotated as “transcriptional regulatory protein AlgP’’ in NCBI; however, (1) no functional domains were identified by InterProScan, (2) no BLAST hits to known proteins (E-value < 1e-5), (3) highly disordered throughout the whole protein (Fig. 5C), and (4) a predicted nuclear localization signal (DeepLoc2, 0.7715 probability), suggesting that this protein may be a Dsup ortholog. The AlphaFold2 structure prediction also implied a lack of globular structure (Fig. 5D). PmDsup was significantly upregulated in females (TPM, female: 280, male: 82, FDR = 1.31 × 10–6, Fig. 5E), implying the importance of this gene in females.


In this study, we focused on gonochoristic tardigrade P. metropolitanus to identify possible factors that affect sexual dimorphism. Cytological studies have not identified definitive sex-linked chromosomes in tardigrades [80, 81], and multiple reports have observed biased sex ratios in tardigrades [23, 82,83,84,85]. These observations suggests that sex determination in tardigrades may not depend on the random distribution of sex chromosomes (or the existence of a sex chromosome). Even in the absence of sex chromosomes, as hypothesized in tardigrades, genomic loci affecting sexual dimorphism could exist, which may be detected by comprehensive omics methods.

Therefore, we aimed to characterize the molecular basis of sexual dimorphism in tardigrades by comparing the transcriptome and the genome between the sexes of P. metropolitanus. We hypothesized that sex-linked genes may be related to sex determination or dimorphism, and if focused on a small genomic region, may imply a sex-determining region, such as the M factor found in many eukaryotes [86]. Transcriptome analysis of both sexes indicated a large number of sex-biased genes, despite the small morphological sex-linked differences in Macrobiotidae, with the exception of their germline [20]. We observed upregulation of genes related to spermatogenesis in males, which reflects the activation of spermatogenesis, and large amounts of sperm are continuously produced in adult males [23, 83]. In contrast to that in males, DNA replication- and meiosis-related genes were highly expressed in females. Females undergo DNA replication not only to produce oocytes through meiosis [23] but also to shed the cuticular exoskeleton during the last stage of the reproductive process (simplex stage) [23, 87]. Mitotic cells are generally observed in the post-simplex stage [88]. Together, the regulation of DNA replication and meiosis is consistent with the production of mitotic cells and extensive replication of the epidermal layer [88, 89].

We identified a small gene set highly biased toward males, but missing in females, which we hypothesized may be related to sexual dimorphism. Genome loci enrichment analysis of this gene set identified approximately 325 bins spanning 29 scaffolds as male-biased. This region was enriched in sperm and ion transport-related genes, which is consistent with the production of sperm at the adult male life stage. To evaluate sex specificity, we produced low-coverage genome sequencing data to evaluate sex-specific regions and observed that most regions were present in the genomes of both sexes. Genome-wide analysis revealed several male-specific regions; however, PCR evaluation produced contradictory results. We used a laboratory-cultured TYO strain of P. metropolitanus for genome and transcriptome sequencing, therefore, we anticipated low levels of heterozygosity within the culture population. However, the results obtained at this stage suggest that the genomic differences we detected as sex-linked can be explained by individual variability. Additionally, during the YGS analysis, we observed a high number of contigs that showed approximately 50% P_VSC_UK, suggesting that there are a large number of contigs that contain sequences specific for both sexes, which we hypothesize that individual variability may have caused this abnormal distribution. Together, the lack of sex-specific regions may indicate that the difference between sexes is due to epigenetic modifications.

One of the key findings of this study is the accumulation of knowledge for sex determination cascade-related genes, particularly the DMRT gene family. The DMRT family is a highly conserved transcription factor that plays an important role in sex differentiation in many animals and has been studied extensively in insects [7]. Several studies have identified DMRT orthologs to be located on the sex chromosomes and regulate the growth of sex-specific tissues [8, 9]. The evolutionary background of this gene family has been extensively analyzed in other lineages [7]; however, such analysis has been overlooked. In our analysis, we identified a Macrobiotidae-specific Dmrt93B subfamily located in a male-biased region, which we termed the 3090/3093 complex in addition to the Dmrt99E, Dmrt93B, and Drmt11E subfamilies. Orthologs of DMRT genes in both gonochoristic and parthenogenetic tardigrades are expressed during several stage of development, thus the genes may be functional. Family-specific orthologs of Dmrt93B subfamily have been found in Macrobiotidae and several Echiniscidae. The lack of two copies in other Macrobiotidae species may be the result of misassembly in their genomes; the analyzed genomes are based on Illumina short reads, and the extremely similar 30–180 aa (corresponding to approximately 450 bp) may have resulted in a misassembly. While conservation in Echiniscidae complicates the evolution of this subfamily, the identification of orthologs in various Macrobiotidae species suggests that this is an important DMRT subfamily. In fact, the two 3090/3093 complex paralogs were expressed higher in males, similar to Daphnia Dsx1 [11, 12], suggesting that these subfamily orthologs may inhibit feminization or progress musculation. Furthermore, we did not find any orthologs of the dsx subfamily in any of the tardigrade genomes analyzed, and we did not identify splicing variants in any of the P. metropolitanus DMRT orthologs, suggesting a sex differentiation cascade different from those that rely on sex-specific dsx splicing variants like those observed in insects. In addition, we observed sex-biased expression in tra2 splicing variants which may be functional in the P. metropolitanus sex determination cascade; however, several factors in this cascade may be lost in this lineage. The Bombyx mori sxl gene induces dimorphism of the sperm, not sex determination [90]; therefore, it is possible that the lack of sxl may imply a different regulatory pathway than is known.

Tardigrades are renowned for their ability to tolerate extreme stress [22], and P. metropolitanus also shows a high tolerance to desiccation stress [18]. Interestingly, we observed sex-biased expression of several anhydrobiosis genes, hypothesized to play protective roles during anhydrobiosis [30,31,32,33, 91]. For instance, CAHS genes are tardigrade-specific proteins that form gel filaments that possibly protect cells [92,93,94]. Recent studies have observed tissue/organelle specificity for these proteins, which further implies the existence of orthologs with sex-specific expression [78]. Therefore, we hypothesized that orthologs of such genes may exhibit sex-specific expression to protect sex-specific organs. Indeed, we identified CAHS, SAHS, and AMNP orthologs with sex-specific expression. Two of the three male-induced CAHS orthologs were highly expressed, but were not induced differentially between active and anhydrobiotic conditions. This may imply that these CAHS orthologs participate in the protection of male-specific tissues or sperm. Furthermore, we identified a P. metropolitanus Dsup ortholog that is highly expressed in females. Coupled with the observation of the enrichment of meiosis-related genes from transcriptome analysis, we suggest that Dsup may actively function to accommodate the production of oocytes/oogenesis rather than spermatozoa/spermatogenesis. In contrast, AMNP, a tardigrade-specific peroxidase, was highly expressed in males, suggesting enhanced protection against oxidative stress. Similar observations have been made in the sperm of many animals [95, 96]. Together, the sex-biased expression of anhydrobiosis genes may provide protection for sex-specific tissues.


In this study, we identified male-biased regions that may harbor potential candidates that regulate sexual dimorphism in the gonochoristic tardigrade P. metropolitanus. Simultaneously, these data denied the sex-chromosome-based sex determination scheme. We also provide evidence for a new DMRT subfamily that may contribute to sex differentiation in this family. The 3090/3093 complex DMRT paralogs may be initial candidates for disruption or gene editing for evaluation their relationships with sex determination [78, 97,98,99]. Future studies utilizing high-quality genomes and careful physiological experiments are required to reveal sex determination cues not only in this species but also in other tardigrades.

Availability of data and materials

The raw reads for genome DNA-Seq were submitted to NCBI SRA under Bioproject PRJNA1063779. The raw reads and processed expression profiles were uploaded to NCBI GEO under the accession ID GSE253242. Other datasets analyzed in this study have been uploaded to figshare (


  1. Hamilton WD, Axelrod R, Tanese R. Sexual reproduction as an adaptation to resist parasites (a review). Proc Natl Acad Sci U S A. 1990;87(9):3566–73.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  2. Crow JF. Advantages of sexual reproduction. Dev Genet. 1994;15(3):205–13.

    Article  CAS  PubMed  Google Scholar 

  3. Hopkins BR, Kopp A. Evolution of sexual development and sexual dimorphism in insects. Curr Opin Genet Dev. 2021;69:129–39.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  4. Marin I, Baker BS. The evolutionary dynamics of sex determination. Science. 1998;281(5385):1990–4.

    Article  CAS  PubMed  Google Scholar 

  5. Zarkower D. Establishing sexual dimorphism: conservation amidst diversity? Nat Rev Genet. 2001;2(3):175–85.

    Article  CAS  PubMed  Google Scholar 

  6. Haag ES, Doty AV. Sex determination across evolution: connecting the dots. PLoS Biol. 2005;3(1):e21.

    Article  PubMed  PubMed Central  Google Scholar 

  7. Matson CK, Zarkower D. Sex and the singular DM domain: insights into sexual regulation, evolution and plasticity. Nat Rev Genet. 2012;13(3):163–74.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Yoshimoto S, Okada E, Umemoto H, Tamura K, Uno Y, Nishida-Umehara C, et al. A W-linked DM-domain gene, DM-W, participates in primary ovary development in Xenopus laevis. Proc Natl Acad Sci U S A. 2008;105(7):2469–74.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Matsuda M, Nagahama Y, Shinomiya A, Sato T, Matsuda C, Kobayashi T, et al. DMY is a Y-specific DM-domain gene required for male development in the medaka fish. Nature. 2002;417(6888):559–63.

    Article  CAS  PubMed  Google Scholar 

  10. Kleiven OT, Larsson P, Hobæk A. Sexual reproduction in Daphnia magna requires three stimuli. Oikos. 1992;65(2):197–206.

    Article  Google Scholar 

  11. Kato Y, Kobayashi K, Watanabe H, Iguchi T. Environmental sex determination in the branchiopod crustacean Daphnia magna: deep conservation of a Doublesex gene in the sex-determining pathway. PLoS Genet. 2011;7(3):e1001345.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  12. Kato Y, Perez CAG, Mohamad Ishak NS, Nong QD, Sudo Y, Matsuura T, et al. A 5’ UTR-overlapping lncRNA activates the male-determining gene Doublesex1 in the crustacean Daphnia magna. Curr Biol. 2018;28(11):1811-7 e4.

    Article  CAS  PubMed  Google Scholar 

  13. Actual checklist of Tardigrada species. 2023. Available from: Cited 24 Jan 2024.

  14. Rahm PG. A new order of tardigrades from the hot springs of Japan (Furu-Yusection, Unzen). Zool Sci. 1937;16(4):345–52.

    Google Scholar 

  15. Grothman GT, Johansson C, Chilton G, Kagoshima H, Tsujimoto M, Suzuki AC. Gilbert Rahm and the status of Mesotardigrada Rahm, 1937. Zoolog Sci. 2017;34(1):5–10.

    Article  PubMed  Google Scholar 

  16. Murai Y, Yagi-Utsumi M, Fujiwara M, Tanaka S, Tomita M, Kato K, et al. Multiomics study of a heterotardigrade, Echinisicus testudo, suggests the possibility of convergent evolution of abundant heat-soluble proteins in Tardigrada. BMC Genomics. 2021;22(1):813.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Arakawa K. Examples of extreme survival: Tardigrade genomics and molecular anhydrobiology. Annu Rev Anim Biosci. 2022;10:17–37.

    Article  CAS  PubMed  Google Scholar 

  18. Hara Y, Shibahara R, Kondo K, Abe W, Kunieda T. Parallel evolution of trehalose production machinery in anhydrobiotic animals via recurrent gene loss and horizontal transfer. Open Biol. 2021;11(7):200413.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Altiero T, Suzuki AC, Rebecchi L. Reproduction, development and life cycles. Water Bears: The Biology of Tardigrades. Zoological Monographs. vol. 2. Barsel (Switzerland). Springer Nature. 2018. p. 211–47.

  20. Nelson DR, Guidetti R, Rebecchi L. Phylum Tardigrada. Ecology and General Biology, Vol I: Thorp and Covich's Freshwater Invertebrates, 4th Edition. 2015:347–80.

  21. Sugiura K, Matsumoto M. Sexual reproductive behaviours of tardigrades: a review. Invertebr Reprod Dev. 2021;65(4):279–87.

    Article  Google Scholar 

  22. Yoshida Y, Tanaka S. Deciphering the biological enigma-genomic evolution underlying anhydrobiosis in the Phylum Tardigrada and the chironomid Polypedilum vanderplanki. Insects. 2022;13(6):577.

  23. Sugiura K, Minato H, Suzuki AC, Arakawa K, Kunieda T, Matsumoto M. Comparison of sexual reproductive behaviors in two species of Macrobiotidae (Tardigrada: Eutardigrada). Zoolog Sci. 2019;36(2):120–7.

    Article  CAS  PubMed  Google Scholar 

  24. Sugiura K, Matsumoto M, Kunieda T. Description of a model tardigrade Paramacrobiotus metropolitanus sp. nov. (Eutardigrada) from Japan with a summary of its life history, reproduction and genomics. Zootaxa. 2022;5134(1):92–112.

    Article  PubMed  Google Scholar 

  25. Sugiura K, Shiba K, Inaba K, Matsumoto M. Morphological differences in tardigrade spermatozoa induce variation in gamete motility. BMC Zool. 2022;7(1):8.

    Article  PubMed  PubMed Central  Google Scholar 

  26. Arakawa K, Yoshida Y, Tomita M. Genome sequencing of a single tardigrade Hypsibius dujardini individual. Sci Data. 2016;3:160063.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Wang J, Chitsaz F, Derbyshire MK, Gonzales NR, Gwadz M, Lu S, et al. The conserved domain database in 2023. Nucleic Acids Res. 2023;51(D1):D384–8.

    Article  CAS  PubMed  Google Scholar 

  28. Thumuluri V, Almagro Armenteros JJ, Johansen AR, Nielsen H, Winther O. DeepLoc 2.0: multi-label subcellular localization prediction using protein language models. Nucleic Acids Res. 2022;50(W1):W228–34.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. Jones P, Binns D, Chang HY, Fraser M, Li W, McAnulla C, et al. InterProScan 5: genome-scale protein function classification. Bioinformatics. 2014;30(9):1236–40.

    Article  CAS  PubMed  Google Scholar 

  30. Yamaguchi A, Tanaka S, Yamaguchi S, Kuwahara H, Takamura C, Imajoh-Ohmi S, et al. Two novel heat-soluble protein families abundantly expressed in an anhydrobiotic tardigrade. PLoS ONE. 2012;7(8):e44209.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Hashimoto T, Horikawa DD, Saito Y, Kuwahara H, Kozuka-Hata H, Shin IT, et al. Extremotolerant tardigrade genome and improved radiotolerance of human cultured cells by tardigrade-unique protein. Nat Commun. 2016;7:12808.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  32. Yoshida Y, Koutsovoulos G, Laetsch DR, Stevens L, Kumar S, Horikawa DD, et al. Comparative genomics of the tardigrades Hypsibius dujardini and Ramazzottius varieornatus. PLoS Biol. 2017;15(7):e2002266.

    Article  PubMed  PubMed Central  Google Scholar 

  33. Yoshida Y, Satoh T, Ota C, Tanaka S, Horikawa DD, Tomita M, et al. Time-series transcriptomic screening of factors contributing to the cross-tolerance to UV radiation and anhydrobiosis in tardigrades. BMC Genomics. 2022;23(1):405.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Pertea G, Pertea M. GFF Utilities: GffRead and GffCompare. F1000Res. 2020;9:304.

  35. Mirdita M, Schutze K, Moriwaki Y, Heo L, Ovchinnikov S, Steinegger M. ColabFold: making protein folding accessible to all. Nat Methods. 2022;19(6):679–82.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Meng EC, Goddard TD, Pettersen EF, Couch GS, Pearson ZJ, Morris JH, et al. UCSF ChimeraX: Tools for structure building and analysis. Protein Sci. 2023;32(11):e4792.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Hoencamp C, Dudchenko O, Elbatsh AMO, Brahmachari S, Raaijmakers JA, van Schaik T, et al. 3D genomics across the tree of life reveals condensin II as a determinant of architecture type. Science. 2021;372(6545):984–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12:323.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.

    Article  PubMed  PubMed Central  Google Scholar 

  40. Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011;29(7):644–52.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Falcon S, Gentleman R. Using GOstats to test gene lists for GO term association. Bioinformatics. 2007;23(2):257–8.

    Article  CAS  PubMed  Google Scholar 

  42. Morgan M, Falcon S, Gentleman R. GSEABase: Gene set enrichment data structures and methods 2023. Available from:

  43. Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 2010;26(6):841–2.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Pertea M, Pertea GM, Antonescu CM, Chang TC, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. 2015;33(3):290–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  45. Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12(4):357–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Haas BJ, Delcher AL, Mount SM, Wortman JR, Smith RK Jr, Hannick LI, et al. Improving the Arabidopsis genome annotation using maximal transcript alignment assemblies. Nucleic Acids Res. 2003;31(19):5654–66.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Haas BJ, Salzberg SL, Zhu W, Pertea M, Allen JE, Orvis J, et al. Automated eukaryotic gene structure annotation using EVidenceModeler and the program to assemble spliced alignments. Genome Biol. 2008;9(1):R7.

    Article  PubMed  PubMed Central  Google Scholar 

  48. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence alignment/map format and SAMtools. Bioinformatics. 2009;25(16):2078–9.

    Article  PubMed  PubMed Central  Google Scholar 

  49. Altschul SF, Madden TL, Schaffer AA, Zhang J, Zhang Z, Miller W, et al. Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. Nucleic Acids Res. 1997;25(17):3389–402.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  50. Toyota K, Kato Y, Sato M, Sugiura N, Miyagawa S, Miyakawa H, et al. Molecular cloning of doublesex genes of four cladocera (water flea) species. BMC Genomics. 2013;14:239.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  51. Katoh K, Standley DM. MAFFT multiple sequence alignment software version 7: improvements in performance and usability. Mol Biol Evol. 2013;30(4):772–80.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. Minh BQ, Schmidt HA, Chernomor O, Schrempf D, Woodhams MD, von Haeseler A, et al. IQ-TREE 2: New models and efficient methods for phylogenetic inference in the genomic era. Mol Biol Evol. 2020;37(5):1530–4.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  53. Yoshida Y, Sugiura K, Tomita M, Matsumoto M, Arakawa K. Comparison of the transcriptomes of two tardigrades with different hatching coordination. BMC Dev Biol. 2019;19(1):24.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. Fleming JF, Pisani D, Arakawa K. The evolution of temperature and desiccation-related protein families in Tardigrada reveals a complex acquisition of extremotolerance. Genome Biol Evol. 2024;16(1): evad217.

  55. Koren S, Walenz BP, Berlin K, Miller JR, Bergman NH, Phillippy AM. Canu: scalable and accurate long-read assembly via adaptive k-mer weighting and repeat separation. Genome Res. 2017;27(5):722–36.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  56. Hu J, Wang Z, Sun Z, Hu B, Ayoola AO, Liang F, et al. An efficient error correction and accurate assembly tool for noisy long reads. bioRxiv. 2023;2023.03.09.531669.

  57. Shafin K, Pesout T, Lorig-Roach R, Haukness M, Olsen HE, Bosworth C, et al. Nanopore sequencing and the Shasta toolkit enable efficient de novo assembly of eleven human genomes. Nat Biotechnol. 2020;38(9):1044–53.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  58. Kolmogorov M, Yuan J, Lin Y, Pevzner PA. Assembly of long, error-prone reads using repeat graphs. Nat Biotechnol. 2019;37(5):540–6.

    Article  CAS  PubMed  Google Scholar 

  59. Ruan J, Li H. Fast and accurate long-read assembly with wtdbg2. Nat Methods. 2020;17(2):155–8.

    Article  CAS  PubMed  Google Scholar 

  60. Wong J, Coombe L, Nikolic V, Zhang E, Nip KM, Sidhu P, et al. Linear time complexity de novo long read genome assembly with GoldRush. Nat Commun. 2023;14(1):2906.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  61. Bankevich A, Nurk S, Antipov D, Gurevich AA, Dvorkin M, Kulikov AS, et al. SPAdes: a new genome assembly algorithm and its applications to single-cell sequencing. J Comput Biol. 2012;19(5):455–77.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  62. Nie F, Huang N, Zhang J, Ni P, Wang Z, Xiao C, et al. de novo diploid genome assembly using long noisy reads via haplotype-aware error correction and inconsistent overlap identification. bioRxiv. 2023;2022.09.25.509436;.

  63. Vaser R, Sikic M. Time- and memory-efficient genome assembly with Raven. Nat Comput Sci. 2021;1(5):332–6.

    Article  PubMed  Google Scholar 

  64. Hu J, Fan J, Sun Z, Liu S. NextPolish: a fast and efficient genome polishing tool for long-read assembly. Bioinformatics. 2020;36(7):2253–5.

    Article  CAS  PubMed  Google Scholar 

  65. Huang N, Li H. compleasm: a faster and more accurate reimplementation of BUSCO. Bioinformatics. 2023;39(10): btad595.

  66. Simao FA, Waterhouse RM, Ioannidis P, Kriventseva EV, Zdobnov EM. BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs. Bioinformatics. 2015;31(19):3210–2.

    Article  CAS  PubMed  Google Scholar 

  67. Vasimuddin M, Misra S, Li H, Aluru S. Efficient architecture-aware acceleration of BWA-MEM for multicore systems. 2019 IEEE International Parallel and Distributed Processing Symposium (IPDPS)2019. p. 314–24.

  68. Carvalho AB, Clark AG. Efficient identification of Y chromosome sequences in the human and Drosophila genomes. Genome Res. 2013;23(11):1894–907.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  69. Marcais G, Kingsford C. A fast, lock-free approach for efficient parallel counting of occurrences of k-mers. Bioinformatics. 2011;27(6):764–70.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  70. Tang H, Bowers JE, Wang X, Ming R, Alam M, Paterson AH. Synteny and collinearity in plant genomes. Science. 2008;320(5875):486–8.

    Article  CAS  PubMed  Google Scholar 

  71. Jones DT, Cozzetto D. DISOPRED3: precise disordered region predictions with annotated protein-binding activity. Bioinformatics. 2015;31(6):857–63.

    Article  CAS  PubMed  Google Scholar 

  72. Erdos G, Pajkos M, Dosztanyi Z. IUPred3: prediction of protein disorder enhanced with unambiguous experimental annotation and visualization of evolutionary conservation. Nucleic Acids Res. 2021;49(W1):W297–303.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  73. Ahringer J. Reverse genetics. WormBook. The C. elegans Research Community. 2006.

  74. Untergasser A, Cutcutache I, Koressaar T, Ye J, Faircloth BC, Remm M, et al. Primer3–new capabilities and interfaces. Nucleic Acids Res. 2012;40(15):e115.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  75. Ammermann D. The cytology of parthenogenesis in the tardigrade Hypsibius dujardini. Chromosoma. 1967;23(2):203–13.

    Article  CAS  PubMed  Google Scholar 

  76. Keilin D. The problem of anabiosis or latent life: history and current concept. Proc R Soc Ser B-Bio. 1959;150(939):149–91.

    CAS  Google Scholar 

  77. Bertolani R, Guidetti R, Jönsson KI, Altiero T, Boschini D, Rebecchi L. Experiences with dormancy in tardigrades. J Limnol. 2004;63(Supplement 1):16–25.

    Article  Google Scholar 

  78. Tanaka S, Aoki K, Arakawa K. In vivo expression vector derived from anhydrobiotic tardigrade genome enables live imaging in Eutardigrada. Proc Natl Acad Sci U S A. 2023;120(5):e2216739120.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  79. Chavez C, Cruz-Becerra G, Fei J, Kassavetis GA, Kadonaga JT. The tardigrade damage suppressor protein binds to nucleosomes and protects DNA from hydroxyl radicals. Elife. 2019;8:e47682.

  80. Rebecchi L, Altiero T, Bertolani R. Banding techniques on tardigrade chromosomes: the karyotype of Macrobiotus richtersi (Eutardigrada, Macrobiotidae). Chromosome Res. 2002;10(6):437–43.

    Article  CAS  PubMed  Google Scholar 

  81. Altiero T, Rebecchi L. First evidence of achiasmatic male meiosis in the water bears Richtersius coronifer and Macrobiotus richtersi (Eutardigrada, Macrobiotidae). Hereditas. 2003;139(2):116–20.

    Article  PubMed  Google Scholar 

  82. Bertolani R, Rebecchi L, Beccaccioli G. Dispersal of Ramazzottius and other tardigrades in relation to type of reproduction. Invertebr Reprod Dev. 1990;18(3):153–7.

    Article  Google Scholar 

  83. Rebecchi L, Bertolani R. Maturative pattern of ovary and testis in eutardigrades of freshwater and terrestrial habitats. Invertebr Reprod Dev. 1994;26(2):107–17.

    Article  Google Scholar 

  84. Rebecchi L, Rossi V, Altiero T, Bertolani R, Menozzi P. Reproductive modes and genetic polymorphism in the tardigrade Richtersius coronifer (Eutardigrada, Macrobiotidae). Invertebr Biol. 2003;122(1):19–27.

    Article  Google Scholar 

  85. Suzuki AC. Appearance of males in a thelytokous strain of Milnesium cf. tardigradum (Tardigrada). Zoolog Sci. 2008;25(8):849–53.

    Article  PubMed  Google Scholar 

  86. Davey J. Mating pheromones of the fission yeast Schizosaccharomyces pombe: purification and structural characterization of M-factor and isolation and analysis of two genes encoding the pheromone. EMBO J. 1992;11(3):951–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  87. Guidetti R. Evolution of egg deposition strategies, exaptations of exuvia, and thanatochresis in tardigrades. Organisms Diversity & Evolution. 2024.

  88. Czernekova M, Jönsson KI. Mitosis in storage cells of the eutardigrade. Zool J Linn Soc. 2016;178(4):888–96.

    Article  Google Scholar 

  89. Gross V, Bahrle R, Mayer G. Detection of cell proliferation in adults of the water bear Hypsibius dujardini (Tardigrada) via incorporation of a thymidine analog. Tissue Cell. 2018;51:77–83.

    Article  CAS  PubMed  Google Scholar 

  90. Sakai H, Oshima H, Yuri K, Gotoh H, Daimon T, Yaginuma T, et al. Dimorphic sperm formation by Sex-lethal. Proc Natl Acad Sci U S A. 2019;116(21):10412–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  91. Tanaka S, Tanaka J, Miwa Y, Horikawa DD, Katayama T, Arakawa K, et al. Novel mitochondria-targeted heat-soluble proteins identified in the anhydrobiotic tardigrade improve osmotic tolerance of human cells. PLoS ONE. 2015;10(2):e0118272.

    Article  PubMed  PubMed Central  Google Scholar 

  92. Tanaka A, Nakano T, Watanabe K, Masuda K, Honda G, Kamata S, et al. Stress-dependent cell stiffening by tardigrade tolerance proteins through reversible formation of cytoskeleton-like filamentous network and gel-transition. bioRxiv. 2022.

  93. Yagi-Utsumi M, Aoki K, Watanabe H, Song C, Nishimura S, Satoh T, et al. Desiccation-induced fibrous condensation of CAHS protein from an anhydrobiotic tardigrade. Sci Rep. 2021;11(1):21328.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  94. Malki A, Teulon JM, Camacho-Zarco AR, Chen SW, Adamski W, Maurin D, et al. Intrinsically disordered tardigrade proteins self-assemble into fibrous gels in response to environmental stress. Angew Chem Int Ed Engl. 2022;61(1):e202109961.

    Article  CAS  PubMed  Google Scholar 

  95. O’Flaherty C, Scarlata E. Oxidative stress and reproductive function: The protection of mammalian spermatozoa against oxidative stress. Reproduction. 2022;164(6):F67–78.

    Article  CAS  PubMed  Google Scholar 

  96. Tsai Y, Lin YC, Lee YH. Octopamine-MAPK-SKN-1 signaling suppresses mating-induced oxidative stress in Caenorhabditis elegans gonads to protect fertility. iScience. 2023;26(3):106162.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  97. Tenlen JR, McCaskill S, Goldstein B. RNA interference can be used to disrupt gene function in tardigrades. Dev Genes Evol. 2013;223(3):171–81.

    Article  CAS  PubMed  Google Scholar 

  98. Kumagai H, Kondo K, Kunieda T. Application of CRISPR/Cas9 system and the preferred no-indel end-joining repair in tardigrades. Biochem Biophys Res Commun. 2022;623:196–201.

    Article  CAS  PubMed  Google Scholar 

  99. Kondo K, Tanaka A, Kunieda T. Single-step generation of homozygous knock-out/knock-in individuals in an extremotolerant parthenogenetic tardigrade using DIPA-CRISPR. bioRxiv. 2024;

Download references


We thank Yuki Takai and Naoko Ishii (Keio University) for experimental support. For giving us many helpful comments, we also thank two anonymous reviewers, Dr. Diane R. Nelson (East Tennessee University), Dr. Hajime Watanabe (University of Osaka), Dr. Yasuhiko Kato (University of Osaka), Dr. Atsushi C. Suzuki (Keio University), Yu Saito (Keio University), and Ryo Ogushi (Keio University). RNA-seq was performed by Chemical Dojin Co. Ltd. We also thank the members of the Japanese Society for Tardigradology for fruitful discussions.


Grant-in-aid KAKENHI (JP18J21345) to KS, (21H05279) to KA, grant for Research Project from Research and Education Center for Natural Science, Keio University for MM and KS. Joint Research of the 13 Exploratory Research Centers on Life and Living Systems (ExCELLS program 19–501, 22EXC601) to KA.

Author information

Authors and Affiliations



KS, KH, and MM prepared specimens. KA performed the genome sequencing. KS, TK, and MM performed RNA-seq. KS and YY analyzed the data. KS and YY drafted the manuscript. AK, TK, and MM improved the manuscript. KS, YY, AK, TK, and MM designed this study.

Corresponding author

Correspondence to Midori Matsumoto.

Ethics declarations

Ethics approval and consent to participate

Tardigrades (invertebrates) were used for this study. No vertebrate or human individuals or tissues were used. The laws and regulations set forth by the Ethics Committees of Keio University and the University of Tokyo were followed.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information


Additional file 1: Figure S1. Gene ontology enrichment analysis of differentially expressed genes between females and males. Gene ontology terms enriched in genes higher expressed in [A] females and [B] males.


Additional file 2: Figure S2. Tra2 genes of P. metropolitanus. Gene structure and RNA-seq-based evidence of intronic regions. From the upper row, (1) the gene structures predicted by Hara et al., (2–4) RNA-Seq read mapping of male samples, (5–8) RNA-Seq read mapping of female samples, and (9) EvidenceModeler and PASA expanded gene structure. These structures were visualized using Jbrowse2 instance. The blue and red arrows indicate the male- and female-biased tra2 splicing variants.


Additional file 3: Figure S3. Expression of DMRT orthologs in H. exemplaris and R. varieornatus. Error bars indicate the standard deviation. On the X-axis, E and B time points indicate #day after oviposition (embryo) and #days after hatching (baby), and adults (active and tun).


Additional file 4: Figure S4. Structures of DMRT orthologs. [AB] Multiple alignment of the P. metropolitanus 3090/3093 complex orthologs [A] Amino acid sequences aligned by ClustalO and visualized by Seaview [B] Nucleotide sequences for the whole gene sequence (exon+intron) aligned by LAST through the MAFFT web site ( The red line indicates the matching regions between the two sequences. [CDE] AlphaFold2 predicted the 3D structure of [C] full-length [D] DM domain, and [E] the CUE-DMA domain. The arrowheads in cyan and magenta indicate the DM and CUE-DMA domains, respectively. Dm indicates D. melanogaster.

Additional file 5: Table S1. Dmrt gene accession IDs and sequences used in the phylogenetic analysis.

Additional file 6: Table S2. Primer sequences for genotyping.

Additional file 7: Table S3. Summary of NGS data.

Additional file 8: Table S4. GO enrichment analysis for genes within female-biased bins.

Additional file 9: Table S5. GO enrichment analysis for genes within male-biased bins.

Additional file 10: Table S6. Statistics of genome reassemblies.

Additional file 11: Table S7. Expression profiles for anhydrobiosis related genes.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Sugiura, K., Yoshida, Y., Hayashi, K. et al. Sexual dimorphism in the tardigrade Paramacrobiotus metropolitanus transcriptome. Zoological Lett 10, 11 (2024).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: