Article The molecular evolution of spermatogenesis across mammals ­­­ https://doi.org/10.1038/s41586-022-05547-7 Florent Murat1,2,19 ✉, Noe Mbengue1,19 ✉, Sofia Boeg Winge3,4,5, Timo Trefzer6, Evgeny Leushkin1, Mari Sepp1, Margarida Cardoso-Moreira7, Julia Schmidt1, Received: 2 November 2021 Celine Schneider1, Katharina Mößinger1, Thoomke Brüning1, Francesco Lamanna1, Accepted: 9 November 2022 Meritxell Riera Belles5, Christian Conrad6, Ivanela Kondova8, Ronald Bontrop8, Rüdiger Behr9,10, Philipp Khaitovich11, Svante Pääbo12, Tomas Marques-Bonet13,14,15,16, Published online: xx xx xxxx Frank Grützner17, Kristian Almstrup3,4,18, Mikkel Heide Schierup5 & Henrik Kaessmann1 ✉ Open access Check for updates The testis produces gametes through spermatogenesis and evolves rapidly at both the morphological and molecular level in mammals1–6, probably owing to the evolutionary pressure on males to be reproductively successful7. However, the molecular evolution of individual spermatogenic cell types across mammals remains largely uncharacterized. Here we report evolutionary analyses of single- nucleus transcriptome data for testes from 11 species that cover the three main mammalian lineages (eutherians, marsupials and monotremes) and birds (the evolutionary outgroup), and include seven primates. We find that the rapid evolution of the testis was driven by accelerated fixation rates of gene expression changes, amino acid substitutions and new genes in late spermatogenic stages, probably facilitated by reduced pleiotropic constraints, haploid selection and transcriptionally permissive chromatin. We identify temporal expression changes of individual genes across species and conserved expression programs controlling ancestral spermatogenic processes. Genes predominantly expressed in spermatogonia (germ cells fuelling spermatogenesis) and Sertoli (somatic support) cells accumulated on X chromosomes during evolution, presumably owing to male-beneficial selective forces. Further work identified transcriptomal differences between X- and Y-bearing spermatids and uncovered that meiotic sex-chromosome inactivation (MSCI) also occurs in monotremes and hence is common to mammalian sex-chromosome systems. Thus, the mechanism of meiotic silencing of unsynapsed chromatin, which underlies MSCI, is an ancestral mammalian feature. Our study illuminates the molecular evolution of spermatogenesis and associated selective forces, and provides a resource for investigating the biology of the testis across mammals. The rapid evolution of the testis across mammals is probably mainly tend to be predominantly expressed in the testis and thus probably explained by positive selection associated with sperm competition, also contribute to its rapid phenotypic evolution3,9. which reflects the evolutionary pressure on males to achieve repro- The testis also shows several other unique molecular features. First, ductive success7. Consequently, testis sizes, sperm production rates, chromatin in spermatogenic cells is massively remodelled during sper- sperm morphologies and other cellular traits substantially vary across matogenesis, a process that culminates in the tight packaging of DNA mammals, even between closely related species such as the great apes, around protamines in the compact sperm head10. This remodelling due to great mating system differences, especially regarding the extent leads to widespread leaky transcription in the genome11, which in turn of female promiscuity7. The rapid evolution of the testis is reflected at probably facilitates the initial transcription and, hence, the frequent the molecular level. Previous gene expression comparisons for various emergence of new testis-expressed genes and alternative exons dur- organs across mammals revealed that rates of evolutionary expres- ing evolution3,9,11,12. Second, the differentiation of sex chromosomes sion change are highest in the testis, probably due to frequent adap- from ancestral autosomes triggered the emergence of MSCI in euthe- tive changes but potentially also widespread relaxation of purifying rians and marsupials13 (therians), which led to the establishment of selection1–6. Consistently, genes with testis-specific expression tend to backup gene copies that substitute for parental genes on the X during be enriched with genes whose coding sequences have been shaped by meiosis3,14. In spite of MSCI, the X chromosome has become enriched positive selection8. In addition, new genes that emerge during evolution with testis-expressed genes during evolution3,15–20, presumably due A list of affiliations appears at the end of the paper. Nature | www.nature.com | 1 Article a Human Chimp. Bonobo Gorilla Gibbon Macaque Marmoset Mouse Opossum Platypus Chicken Other somatic Sertoli 45 90 Lept. SC Zyg. SC Dif. SG 180 Undif. SG Pach./Dipl. SC Other/mixture Late round SD Elongated SD Early round SD 310 UMAP Million years ago b c Cell type 50 Sertoli chicken.1 Other somatic chicken.2 Chicken platypus.1 SG platypus.2 SC opossum.1 PC2 (8% variance explained) rSD opossum.2 eSD opossum.3 25 Platypus mouse.1 Species mouse.2 marmoset.1 Human marmoset.2 Opossum Chimp. gibbon.1 macaque.1 Mouse Bonobo macaque.2 0 Gorilla human.1 Gibbon human.2 gorilla.1 Macaque gorilla.2 Marmoset <0.9 chimp.1 Mouse <0.5 chimp.3 Opossum chimp.2 Primates bonobo.1 –25 Platypus 0.1 bonobo.2 Chicken –20 0 20 40 PC1 (11% variance explained) Fig. 1 | snRNA profiling across ten mammals and a bird. a, Species sampled Chimp., chimpanzee. b, Principal component (PC) analysis of cell-type pseudo- and uniform manifold approximation and projection (UMAP) of snRNA-seq bulks. Species and lineages are encircled by a dashed line. Each symbol datasets. UMAP of the integrated primate dataset (dashed box), showing represents an individual. c, Gene expression phylogeny based on pseudo-bulk undifferentiated and differentiated SG (undif. SG and dif. SG, respectively), transcriptomes for whole testes. Bootstrap values (4,498 1:1 orthologous amniote leptotene, zygotene, pachytene and diplotene SCs (lept. SC, zyg. SC, pach. genes were randomly sampled with replacement 1,000 times) are indicated SC and dipl. SC, respectively), spermatids (SD) and somatic cell types. by circles, ≥0.9 (white fill). to sexually antagonistic selective forces favouring the fixation of representing all simian (anthropoid) lineages (Fig. 1a): eutherian mam- male-beneficial mutations on this chromosome21. Finally, translational mals (representatives for five of the six extant ape lineages, including regulation of transcriptomes is widespread across spermatogenesis6. humans; rhesus macaque, an Old World monkey; common marmoset, Previous large-scale transcriptomic investigations of testis evolution a New World monkey; and mouse), marsupials (grey short-tailed opos- were largely limited to bulk-organ samples1–6,19. Recent high-throughput sum) and egg-laying monotremes (platypus). Corresponding data single-cell (sc) or single-nucleus (sn) RNA-sequencing (RNA-seq) tech- were generated for a bird (red jungle fowl, the progenitor of domestic nologies enable detailed investigations of the cellular and molecular chicken; hereafter referred to as ‘chicken’), to be used as an evolution- evolution of the testis, as exemplified by two scRNA-seq comparisons ary outgroup. The dataset consists of 27 libraries, with one to three between human, macaque and mouse22,23, but a comprehensive investi- biological replicates per species and a median of roughly 275 million gation of the evolution of spermatogenesis across all main mammalian snRNA-seq reads per library (Supplementary Table 1). We refined and lineages is lacking. extended existing genome annotations across all species on the basis of Here we provide an extensive snRNA-seq resource covering testes bulk-testis RNA-seq data (seven libraries) (Supplementary Tables 1 and 2 from ten representative mammals and a bird, enabling detailed com- and Methods), to ensure optimal read-mapping and prevent biases parisons of spermatogenic cells and underlying gene expression in cross-species analyses. After quality controls and filtering steps programs within and across mammals (https://apps.kaessmannlab. (Methods), we obtained transcriptomes for a total of 97,521 high-quality org/SpermEvol/). Our evolutionary analyses of these data unveiled nuclei for the 11 species, with a mean of 8,866 cells per species, a median ancestral as well as species- and lineage-specific cellular and molecular of 1,856 RNA molecules (unique molecular identifiers (UMIs)) detected characteristics of mammalian spermatogenesis. per cell and low percentages of mitochondrial UMIs (Supplementary Fig. 1 and Supplementary Table 1). We identified the main germ cell types11 along the continuous cellular Spermatogenesis across 11 species proliferation and differentiation path of spermatogenesis across all We generated snRNA-seq data for testes from ten species that cover species: spermatogonia (SG), the mitotic cells fuelling spermatogen- the three main mammalian lineages and include key primate species, esis, including spermatogonial stem cells; spermatocytes (SC), where 2 | Nature | www.nature.com meiosis takes place; and the haploid round spermatids (rSD) and elon- cellular level the notion that the conservation of core spermatogenic gated spermatids (eSD), which together reflect spermiogenesis (Fig. 1a, functions restricts transcriptome divergence. Extended Data Figs. 1a–f, Supplementary Table 3 and Methods). We also identified separate clusters corresponding to somatic testicular cells, in particular Sertoli cells (except for platypus, in which Sertoli cells Evolutionary forces could not be unambiguously distinguished), the main spermatogenesis We sought to trace the evolutionary forces underlying the rapid evolu- support cells, but also other cell types, such as Leydig cells, peritubular tion of late spermatogenesis. Two non-mutually exclusive patterns of cells, endothelial cells and macrophages (Extended Data Figs. 1a–e). natural selection may account for this observation. First, later stages The close evolutionary relationship of the seven primates in our study of spermatogenesis might evolve under weaker purifying selection enabled the direct integration of datasets across these species and thus (that is, reduced functional constraints) and hence be less refractory the identification of sub-cell types that correspond to intermediate to change. Second, the greater divergence in later stages might result events during spermatogenesis (Fig. 1a). For all species, we traced the from stronger positive selection, increasing the rate of fixation of adap- dynamic gene expression programs underlying spermatogenic dif- tive changes. To investigate patterns of functional constraint during ferentiation and key molecular events, thus also identifying a host of spermatogenesis, we assessed the tolerance to functional mutations new marker genes (Extended Data Fig. 1g and Supplementary Table 4). of genes24 used in different spermatogenic stages in humans, which To obtain an overview of cell-type relationships across species, we showed a progressive increase of mutational tolerance starting during performed a principal component analysis (PCA) based on pseudo-bulk meiosis and culminating in early spermiogenesis (Fig. 2c and Extended cell-type transcriptomes (Fig. 1b). The first principal component (PC1) Data Fig. 3a). Consistently, on the basis of a set of neutrally ascertained orders the spermatogenic cell types according to the progression of mouse knockouts25, we found that the percentage of expressed genes spermatogenesis for all species (Fig. 1b, from left to right). This observa- associated with lethality decreases during spermatogenesis (Fig. 2d). tion suggests that our data capture ancestral aspects of spermatogenic Also, in agreement with a progressive reduction in functional con- gene expression programs that are shared across mammals or amni- straints towards later spermatogenic stages, we find that the normal- otes, despite the rapid evolution of the testis1–6 and the long divergence ized rate of amino acid altering substitutions in coding sequences times of 310 million years (Fig. 1a). PC2 separates the data by species or across primates is higher in late spermatogenesis (Fig. 2e and Extended lineages, reflecting diverged aspects of spermatogenesis, whereas PC1 Data Fig. 3b), although this increase might additionally reflect a higher and PC2 together separate somatic and spermatogenic cell types for proportion of genes under positive selection. Indeed, an examination each species. The close clustering of biological replicates is a further of the temporal expression pattern of genes whose encoded protein indicator of the high data quality. sequences have been shaped by positive selection revealed a notable increase in percentages of positively selected genes used during sper- matogenesis, with a peak in rSD (Fig. 2f and Extended Data Fig. 3c). Rates of evolution along spermatogenesis Because new genes also contribute to evolutionary innovations, we To investigate rates of gene expression evolution across cell types, investigated the temporal contribution of recently emerged genes we reconstructed gene expression trees (Methods). A tree based on to gene expression programs in germ cells, using an index that com- pseudo-bulk transcriptomes for the whole testis (Fig. 1c) recapitu- bines the phylogenetic age of genes with their expression2 (Methods). lates the known mammalian phylogeny (Fig. 1a, except for the gib- This analysis revealed that transcriptomes become younger during bon–macaque grouping), akin to previous trees based on bulk-tissue spermatogenesis (Fig. 2g and Extended Data Fig. 3d), indicating that RNA-seq data across mammalian organs1,2. This observation is consist- new genes have increasingly more prominent roles in later stages, ent with the view that regulatory changes steadily accumulated over particularly in rSD, consistent with previous observations22,23. Pre- evolutionary time1, with present-day RNA abundances reflecting the vious work based on bulk cell-type analyses in mouse11 uncovered a evolution of mammalian lineages and species. transcriptionally permissive chromatin environment during spermato- To trace the cellular source of the rapid evolution of the testis, we genesis, in particular in rSD, which was suggested to have facilitated built expression trees for the different cell types, which also recapitu- the emergence of new genes during evolution3,9,11. Consistently, we late the known species relationships (Extended Data Fig. 2a). Notably, detect in all species considerably increased contributions of inter- the total branch lengths of the trees, which reflect the amount of genic transcripts after meiosis and a concomitant decrease in the evolutionary expression change, vary substantially between cell contributions of protein-coding genes (Fig. 2h and Extended Data types (Fig. 2a). Whereas the rate of expression evolution is similar Fig. 3e). Notably, an analysis of translatome data6 revealed a decline in Sertoli cells and diploid spermatogenic cells (and lower than that of translational efficiencies of transcripts during spermatogenesis in other somatic cell types), it is substantially higher in the post- (reaching a minimum in rSD) in all species (Fig. 2i and Extended Data meiotic haploid cell types (rSD and eSD), consistent with a recent Fig. 3f). This decline is consistent with observations from mouse bulk inference based on data for three eutherians23. The higher resolu- data for a restricted number of cell types6 and probably mitigates the tion afforded by a primate-specific analysis provides further details functional consequences of the concurrent increase in transcriptional (Fig. 2a). Starting in late meiosis (pachytene SC), evolutionary rates promiscuity of the genome. progressively increase until the end of spermiogenesis (late eSD). We next explored the reasons underlying the dynamic changes of Thus, late spermatogenic stages drive the previously observed rapid selective forces and patterns of innovation during spermatogenesis. evolution of the testis1–6. The breadth of expression across tissues and developmental processes Pairwise species comparisons, including downsampling analyses (here referred to as expression pleiotropy) was proposed to represent (Extended Data Fig. 2b), confirm the rapid expression evolution of a key determinant of the types of mutation that are permissible under postmeiotic cell types across amniotes and that gene expression diver- selection26. We, therefore, assessed patterns of expression pleiotropy gence increases with evolutionary time (Fig. 2b), in accord with the across spermatogenesis using spatiotemporal transcriptome data for expression phylogeny results (Extended Data Fig. 2a). It is notewor- several mammalian organs2, which revealed that genes used later in thy, however, that expression divergence levels are roughly as similar spermatogenesis, in particular those in rSD, have substantially more between human and chicken as they are between human and platypus, specific spatiotemporal profiles than genes used earlier in spermato- although the bird lineage diverged 110 million years before the separa- genesis and in somatic cells (Fig. 2j and Extended Data Fig. 3g,h). Given tion of monotremes and therians (that is, eutherians and marsupials). that a decrease in expression pleiotropy can explain both a decrease This observation, previously made for whole organs1, supports at the in functional constraints and an increase in adaptation2,26, we suggest Nature | www.nature.com | 3 Article a b ST 1.8 Amniotes Primates 0.8 OS SG 0.8 SC 1.6 0.7 rSD Total tree length Total tree length eSD Spearman’s U 1.4 0.6 0.6 1.2 0.5 0.4 1.0 0.4 Divergence time (Ma) 0 100 200 300 OS ST SG SC rSD D OS ST G G SC Pa SC Ea . SC La SD rSD D eS eS .S .S Primates Mouse Opossum Platypus Chicken r pt. g. ch dif Dif rly te Zy Le Un c d e 50 0.300 45 40 0.275 Lethals (%) pLI percentile dN/dS 30 40 0.250 20 0.225 35 10 Spermatogenesis Spermatogenesis Spermatogenesis OS OS OS ST ST ST f g h Transcriptome phylogenetic age 1.75 5 Percentage of the transcripome Positive selection (%) 1.50 75 4 3 1.25 Protein coding genes 50 1.00 Intergenic elements 2 0.75 25 1 0 0.50 Spermatogenesis Spermatogenesis Spermatogenesis OS OS OS ST ST ST i j k 5.0 0.50 20 Translational efficiency Tissue specificity 2.5 0.45 15 Infertile (%) 0.40 0 10 0.35 −2.5 5 0.30 −5.0 0 Spermatogenesis Spermatogenesis Spermatogenesis OS OS OS ST ST ST Fig. 2 | Gene expression divergence and evolutionary forces. a, Total branch f, Percentage of expressed genes under positive selection (out of 11,170 genes lengths of expression trees among testicular cell types for amniotes and tested for positive selection) in chimpanzee. g, Mean phylogenetic age of primates. Box plots show the median (central value); upper and lower quartile expressed genes in mouse. h, Percentage of UMIs mapping to protein-coding (box limits) and 95% confidence intervals (whiskers) for 1,000 bootstrap genes (top) or intergenic elements (bottom) in gorilla. i, Translational efficiency replicates. b, Spearman’s correlations between humans and other species values (data from ref. 6) are plotted for all genes with predominant expression in from 100 bootstrap replicates (dots). Lines correspond to linear regression a given cell type in human. j, Mean of tissue-specificity values (data from ref. 2) trends (after log transformation of the time axis. Regression R 2 values range in opossum. k, Percentage of expressed genes associated with infertility from 0.86 to 0.97. Ma, millions of years ago. c, Mean pLI value of expressed (out of 3,552 knockouts) in mouse. c–h,j,k, Plotted is the mean value per cell. genes in human (the pLI score reflects the tolerance of a gene to a loss-of- c–k, Superimposed thick black dots indicate medians from biological replicates. function mutation; lower values mean less tolerance). d, Percentage of Box plots depict the median (centre value); upper and lower quartile (box limits) expressed genes leading to a lethal phenotype when knocked out in mouse with whiskers at 1.5 times the interquartile range. Red lines separate somatic (out of 4,742 knockouts25). e, Mean normalized ratio of nonsynonymous (dN) (OS, other somatic; ST, Sertoli cells) and germ cells. Data for other studied over synonymous (dS) nucleotide substitutions of expressed genes in macaque. species are shown in Extended Data Fig. 3. that it is probably a main contributor to the accelerated molecular associated with infertility25 is relatively high in SC and spermatids evolution in late spermatogenesis. In addition, the specific type of selec- (especially rSD): higher than in SG and somatic cells (Fig. 2k). tion acting on haploid cells27 (haploid selection), in which expressed alleles are directly exposed to selection, may have contributed to the exceptionally rapid evolution of rSD. Gene expression conservation and innovation Whereas the tissue- and time-specific late spermatogenic genes, We next sought to trace the individual genes underlying conserved in general, are not essential for viability (Fig. 2d,j and Extended Data (ancestral) and diverged aspects of germ cells by comparing expression Fig. 3g,h, above), we proposed that the specific aforementioned evolu- trajectories along spermatogenesis of one-to-one (1:1) orthologous tionary forces indicate that many of these genes evolved crucial roles genes across species2 (Supplementary Fig. 2 and Methods). For the pri- in spermatogenesis. Indeed, we find that the proportion of genes mates, this analysis revealed roughly 1,700–2,900 genes with conserved 4 | Nature | www.nature.com a TEX11 PI4KB b DMRT1 IP6K1 63 101 * 23 30 512 2,604 * 22 46 153 39 2,939 46 16 2,479 26 472 Expression levels Expression levels 34 62 75 1,881 416 207 2,001 92 389 106 1,687 SG SC D D SG SC D D SG SC D D SG SC D D rS eS rS eS rS eS rS eS c snRNA-seq (RUBCNL) d smISH (RUBCNL) SD Normlaized staining 0.5 0.4 * SC 0.3 0.2 SG SC eSD D rS 0.1 SG SG SC SD Expression levels Normlaized staining SD 0.8 SG SC eSD D rS 0.6 SC 0.4 0.2 SG SC eSD D rS SG SG SC SD SD 1.00 Normlaized staining SG SC eSD D rS 0.75 snRNA-seq data 0.50 SC not available 0.25 SG 0 SG SC SD 50 μm 10 μm Fig. 3 | Evolution of gene expression trajectories along spermatogenesis. 2,927 1:1 orthologues. DMRT1 and IP6K1 are examples of conserved and a,b, The numbers of changed trajectories (in purple) and conserved changed trajectories, respectively. c, RUBCNL expression trajectory from trajectories (in olive) are indicated. In a–c, red asterisks indicate the branch snRNA-seq data along spermatogenesis for human, bonobo, chimpanzee for which a trajectory change has been called. For each replicate, mean and gorilla (from top to bottom). Available orangutan samples did not allow expression levels across cells of a given cell type were calculated. Marks for the generation of high-quality snRNA-seq data. d, Detection of RUBCNL indicate values for the replicates with the highest and lowest mean expression (each red dot reflects a single transcript) expression in human, chimpanzee levels, dots indicate the median of mean expression values of the replicates and orangutan testis by smISH using RNAScope. Left, seminiferous tubule (three replicates for chimpanzee and opossum; two for human, bonobo, gorilla, cross sections counterstained with haematoxylin, and closeups on areas macaque, marmoset, mouse, platypus and chicken; one for gibbon). In a–c, containing SG, SC and spermatids (SD). Right, quantification of RUBCNL dashed vertical lines separate SG, SCs, rSD and eSD. a, Primate trajectories expression levels in ten tubules per section (n = 3 for human, n = 1 for (human, bonobo, chimpanzee, gorilla, gibbon, macaque and marmoset; from chimpanzee and orangutan). Box plots represent the mean (central value) top to bottom) based on 4,459 1:1 orthologues. TEX11 and PI4KB are examples distribution of staining (dots); upper and lower quartile (box limits) with of conserved and changed trajectories, respectively. b, Amniote trajectories whiskers at 1.5 times the interquartile range. (human, mouse, opossum, platypus and chicken; from top to bottom) based on expression trajectories across different lineages or species (Fig. 3a and expression of the gene RUBCNL, encoding a regulator of autophagy28, Supplementary Table 5). For example, the temporal expression of 1,687 changed in the human lineage, with a relative reduction of expression genes is conserved across the seven primates and probably reflects the levels in spermatids (Fig. 3c,d). Thus, the role of RUBCNL in autophagy core ancestral gene expression program of the simian testis. during spermatogenesis29 is potentially different in humans compared By contrast, we detected 635 trajectory changes during the evo- to other Old World anthropoids. We also confirmed the inversion of lution of apes and Old World monkeys (Fig. 3a and Supplementary the expression trajectory of the myosin-encoding gene MYO3B in SG Table 5). For example, 63 and 94 trajectory changes occurred in the and rSD on the chimpanzee–bonobo lineage (Extended Data Fig. 4a), human and chimpanzee–bonobo lineages, respectively, since their as well as the expression increase of ADAMTS17, a family-member of divergence roughly 7 million years ago. We spatially validated three of proteases with key spermatogenic functions30, in rSD relative to SG these changes using single-molecule RNA in situ hybridization (smISH) in humans (Extended Data Fig. 4b). We note that observed quantita- for three great apes (human, chimpanzee and orangutan as the out- tive and partly qualitative expression pattern differences between group) (Fig. 3c, Extended Data Fig. 4, Supplementary Fig. 3, Supplemen- the two complementary data types are expected because of technical tary Table 6 and Methods). Our smISH experiments confirm that the differences. Nature | www.nature.com | 5 Article We also assessed the conservation of expression trajectories across the three main mammalian lineages and amniotes (Fig. 3b and Supple- Sex chromosomes mentary Table 7). The 416 genes with conserved expression across all Sex chromosomes emerged twice in parallel in mammals from differ- mammals probably trace back to ancestral gene expression programs of ent sets of ancestral autosomes. The therian XY chromosome system the ancestral mammalian testis, whereas 389 genes may represent the originated just before the split of eutherians and marsupials (Fig. 4a) core ancestral spermatogenic program of amniotes. A notable example and hence has evolved largely independently in these two lineages. is DMRT1, which is highly expressed specifically in SG across amniotes Around the same time, the monotreme sex-chromosome system arose and is required for mouse spermatogonial stem cell maintenance and from a different pair of autosomes and subsequently expanded to five replenishment31. XY pairs40. These sex-chromosome formation events entailed sub- In agreement with these highly conserved sets of genes having key stantial remodelling of gene contents and expression patterns due to roles in spermatogenesis, our analyses of fertility phenotypes25 unveiled structural changes and sex-related selective forces3. We used our data that genes involved in fertility are significantly more conserved in to systematically assess testicular expression patterns of sex chromo- their expression trajectories than genes not associated with fertility somal genes and their evolution across mammals. (Extended Data Fig. 5b). Consistently, a Gene Ontology32 enrichment We first contrasted cell-type specificities of X-linked and autosomal analysis indicates an involvement of conserved genes in fundamen- genes, which showed a notable excess (12–60%) of X-linked genes with tal spermatogenic processes that are typical of the cell type in which predominant expression in SG across all eutherians, in agreement they show peak expression (Extended Data Fig. 6a,b). Thus, genes with with a previous mouse study17 and includes conserved genes with key conserved trajectories for which spermatogenesis functions remain spermatogenic functions such as TEX11 (Figs. 3a and 4b, Extended uncharacterized represent promising candidates for the exploration of Data Fig. 9a and Supplementary Table 10). The higher-resolution fertility phenotypes (Supplementary Tables 5, 7 and 8). Notably, genes primate data showed that the X chromosome is also enriched for with lineage-specific trajectory changes are enriched with broader, genes expressed in leptotene SC (Fig. 4b), presumably reflecting tran- typically metabolic, processes (Extended Data Fig. 6c,d). They are also script carry-over from SG, given the global transcriptional silencing significantly more tissue- and time-specific than genes with conserved of the genome during the leptotene stage41 (Supplementary Fig. 1b expression (Extended Data Fig. 6e), which may have facilitated their and Supplementary Table 3) and the observation that most X-linked expression change during evolution because of reduced pleiotropic genes expressed in leptotene SC are also expressed in differentiated constraints. However, genes with changed trajectories neverthe- SG (Supplementary Fig. 4). Notably, we detected enrichments of genes less include many genes for which key fertility functions have been with SG-specific expression also on the opossum and platypus X chro- described (Supplementary Tables 7 and 8). For example, IP6K1, which mosomes, thus revealing this to be a shared pattern across mammals elicits infertility when knocked out in the mouse33, shows strongly (Fig. 4b and Extended Data Fig. 9a,b). increasing expression towards the end of spermiogenesis in all amni- We also uncovered an enrichment of genes with predominant otes except the primates, for which expression is high in SG and then expression in Sertoli cells across X chromosomes (except for platypus, declines (Fig. 3b and Extended Data Fig. 5c). Thus, the primary function for which Sertoli cells could not be unambiguously distinguished) of IP6K1 probably shifted from late to early spermatogenesis during (Fig. 4b). Altogether, our observations suggest that sex-related forces primate evolution. favoured the independent accumulation of SG- and Sertoli-specific In conjunction with mouse fertility data, we used our data to explore genes on the X during the evolution of the different sex-chromosome the contribution of new genes (mostly arising from gene duplica- systems. Consistently, autosomes in outgroup species corresponding tions) to the evolution of new spermatogenic functions, focusing to the different mammalian X chromosomes (for example, platy- on the rodent lineage leading to mouse. This analysis revealed the pus chromosome 6, which is homologous to the therian X) do not emergence of key spermatogenic genes at different time points during show any excess of SG- or Sertoli-expressed genes (Extended Data evolution (Extended Data Fig. 7a and Supplementary Table 8), such as Fig. 9a,c), which means that the ancestral autosomes that gave rise two rodent-specific retrogenes—D1Pas1 and H2al2a—that originated to present-day sex chromosomes were not enriched for such genes. through RNA-based duplication from parental genes on the X chromo- Overall, we propose that the accumulation of SG- and Sertoli-specific some (Extended Data Fig. 7c). D1Pas1 is essential for meiosis34, whereas genes was facilitated by the specific selective environment on the H2al2a is essential for genome compaction in late spermatogenesis35. X chromosome, in which male-beneficial mutations are always visible Both genes show strongly increasing expression levels in late sper- to selection because of the single-copy (hemizygous) status of the matogenesis (Extended Data Fig. 7b) and are thus in agreement with X chromosome in males3,17,21. new genes contributing to functional roles predominantly during late We next sought to separate X- and Y-bearing spermatids to inves- spermatogenesis (above, Fig. 2g and Extended Data Fig. 3d). tigate their distinct transcriptomal properties during spermiogen- Finally, we used our data to investigate ligand–receptor interac- esis. Such an analysis is probably not possible using single-whole-cell tions underlying the communication between Sertoli cells, which have transcriptomic data because X and Y spermatids remain connected a central role in supporting and controlling spermatogenesis11, and by cytoplasmic bridges and hence are thought to contain similar cyto- germ cells across species (Methods). Our cross-species comparisons plasmic transcript pools42. However, our single-nucleus data should revealed various conserved known and new ligand–receptor interac- afford the separation of X and Y spermatids. Indeed, on the basis of tions (Extended Data Fig. 8a,b and Supplementary Table 9). For exam- differential X and Y transcript contents (Methods), we were able to ple, our data provide evidence that communication of Sertoli cells separate spermatids into distinct X and Y lineages across mammals with SG, SC and rSD occurs in all amniotes through interactions of (Fig. 4c and Supplementary Table 3). As expected, our approach failed the cell adhesion molecule CADM1 (ref. 36) or between CADM1 (Sertoli to separate X and Y spermatids in available human23 and mouse43 cells) and NECTIN3 (SC and rSD) (Supplementary Table 9). We note scRNA-seq datasets (Extended Data Fig. 10), supporting the notion that CADM1 was previously thought to not to be expressed in Sertoli of substantial transcript exchange across X- and Y-spermatid cells cells37, but our data show high expression in this cell type across amni- through cytoplasmic bridges42, although this equilibration may not otes, in agreement with human protein atlas data (proteinatlas.org/ be complete44. ENSG00000182985-CADM1/tissue/testis#). Our work also supports A differential expression analysis between X and Y spermatids the notion that the NECTIN2–NECTIN3 complex mediates the com- identified, as expected, most sex-chromosome genes, including munication of Sertoli cells with spermatids not only in mouse38 but gametologues (that is, genes with homologous counterparts on also in humans39 (Supplementary Table 9). X and Y chromosomes), such as the translational regulatory genes 6 | Nature | www.nature.com a c **** **** d n = 1,005 2.0 Percentage of Percentage of –1.6 Y transcripts X transcripts log10(X/A) 1.5 X 2 n = 990 1.0 –2.0 1 0.5 –2.4 Y 0 90 Ma **** n = 1,237 **** –2 Percentage of 3 7.5 X transcripts Percentage of –3 Y transcripts log10(X/A) 2 n = 1,057 5.0 –4 XY 1 2.5 –5 159 Ma 0 –6 **** **** –1 n = 1,848 4 Percentage of Percentage of Y transcripts log10(X/A) X transcripts 2 3 –2 n = 1,245 2 177 Ma 1 –3 1 0 0 –4 X1-5 **** n = 2,270 **** Y1-5 17.5 Percentage of Percentage of 1.5 –0.8 X transcripts Y transcripts log10(X/A) 15.0 X5 n = 2,096 1.0 12.5 –1.0 Y5 10.0 0.5 X1 X2 X3 X4 7.5 0 –1.2 Y3 Y4 **** **** –0.8 12.5 9 Percentage of Percentage of Y1 X transcripts X transcripts Y2 from SDRs from PARs 10.0 8 log10(X/A) UMAP –1.2 7.5 7 X-bearing spermatids 6 –1.6 5.0 Y-bearing spermatids 5 PARs Non-spermatids or 2.5 SDRs –2.0 unassigned g g g g in rin in rin ar ar ea ea be be b b X- Y- X- Y- b 10 12 Percentage X-linked genes Percentage X-linked genes Percentage X-linked genes Percentage X-linked genes Percentage X-linked genes 9 4 6 25 *** 200 40 7 * **** 20 *** 60 *** 13 **** *** 40 20 5 30 ***** 15 5 30 40 15 20 10 * 6 20 10 3 2 13 14 10 10 1 20 5 5 10 10 18 5 1 5 2 1 1 1 20 2 0 0 0 0 0 0 0 0 0 0 0 0 0 NA S if. T Di SG p G Payg. C Each. SC La y rSC te D D D nd SS if. T Le if. SG p G Pa g. C Each. C La y rSC te D D D S ST SG SC D D S ST SG SC D D S SG SC D D nd S O O O O O rS eS rS eS rS eS rS eS rS eS Z t. S rl S Zy t. S S rl S Le f. S D S U U Fig. 4 | Mammalian sex-chromosome evolution. a, Phylogenetic tree for enriched in each cell type is indicated above each bar. c, UMAP representation human, mouse, opossum and platypus (from top to bottom). Arrows show sex of germ cells (left). Progression of spermatogenesis (grey arrow) and meiotic systems origination40. Illustration of therian (human, top) and monotreme divisions (black arrows) are indicated. Spermatids identified as X- and Y-bearing (platypus, bottom) sex chromosomes. Recombining (crosses) PARs in turquoise are coloured in orange and purple, respectively. Box plots show the median and SDRs in red. b, Percentages of X-linked genes among testis-specific genes (central value) percentages of X and Y transcripts in X- and Y-bearing spermatids; with predominant expression in a given cell type for human, macaque, mouse, upper and lower quartile (box limits) with whiskers at 1.5 times the interquartile opossum and platypus (from left to right; Supplementary Table 10). The red range. Two-sided Wilcoxon rank-sum tests were performed for statistical dashed line represents the expected percentages of X-linked genes, if testis- comparisons (Benjamini–Hochberg corrected ****P < 2.2 × 10 −16). d, X-to- specific genes with predominant expression in the different cell types were autosome transcript ratios in individual germ cells across spermatogenesis randomly distributed across the genome. Asterisks indicate significance after (from left to right). For spermatids, only X-bearing cells are considered. Lines two-sided exact binomial testing (Benjamini–Hochberg corrected P values from depict generative additive model trend. Red arrows indicate MSCI. c,d, For left to right: 0.000353, 1.18 × 10 −7, 0.000687, 0.0142, 1.08 × 10 −7, 3.45 × 10 −8, platypus, X transcripts are dissected according to their location on the 0.000792, 0.0355 and 0.00043). The number of testis-specific X-linked genes X chromosomes (Methods). DDX3X/DDX3Y (Extended Data Fig. 11 and Supplementary Table 11). differentiated regions (SDRs) versus pseudoautosomal regions (PARs), However, we also found autosomal genes, especially in some species which are large in monotremes (Fig. 4a). We proposed that expression (for example, human), which might reflect trans-regulatory effects signals from the large PARs, which are expected to synapse and hence associated with the X and Y chromosomes (Extended Data Fig. 11, Sup- not to be affected by MSUC, might have prevented the detection of plementary Table 11 and Methods). MSCI in previous studies. We then assessed gene expression across spermatogenesis separately Indeed, whereas the joint analysis of all platypus X-linked genes only for X- and Y-linked genes. We traced a substantial dip in X transcript abun- shows a small expression dip around the pachytene stage (Fig. 4d, dances around the pachytene stage of meiosis across therians (Fig. 4d), upper platypus graph), an analysis only of SDR genes reveals a strong which reflects the process of MSCI13, a sex-chromosome-specific reduction of X transcript levels. By contrast, PAR genes show stable instance of the general epigenetic phenomenon of meiotic silencing expression levels across spermatogenesis (Fig. 4d, lower graph). of unsynapsed chromatin13 (MSUC). An analysis of Y transcripts pro- Moreover, the difference in transcript abundances between SDRs vides consistent results but at lower resolution owing to the very small and PARs owing to MSCI is visible for all five platypus X chromosomes number of Y-linked genes (Extended Data Fig. 12a). Previous work did (Extended Data Fig. 12b). Notably, our assessment of the completeness not find evidence for MSCI in monotremes45, suggesting that MSCI of MSCI across species reveals that platypus MSCI is as complete as originated in the therian ancestor after the separation from the mono- that in other species for which there is little or no MSCI escape23,47,48 treme lineage45. We revisited this question in platypus using a new (Extended Data Fig. 12c). The presence of MSCI at the SDRs in mono- assembly and annotation46 that includes a detailed definition of sexually tremes is consistent with the partial association of platypus sex Nature | www.nature.com | 7 Article chromosomes with perinucleolar repressive histone modifications at the pachytene stage45. Altogether, our data reveal that efficient Online content MSCI is common to all mammalian sex-chromosome systems, which Any methods, additional references, Nature Portfolio reporting summa- indicates that the general mechanism of MSUC is an ancestral mam- ries, source data, extended data, supplementary information, acknowl- malian feature. edgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41586-022-05547-7. Discussion Our analyses uncovered that the previously observed rapid evolution 1. Brawand, D. et al. The evolution of gene expression levels in mammalian organs. Nature of the testis1–7 is driven by an accelerated rate of molecular innova- 478, 343–348 (2011). tion in late spermatogenesis, in particular in rSD. Our findings suggest 2. Cardoso-Moreira, M. et al. Gene expression across mammalian organ development. a scenario in which the accelerated fixation of regulatory changes, Nature 571, 505–509 (2019). 3. Necsulea, A. & Kaessmann, H. Evolutionary dynamics of coding and non-coding amino acid altering substitutions and new genes during evolution in transcriptomes. Nat. Rev. Genet. 15, 734–748 (2014). late spermatogenesis (presumably driven by sperm competition) was 4. Necsulea, A. et al. The evolution of lncRNA repertoires and expression patterns in facilitated by reduced pleiotropic constraints, a transcriptionally per- tetrapods. Nature 505, 635–640 (2014). 5. Sarropoulos, i, Marin, R., Cardoso-Moreira, M. & Kaessmann, H. Developmental dynamics missive chromatin environment and potentially haploid selection for of lncRNAs across mammalian organs and species. Nature 571, 510–514 (2019). spermatid genes that are not rendered effectively diploid through tran- 6. Wang, Z. Y. et al. Transcriptome and translatome co-evolution in mammals. Nature 588, script exchange through cytoplasmic bridges. In agreement with late 642–647 (2020). 7. Ramm, S. A., Scharer, L., Ehmcke, J. & Wistuba, J. Sperm competition and the evolution of spermatogenic stages explaining the rapid evolution of the testis, early spermatogenesis. Mol. Hum. Reprod. 20, 1169–1179 (2014). spermatogenic cell types and somatic cells show patterns of constraint 8. Nielsen, R. et al. A scan for positively selected genes in the genomes of humans and and innovation similar to those of cell types in the brain (Extended Data chimpanzees. PLoS Biol. 3, e170 (2005). 9. Kaessmann, H. Origins, evolution, and phenotypic impact of new genes. Genome Res. Fig. 3a,b,d), an organ that evolves slowly at the molecular level1–3. We 20, 1313–1326 (2010). note that differences in cell-type abundances across mammals, which 10. Kimmins, S. & Sassone-Corsi, P. Chromatin remodelling and epigenetic features of germ are pronounced for the testis7 (Extended Data Fig. 1f), presumably also cells. Nature 434, 583–589 (2005). 11. Soumillon, M. et al. Cellular source and mechanisms of high transcriptome complexity in contributed to the rapid gene expression divergence of this organ6. the mammalian testis. Cell Rep. 3, 2179–2190 (2013). Our cross-species comparisons of individual genes revealed tempo- 12. Mazin, P. V., Khaitovich, P., Cardoso-Moreira, M. & Kaessmann, H. Alternative splicing ral expression differences across species, which were probably facili- during mammalian organ development. Nat. Genet. 53, 925–934 (2021). 13. Turner, J. M. Meiotic silencing in mammals. Annu. Rev. Genet. 49, 395–412 (2015). tated by reduced pleiotropic constraints. Our results thus provide an 14. Potrzebowski, L. et al. Chromosomal gene movements reflect the recent origin and extensive list of candidates whose contributions to the evolution of biology of therian sex chromosomes. PLoS Biol. 6, e80 (2008). species-specific spermatogenesis phenotypes can be experimentally 15. Mueller, J. L. et al. The mouse X chromosome is enriched for multicopy testis genes showing postmeiotic expression. Nat. Genet. 40, 794–799 (2008). scrutinized. We also uncovered conserved expression programs under- 16. Mueller, J. L. et al. Independent specialization of the human and mouse X chromosomes lying spermatogenic processes ancestral to individual mammalian for the male germ line. Nat. Genet. 45, 1083–1087 (2013). 17. Wang, P. J., McCarrey, J. R., Yang, F. & Page, D. C. An abundance of X-linked genes lineages and mammals as a whole. expressed in spermatogonia. Nat. Genet. 27, 422–426 (2001). Further analyses illuminated the role of sex chromosomal genes in 18. Julien, P. et al. Mechanisms and evolutionary patterns of mammalian and avian dosage spermatogenesis. We found that genes predominantly expressed in compensation. PLoS Biol. 10, e1001328 (2012). 19. Meunier, J. et al. Birth and expression evolution of mammalian microRNA genes. Genome SG and Sertoli cells independently accumulated on X chromosomes Res. 23, 34–45 (2013). across mammals and their two sex-chromosome systems during evolu- 20. Zhang, Y. E., Vibranovski, M. D., Landback, P., Marais, G. A. & Long, M. Chromosomal tion. This suggests that X chromosomes have been shaped by strong redistribution of male-biased genes in mammalian evolution with two bursts of gene gain on the X chromosome. PLoS Biol. 8, e1000494 (2010). male-related selective forces21, leading to the emergence of X-linked 21. Rice, W. R. Sex chromosomes and the evolution of sexual dimorphism. Evolution 38, genes with functional roles in testis cell types in which active transcrip- 735–742 (1984). tion is possible. Indeed, in addition to the SG and Sertoli-cell enrichment 22. Lau, X., Munusamy, P., Ng, M. J. & Sangrithi, M. Single-cell RNA sequencing of the cynomolgus macaque testis reveals conserved transcriptional profiles during found here, previous work showed that these selective pressures led mammalian spermatogenesis. Dev. Cell 54, 548–566.e7 (2020). to the repeated duplication of genes on the X that facilitated expres- 23. Shami, A. et al. Single-cell RNA sequencing of human, macaque, and mouse testes sion after meiosis15. uncovers conserved and divergent features of mammalian spermatogenesis. Dev. Cell 54, 529–547 e12 (2020). Our ability to separate X- and Y-bearing spermatids and the avail- 24. Dickinson, M. E. et al. High-throughput discovery of novel developmental phenotypes. ability of a new platypus genome assembly46 unveiled that MSCI is a Nature 537, 508–514 (2016). general feature of mammalian sex-chromosome systems, indicating 25. Koscielny, G. et al. The International Mouse Phenotyping Consortium Web Portal, a unified point of access for knockout mice and related phenotyping data. Nucleic Acids that MSUC was already present in the common mammalian ancestor. Res. 42, 802–809 (2014). Previous work did not find evidence for MSCI in birds49, which raises 26. Carroll, S. B. Evolution at two levels: on genes and form. PLoS Biol. 3, 1159–1166 the question whether MSUC arose in the common mammalian ances- (2005). 27. Immler, S. Haploid selection in “diploid” organisms. Annu. Rev. Ecol. Evol. Syst. 50, tor after the split from the reptile lineage roughly 310 million years 219–236 (2019). ago or whether this mechanism was lost in the avian lineage and arose 28. Cheng, X. & Sun, Q. RUBCNL/Pacer and RUBCN/Rubicon in regulation of autolysosome earlier in evolution. The latter scenario would be consistent with the formation and lipid metabolism. Autophagy 15, 1120–1121 (2019). 29. Zhu, Y. et al. Autophagy in male reproduction. Syst. Biol. Reprod. Med. 65, 265–272 observation of MSCI in invertebrate species13. (2019). Our data and results, together with the accompanying online resource 30. Russell, D. L., Brown, H. M. & Dunning, K. R. ADAMTS proteases in fertility. Matrix Biol. we developed (https://apps.kaessmannlab.org/SpermEvol/) to facili- 44–46, 54–63 (2015). 31. Zhang, T., Oatley, J., Bardwell, V. J. & Zarkower, D. DMRT1 is required for mouse tate the exploration of our data, provide an extensive resource for inves- spermatogonial stem cell maintenance and replenishment. PLoS Genet. 12, e1006293 tigating the biology of the testis and associated fertility disorders across (2016). mammals. Future studies should seek to complement our snRNA-seq 32. Ashburner, M. et al. Gene Ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat. Genet. 25, 25–29 (2000). data to overcome its limitations. scRNA-seq data will be valuable for 33. Fu, C. et al. Multiple aspects of male germ cell development and interactions with Sertoli inferring transcriptome patterns unique to the cytoplasm, single-cell cells require inositol hexakisphosphate kinase-1. Sci. Rep. 8, 7039 (2018). full-length transcript data are needed to assess the pronounced isoform 34. Inoue, H. et al. Mouse D1Pas1, a DEAD-box RNA helicase, is required for the completion of first meiotic prophase in male germ cells. Biochem. Biophys. Res. Commun. 478, diversity of the testis11, and single-cell translatome data50 are required 592–598 (2016). to understand the contribution of post-transcriptional changes6 to the 35. Barral, S. et al. Histone variant H2A.L.2 guides transition protein-dependent protamine evolution of spermatogenesis. assembly in male germ cells. Mol. Cell 66, 89–101 e8 (2017). 8 | Nature | www.nature.com 36. Yamada, D. et al. Disruption of spermatogenic cell adhesion and male infertility in mice lacking TSLC1/IGSF4, an immunoglobulin superfamily cell adhesion molecule. Mol. Cell Open Access This article is licensed under a Creative Commons Attribution Biol. 26, 3610–3624 (2006). 4.0 International License, which permits use, sharing, adaptation, distribution 37. Wakayama, T. & Iseki, S. Role of the spermatogenic-Sertoli cell interaction through and reproduction in any medium or format, as long as you give appropriate cell adhesion molecule-1 (CADM1) in spermatogenesis. Anat. Sci. Int. 84, 112–121 credit to the original author(s) and the source, provide a link to the Creative Commons licence, (2009). and indicate if changes were made. The images or other third party material in this article are 38. Ozaki-Kuroda, K. et al. Nectin couples cell-cell adhesion and the actin scaffold at included in the article’s Creative Commons licence, unless indicated otherwise in a credit line heterotypic testicular junctions. Curr. Biol. 12, 1145–1150 (2002). to the material. If material is not included in the article’s Creative Commons licence and your 39. Bronson, R., Mikhailik, A., Schwedes, J., Gnatenko, D. & Hatchwell, E. Detection of intended use is not permitted by statutory regulation or exceeds the permitted use, you will candidate nectin gene mutations in infertile men with severe teratospermia. J. Assist. need to obtain permission directly from the copyright holder. To view a copy of this licence, Reprod. Genet. 34, 1295–1302 (2017). visit http://creativecommons.org/licenses/by/4.0/. 40. Cortez, D. et al. Origins and functional evolution of Y chromosomes across mammals. Nature 508, 488–493 (2014). © The Author(s) 2022 41. Chen, Y. et al. Single-cell RNA-seq uncovers dynamic processes and critical regulators in mouse spermatogenesis. Cell Res. 28, 879–896 (2018). 1 Center for Molecular Biology (ZMBH), DKFZ-ZMBH Alliance, Heidelberg University, 42. Braun, R. E., Behringer, R. R., Peschon, J. J., Brinster, R. L. & Palmiter, R. D. Genetically Heidelberg, Germany. 2INRAE, LPGP, Rennes, France. 3Department of Growth and haploid spermatids are phenotypically diploid. Nature 337, 373–376 (1989). Reproduction, Rigshospitalet, Copenhagen University Hospital, Copenhagen, Denmark. 43. Ernst, C., Eling, N., Martinez-Jimenez, C. P., Marioni, J. C. & Odom, D. T. Staged 4 International Center for Research and Research Training in Endocrine Disruption of Male developmental mapping and X chromosome transcriptional dynamics during mouse spermatogenesis. Nat. Commun. 10, 1251 (2019). Reproduction and Child Health, Rigshospitalet, Copenhagen University Hospital, 44. Bhutani, K. et al. Widespread haploid-biased gene expression enables sperm-level Copenhagen, Denmark. 5Bioinformatics Research Centre, Aarhus University, Aarhus, natural selection. Science 371, eabb1723 (2021). Denmark. 6Berlin Institute of Health at Charité, University of Medicine Berlin, Corporate 45. Daish, T. J., Casey, A. E. & Grutzner, F. Lack of sex chromosome specific meiotic Member of the Free University of Berlin, Humboldt-University of Berlin, Berlin, Germany. silencing in platypus reveals origin of MSCI in therian mammals. BMC Biol. 13, 106 7 Evolutionary Developmental Biology Laboratory, Francis Crick Institute, London, UK. (2015). 8 Biomedical Primate Research Center (BPRC), Rijswijk, the Netherlands. 9German Primate 46. Zhou, Y. et al. Platypus and echidna genomes reveal mammalian biology and evolution. Center (DPZ), Platform Degenerative Diseases, Göttingen, Germany. 10German Center for Nature 592, 756–762 (2021). Cardiovascular Research (DZHK), Partner Site Göttingen, Göttingen, Germany. 11Center for 47. Mahadevaiah, S. K. et al. Key features of the X inactivation process are conserved between marsupials and eutherians. Curr. Biol. 19, 1478–1484 (2009). Neurobiology and Brain Restoration, Skolkovo Institute of Science and Technology, Moscow, 48. Sin, H. S. & Namekawa, S. H. The great escape: active genes on inactive sex Russia. 12Max Planck Institute for Evolutionary Anthropology, Leipzig, Germany. 13Institute of chromosomes and their evolutionary implications. Epigenetics 8, 887–892 (2013). Evolutionary Biology (UPF-CSIC), Barcelona, Spain. 14Catalan Institution of Research and 49. Guioli, S., Lovell-Badge, R. & Turner, J. M. Error-prone ZW pairing and no evidence for Advanced Studies (ICREA), Barcelona, Spain. 15CNAG-CRG, Centre for Genomic Regulation meiotic sex chromosome inactivation in the chicken germ line. PLoS Genet. 8, e1002560 (CRG), Barcelona Institute of Science and Technology (BIST), Barcelona, Spain. 16Miquel (2012). Crusafont Catalan Institute of Paleontology, Autonomous University of Barcelona, Barcelona, 50. VanInsberghe, M., van den Berg, J., Andersson-Rolf, A., Clevers, H. & van Oudenaarden, Spain. 17The Robinson Research Institute, School of Biological Science, University of Adelaide, A. Single-cell Ribo-seq reveals cell cycle-dependent translational pausing. Nature 597, Adelaide, South Australia, Australia. 18Department of Cellular and Molecular Medicine, Faculty 561–565 (2021). of Health and Medical Sciences, University of Copenhagen, Copenhagen, Denmark. 19These Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in authors contributed equally: Florent Murat, Noe Mbengue. ✉e-mail: florent.murat@inrae.fr; published maps and institutional affiliations. n.mbengue@zmbh.uni-heidelberg.de; h.kaessmann@zmbh.uni-heidelberg.de Nature | www.nature.com | 9 Article Methods Tris-HCl to a final concentration of 20 mM. The fixed nuclei were pel- leted at 500g for 5 min at 4 °C. The supernatant was removed and the Data reporting pellet resuspended in 0.67 vol. of Wash Buffer (250 mM sucrose, 25 mM No statistical methods were used to predetermine sample size. The KCl, 5 mM MgCl2, 10 mM Tris-HCl pH 8, 1% BSA and freshly added 1 µM experiments were not randomized and investigators were not blinded DTT, 0.4 U µl−1 RNase Inhibitor, 0.2 U µl−1 SUPERasIn). This was centri- to allocation during experiments and outcome assessment. fuged at 500g for 5 min at 4 °C. The supernatant was removed and the pellet resuspended in 0.5 vol. of PBS. Then nuclei were strained using Biological samples and ethics statement 40 μm Flowmi strainers (Sigma). For estimation of nuclei concentra- We generated snRNA-seq data for adult testis samples from human tion, Hoechst DNA dye was added and the nuclei were counted using (Homo sapiens), chimpanzee (Pan troglodytes), bonobo (Pan paniscus), Countess II FL Automated Cell Counter (ThermoFischer Scientific). gorilla (Gorilla gorilla), lar gibbon (Hylobates lar), rhesus macaque For platypus and chicken, a similar preparation method was used, but (Macaca mulatta), common marmoset (Callithrix jacchus), mouse the nuclei were not fixed, given that this protocol gave optimal results (Mus musculus, strain: RjOrl:SWISS; Janvier Laboratories), grey short- for these species (the fixation protocol failed to yield data of adequate tailed opossum (Monodelphis domestica), platypus (Ornithorhynchus quality). In brief, tissue pieces weighing roughly 5 mg were homog- anatinus) and chicken (red jungle fowl, Gallus gallus) (Supplemen- enized in 100–150 µl mg−1 of prechilled lysis buffer (250 mM sucrose, tary Table 2). In addition, we produced bulk RNA-seq data for chim­ 25 mM KCl, 5 mM MgCl2, 10 mM Tris-HCl pH 8, 1% BSA, 0.1% IGEPAL and panzee, gorilla, gibbon and marmoset from the same individuals. freshly added 1 µM DTT, 0.4 U µl−1 RNase Inhibitor, 0.2 U µl−1 SUPERasin) Adult human testis samples used for in situ hybridization experiments and lysed for 5 min on ice. The lysate was centrifuged at 100g for 1 min were obtained from orchiectomy specimens from three individuals at 4 °C. The supernatant was transferred to a new reaction tube and with testicular cancer. Tissue adjacent to the tumour that was devoid centrifuged at 500g for 5 min at 4 °C. The supernatant was removed of cancer cells and germ cell neoplasia in situ, and tubules with nor- and the pellet resuspended in 0.67 vol. of Wash Buffer. This was centri- mal spermatogenesis were used. Other adult primate testis tissue fuged at 500g for 5 min at 4 °C. The supernatant was removed and the was obtained from a western chimpanzee and a Bornean orangutan pellet resuspended in 0.5 vol. of PBS. Then nuclei were strained using (Pongo pygmaeus). 40-μm Flowmi strainers (Sigma). For estimation of nuclei concentra- Our study complies with all relevant ethical regulations with respect tion, Hoechst DNA dye was added and the nuclei were counted using to both human and other species’ samples. Human samples underlying Countess II FL Automated Cell Counter (ThermoFischer Scientific). the snRNA-seq data were obtained from scientific tissue banks (http:// medschool.umaryland.edu/btbank/) or dedicated companies (https:// Preparation of sequencing libraries www.tissue-solutions.com/); informed consent was obtained by these For the construction of snRNA-seq libraries, Chromium Single Cell 3′ sources from donors before death or from next of kin. The samples Reagent Kits (10X Genomics; v.2 chemistry for human, chimpanzee, used for the RNA in situ hybridization experiments were obtained from bonobo, gorilla, gibbon, macaque, marmoset, mouse and opossum; v.3 the tissue biobank at the Department of Growth and Reproduction chemistry for platypus and chicken) were used according to the manu- (Rigshospitalet, Copenhagen, Denmark) containing orchiectomy speci- facturer’s instructions. Then 15,000 to 20,000 nuclei were loaded per mens from individuals with testicular cancer (The Danish Data Protec- lane in the Chromium microfluidic chips and complementary DNA was tion Agency, permit number J. no. 2001-54-0906). All patients have amplified in 12 PCR cycles. Sequencing was performed with NextSeq550 given informed consent for donating the residual tissues for research. (Illumina) according to the manufacturer’s instructions using the The use of all human samples for the type of work described in this study NextSeq 500/550 High Output Kit v.2.5 (75 cycles) with paired-end was approved by an ethics screening panel from the European Research sequencing (read lengths of read1 26 bp, read2 57 bp; index1 8 bp, Council (ERC) and local ethics committees: from the Cantonal Ethics roughly 170 to 380 million reads per library for v.2 chemistry; read Commission in Lausanne (authorization 504/12); from the Ethics Com- lengths of read1 28 bp, read2 56 bp Index1 8 bp, roughly 247 to 306 million mission of the Medical Faculty of Heidelberg University (authorization reads per library for v.3 chemistry) (Supplementary Table 2). S-220/2017) and from the regional medical research ethics committee For bulk RNA-seq data generation, RNA was extracted using the of the capital region of Copenhagen (H-16019637). All primates used in RNeasy Micro kit (QIAGEN). The tissues were homogenized in RLT this study suffered sudden deaths for reasons other than their partici- buffer supplemented with 40 mM DTT. The RNA-seq libraries were pation in this study and without any relation to the organ sampled. The constructed using the TruSeq Stranded messenger RNA LT Sample use of all other mammalian samples for the type of work described in Prep Kit (Illumina) as described in ref. 2. Libraries were sequenced on this study was approved by ERC ethics screening panels. Illumina NextSeq550 using a single-end run (read1 159 bp; index1 7 bp) with roughly 24–60 million reads per library (Supplementary Table 2). Nuclei isolation For the samples from therian species we developed a nuclei preparation Genome and transcript isoform annotation method that includes fixation with dithio-bis(succinimidyl propionate) Given that the quality of genome annotation differs substantially (DSP; or Lomant’s Reagent), a reversible cross-linker that stabilizes the between the studied species and given the specific and widespread isolated nuclei. The method was adapted from protocols used for the transcription of the genome in the testis11, we refined and extended fixation of single-cell suspensions51 and for the isolation of single nuclei previous annotations from Ensembl53 on the basis of testis RNA-seq from archived frozen brain samples52. Tissue pieces weighing roughly data. Specifically, akin to the procedure described in refs. 2,6, we used 5 mg were homogenized in 100–150 µl mg−1 of prechilled lysis buffer previous extensive stranded poly(A)-selected RNA-seq datasets2,54 (250 mM sucrose, 25 mM KCl, 5 mM MgCl2, 10 mM HEPES pH 8, 1% BSA, (100 nt, single-end) for human, macaque, mouse, opossum, platy- 0.1% IGEPAL and freshly added 1 µM DTT, 0.4 U µl−1 RNase Inhibitor (New pus and chicken, and generated and used stranded poly(A)-selected England BioLabs), 0.2 U µl−1 SUPERasIn (ThermoFischer Scientific)) RNA-seq datasets (159 nt, single-end) for chimpanzee, gorilla, gibbon and lysed for 5 min on ice. The lysate was centrifuged at 100g for 1 min and marmoset. For each species, we downloaded the reference genome at 4 °C. The supernatant was transferred to a new reaction tube and from Ensembl release 87 (ref. 53): hg38 (human), CHIMP2.1.4 (chim- centrifuged at 500g for 5 min at 4 °C. The supernatant was removed panzee), rheMac8 (rhesus macaque), C_jacchus3.2.1 (marmoset), and the pellet resuspended in 0.67 vol. (of volume lysis buffer used) of mm10 (mouse), monDom5 (opossum), and galGal5 (chicken); from freshly made fixation solution (1 mg ml−1 DSP in PBS) and incubated for Ensembl release 96 (ref. 55): gorGor4 (gorilla) and Nleu_3.0 (gibbon); 30 min at room temperature. The fixation was quenched by addition of and from Ensembl release 100 (ref. 56): mOrnAna1.p.v1 (platypus). Raw reads were first trimmed with cutadapt v.1.8.3 (ref. 57) to remove we removed potential doublets using doubletFinder_v3 function of adapter sequences and low-quality (Phred score <20) nucleotides, DoubletFinder63 v.2.0.1 (parameters PCs = 1:20, pN = 0.25, nExp = 5% of then reads shorter than 50 nt were filtered out (parameters: --adapter= the total number of cells, identifying pk using paramSweep_v3, sum- AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC --match-read- marizeSweep and find.pK functions). wildcards --minimum-length=50 -q 20). Processed reads were then mapped to the reference transcriptome and genome using Tophat2 Integration of datasets v.2.1.1 (ref. 58) (parameters: --bowtie1 --read-mismatches 6 --read- From the previously filtered UMI count tables, we created Seurat gap-length 6 --read-edit-dist 6 --read-realign-edit-dist 0 --segment- objects for every sample independently, normalized the data and length 50 --min-intron-length 50 --library-type fr-firststrand --max- identified the top 10,000 most variable genes. Next, for each spe- insertion-length 6 --max-deletion-length 6). Next, we assembled models cies independently, we applied the Seurat62 anchoring approach of transcripts expressed using StringTie v.1.3.3 (ref. 59) (parameters: using FindIntegrationAnchors and IntegrateData functions with -f 0.1 -m 200 -a 10 -j 3 -c 0.1 -v -g 10 -M 0.5). Stringent requirements on 20 principal components to integrate all datasets together into a the number of reads supporting a junction (-j 3), minimum gap between single Seurat object correcting for the batch effect. For each inte- alignments to be considered as a new transcript (-g 10) and fraction grated species-specific Seurat object, we normalized (NormalizeData covered by multi-hit reads (-M 0.5) were used to avoid merging inde- function) and scaled the data (ScaleData function) and performed a pendent transcripts and to reduce the noise caused by unspliced or PCA (RunPCA function). Louvain clusters were calculated using Find- incompletely spliced transcripts. We compared the assembled tran- Neighbors and FindClusters functions (parameters dims = 1:20, 1:20, script models to the corresponding reference Ensembl annotations 1:20, 1:20, 1:20, 1:20, 1:20, 1:17, 1:8, 1:10 and 1:10, and resolution = 0.5, using the cuffcompare program v.2.2.1 from the cufflinks package60. 0.5, 0.2, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.3 and 0.5, respectively, for human, Finally, we combined the newly identified transcripts with the respec- chimpanzee, bonobo, gorilla, gibbon, macaque, marmoset, mouse, tive Ensembl gene annotation into a single gtf file. We extended the opossum, platypus and chicken). The uniform manifold approxima- original Ensembl annotations by 2–61 Mbp with new transcripts and by tion and projection (UMAP) embedding coordinates were calculated 23–49 Mbp with new splice isoforms (Supplementary Table 1). using the RunUMAP function (parameters dims = 1:20, 1:20, 1:20, 1:20, 1:20, 1:20, 1:20, 1:17, 1:10, 1:10 and 1:10, and min_dist = 0.3, 0.3, Raw reads processing 0.1, 0.1, 0.3, 0.3, 0.3, 0.1, 0.2, 0.3 and 0.6, respectively, for human, CellRanger v.3.0.2 was used for platypus and chicken, and CellRanger chimpanzee, bonobo, gorilla, gibbon, macaque, marmoset, mouse, v.2.1.1 for the other species in line with the used Chromium chemis- opossum, platypus and chicken). We note that—consistent with the try. The CellRanger mkref function was used with default settings to high correlation between biological replicates (Supplementary build each species reference from genomic sequences and custom- Fig. 5a,b)—the data already integrate well before the batch correc- ized extended annotation files (Supplementary Table 1). Given that tion (Supplementary Fig. 18c,d). We also note that key marker genes pre-mRNA transcripts are abundant in nuclei61, exons and introns are expressed in the same integrated areas across replicates when features were concatenated as described in the CellRanger v.2.1.1 assessing their expression in the different replicates using the inte- documentation for considering intronic and exonic reads for gene grated object coordinates (Supplementary Fig. 18e), which supports expression quantification. The CellRanger count function was used that the integration is correct. with default settings to correct droplet barcodes for sequencing errors, All primate datasets were merged using the LIGER64 (v.0.5.0) integra- align reads to the genome and count the number of UMIs for every gene tion tool. A LIGER object was created using the createLiger function and barcode combination. based on primate 1:1 orthologues from Ensembl release 87, and nor- malized with normalize, selectGenes and scaleNotCenter functions Identification of usable nuclei with default settings. Then, the joint matrix was factorized using the We used a combined approach for detection of usable nuclei. This optimizeALS function (k = 20) and the quantile normalization was was done to optimally account for the lower RNA content of nuclei performed with the quantile_norm and default settings. The Louvain compared to the whole cells. Specifically, to identify usable nuclei, clusters were calculated with the louvainCluster function and default we used a knee point-based approach combined with the fraction of settings as well as UMAP coordinates with the runUMAP function intronic reads as a marker of pre-mRNA transcripts (abundant in the (n_neighbors = 100, min_dist = 0.2). nucleus) (Supplementary Fig. 1f) and MALAT1 (nuclear-enriched long non-coding (lnc)RNA) expression as a marker of nuclei (when present Estimation of expression levels and normalization in the genome of a given species). The gene UMI counts per cell were normalized using the Seurat R pack- age and its NormalizeData function. Therefore, the UMI counts of each Quality control of filtered cells gene in each cell are divided by the total UMI counts of each cell, mul- For each sample independently, high-quality nuclei were selected tiplied by 10,000 and log transformed. removing outliers on the basis of the number of UMIs and the percent- age of mitochondrial RNA (Supplementary Fig. 1b,c,e). We created a Cell-type assignment Seurat62 object using the Seurat R package v.3.1.4 from the subset raw We identified the main cell-type populations from the primate inte- UMI count table generated by CellRanger corresponding to the usable grated, mouse, opossum, platypus and chicken objects indepen- droplets identified upstream, normalized the data using the Normalize­ dently using known marker genes65,66 mostly from human and mouse Data function, identified the top 10,000 most variables genes using and their respective 1:1 orthologues in the other species. CLU marks the FindVariableFeatures function, scaled the data using the ScaleData Sertoli cells; TAGLN and ACTA2 peritubular and smooth muscle cells; function, performed the PCA using the RunPCA function and calculated CD34 and TM4SF1 endothelial cells; APOE and CD74 macrophages; the Louvain clusters using the FindNeighbors (parameters dims = 1:20) STAR and CYP11A1 Leydig cells; GFRA1, PIWIL4 (undifferentiated), and FindClusters (dims = 1:20, resolution = 0.5) functions. To opti- DMRT1(differentiated) and STRA8 SG; SYCE1 (leptotene), SYCP1 (zygo- mally account for the fact that testis cell types have diverse transcrip- tene), PIWIL1 (pachytene), SYCP2, TANK and AURKA SC; LRRIQ1 (early), tome characteristics11, we filtered out outlier droplets for each cluster ACRV1 and SPACA1 (late) rSD; and SPATA3, NRBP1, PRM1 and GABBR2 independently with values lower than the first quartile (Q1) − 1.5 × IQR eSD. Cell-type assignment was robustly reinforced by complementary (interquartile range) and higher than the third quartile (Q3) + 1.5 × IQR analyses and metrics such as UMAP coordinates, pseudotime trajecto- for both the UMI content and the fraction of mitochondrial RNA. Then, ries, transcriptional activities (UMI counts) and previous knowledge. Article the branching pattern of the majority-rule consensus tree shown in Pseudotime the figures (Fig. 1c and Extended Data Fig. 2a). The total tree length Pseudotime trajectories were calculated using the slingshot v.1.2.0 was calculated by removing the intra-species variability between indi- R package67. We applied the getLineages function with the upstream viduals (Fig. 2a). calculated clusters and UMAP embedding coordinates of the germ cells to obtain connections between adjacent clusters using a minimum Evolutionary forces spanning tree. We provided the starting and ending clusters on the In Fig. 2c, we plotted the median pLI score69 across expressed genes (≥ basis of the previous cell-type assignment with known marker gene 1 UMI) in each nucleus. We obtained the pLI scores from ref. 70. For expression. Then we applied the getCurves function to the obtained Fig. 2d, we used a set of neutrally ascertained knockouts consisting in lineages to construct smooth curves and order the cells along a pseu- 4,742 protein-coding genes, 1,139 of which are classified as lethal. For dotime trajectory. Pseudotime values are highly variable depending each cell, the denominator is the number of genes expressed that were on used tools, thus we ordered the cells one by one on the basis of their tested for lethality and the numerator the genes among those that pseudotime values and divided their rank by the total number of cells, resulted in a lethal phenotype. Tested genes for viability and associated to obtain evenly distributed values between 0 and 1. Finally, we vali- phenotype information were downloaded from the International dated the obtained pseudotime trajectories on the basis of previous Mouse Phenotyping Consortium25. For Fig. 2e (and Extended Data cell-type assignments, expression patterns of marker genes and UMAP Fig. 3a), we used the average dN/dS values across 1:1 orthologues in embedding coordinates. primates. For each cell, the mean dN/dS value is plotted. Conserved 1:1 orthologues across six primates (human, chimpanzee, gorilla, gibbon, Marker gene identification macaque and marmoset) as well as their coding and protein sequences To precisely identify marker genes along spermatogenesis, we grouped were extracted from Ensembl53, providing a set of 11,791 protein-coding the germ cells into 20 evenly distributed bins along the pseudotime tra- genes. For each species and orthologue the longest transcript was jectory for each species. Then, we applied the FindAllMarkers function extracted. Orthologous protein sequences were aligned using clustalo from the Seurat62 R package (parameter only.pos = T, min.pct = 0.25, v.1.2.4; then pal2nal v14 was used (with protein sequences alignments logfc.threshold = 0.25, return.thresh = 0.05) of the Seurat v.3.1.4 R and coding sequences as input) to produce codon-based alignments. package to the 22 groups (20 germline groups, the Sertoli and other The codeml software from the PAML package71 v.4.9 was used to esti- somatic cell groups) (Extended Data Fig. 1 and Supplementary Table 4). mate dN/dS ratios. The M0 site model was applied to the orthologue alignments to estimate one average dN/dS ratio per orthologous gene Phylogenetic trees set across species (parameter NSites = 0, model = 0). In Fig. 2f (and Phylogenetic trees and indicated divergence times (Figs. 1a and 4a) are Extended Data Fig. 3b), we plotted the percentage of positively selected based on TimeTree68 (v.5) (http://www.timetree.org/). genes expressed across nuclei. For each nucleus, the denominator is the number of expressed genes that were tested for signatures of Orthologous gene sets positive selection, and the numerator is the number of genes among We used four different sets of orthologous genes in our study: (1) com- those with evidence for positive selection. We used sets of genes previ- parative analyses involving all 11 amniote species were performed using ously identified as carrying evidence for coding-sequence adaptation 4,498 1:1 orthologue genes that are expressed (that is, one UMI in at least in primates72 (331 positively selected genes out of 11,170 genes tested) three cells of any cell type) across all species (among a total of 8,045 1:1 and mammals73 (544 positively selected genes out of 16,419 genes orthologues). (2) Comparative analyses involving the seven primate tested). species were performed using 8,451 1:1 orthologue genes expressed In Fig. 2g (and Extended Data Fig. 3c), we plotted the average phy- across all primate species (among 11,948 1:1 orthologues). (3) The logenetic age of expressed genes across somatic and germ cells. The comparative Sertoli-germ cell communication analysis was based on phylogenetic age of genes is an index that gives greater weight to mapping 35,186 human testis-expressed genes to 1:1 orthologous genes young new genes (as described in ref. 2,74). The range of the score dif- in the other species (macaque 13,090; mouse 14,302; opossum 10,865 fers between species depending on the number of outgroup lineages and chicken 10,515). (4) Species-specific analyses were performed using available (more lineages allowed for more details in the phylogeny) all genes expressed in a given species (roughly 15,000 genes per cell and therefore this index cannot be compared across species, only type; Supplementary Fig. 1d). Orthologous gene sets were extracted within species (that is, across cells and cell types). The phylogenetic from Ensembl53 using the biomaRt R package v.2.40.5. age of genes was obtained from GenTree (http://gentree.ioz.ac.cn/) with Ensembl release 69 (ref. 74). In Fig. 2h (and Extended Data Fig. 3d), Global patterns of gene expression differences across mammals we plotted the percentage of the cell transcripts originating from Pseudo-bulk samples were generated using the AverageExpression protein-coding genes and intergenic elements. Gene biotypes were function of the Seurat R package with various groups of cells depending obtained from Ensembl. Intergenic elements are all elements that are on the pseudo-bulk samples produced in the study. For the analyses not protein-coding genes (lncRNAs, pseudogenes, pseudogenes and presented in Fig. 1b, we performed the PCA of normalized expression in other sequences). In Fig. 2i (and Extended Data Fig. 3e), normalized amniote testicular cell types (pseudo-bulks) for each individual based log2-transformed median expression values across replicates at the on 4,498 1:1 amniote orthologues. PCA was performed using the prcomp transcriptome (etr) and translatome (etl) layers were used to calculate function of the stats R package. For Fig. 1c, we constructed gene expres- translation efficiency (TE = log2(etr) − log2(etl)) in testis (as described sion trees (as described in ref. 1) using the neighbour-joining approach, in ref. 6) from RNA-seq and Ribo-seq data6, respectively. Translation on the basis of pairwise expression distance matrices between cor- efficiency values were calculated at the whole testis level, thus only responding pseudo-bulk samples for the different cell types across cell-type-specific genes (for which 60% of their transcripts at the species. The distance between samples was computed as 1 − ρ, where whole testis level are concentrated in a single cell type) were used. ρ is Spearman’s correlation coefficient and was computed using the Higher values show a more efficient translation of transcripts, whereas cor function of the stats R package. The neighbour-joining trees were lower values indicate relative translational repression. For Fig. 2j (and constructed using the ape R package v.5.3. The reliability of branching Extended Data Fig. 3f,g), we used time- and tissue-specificity indexes of patterns was assessed with bootstrap analyses (the 4,498 1:1 amniote expressed genes across somatic and germ cells in testis. As described in orthologues were randomly sampled with replacement 1,000 times). ref. 2, tissue and time specificity indexes are calculated from RNA-seq The bootstrap values are the proportions of replicate trees that share data across organs and developmental stages. Both indexes range from 0 (broad expression) to 1 (restricted expression). The indexes The sections were counterstained with Mayer’s haematoxylin and were obtained from ref. 2. For each nucleus, we plotted the median mounted with Vectamount Permanent Mounting Medium (Vector Labora- index across expressed genes. In Fig. 2k, we plotted the percentage of tories). The negative control probe DapB (a bacterial RNA) was run in para­ genes causing infertility when knocked out (out of 3,252 knockouts, llel with the target probes and showed ≤5% positive cells in each section. 173 of which caused infertility). Tested genes for infertility and associ- For ADAMTS17 and MYO3B, positive (that is, with red dots) rSD and ated phenotype information were downloaded from the International SG were counted. For each section (human n = 3, chimpanzee n = 1, Mouse Phenotyping Consortium database25. For each nucleus, the orangutan n = 1), ten tubules were counted using the NDP.viewPlus denominator corresponds to the number of genes expressed that were software (Hamamatsu Photonics). Two independent observers (S.B.W. tested for infertility and the numerator to the genes among those that and K.A.) counted positive and negative rSD and SG. No discrimination resulted in an infertility phenotype. in intensity of the dots or the number of dots per cell was performed. Cell-type identification was performed on the basis of nucleus mor- Gene expression trajectories along spermatogenesis phology and localization in the tubule. Only SG lining the edge of the We compared gene expression trajectories along spermatogenesis tubules were counted (Supplementary Fig. 3a). The inter-observer across primates using human, chimpanzee, bonobo, gorilla, gibbon, variance was found to be 7 and 9% for rSD and SG, respectively. For macaque and marmoset (based on Ensembl 87 orthologues), and across RUBCNL, quantification of staining intensity was performed with R amniotes using human (as a representative of primates), mouse, opos- v.3.6.1 using countcolors76 (v.0.9.1) and colordistance77 (v.1.1.1) pack- sum, platypus and chicken (based on Ensembl 100 orthologues). To ages. For each section (human n = 3, chimpanzee n = 1, orangutan n = 1), compare robustly expressed genes, we used genes that are expressed ten tubules were divided into three parts: area dominated by SG, area in at least 5% of the cells in at least one cluster in all considered species. dominated by SCs (also containing Sertoli cells) and area dominated by We used the mfuzz package75 (v.2.44.0), an unsupervised soft clustering spermatids (no distinction between the different types of spermatid) method, to cluster gene expression trajectories along spermatogenesis (Supplementary Fig. 3b). In each tubular area, the number of cells was (eight cell types in primates; four cell types in amniotes) across species counted manually using the NDP.view2Plus software (v.2.8.24). Then, using Dmin and mestimate functions to estimate the number of clusters the pixels occupied by red staining were quantified and the expression and the fuzzification parameter (Supplementary Fig. 2). As described level for each cell type was calculated by dividing the stained pixels by in ref. 2, we inferred within a phylogenetic framework the probability the number of cells. For each picture, the stained pixels for each cell that there were changes in trajectories along spermatogenesis, that type were normalized by the total amount of stained pixels. is, that genes changed their cluster assignment in specific branches, using a 5% probability cut-off. Gene Ontology analysis Enriched terms in the Gene Ontology32 analyses of genes with conserved Trajectory conservation score and diverged expression trajectories were identified using the goana We calculated a global trajectory conservation score across species for function of the limma R package, v.3.40.6 (default parameters). each 1:1 orthologous gene set. For a given orthologous gene set, this score corresponds to the log-transformed probability that all members Sertoli-germ cell communication analysis fall into the same mfuzz trajectory cluster as: We identified ligand–receptor interactions underlying Sertoli-germ cell communications for human, macaque, mouse, opossum and  Conservation_scoreg = log 2  ∑i∈c ∏j∈s Pg , i, j   chicken, respectively, using the CellPhoneDB78 (v.2) approach and recommended parameter settings (parameters method statistical_anal- ysis). To apply CellPhoneDB, which uses a database of human ligand– where g corresponds to a given orthologous gene set, c to all mfuzz receptor interactions, to the data for the other species, we mapped trajectory clusters (1–9 for primates, 1–12 for amniotes), s to all species human testis-expressed genes to their corresponding 1:1 orthologues in (human, chimpanzee, bonobo, gorilla, gibbon, macaque and marmo- each of these species. Enriched receptor–ligand interactions between set for primates; human, mouse, opossum, platypus and chicken for two cell types are predicted on the basis of expression of a receptor amniotes) and Pg,i,j to the probability that the gene g of the species j by one cell type and a ligand by another cell type78. Only receptors falls into the cluster i. A higher conservation score indicates a greater and ligands expressed in more than 10% of the cells in each cell type global trajectory conservation. As a proof of concept, we plotted the were considered for pairwise comparisons between all cell types in trajectory conservation score for conserved and changed trajectories the dataset. CellPhoneDB uses empirical shuffling to calculate which that revealed a significant higher conservation score for conserved ligand–receptor pairs show significant (P < 0.05) cell-type specificity78. trajectories (Extended Data Fig. 5a). Significant interactions across species were illustrated using the R package UpSetR (v.1.4.0). Finally, given potential false positive (and RNA in situ hybridization and expression quantification negative) predictions of CellPhoneDB and similar approaches79, we Fresh testicular tissue was fixed in GR-fixative (7.4% formaldehyde, 4% consider interactions that are predicted for several species and prob- acetic acid, 2% methanol, 0.57% sodium phosphate, dibasic and 0.11% ably reflect evolutionary conservation (Extended Data Fig. 8 and potassium phosphate, monobasic) overnight (for at least 16 h) at 4 °C, Supplementary Table 9), such as those reported in the main text, to be dehydrated and embedded in paraffin. The in situ hybridization experi- more reliable than species-specific predictions. ments were carried out on 4 μm sections mounted on SuperFrost Plus Slides (ThermoFisher Scientific) using the RNAScope 2.5 HD Detection Cell-type and testis-specific genes per chromosome Reagent RED kit according to the manufacturer’s recommendations Testis-specific genes were obtained from previously generated RNA-seq (Advanced Cell Diagnostics). Briefly, testicular tissue sections were data2 of adult organs (RPKM ≥ 1 in testis and RPKM < 1 in brain, cerebel- dewaxed in xylene and washed in 100% ethanol followed by treatment lum, heart, kidney and liver). Among these, cell-type-specific genes with hydrogen peroxide for 10 min. Target retrieval was performed for were studied for each chromosome. Genes with predominant expres- 15 or 30 min (see Supplementary Table 6 for specifications for each sion in specific somatic cell types were identified using the FindAllMark- probe and species) using a steamer, followed by treatment with pro- ers function (parameter only.pos = TRUE, min.pct = 0.05, logfc. tease plus for 30 min at 40 °C. The slides were hybridized with the target threshold = 0.25, return.thresh = 0.05). Predominant expression of probe (Supplementary Table 6) for 2 h at 40 °C followed by a series of genes in specific germ cell types was assigned on the basis of the trajec- signal amplifications (with amplification around 5 for 30 or 60 min). tory analyses (above); that is, predominant expression was assigned Article on the basis of the cell type in which the expression level of the gene indeed actively transcribed in SC (that is, they do not merely represent peaks in the trajectory analysis. We then first calculated the percentage genes expressed in SG with stable transcripts still detectable in SC), akin of genes located on a given chromosome among all genes in the genome to previous work23. Specifically, we identified differentially expressed (x axis of plots in Extended Data Fig. 9a,c; red horizontal line for X-linked genes between SC and SG using the FindMarkers function from Seurat62 genes in Fig. 4b). We then contrasted this with the percentage of R package (parameters, default). A Wilcoxon rank-sum test was used testis-specific genes with predominant expression in a given cell type to calculate P values, which were then adjusted using a Bonferroni cor- (y axis of plots in Extended Data Fig. 9a,c; y axis of plots in Fig. 4b for rection for several tests. Only genes that were detected as expressed in X-linked genes). Finally, the percentage of testis-specific genes per cell at least 10% of cells from either of the two cell-type populations were type and chromosome was statistically compared to the percentage tested. Genes showing, on average, at least a 0.25-fold expression dif- of genes per chromosome in the genome using exact binomial tests. ference (loge-scale) between the two groups, and an adjusted P value below 0.05 were considered to be differentially expressed. X-linked Classification of X- and Y-bearing spermatids genes in SDRs with significantly higher expression in SC than SG were The Y chromosome carries a low number of genes and is missing in considered to be potential escapees (Extended Data Fig. 12c). some genome assemblies. Thus, we focused on the fraction of X tran- scripts in spermatids to classify them as X- or Y-bearing cells. For this, General statistics and plots we fitted a Gaussian Mixture Model to the data with two components Unless otherwise stated, all statistical analyses and plots were done (bimodal distribution) independently for each replicate, using the in R v.3.6.2 (ref. 83). Plots were created using ggplot2 v.3.2.1, tidyverse function normalmixEM of the mixtools (v.1.2.0) R package. The two v.1.3.0, dplyr v.0.8.5, cowplot v.1.0.0 and pheatmap v.1.0.12. obtained normal distributions were used to classify X- (higher levels of X transcripts) and Y-bearing (lower levels of X transcripts) spermatids Reporting summary using 95% confidence intervals. Outlier and overlapping cells were not Further information on research design is available in the Nature Port- assigned to either category. Finally, we checked that the fraction of Y folio Reporting Summary linked to this article. transcripts was significantly higher in Y-bearing spermatids (Fig. 4c). Bifurcating UMAPs (Fig. 4c) were obtained using X- and Y-linked genes in addition to previously identified highly variable genes to perform Data availability the PCA associated with the UMAP coordinate calculation. For platy- Raw and processed bulk and snRNA-seq data have been deposited pus, X transcripts are separated according to their location on the in ArrayExpress with the accession codes E-MTAB-11063 (human X chromosomes (that is, PARs and SDRs, respectively, as annotated in a snRNA-seq), E-MTAB-11064 (chimpanzee snRNA-seq), E-MTAB-11067 previous study46). In the platypus genome assembly used46, X and Y PARs (bonobo snRNA-seq), E-MTAB-11065 (gorilla snRNA-seq), E-MTAB-11066 are both assigned to the X chromosome. Thus, reported X transcripts (gibbon snRNA-seq), E-MTAB-11068 (macaque snRNA-seq), E-MTAB- may stem from X SDRs, X PARs or Y PARs, whereas reported Y tran- 11069 (marmoset snRNA-seq), E-MTAB-11071 (mouse snRNA-seq), scripts only stem from Y SDRs. Illustrations of human and platypus sex E-MTAB-11072 (opossum snRNA-seq), E-MTAB-11070 (platypus snRNA- chromosomes (with their respective PARs/SDRs) in Fig. 4a are based seq), E-MTAB-11073 (chicken snRNA-seq) and E-MTAB-11074 (chimpanzee, on previous work46,80. gorilla, gibbon and marmoset bulk RNA-seq) (https://www.ebi.ac.uk/ arrayexpress/). All other data are available as Supplementary Infor- Transcriptomal differences between X and Y spermatids mation or available upon request. The testis gene expression at the We identified differentially expressed genes between X- and Y-spermatid single-nucleus level across the 11 studied species can be visualized populations using the FindMarkers function from Seurat62 R package using the shiny app we developed: https://apps.kaessmannlab.org/ (parameters, default). A Wilcoxon rank-sum test was used to calculate SpermEvol/. P values that were adjusted using Bonferroni corrections for several tests. Only genes that were detected as expressed in at least 10% of cells from either of the two populations were tested. Genes that show, on Code availability average, at least a 0.25-fold higher expression (log2-scale) in one of the Custom scripts used to generate the results reported in the paper and populations, that are at the same time expressed in twice the number processed data are available at https://github.com/evo-bio/Spermato- of cells in that population and that have an adjusted P value below 0.01 genesis. were considered to be differentially expressed. We note that several of the differentially expressed genes, including the most significant cases 51. Attar, M. et al. A practical solution for preserving single cells for RNA sequencing. Sci. Rep. in human (Extended Data Fig. 11a), are putatively non-coding, which 8, 2151 (2018). 52. Krishnaswami, S. R. et al. Using single nuclei for RNA-seq to capture the transcriptome of is noteworthy because lncRNAs are typically nuclear81 and hence their postmortem neurons. Nat. Protoc. 11, 499–524 (2016). differential expression levels are unlikely to be offset by transcript 53. Yates, A. et al. Ensembl 2016. Nucleic Acids Res. 44, D710–D716 (2016). exchange between spermatid cells through cytoplasmic bridges. The 54. Marin, R. et al. Convergent origination of a Drosophila-like dosage compensation mechanism in a reptile lineage. Genome Res. 27, 1974–1987 (2017). three most Y-spermatid-specific transcripts are lncRNAs emanating 55. Cunningham, F. et al. Ensembl 2019. Nucleic Acids Res. 47, D745–D751 (2019). from homologous low copy repeats (that is, segmental duplications) 56. Yates, A. D. et al. Ensembl 2020. Nucleic Acids Res. 48, D682–D688 (2019). on chromosomes 13 (FAM230C), 21 (XLOC-095504) and 22 (FAM230F) 57. Martin, M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet.J. 17, 10–12 (2011). that cause genomic disorders by triggering non-allelic homologous 58. Kim, D. et al. TopHat2: accurate alignment of transcriptomes in the presence of insertions, recombination events. They include the FAM230F lncRNA in the q11.2 deletions and gene fusions. Genome Biol. 14, R36 (2013). low copy repeat region on chromosome 22 (22q11.2) that is particularly 59. Pertea, M. et al. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat. Biotechnol. 33, 290–295 (2015). susceptible to non-allelic homologous recombination-generated dele- 60. Trapnell, C. et al. Transcript assembly and quantification by RNA-Seq reveals unannotated tions that lead to various congenital malformation disorders, including transcripts and isoform switching during cell differentiation. Nat. Biotechnol. 28, 511–515 the DiGeorge syndrome, the most frequent microdeletion disorder82. (2010). 61. Habib, N. et al. Massively parallel single-nucleus RNA-seq with DroNc-seq. Nat. Methods 14, 955–958 (2017). MSCI completeness analysis 62. Stuart, T. et al. Comprehensive integration of single-cell data. Cell 177, 1888–1902 (2019). To identify potential MSCI escapee genes, we screened for X-linked 63. McGinnis, C. S., Murrow, L. M. & Gartner, Z. J. DoubletFinder: doublet detection in single- cell RNA sequencing data using artificial nearest neighbors. Cell Syst. 8, 329–337 (2019). genes with a significant increase in transcript abundance from SG to SC 64. Welch, J. D. et al. Single-cell multi-omic integration compares and contrasts features of stages subject to MSCI, which ensures that potential escape genes are brain cell identity. Cell 177, 1873–1887 (2019). 65. Green, C. D. et al. A comprehensive roadmap of murine spermatogenesis defined by 83. R Core Team. R: A Language and Environment for Statistical Computing, https://www. single-cell RNA-seq data resources. Dev. Cell 46, 651–667 (2018). R-project.org/ (R Foundation for Statistical Computing, 2018). 66. Lukassen, S., Bosch, E., Ekici, A. B. & Winterpacht, A. Characterization of germ cell 84. Bakken, T. E. et al. Comparative cellular analysis of motor cortex in human, marmoset and differentiation in the male mouse through single-cell RNA sequencing. Sci. Rep. 8, 6521 mouse. Nature 598, 111–119 (2021). (2018). 67. Street, K. et al. Slingshot: cell lineage and pseudotime inference for single-cell transcriptomics. BMC Genomics 19, 477 (2018). Acknowledgements We thank all members of the Kaessmann group, S. Tirier for discussions 68. Kumar, S. et al. TimeTree 5: an expanded resource for species divergence times. Mol. Biol. and N. Trost for the administration of the Kaessmann laboratory server. Computations were Evol. 39, msac174 (2022). performed on the Kaessmann laboratory server and the bwForCluster from the Heidelberg 69. Lek, M. et al. Analysis of protein-coding genetic variation in 60,706 humans. Nature 536, University Computational Center (supported by the state of Baden-Württemberg through 285–291 (2016). bwHPC and the German Research Foundation grant no. INST 35/1134-1 FUGG). This research 70. Cassa, C. A. et al. Estimating the selective effects of heterozygous protein-truncating was supported by grants from the ERC (grant no. 615253, OntoTransEvol) and German Research variants from human exome data. Nat. Genet. 49, 806–810 (2017). Council (DFG, grant nos. SFB 873 and KA 1710/4-1) to H.K., by the CellNetworks Postdoc 71. Yang, Z. PAML 4: phylogenetic analysis by maximum likelihood. Mol. Biol. Evol. 24, Fellowship and EMBO Long-Term Fellowship to F.M. (grant no. ALTF 591-2017), and by the 1586–1591 (2007). Australian Research Council (grant no. FT160100267) to F.G. and by the Novo Nordisk 72. Van Der Lee, R., Wiel, L., Van Dam, T. J. P. & Huynen, M. A. Genome-scale detection of Foundation (grant no. NNF21OC0069913) to K.A. and (grant no. NNF18OC0031004) to M.H.S. positive selection in nine primates predicts human-virus evolutionary conflicts. Nucleic The use of all other mammalian samples for the type of work described in this study was Acids Res. 45, 10634–10648 (2017). approved by ERC ethics screening panels (ERC starting grant no. 242597, SexGenTransEvolution 73. Kosiol, C. et al. Patterns of positive selection in six mammalian genomes. PLoS Genet. 4, and ERC consolidator grant no. 615253, OntoTransEvol). e1000144 (2008). 74. Shao, Y. et al. GenTree, an integrated resource for analyzing the evolution and function of Author contributions F.M., N.M. and H.K. conceived and organized the study based on H.K.’s primate-specific coding genes. Genome Res. 29, 682–696 (2019). original design. F.M., N.M. and H.K. wrote the manuscript, with input from all authors. F.M. 75. Futschik, M. E. & Carlisle, B. Noise-robust soft clustering of gene expression time-course performed all analyses, developed the shiny app and drew the species icons. K.M. and T.B. data. J. Bioinform. Comput. Biol. 3, 965–988 (2005). collected the samples. R. Behr, P.K, F.G, S.P., I.K., R. Bontrop, T.M.-B., K.A. and M.H.S. provided 76. Hooper, S. E., Weller, H. & Amelon, S. K. Countcolors, an R package for quantification of the samples. F.M., N.M., T.T. and M.S. established snRNA-seq methods. N.M. performed the fluorescence emitted by Pseudogymnoascus destructans lesions on the wing membranes snRNA-seq experiments with support from J.S. C.S., K.M. and T.B. S.B.W. performed the smISH of hibernating bats. J. Wildl. Dis. 56, 759–767 (2020). experiments, and analysed the results together with K.A. E.L. performed the testis genomic 77. Weller, H. I. & Westneat, M. W. Quantitative color profiling of digital images with earth annotations. M.S., M.C.-M, F.L., M.R.B., C.C., K.A. and M.H.S. provided useful feedback and mover’s distance using the R package colordistance. PeerJ 7, e6398 (2019). discussions. 78. Efremova, M., Vento-Tormo, M., Teichmann, S. A. & Vento-Tormo, R. CellPhoneDB: inferring cell-cell communication from combined expression of multi-subunit ligand-receptor Competing interests The authors declare no competing interests. complexes. Nat. Protoc. 15, 1484–1506 (2020). 79. Dimitrov, D. et al. Comparison of methods and resources for cell-cell communication Additional information inference from single-cell RNA-Seq data. Nat. Commun. 13, 3224 (2022). Supplementary information The online version contains supplementary material available at 80. Flaquer, A., Rappold, G. A., Wienker, T. F. & Fischer, C. The human pseudoautosomal https://doi.org/10.1038/s41586-022-05547-7. regions: a review for genetic epidemiologists. Eur. J. Hum. Genet. 16, 771–779 (2008). Correspondence and requests for materials should be addressed to Florent Murat, 81. Derrien, T. et al. The GENCODE v7 catalog of human long noncoding RNAs: analysis of Noe Mbengue or Henrik Kaessmann. their gene structure, evolution, and expression. Genome Res. 22, 1775–1789 (2012). Peer review information Nature thanks Camille Berthelot, Sue Hammoud and Susana Chuva 82. Babcock, M., Yatsenko, S., Stankiewicz, P., Lupski, J. R. & Morrow, B. E. AT-rich repeats de Sousa Lopes for their contribution to the peer review of this work. Peer reviewer reports are associated with chromosome 22q11.2 rearrangement disorders shape human genome available. architecture on Yq12. Genome Res. 17, 451–460 (2007). Reprints and permissions information is available at http://www.nature.com/reprints. Article Extended Data Fig. 1 | Cellular composition and marker gene assessment. (differentiated), and STRA8 spermatogonia; SYCE1 (leptotene), SYCP1 a–e, Clustering and UMAP representations (left), marker gene-based cell type (zygotene), PIWIL1 (pachytene), SYCP2, TANK and AURKA spermatocytes; assignment (right), for a, primates, b, mouse, c, opossum, d, platypus and LRRIQ1 (early), ACRV1 (late) and SPACA1 (late) round spermatids; SPATA3, e, chicken. CLU marks Sertoli cells; TAGLN and ACTA2 peritubular and smooth NRBP1, PRM1 and GABBR2 elongated spermatids. f, Cell-type composition. muscle cells; CD34 and TM4SF1 endothelial cells; APOE and CD74 macrophages; g, Scaled expression of specific genes (y-axis) for all species and cell types STAR and CYP11A1 Leydig cells; GFRA1, PIWIL4 (undifferentiated), DMRT1 (x-axis). Extended Data Fig. 2 | Mammalian testicular gene expression divergence. trees that have the branching pattern as in the majority-rule consensus tree a, Gene expression phylogenies. Neighbor-joining trees based on pairwise shown) are indicated by circles at the corresponding nodes: ≥ 0.9 (white fill). expression level distances (1−rho, Spearman’s correlation coefficient) for the b, Human-macaque (left, n = 304) and human-mouse (right, n = 257) pairwise main testicular cell types, respectively (“Other somatic” refers to peritubular, Spearman correlations after downsampling to the same number of nuclei and smooth muscle, endothelial, macrophage, and Leydig cells as detailed in the UMI across cell types for each species. Boxplots show the median (central legend of Extended Data Fig. 1 and Methods). Trees are drawn to the same scale value); upper and lower quartile (box limits) and 95% confidence intervals (indicated by the scale bar). Bootstrap values (i.e., proportions of replicate (whiskers) for 100 bootstrap replicates. Article Extended Data Fig. 3 | Evolutionary forces. a, Tolerance to functional map to protein coding genes (top) or to intergenic elements (bottom). XLOC_ variants. b, Average normalized ratio of nonsynonymous over synonymous 030429 drives a different mouse pattern. f, Translational efficiency (TE; nucleotide substitutions (dN/dS) of expressed genes across primates. normalized log 2-transformed values) of expressed genes. g, Tissue specificity. c, Percentage of positively selected genes along spermatogenesis. Opossum, h, Time specificity (development). In a, b, and d, human middle temporal platypus, and chicken were not used in the original studies72,73 and so the gyrus (MTG) snRNA-seq data84 was used for non-gonadal (brain) comparisons. 1:1 orthologues in these species were used. d, Phylogenetic age of cellular See legend of main Fig. 2 for further details regarding the plotted measures. transcriptomes. Higher values reflect increased lineage-specific gene Box plots depict the median (center value); upper and lower quartile contributions (that is, younger transcriptomes). e, Percentages of UMIs that (box limits) with whiskers at 1.5 times the interquartile range. Extended Data Fig. 4 | See next page for caption. Article Extended Data Fig. 4 | Examples of gene expression trajectory changes cross sections from human, chimp., and orangutan (from top to bottom) along primate spermatogenesis validated by smISH. Left: Median gene stained with smISH for MYO3B (a) and ADAMTS17 (b) using RNAScope. Red expression trajectory from snRNA-seq along spermatogenesis for human, asterisks indicate on which branch a trajectory change was called. Examples of bonobo, chimp., and gorilla (from top to bottom). For each replicate, mean spermatogonia (SG) and round spermatids (rSD) closeups are provided for expression levels across cells of a given cell type were calculated. Marks visualization purposes. Right: percentage of positive cells in SG and rSD from indicate values for the replicates with the highest and lowest mean expression 10 tubules per section (n = 3 for human, n = 1 for chimp. and orangutan). Box levels, dots indicate the median of mean expression values of the replicates. plots represent the mean (central value) distribution of staining; upper and There are two biological replicates for human, bonobo, and gorilla and three lower quartile (box limits) with whiskers at 1.5 times the interquartile range. biological replicates for chimpanzee. Center: images of seminiferous tubule Extended Data Fig. 5 | Conservation of gene expression trajectories along spermatogenesis (spermatogenesis, from left to right). For each replicate, spermatogenesis. a, Global trajectory conservation score across amniotes for mean expression levels across cells of a given cell type were calculated. Marks changed (n = 237) and conserved (n = 389) trajectories (P < 2.2x10 −16, two-sided indicate values for the replicates with the highest and lowest mean expression Wilcoxon rank-sum test). b, Trajectory conservation score between human and levels, dots indicate the median of mean expression values of the replicates. mouse of genes that lead to either a fertile (n = 1077) or an infertile (n = 108) There are three replicates for chimpanzee; two replicates for human, bonobo, phenotype when knocked out in mouse (P = 0.007, two-sided Wilcoxon rank- gorilla, macaque, and marmoset; one replicate for gibbon. The dashed vertical sum test). a,b, Box plots depict the median (center value); upper and lower lines separate spermatogonia (SG), spermatocytes (SC), round (rSD) and quartile (box limits) with whiskers at 1.5 times the interquartile range. elongated spermatids (eSD). c, Gene expression trajectories (DMRT1, top; IP6K1, bottom) along primate Article Extended Data Fig. 6 | Functions of conserved and changed trajectories. d, The top 20 enriched biological process GO terms in genes showing trajectory a, Top 20 enriched biological process GO terms of genes showing conserved changes in amniotes. e, Human tissue (left) and time (right) specificity for expression trajectories across primates and peak expression in the different changed (n = 237) and conserved (n = 389) trajectories across amniotes, cell types, respectively. b, Top 20 enriched biological process GO terms of respectively (two-sided Wilcoxon rank-sum tests were performed for statistical genes showing conserved trajectories between human and mouse and peak comparisons). Box plots depict the median (center value); upper and lower expression in the different cell types, respectively. c, The top 20 enriched quartile (box limits) with whiskers at 1.5 times the interquartile range. biological process GO terms of genes with trajectory changes in primates. Extended Data Fig. 7 | New genes with crucial spermatogenic functions. For each of the two replicates, mean expression levels across cells of a given cell a, percentages of genes among new genes of different ages (i.e., emergence in type were calculated. Marks indicate values for the replicates and dots indicate last common vertebrate, tetrapod, mammalian, eutherian, or rodent ancestors) the median of the values of the replicates. c, Trees of H2al2a and D1Pas1 and that show infertility phenotypes when knocked out in the mouse (numbers of paralogs from which they originated through the mechanism of RNA-based genes with infertility phenotypes and total number of genes considered for an gene duplication (i.e., D1Pas1 stems from Ddx3y, whereas the precise ancestral age category are provided above each bar). b, Expression level trajectories of parental gene of H2al2a is not easily discernable due to a complex history of a the retrogenes H2al2a and D1Pas1 in spermatogenic and somatic cell types. number of RNA- and DNA-based duplication events). Article Extended Data Fig. 8 | Sertoli-germ cell communications mediated by (details in Supplementary Table 9). b, Overview of all distinct significant ligand-receptor interactions in testis across mammals. a, Significant ligand-receptor interactions between Sertoli and the four principal germ cell ligand-receptor interactions (as assessed by CellPhoneDB78; see Methods) types across species. The number of detected significant interactions for each across species for Sertoli-spermatogonia, Sertoli-spermatocytes, Sertoli- species in both panels is indicated to the bottom left of each plot (Set size). round spermatids and Sertoli-elongated Spermatids communications Extended Data Fig. 9 | Per chromosome testis-specific genes enriched in genome-wide”; i.e., the number of genes on a given chromosome divided by the spermatogonia and Sertoli cells. a,c, Per chromosome percentage of total number of genes in the genome). X chromosomes are colored in red. The expressed testis-specific genes (y-axis) with predominant expression in diagonal shows the numbers expected if testis-specific genes were randomly differentiated spermatogonia (human and macaque) or spermatogonia distributed across the genome. b, Location of platypus testis-specific genes (mouse, opossum and platypus) (a) and Sertoli cells (c) versus the percentage enriched in spermatogonia on chromosomes X1, X3 and X5. of all genes in the genome located on a given chromosome (x-axis: “% genes Article Extended Data Fig. 10 | Classification of X and Y bearing spermatids. The single-nucleus (left) and publicly available single-cell data in mouse43 (a) and scatter and bar plots show the percentages of X and Y transcripts, and the human23 (b). The red and green lines depict the fitted curves for the bimodal distribution of X transcripts (%), respectively, across nuclei or cells in our distributions. Extended Data Fig. 11 | Transcriptome differences between X and Y bearing in Supplementary Table 11). Genes marked with an asterisk are long non-coding spermatids. The heatmaps indicate transcript abundances of significantly RNAs. Genes in red are located on autosomes, those in black on sex chromosomes. enriched genes in X- and Y-bearing spermatids in the different species (a: human, “XLOC” denotes loci annotated in our previous work6; notably, only few XLOC b: mouse, c: opossum, d: platypus). The corresponding bar plots (genes in the loci (opossum: XLOC-076551, XLOC-076500, XLOC-076152, XLOC-076546) same order) show the genes with significantly enriched expression in X- and were found to have coding potential (as assessed by analyses based on Y-bearing cells, respectively using two-sided Wilcoxon rank-sum tests with ribosome profiling data generated in our previous work6), whereas all others Bonferroni correction (ordered according to increasing P-values) (see also lists likely represent lncRNAs. Article Extended Data Fig. 12 | Mammalian sex chromosome evolution. not merely represent genes expressed in SG with stable transcripts still a, Y-to-autosome transcript ratios in individual germ cells along spermatogenesis. detectable in SC) (Methods). Expression differences (loge-fold change, x-axis) For spermatids, only Y-bearing cells are considered. Whiskers depict the of X and autosomal genes between SC and SG were assessed using differential median ± 25th and 75th percentiles. b, Average expression along platypus expression analysis (two-sided Wilcoxon rank-sum tests; p-values on the Y-axis X chromosomes in spermatogonia and pachytene spermatocytes. Genes are are Bonferroni-corrected). Vertical red lines indicate a 0.25-fold expression grouped into 1,000,000 bp bins along the chromosomes. c, To identify potential change (loge-scale); horizontal red lines indicate corrected P < 0.05. Potential MSCI escapees, we screened for X-linked genes in SDRs with a significant increase escapees (i.e., X-linked genes in SDRs with significantly higher expression in in transcript abundance from SG to SC stages subject to MSCI, which ensures SC than SG) are indicated by red arrows. that potential escape genes are indeed actively transcribed in SC (i.e., they do