Abstract
Summary
Male germ cells share a common origin across animal species, therefore they likely retain a conserved genetic program that defines their cellular identity. However, the unique evolutionary dynamics of male germ cells coupled with their widespread leaky transcription pose significant obstacles to the identification of the core spermatogenic program. Through network analysis of the spermatocyte transcriptome of vertebrate and invertebrate species, we describe the conserved evolutionary origin of metazoan male germ cells at the molecular level. We estimate the average functional requirement of a metazoan male germ cell to correspond to the expression of approximately 10,000 protein-coding genes, a third of which defines a genetic scaffold of deeply conserved genes that has been retained throughout evolution. Such scaffold contains a set of 79 functional associations between 104 gene expression regulators that represent a core component of the conserved genetic program of metazoan spermatogenesis. By genetically interfering with the acquisition and maintenance of male germ cell identity, we uncover 161 previously unknown spermatogenesis genes and three new potential genetic causes of human infertility. These findings emphasize the importance of evolutionary history on human reproductive disease and establish a cross-species analytical pipeline that can be repurposed to other cell types and pathologies.
Main Text
Understanding what defines the uniqueness of a given cell type out of the 843 predicted cellular fates in the human body is a complex and fascinating problem1. Through Conrad Waddington’s foundational work, we have come to appreciate that developmental trajectories ultimately dictate cell type identity via the establishment of specific transcriptional programs2. The fact that transcriptomes tend to cluster by tissue type rather than by species3 clearly indicates that gene expression identity can be maintained across many million years of evolutionary divergence. This echoes the modular nature of eukaryotic biological processes, whose intervening macromolecular complexes are typically built by the addition of younger components to a core block of ancient subunits4.
The emergence of germ cells is considered one of the first cell type specializations in metazoan history5. Since the capability to undergo both sexual reproduction and gametogenesis were already present in the unicellular ancestor of all Metazoa6, the split between germ line and soma presumably provided early multicellular organisms increased robustness against mutations while minimizing genetic conflict between different cell lineages7. Yet, reproductive protein genes are some of the most rapidly evolving across the tree of life8, leading to substantial variation in male reproductive tissues even between closely related species9,10. At the molecular level, this variation can be mainly traced back to rapid evolutionary changes in gene expression in the testis11 - likely facilitated by widespread transcriptional leakage in male germ cells12 - and to the preferential emergence of new genes in this organ13, where they often acquire functions in spermatogenesis14. Although the rapid divergence of genes directly involved in reproduction has traditionally been interpreted from an adaptationist stance, relaxed selection and drift have been recently proposed to account for this pattern15. Central to this debate is the often-overlooked contribution of old genes for germ cell development and function.
We have recently observed, in a wide range of plant species, a substantial contribution of old genes to the pollen transcriptome, suggestive of an ancient transcriptional program common to plant male gametes16. Such observation finds parallel in the concept of a metazoan germ line multipotency program17, and is supported by transcriptional similarities between equivalent spermatogenic cell types across mammalian species18–20. This led us to explore the possibility that, despite the rapid evolution of reproductive protein genes, the functional basis of metazoan male germ cell identity relies on an old, evolutionarily conserved genetic program that can provide relevant insight into spermatogenesis and human infertility.
To test this hypothesis, we devised an interdisciplinary research platform based on four combined approaches (Fig. 1a). First, we determined, through comparison of genome sequences, the age of the gene expression program of male germ cells from three evolutionarily distant metazoan species: humans (Homo sapiens), mice (Mus musculus) and fruit flies (Drosophila melanogaster). Then, we used network science to infer the significance of deeply conserved genes within the context of the complex male germ line transcriptome. Subsequently, through developmental biology (in vivo RNAi in fruit fly testes), we defined the role of a key subset of the conserved germ cell transcriptome in male reproductive fitness. Finally, we merged this information with clinical genetics to identify new potential causes of human infertility. Overall, we show that deeply conserved genes play a prominent role in male germ cell regulation and that the disruption of this ancient genetic program leads to human reproductive disease.
A more conservative estimate of the complexity of the male germ cell transcriptome
Male germ cell development is divided into three conserved stages21. The first is the pre-meiotic stage and corresponds to the mitotic expansion of committed precursors (spermatogonia). Meiosis defines the second stage, with the newly differentiated spermatocytes undergoing reductive division. In the third stage, the post-meiotic cells (spermatids) embark on a cytodifferentiation program that culminates in the formation of mature male gametes.
To understand to what extent male germ cell transcription quantitatively differs from that of somatic lineages, we collected previously published high-quality RNA-Seq datasets from pre- meiotic, meiotic and post-meiotic germ cell populations, and compared them with representative somatic cell types of the primary embryonic layers: neurons (ectoderm), muscle (mesoderm) and enterocytes (endoderm; Fig. 1b and Sup. Table 1). We centered our analysis on three evolutionarily distant gonochoric species with excellent genome annotations - humans, mice and fruit flies (the last two also being well-established animal models for the functional study of human spermatogenesis) - and observed that, as previously reported12,22, male germ cells have a generally more diversified transcriptome than their somatic counterparts (as measured by the percentage of the entire genome that each cell type expresses). Yet, this increased complexity was mainly restricted to extremely permissive minimum average expression cut-offs ranging from >0.01 to >0.5 transcripts per million (TPM), suggestive of a significant degree of transcriptional noise or other previously hypothesized translation-independent roles for lowly expressed transcripts23. Indeed, the use of a more suitable threshold for protein-coding transcripts (TPM >1)24 revealed that, on average, metazoan male germ cells tend to express approximately 10,000 genes (9,468 ± 498 to 12,195 ± 1,144, depending on the species), corresponding to half of the entire protein-coding genome (48.0 - 67.9%; Fig. 1b and Sup. Table 1). In both humans and mice, this range is comparable to that observed in somatic cell types such as neurons and enterocytes (53.4 - 53.7% and 37.3 - 45.6% of the genome, respectively). Only in fruit flies, a species with a smaller genome (13,947 protein-coding genes vs. 22,802 and 22,287 in humans and mice, respectively), the diversity of the male germ cell transcriptome diverged substantially from that found in representative somatic tissues (with an average of 67.9 ± 3.6% of the genome expressed in germ cells vs. 33.3 ± 12.5% in the soma). It should be noted that this percentage may be higher in fruit flies due to the presence of transcripts of somatic origin in the germ cell dataset (see Sup. Table 1). Collectively, these observations suggest that the average functional requirement of a metazoan male germ cell corresponds to the expression of approximately 10,000 protein-coding genes. Such a number brings a significant but largely overlooked functional constraint to the general notion that almost all protein-coding genes are transcribed in the male germ line.
The male germ cell transcriptome has an old evolutionary origin
We next assessed the contribution of deeply conserved genes to the male germ cell transcriptome. For this we used phylostratigraphy, a technique that determines the evolutionary age of groups of homologous genes (defined as orthogroups in our analysis, including both orthologs and paralogs) by mapping their last common ancestor in a species tree25. We assembled our tree based on the proteomes of 25 phylogenetically representative eukaryotic species, and assigned each orthogroup to the oldest phylogenetic node (phylostratum) it could be traced back to (Fig. 1c and Sup. Fig. 1). This way, phylostrata ranked 1 to 5 were the oldest and contained orthogroups common to all Metazoa, while phylostratum 16 contained the youngest, species-specific orthogroups. A total of 113,757 orthogroups were identified in the 25 representative eukaryotes, 85.5% of which (97,270) were species-specific when adding all phylostrata. By using TPM >1 as minimum expression cut-off (to minimize the effects of widespread, low-level transcription in the male germ line), we observed that 48.0 - 67.9% of the entire protein-coding genome was expressed in male germ cells, depending on the species. Despite the rapid divergence typically associated with reproduction-related genes, we found that the majority (65.2 to 70.3%) of genes expressed in male germ cells mapped to the oldest-ranking phylostrata containing orthogroups common to all Metazoa, henceforth referred to as deeply conserved genes (Fig. 1d). Indeed, in all three tested species, this percentage was not significantly different from that recorded in the transcriptome of the somatic cell types (also at a threshold of TPM >1). The fraction of deeply conserved genes that were ubiquitously expressed (proxy to their involvement in cellular housekeeping processes) varied from 49.7 to 63.4% (Sup. Fig. 2a-b). This strongly suggests that almost half of all deeply conserved genes expressed by male germ cells are involved in more specific roles than just the maintenance of basal cellular functions.
By summing the products of the age of all expressed genes and their expression levels at a given developmental stage (a metric known as the transcriptome age index - TAI26) we determined the transcriptome age of the different male germ cell stages (Fig. 1e). As recently reported in mammalian spermatogenesis20, meiotic and pre-meiotic cells across all three tested species had lower TAIs than post-meiotic cells, indicative of older transcriptomes. This trend was less obvious in the fruit fly as in this species post-meiotic transcription is largely residual27. Collectively, we observed that, both in vertebrates and invertebrates, the male germ cell transcriptome has an old evolutionary origin that is tempered by the increased expression of younger genes at later developmental stages. Our data, indicative of a core spermatogenic program spanning more than 600 million years of evolution, significantly expand the breadth of previous observations in the mammalian lineage18–20, while emphasizing the potential relevance of using conserved expression between distant animal species to gain insight into the molecular basis of human infertility.
Deeply conserved genes are central components of the male germ cell transcriptome
We next addressed the possible significance of the abundant expression of deeply conserved genes in the male germ line. For this, we focused on meiosis, a fundamental process conserved through eukaryotes. Besides the existence of an ancient meiotic toolkit required for recombination and reduction division across sexually reproducing species28, the mitosis-to- meiosis transition in the male germ line is also home to a transcriptional burst associated with the acquisition of male germ cell identity29. The transcripts produced in this burst are mainly required for post-meiotic development, with their translation being shifted towards the later spermatogenic stages30. Accordingly, when analyzing the expression pattern of all human and mouse genes annotated to the “spermatid differentiation” Gene Ontology (GO) term (GO:0048515), we observed that only a small fraction of these was expressed primarily/exclusively after meiosis [5.7% (13 out of 228) and 8.7% (25 out of 289) in humans and mice, respectively]. Similar to previous reports in fruit flies27,29,31, most of the mammalian spermatid differentiation genes were found to be expressed already at the meiotic stage, with 82.9% (189 out of 228, in humans) and 77.9% (225 out of 289, in mice) being expressed in spermatocytes or in spermatocytes and spermatogonia. These results show that the meiotic transcriptome provides a particularly suitable entry-point into the genetic basis of the core spermatogenic program. Hence, we assembled protein-protein interaction (PPI) networks based on the transcriptome data of human, mouse, and fruit fly spermatocytes, where nodes represent all expressed genes and edge weights indicate the probability of the connected genes contributing to a specific function according to the STRING database32 (Fig. 1f). The structure of these networks reflects the multiple genetic interdependencies responsible for cellular function, as illustrated by the characteristic clustering of functionally related genes into topologically defined modules33. Network edges were filtered to remove low confidence interactions (those with a combined confidence score <0.5, Sup. Fig. 3a-c), and, as before, only genes expressed at TPM >1 were included to minimize the confounding effects of widespread low-level transcription. The resulting networks contained between 7,961 and 11,322 genes, depending on the species, and an average of 301,036 ± 66,416 edges.
Consistent with phylostratigraphy, the spermatocyte PPI networks included a substantial fraction of genes conserved across all Metazoa (45.1 to 58.1% of all genes; Fig. 1g). As a starting point for our analyses, we set out to determine to what extent the previously reported increased connectivity of conserved genes in PPI networks34,35 is maintained in a cell type characterized by rapid rates of evolution. For comparison purposes, we assembled, as before, PPI networks based on the transcriptome data of human, mouse, and fruit fly enterocytes (a somatic cell type with a comparably young transcriptome; Sup. Fig. 2b). These somatic cell networks, albeit smaller due to lower transcriptome diversity in enterocytes (Fig. 1b), were comparable in size to those of the spermatocytes, containing between 3,310 and 9,837 genes, depending on the species, and an average of 203,296 ± 117,545 edges (Sup. Fig. 3a-d). Similarly to what we observed in enterocytes, deeply conserved genes in the spermatocyte PPI networks were significantly more connected than their non-conserved counterparts (higher degree centrality), and their interactors were themselves more connected than those of non-conserved genes (higher page rank, Kolmogorov-Smirnov test for the two analyses; Fig. 1h; Sup. Fig. 4). Although the connectivity properties of the spermatocyte conserved genes were lower than those expressed in enterocytes, they were still salient enough for machine learning classifiers to predict reliably if a gene was conserved just based on spermatocyte network features. More specifically, using a Random Forest classifier, the AUC (area under the receiver operating characteristic curve) score was 0.74, 0.75, and 0.82 in the human, mouse, and fruit fly spermatocyte datasets, respectively (Fig. 1i). Classification performance was equally high when using a different supervised learning algorithm (linear Support Vector Machine) and an alternative performance metric (precision and recall; Sup. Fig. 5).
To address a possible ascertainment bias associated with more available information on conserved genes in the STRING database, we tested to what extent degree centrality and page rank were affected by network rewiring. In this approach, a variable percentage (20 to 100%) of all edges is randomly shuffled across the network, thus diluting any latent biases in the datasets. We observed that both network properties remained higher in conserved genes even when 80% of all edges were rewired in the spermatocyte PPI networks (Sup. Fig. 6). Such result suggests that the increased connectivity of conserved genes is mostly driven by their intrinsic properties, rather than by differences stemming from the amount of source data available. Based on these analyses, we conclude that despite the testis being a rapidly evolving organ and a preferential birthplace for new genes, deeply conserved genes remain central components in the male germ cell transcriptome of evolutionarily distant species. Accordingly, we decided to explore this increased connectivity as a means of gaining insight into the conserved molecular basis of metazoan spermatogenesis.
A core component of the ancient genetic program of spermatogenesis
We have previously shown that it is possible to simplify the complexity of biological networks by removing redundancy and only retaining the most relevant interactions that sustain the dynamics of the system36. To do so, one can extract from a network its distance backbone, i.e., the small subset of metric edges that is sufficient to compute all shortest paths in the original network37.
This technique retains the minimum number of required edges to preserve all network nodes, as well as all of the natural hierarchy of connections and community structures38. In the context of biological networks, such distance backbone uncovers the subgraph of gene interactions more likely to convey the key functional characteristics of the gene expression program39.
By extending orthology to the backbone, we have developed a new type of multilayer network: the orthoBackbone. It expands on the concept of interologs40,41, as it consists of only the edges that connect the same pair of orthologs in the network backbone of different species (Fig. 2a). Since this technique excludes all backbone edges (and corresponding nodes) that do not involve the same orthologs across multiple species, it ultimately identifies a subset of conserved genes whose functional interactions at the protein level are part of the distance backbone across evolution. In this regard, the orthoBackbone offers access to what can be considered core features of the genetic program of a given cellular state. Importantly, the use of fruit flies (an animal without male crossing over and recombination) in this inter-species approach significantly reduces the contribution of the ancient meiotic toolkit to the orthoBackbone, thus emphasizing functional interactions associated with spermatogenesis but not with meiosis per se.
We observed that the orthoBackbone drastically reduced the size of the spermatocyte PPI networks, retaining only 1.7 to 2.7% of all edges, depending on the species (Fig. 2b). This almost residual fraction of edges nevertheless contained 70% of all deeply conserved genes in the network (Fig. 2c and Sup. Table 2), further illustrating that strong functional relationships between deeply conserved genes are also maintained in the male germ line. GO term enrichment analysis for biological processes revealed that orthoBackbone genes were preferentially involved in gene expression and protein regulation, in contrast with the other (not retained) conserved genes that were mainly associated with cell signaling pathways (Fig. 2d and Sup. Fig. 7). To better understand to what extent the orthoBackbone can provide insight into cell-type specific processes, we determined the degree of similarity between the spermatocyte and enterocyte orthoBackbones. We noted a modest overlap between both subnetworks, with just 16.2% to 18.0% of common edges depending on the species (against background expectations of 33.6 to 61.5%, based on Jaccard similarity index scores between the PPI networks of both cell types; Sup. Fig. 3d).
Using GO annotations, we further selected spermatocyte orthoBackbone interactions involving gene expression regulators previously linked to spermatogenesis. By doing so, we obtained a subset of 79 functional interactions between 104 human genes (out of the total of 3,596 in the orthoBackbone; Fig. 2e). These were ascribed to a wide gamut of regulatory processes and involved noteworthy spermatogenesis genes such as RFX2 (transcription), CDYL (chromatin remodeling) and BOLL (translation), among others42–44. Of the 104 genes, 34 were testis- specific/enriched while 60 had low tissue specificity, consistent with previous observations that cell identity relies on the integration of cell type-specific genes with more broadly expressed regulators45. Based on its conservation in the orthoBackbone of evolutionarily distant species and on the significant spermatogenic role of the intervening genes, we posit that this subnetwork of 79 functional interactions involving 104 genes likely represents a core component of the conserved genetic program of spermatogenesis. This number fits well with previous estimates, based on expression trajectories of 1:1 orthologous genes, of an ancestral program of amniote spermatogenesis constituted by 389 genes20.
Male germ cell identity has a broad functional basis
We next set out to determine the functional consequences of disrupting the spermatocyte orthoBackbone. Since the latter represents a still sizeable interaction subnetwork consisting of more than 3,000 genes, we decided to select candidate genes likely to affect the acquisition and maintenance of male germ cell identity. For this, we again took into account how the mitosis-to- meiosis transcriptional burst activates the expression of spermatogenesis genes - a process that has been likened to that of cellular reprogramming46,47. We reasoned that genes that are upregulated at meiotic entry and/or downregulated at meiotic exit (thus part of the transcriptional burst) likely represent instructive elements for the male germ cell state, hence preferential routes to tamper with this cellular identity.
Through differential gene expression analysis, we identified, in human and mouse spermatocytes, all expressed genes that shared a similar upregulation at meiotic entry and/or downregulation at meiotic exit. Of these 970 mammalian differentially expressed genes (DEGs), we selected as candidates for functional assessment only those that were also expressed in fruit fly germ cells (Fig. 2f). We did not take into account differential expression in the latter species, since the largely residual levels of post-meiotic transcription in fruit fly spermatogenesis thwarts direct comparisons with the mammalian system27. The resulting 920 fruit fly genes (homologous to 797 and 850 in humans and mice, respectively) were silenced specifically at the mitosis-to-meiosis transition by Drosophila in vivo RNAi, using the well-established bam-GAL4 driver48. Of these, a total of 250 (27.2%) were essential for male fertility, as their silencing resulted, upon mating with wildtype females, in egg hatching rates below the cut-off of 75% (>2 standard deviations of the mean observed in negative controls: 91.6 ± 8.5%; Fig. 2g). We observed that 76.0% of the hits (190 out of 250) were part of the orthoBackbone, compared with 60.0% (402 out of 670) in genes without any obvious effect on male reproductive fitness upon silencing, with this difference being statistically significant (two-sided Fisher’s exact test P<0.0001). Cytological analysis of all 250 genes by testicular phase-contrast microscopy in the RNAi males revealed diverse origins for the infertility phenotype, with the earliest manifestation of cellular defects ranging from the pre-meiotic to the mature gamete stage (Fig. 2h). By exploring publicly available information (see Methods), we determined that 161 (64.4%) of all hits had never been previously associated with male reproduction in any species (Fig. 2i). Accordingly, these 161 new deeply conserved spermatogenesis genes (homologous to 179 and 187 in humans and mice, respectively) represent both a significant advance in our understanding of the genetic basis of male germ cell development, and a valuable resource to explore from a precision medicine perspective. To facilitate open access to this information, we made all data available in the form of a user-friendly gene browser (https://pages.igc.pt/meionav; Fig. 2j). Overall, by merging our results with previously published data we conclude that at least 43.5% of the mitosis-to-meiosis transcriptional burst (our 250 hits + 150 previously reported male fertility genes that gave a negative result in our assay) is required for a diverse range of germ cell functions. This illustrates the pervasive influence of the mitosis-to-meiosis transition on multiple facets of the spermatogenic program and argues for a broad functional basis underlying male germ cell identity.
An ancient regulator of human spermatogenesis
One of our 161 newly identified spermatogenesis genes - the Drosophila RING finger protein 113 (dRNF113, a spliceosomal component also known as mdlc49) - emerged as a particularly interesting germ cell regulator from a clinical perspective. More specifically, by analyzing our in-house whole exome database containing sequencing data from 74 cases of human male meiotic arrest, we identified an infertile man harboring a homozygous loss of function (LoF) variant in a human paralog (RNF113B). This frameshift variant c.556_565del;p.(Thr186GlyfsTer119) leads to the abrogation of the protein’s two functional domains (a C3H1-type zinc finger and a RING finger) and to its truncation (Fig. 3a and Sup. Fig. 8a). The identified man (M1911) is part of a consanguineous family of Middle-Eastern ancestry and shared the homozygous state of the RNF113B LoF variant with his equally infertile brother, but not with his fertile heterozygous brother (Sup. Fig. 9a-c and Sup. Information). We excluded the possibility that this genetic variant is an indicator of another, nearby variant that is homozygous in the affected individuals and represents the actual cause of infertility in both men (Sup. Information). Remarkably, the results from the testicular biopsy of M1911 revealed an equivalent meiotic arrest phenotype to that observed in dRNF113-silenced fruit flies. Indeed, spermatocytes were, by far, the predominant cell type in the male gonads of both species: they were the most advanced cell stage observed in 89.0% of all assessed human seminiferous tubules (vs. 9% in controls) and occupied an average of 64.1% of the entire testis area in fruit flies (vs. 12.3% in controls; Fig. 3b-c and Sup. Fig. 8b-c). Early (round) spermatids were practically absent in both species, with cellular debris accumulating in the post-meiotic region of the fruit fly gonad. The similar reproductive impairment phenotypes recorded in both species highlights the often-complex relationship between tissue specificity and functional significance. In insects, dRNF113 is a broadly expressed gene essential for organismal viability49, whose particular spermatogenic function was uncovered by germ cell-specific silencing. Both roles were functionally uncoupled in primates via the existence of two RNF113 copies, with the correctly spliced form of RNF113B being mainly testis-specific50. The existence of an independent duplication event in the rodent lineage offered us the opportunity to further explore this concept. Compared with its human ortholog (RNF113B), the mouse duplication (Rnf113a2) is broadly expressed, similarly to Rnf113a1, human RNF113A and fruit fly dRNF113. Yet, the generation of a whole-body Rnf113a2 knockout (KO) mouse revealed a comparable phenotype to that of the human RNF113B LoF variant: severe spermatogenic impairment without any overt somatic defects. When compared to wildtype littermates, male homozygous Rnf113a2KO/KO mice were sterile and had visibly smaller testes that were largely devoid of germ cells (SCO: Sertoli cell-only phenotype) except for the rare presence of spermatogonia and spermatocytes (Fig. 3d and Sup. Fig. 8d-f).
Strengthening our interest in the RNF113 genes was their possible involvement with male germ cell identity. By analyzing our recently published single-cell RNA-Seq dataset of normal human spermatogenesis51, we observed that RNF113B was predominantly expressed at the mitosis-to- meiosis transcriptional burst, peaking at the diplotene stage of prophase I (Fig. 3e). This mirrored the protein localization pattern of fruit fly dRNF113, characterized by a substantial nuclear accumulation in primary spermatocytes (Fig. 3f). Also in mice, Rnf113a2 was expressed at the transcriptional burst, with an earlier expression in undifferentiated spermatogonia being the likely cause of the slightly more severe histological phenotype of the Rnf113a2KO/KO mice (Sup. Fig. 8g). Consistent with their role in gene expression, the analysis of M1911’s testicular transcriptome via RNA-Seq revealed the deregulation of 21.7% of all expressed genes (Fig. 3g). Notably, the deregulation detected in M1911’s testicular transcriptome had a clear impact on the spermatocyte orthoBackbone, with 30.0% of its network edges being disrupted (i.e., containing at least one deregulated gene). A similar effect was patent in the dRNF113 RNAi, with 20.3% of the fruit fly testicular transcriptome being deregulated, resulting in the disruption of 26.3% of the orthoBackbone (Fig. 3h). For comparison, the silencing of Prp19, another transcriptional burst gene and spliceosomal component that RNF113 proteins associate to52, had a lower effect on the orthoBackbone with 15.2% of disrupted edges (significantly lower than in the dRNF113 RNAi; two-sided Fisher’s exact test P<0.0001), despite a similar meiotic arrest phenotype (Sup. Fig. 10a-d). Based on the above, we conclude that the RNF113 genes have retained a key role in spermatogenesis for more than 600 million years of evolution. Furthermore, they provided us the means to interfere with a sizable fraction of the male germ cell orthoBackbone in humans and fruit flies - an advantage we next explored from a clinical perspective.
The orthoBackbone as an ancillary tool in clinical genetics
Male infertility is a complex human disease with a poorly defined genetic component. This contributes to a low causative diagnostic yield (typically just 30%), and to a paucity of clinically validated genes [currently just 104, in contrast with the more than 700 already associated with other clinical disorders such as intellectual disability53,54]. Since male infertility affects up to 12% of all men55, addressing such knowledge gap is an issue of clear medical importance. To narrow the gap, we explored the possibility that spermatogenesis is particularly sensitive to disturbances in the male germ cell orthoBackbone. Thus, we harnessed the sizeable effect of RNF113B on this network as a means of identifying additional genetic causes of human infertility. Indeed, we noted that 38.0% (30 out of 79) of the functional interactions that define our previously identified core component of the conserved genetic program of spermatogenesis were disrupted in M1911 (i.e., these interactions contained at least one differentially expressed gene; Sup. Fig. 11). In particular, we observed that 19.2% (20 out of 104) of this gene subset was found to be downregulated, compared with 10.0% (1,685 out of 16,824) of downregulated genes in the testicular transcriptome of M1911, with this difference being statistically significant (two-sided Fisher’s exact test P=0.0046; Fig. 3g). Therefore, we posited that the similar spermatogenic impairment recorded in the human RNF113B LoF variant and on the fruit fly dRNF113 RNAi ultimately reflected the downregulation, in both species, of a common set of orthoBackbone genes.
By defining the overlap between differentially expressed orthoBackbone genes in the testicular transcriptome of M1911 and of the dRNF113 RNAi, we identified 61 conserved human genes that were similarly downregulated in both species. These formed a connected network (based on STRING data) suggestive of their functional co-involvement in related biological processes (Fig. 3i). Of the 61 genes, 27 had already been linked to male germ cell defects in different species: four were clinically validated male infertility genes (CDC14A, CFAP91, DNAI1 and DNAI2), and the other 23 had been previously reported in animal models (see Methods). The fact that among the latter was BOLL, one of the oldest known metazoan gametogenesis genes44, was a particularly noteworthy observation. We next tested if these 27 genes could be used to identify new genetic causes of human infertility. For this, we analyzed whole-exome sequencing data of 1,021 azoospermic men from the MERGE cohort56. Filtering these exomes for LoF variants in the 27 selected orthoBackbone genes revealed two new human male infertility candidate genes: HSPA2 and KPNA2.
HSPA2 - heat shock protein family A member 2 - is a molecular chaperone of the highly conserved 70-kDa heat shock protein (HSP70) gene family. HSP70 members are involved in cellular proteostasis from bacteria to human, with the mouse and fruit fly HSPA2 homologs (Hspa2 and Hsc70-1 to Hsc70-5, respectively) being required for meiotic progression past the primary spermatocyte stage57,58. Indeed, we observed that the silencing of Hsc70-1 in the fruit fly testis, also with the bam-GAL4 driver, resulted in an equivalent meiotic arrest to that of dRNF113, strongly suggesting that the misexpression of this chaperone is a central element of the latter phenotype (Sup. Fig. 12a-b). Two different heterozygous LoF variants were detected in our male infertility cohort, with this gene having a predicted autosomal-dominant inheritance (Sup. Fig. 12c-f, Sup. Information and Methods). The first variant, detected in individual M1678, is the early stop-gain variant c.175C>T;p.(Gln59Ter) that truncates more than 90% of the protein, likely leading to nonsense-mediated decay. There is only one individual listed in the gnomAD database with this variant (minor allele frequency = 0.0004%). The second, identified in individual M2190, is the frameshift variant c.1806dup;p.(Glu603ArgfsTer81) which affects the distal end of the protein’s nucleotide-binding domain. This variant is not listed in the gnomAD database. Histopathological analysis of M2190’s testicular tissue revealed a SCO phenotype in 259 assessed seminiferous tubules, a likely aggravation of the extensive cell death accompanied by tubule vacuolization reported in Hspa2 mutant mice59 (Fig. 3j).
KPNA2 - karyopherin subunit alpha 2 - is a nuclear importin that functions as an adapter protein in nucleocytoplasmic transport. Its mouse ortholog (Kpna2) is required for the nuclear accumulation of Hop260, a conserved regulator of meiotic progression from yeast to mammals61,62. The main function of Hop2 is to repair meiotic DNA double-strand breaks (DSBs), with the corresponding mouse mutant being characterized by a primary spermatocyte arrest coupled to extensive cell death. Since male fruit flies dispense with meiotic DSBs63, the inclusion of KPNA2 in the male germ cell orthoBackbone seems counter-intuitive. Yet, Drosophila importin alpha 2 (dKPNA2, also known as Pendulin) is required for post-meiotic development64, with its silencing specifically at meiotic entry being associated with aborted spermiogenesis (Sup. Fig. 13a-b). This post-meiotic function might also be present in its mammalian orthologs, as suggested by the nuclear localization of Kpna2 in elongating mouse spermatids60. We detected two heterozygous KPNA2 LoF variants in our cohort of infertile men, with this gene also having a predicted autosomal-dominant inheritance (Sup. Fig. 13c-g, Sup. Information and Methods). The splice- site variant (c.667-2A>G;p.?) in individual M1645 is predicted to disrupt the correct splicing of intron 6 due to the loss of an acceptor site. Segregation analysis revealed the de novo occurrence of this variant, strongly supporting both its pathogenicity and our recent results, pointing to the important contribution of heterozygous de novo mutations to male infertility65,66. The other KPNA2 LoF variant, in individual M2098, is the frameshift variant c.1158_1161del;p.(Ser387ArgfsTer14) which affects 27% of the protein, including its minor nuclear localization signal binding site. Both KPNA2 variants are not listed in the gnomAD database. The available testicular histopathology report for the latter individual lists a SCO phenotype as the cause of the azoospermia, again suggestive of the possible deterioration of an initially meiotic phenotype (Fig. 3j).
Based on three in silico prediction tools (Sup. Information and Methods), both HSPA2 and KPNA2 can represent haploinsufficient genes and/or have autosomal dominant inheritance, with the evidence being stronger for KPNA2. Together with the identified de novo occurrence of one of its LoF variants, this makes KPNA2 a strong novel autosomal-dominant candidate gene while more data is still needed before drawing similar conclusions on HSPA2.
In summary, by experimentally disrupting the male germ cell orthoBackbone across evolutionarily distant species, we were able to uncover two new candidate genes for human infertility (HSPA2 and KPNA2) in addition to RNF113B, affecting a combined total of 5 individuals. This successful merger between basic and clinical research highlights the advantages of comparative biology as means of dampening inter-species differences in reproductive physiology, while providing a conceptual framework for a more efficient prioritization of clinically-relevant genetic variants in human reproductive disease.
Discussion
In summary, we have identified a conserved male germ cell gene expression program, spanning more than 600 million years of evolution, that can provide clinical insight into the molecular basis of human infertility. We estimate the average functional requirement of a metazoan male germ cell to correspond to the expression of approximately 10,000 protein-coding genes, the majority of which are deeply conserved. Indeed, despite the testis being a rapidly evolving organ and frequent birthplace of genetic novelty, deeply conserved genes are central components of the male germ cell transcriptome, with younger genes being preferentially expressed in late spermatogenesis. We further pinpoint 79 functional interactions between 104 transcriptional regulators that, owing to their salient network properties across different species, likely define a core component of the ancient genetic program of spermatogenesis. By functionally interfering with 920 deeply conserved genes associated with male germ cell identity, we uncover 161 previously unknown spermatogenesis genes required for a broad range of germ cell functions, as well as three new potential genetic causes of human infertility.
We propose that the transcriptional identity of metazoan male germ cells is built around a relatively small network of deeply conserved gene interactions with an overarching functional impact. These observations may provide a unifying genetic basis for the deep conservation of fundamental germ cell biological processes such as sperm motility and gamete fusion67,68. The existence of a shared genetic identity in metazoan spermatocytes can be regarded as a ramification of an ancestral multipotency program already present in germ line precursor cells17. Indeed, even in species with divergent germ line segregation strategies, conserved functional interactions at the post-translational level are required for primordial germ cell specification69,70. It should be noted that our conclusions were based on a small fraction of the diversity of all metazoan species, as a result of data availability limitations and our focus on well-characterized gonochoric animal models of human spermatogenic impairment. The advent of more comprehensive gene expression datasets from less-studied organisms will surely provide a more encompassing context to cross-species comparisons of this nature.
The actual benefit of comparative biology for the identification of new genetic causes of human disease is often a contentious topic71. Central to this debate is the moderate success in translating animal data to the clinical setting, coupled with the fact that the vast majority of human genetic variants are not shared with other species72. By focusing on the deep evolutionary past of human spermatogenesis, we have identified 79 functional interactions that likely represent a core component of the conserved genetic program of male germ cells, and 179 novel functionally validated candidate genes. Exploring to what extent these core functional interactions reflect direct protein-protein contacts is an important area for future research. In this regard, the recent development of robust in silico predictive systems, such as AlphaFold-Multimer, will greatly facilitate the implementation of such studies73.
The translational application of our cross-species platform in human reproductive healthcare has so far uncovered three new genes associated with human male infertility (RNF113B, HSPA2 and KPNA2). The assumed autosomal dominant inheritance of HSPA2 and KPNA2 awaits further confirmation by additional data from independent cohorts and segregation analyses, but stands in line with other well-established male infertility genes such as DMRT1, NR5A1 and SYCP253. Also, this information can potentially be harnessed in other clinically relevant contexts: the spermatocyte orthoBackbone can be further explored to pinpoint key pathways that may facilitate the reconstitution of meiotic entry in vitro, as well as to maximize the effectiveness of exome sequencing approaches in infertile men through the inclusion of the orthoBackbone genes in prioritization lists for disease-causing variants. Conversely, future iterations of our analytical pipeline can be specifically tailored for the identification of newly emerged genes that may have acquired key roles in reproductive processes, such as those involved in species-specific gamete recognition74.
Taken together, our results emphasize the often-overlooked contribution of evolutionary history to human disease and illustrate how interdisciplinary research can significantly expand our knowledge of fundamental cellular processes. Accordingly, all the code required for repurposing our analytical pipeline to other cell types and pathologies is available as an open-access resource at https://github.com/rionbr/meionav. These resources will likely contribute to a renewed appreciation of comparative biology in the medical field.
Acknowledgements
We gratefully acknowledge all men that gave their informed consent to be included in the Male Reproductive Genomics (MERGE) study cohort. The authors wish to thank Gabriel Martins and José Marques from the Imaging Unit of Instituto Gulbenkian de Ciência (IGC, Portugal) for assistance in acquiring and processing microscopy data. We acknowledge the IGC’s Histology Facility and Genomics Unit, respectively, for the analysis of mouse testicular tissue and preparation of the fruit fly RNA-Seq libraries, Nicole Terwort (University of Münster) for preparing the RNA extractions from human testes, and the Core Facility Genomics of the Medical Faculty of the University of Münster for conducting library preparation and RNA sequencing in samples from individual M1911. Renate Renkawitz-Pohl (Philipps-Universität Marburg, Germany) and Chris Doe (University of Oregon, USA) kindly provided us the bam-GAL4 line and the anti- dRNF113 antibody, respectively. We thank Élio Sucena, Raquel Oliveira and Mónica Bettencourt- Dias (all from the IGC), Patrícia Beldade (Faculty of Sciences, University of Lisbon), and Peter Ellis (University of Kent, UK) for insightful discussions.
Funding
Fundação para a Ciência e a Tecnologia grants PTDC/MEC-AND/30221/2017 and EXPL/MEC- AND/0676/2021 (PNC); CEECIND/03345/2018 (JDB); and PD/BD/114362/2016 (CSM).
National Institutes of Health, National Library of Medicine grant 1R01LM012832 (LMR).
National Science Foundation Research Traineeship grant 1735095 (LMR).
Deutsche Forschungsgemeinschaft (DFG, German Research Foundation), Clinical Research Unit ‘Male Germ Cells’ CRU326, project number 329621271 (FT).
Singaporean Ministry of Education grant MOE2018-T2-2-053 (IJ, MM).
Declaration of interests
The authors declare no competing interests.
Data and code availability
An interactive online application compiling experimental results on 920 evolutionarily conserved spermatocyte genes and associated microscopy images (Meiotic Navigator) is available at https://pages.igc.pt/meionav. Computational data and custom R and Python scripts used in the analysis are available from https://github.com/rionbr/meionav. All testicular RNA-Seq data generated in this study (Sup. Table 1) were deposited in the European Genome-Phenome Archive (human data) and in the Sequence Read Archive (fruit fly data).
Supplementary materials
Methods
Phylostratigraphic analysis
For humans, mice and fruit flies, published RNA-Seq data of different stages of spermatogenesis and of representative cell types of the three primary embryonic layers (Sup. Table 1) were downloaded from the Sequence Read Archive. Data were checked for quality control and preprocessed by trimming adaptor sequences (Trim galore!; v0.5.0; Babraham Bioinformatics). Gene expression was quantified as transcripts per million (TPM) using Salmon v0.14.11. Transcript-level information was aggregated at the gene level, thus minimizing the potential issue of transcript diversity. Orthogroups (set of homologous genes derived from the last common ancestor) were defined based on the proteomes of 25 species representing key phylogenetic positions (Sup. Fig. 1), and were assembled using OrthoFinder v2.4.02, with DIAMOND v0.9.24.125 as aligner under default settings3. A species tree reflecting the current consensus of the eukaryotic phylogeny4–6 served as the basis for the phylostratigraphic analysis (Sup. Fig. 1). Each orthogroup was assigned to a phylostratum (node) by identifying the oldest clade found in the orthogroup7, using ETE v3.0 with the “get_common_ancestor” option. Phylostrata were assigned a node number, ranging from 1 (oldest) to 16 (youngest, species-specific). Only genes with a minimum average TPM >1 were considered expressed. After mapping all expressed genes to the phylostratum containing their corresponding orthogroup, the distribution of the transcriptome allocated to phylostrata 1 to 5 (orthogroups common to all metazoa) was compared, in a pair-wise manner, between germ cells and soma using the Mann–Whitney U test as implemented in the SciPy python package8. The transcriptome age index (TAI) of the germ cell samples was calculated by dividing the product of each gene’s TPM value and node number by the sum of all TPM values9. Higher TAI values represent younger transcriptomes.
Protein-protein interaction (PPI) network construction
To assemble the spermatocyte PPI networks based on transcriptome data, previously published RNA-Seq datasets of purified human and mouse spermatocytes and of the spermatocyte-enriched median region of the fruit fly testes were used to identify all expressed genes (see Sup. Table 1 for the list of all selected datasets). Salmon v0.14.1 was used to quantify gene expression levels, and genes were considered expressed based on a minimum average expression level of TPM >1. All potential interactions between expressed genes were retrieved from the STRING v11 database10. For this, we used the combined confidence score, which corresponds to the estimated likelihood of a given association being true, given the underlying evidence. To remove low-confidence interactions, only PPIs with a STRING combined confidence score ≥0.5 were included (Sup. Fig. 3a-c). A similar approach was employed to assemble enterocyte PPI networks. A multilayer network representation was employed to define orthologous relationships, with each layer corresponding to a particular species and nodes belonging to the same orthogroup being connected across layers. Orthologous connections were established using the EggNOG v5 database at the metazoan level11 and a one-to-many relationship was employed when identifying orthologs between network layers. A variable percentage of all spermatocyte network edges was randomly shuffled to test for possible ascertainment bias in conserved genes (Sup. Fig. 6).
Backbone and orthoBackbone computation
The initial step to define the male germ cell orthoBackbone was the extraction of the metric backbone of each species’ spermatocyte PPI network (i.e., each layer of the multilayer network). The metric backbone is the subgraph that is sufficient to compute all shortest paths in the network, thus removing edges that break the triangle inequality (and are therefore redundant in regards to the shortest paths). Importantly, the metric backbone retains all metric edges while preserving all network nodes, the natural hierarchy of connections and community structures12,13. Considering that a set of the PPIs in each layer are evolutionarily conserved, we computed the orthoBackbone: a single shared subnetwork obtained by collapsing the metric backbone of all layers according to the orthologous relationships. More precisely, the orthoBackbone corresponds to a subgraph of the metric backbone of every layer, where each edge has an analogous edge connecting the same orthologous genes in all other species’ layers. For cases where there was not a one-to-one conserved edge relationship (i.e., an edge of one layer’s backbone mapped to several edges in other species’ layers due to the existence of paralogs), the inclusion criterion was for at least one of these homologous edges to be part of the backbone (Fig. 2a). Note that because only edges between orthologous genes can be part of the orthoBackbone, nodes (including orthologous nodes) that are left with no edges are removed. Thus, the orthoBackbone typically has fewer nodes than its original network layer and it is not necessarily a metric backbone. Next, annotations retrieved from the Gene Ontology (GO) resource database (release 2020-09-10) were used to select all human orthoBackbone edges connecting at least one spermatogenesis-related gene linked to gene expression regulation. This led to the identification of the 79 functional interactions in spermatocytes that form a core component of the conserved genetic program of spermatogenesis (Fig. 2e). To determine the impact of using three evolutionarily distant species to define the orthoBackbone, we calculated the size of this subgraph when calculated based on just two species. We observed that the orthoBackbone of a human and mouse spermatocyte contains 19,333 edges (5.5% of all functional interactions in the human spermatocyte network), and that of a human and fruit fly spermatocyte contains 9,285 edges (2.6% of all functional interactions).
Differential gene expression analysis of the mitosis-to-meiosis transcriptional burst
Reads were aligned to their respective genomes (GRCh38, GRCm38 and BDGP6.22) using HISAT2 v2.1.0 under default parameters14. Uniquely mapped read counts were generated using FeatureCounts v1.5.0-p1 and ENSEMBL GTF annotations15. Differential gene expression analysis was performed using a likelihood test16, as implemented in the edgeR package17. Genes that were significantly upregulated at meiotic entry (spermatocyte vs. spermatogonia) and/or downregulated at meiotic exit (spermatid vs. spermatocyte) were selected for further analysis as part of the mitosis-to-meiosis transcriptional burst that activates the spermatogenic program. Three criteria were employed to define differentially expressed genes (DEGs): false discovery rate ≤0.05; abs[log2(fold change)] ≥1; and average log2(normalized counts per million) ≥1. Subsequently, EggNOG orthogroups defined at the metazoan level11 were used to establish which genes had a similar up/downregulation behavior in both species. This list was further trimmed by only retaining those whose fruit fly orthologs were also expressed at average log2(normalized counts per million) >1 in the corresponding testes, leading to 920 fruit fly genes (homologous to 797 and 850 in humans and mice, respectively).
Drosophila in vivo RNAi
An in vivo UAS-GAL4 system was used to silence the 920 conserved transcriptional burst genes specifically at the mitosis-to-meiosis transition18. Silencing was induced using the bam-GAL4 driver (kindly provided by Renate Renkawitz-Pohl, Philipps Universität Marburg, Germany), which promotes high levels of shRNA expression in late spermatogonia and early primary spermatocytes19. UAS-hairpin lines targeting the selected genes were purchased from the Bloomington Drosophila Stock Centre (BDSC) and the Vienna Drosophila Resource Centre (VDRC). Lines previously associated with a phenotype in the literature, regardless of the tissue, were preferentially chosen for this experiment. In the ten cases where no lines were available, these were generated in-house following standard procedures20. Briefly, shRNA were designed using DSIR21, with the corresponding sequences being cloned into the pWalium 20 vector. Constructs were injected into fruit flies carrying an attp site on the third chromosome (BDSC stock #24749) and transgene-carrying progeny were selected to establish the UAS-hairpin lines. A similar strategy was employed to generate a second independent RNAi reagent for dRNF113 (dRNF113 #2). The antisense sequences selected for each locus are listed in Sup. Table 3, and information on all tested lines is available at the Meiotic Navigator gene browser.
For assessing fertility, gene-silenced males were mated with wildtype Oregon-R virgin females (2 males:4 females) for 12 hours at 25°C. Laid eggs were left to develop for 24 hours at 25°C before the percentage of egg hatching was determined. This percentage (fertility rate) served as a measure of the male reproductive fitness associated with each tested gene. All fruit flies were 3 to 7 days post-eclosion and fertility rates correspond to an average of four independent experiments, with a minimum of 25 eggs scored per replicate. Every batch of experiments included a negative control [RNAi against the mCherry fluorophore (BDSC stock #35785) - a sequence absent in the fruit fly genome], and a positive control [RNAi against Ribosomal protein L3 (BDSC stock #36596) - an essential unit of the ribosome]. A cut-off of <75% fertility rate was established to define impaired reproductive fitness based on the rate observed in the negative controls (>2 standard deviations of 91.6 ± 8.5%). For dRNF113, Hsc70-1 (VDRC stock #106510) and dKPNA2 (VDRC stock #34265) extended fertility tests were conducted, with ∼100 eggs scored per replicate. To this end, egg-laying cages with apple juice agar plates as substrate were set up with a 1:2 male to female ratio. All Drosophila lines were maintained at 25°C in polypropylene bottles containing enriched medium (cornmeal, molasses, yeast, beet syrup and soy flour).
The novelty of the hits was determined based on the following exclusion criteria (applied to all 250 hits and their orthologs): i- clinically validated human infertility genes22; ii- genes associated with a mouse male infertility phenotype (MP:0001925) in the Mouse Genome Informatics database23; and iii- fruit fly genes with an associated male germ cell or male reproductive phenotype in the FlyBase database24. For the latter two categories, information reflects data available in both systems as of April 2023.
Human and mouse testicular imaging
Testicular biopsies of the azoospermic individuals M1911 (RNF113B variant), M2098 (KPNA2 variant) and M2190 (HSPA2 variant) were collected at the Centre for Reproductive Medicine and Andrology (CeRA) in Münster, Germany or at the University Hospital Giessen, Germany. As a control for testicular imaging, tissue of individual M2951 (who underwent a vasectomy reversal at the CeRA) was used. All testicular tissues were fixed in Bouin’s solution and embedded in paraffin using an automatic ethanol and paraffin row (Bavimed Laborgeräte GmbH). Serial sections (5μm- thick) were stained with periodic acid-Schiff according to standard procedures. Testicular histopathology was performed as part of the routine clinical work up. For quantification, at least 100 testicular tubules per testis were screened for the most advanced germ cell type present (categories: spermatogonia, spermatocytes, round spermatids and elongating spermatids), or in their absence, for Sertoli cells or tubular shadows.
Mouse testes were dissected from 10 week-old animals and fixed in modified Davidson’s fluid for 48h at room temperature before being paraffin-embedded according to standard procedures25. Serial sections (3μm-thick) were stained with haematoxylin-eosin. Slides were analyzed with the NanoZoomer-SQ Digital slide scanner (Hamamatsu Photonics).
Drosophila testicular imaging
Squash preparations of freshly dissected Drosophila testes were performed as previously described26 and examined using a phase contrast microscope (Nikon Eclipse E400). Phenotypes of all silenced genes associated with decreased reproductive fitness (250 in total) were assigned to one of five classes, based on the earliest stage in which the cytological defects were detected: 1- pre-meiotic (evidence of increased cell death and/or clear morphological defects in late spermatogonia); 2- meiotic (failure to progress successfully through meiosis, evidence of increased cell death, and/or clear morphological defects in spermatocytes); 3- post-meiotic (abnormal spermiogenesis characterized by sperm individualization defects and/or significant spermatid morphological defects); 4- gamete (decreased mature sperm numbers and/or motility compared with controls after qualitative assessment); 5- undetectable (no cytological defects). Images from at least 2 pairs of testes were acquired per genotype and were corrected for background illumination as previously described27. These are all available in the Meiotic Navigator gene browser.
For the quantification of the meiotic area, phase-contrast images of different testicular regions were acquired at 40x and stitched together, using the MosaicJ tool28 in the ImageJ software (v1.8, National Institutes of Health), to reconstruct a high-resolution image of the entire testis. On average, 8 individual images were acquired to assemble each complete testis. The meiotic area corresponds to the ratio between the area (in pixels) occupied by primary spermatocytes and that of the entire gonad. A total of 15 different testes were quantified per genotype across 3 independent experiments, and groups were compared using unpaired t-tests.
For the dRNF113 immunofluorescence assay, a modified whole-mount protocol specifically designed for the analysis of Drosophila spermatocytes was employed29. Briefly, testes were dissected in testis buffer (183mM KCl, 47mM NaCl, 10mM Tris-HCl, 1mM EDTA and 1mM PMSF), transferred to a pre-fix solution containing 4% formaldehyde (PolySciences) in PBS, and then fixed for 20min using a heptane-fixative mix at 3:1. The fixative consisted of 4% formaldehyde in a PBS + 0.5% NP-40 (Merck) solution. After washing the fixative, samples were incubated for 1h in PBS supplemented with 0.3% Triton X-100 (Sigma-Aldrich), 1% (w/v) bovine serum albumin (BLIRT) and 1% (w/v) donkey serum (Thomas Scientific). Primary antibody incubation (anti-dRNF113 at 1:250; kindly provided by Chris Doe, University of Oregon, USA) was performed overnight at 4°C in PBS supplemented with 1% (w/v) bovine serum albumin and 1% (w/v) donkey serum. After washing the primary antibody solution, samples were incubated for 1h with a goat anti-guinea pig secondary antibody (Alexa Fluor 488 conjugate; 1:1000; Invitrogen) and wheat germ agglutinin to label the nuclear envelope (Alexa Fluor 647 conjugate; 1:500; Invitrogen). Testes were mounted in Vectashield mounting medium with DAPI (Vector Laboratories) and the entire nuclear volume of individual cells was acquired as 1μm-thick slices using a 63x oil immersion objective and 10x digital zoom in a Leica SP5 confocal microscope. Slices were stacked into maximum intensity Z-projections and the relative intensity of the dRNF133 signal was measured as previously described30. A total of 30 late spermatogonia and 30 mature primary spermatocytes (late prophase I) were quantified across 3 independent experiments, and the two groups were compared using an unpaired t-test.
Generation of genetically modified mice
A mutant allele for the mouse Rnf113a2 gene was generated using the CRISPR/Cas9 system in the FVB/N background. gRNA was assembled by in vitro hybridization of an Alt-R scRNA targeting the sequence ATGATCCAGAAGACGAGTGG with the Alt-R tracrRNA (both from Integrated DNA Technologies, Inc.). An active ribonucleoprotein (RNP) complex was obtained by incubation of the gRNA with the Cas9 protein. The RNP containing 1μM of the hybridized gRNA and 100ng/μl of Cas9 were microinjected into the pronucleus of fertilized oocytes and transferred into pseudopregnant females according to standard procedures31. The resulting pups were genotyped by PCR using primers Rnf113-F and Rnf113-R (Sup. Table 3). PCR products were run on a native 15% polyacrylamide gel in TBE buffer, from which 3 pups were identified with a PCR pattern consisting of a mix between a fragment matching the expected size for the wildtype allele and a smaller fragment. The smaller fragments were recovered, amplified, and sequenced, revealing deletions of 3, 9, or 14 nucleotides. Only the allele containing the 14-nucleotide deletion is expected to disrupt Rnf113a2 by generating a frameshift leading to an immediate stop codon after amino acid 64 of the 338 in the full-length. A mouse line was generated from this founder protein (Rnf113a2em1Mllo, here forth referred to as Rnf113a2KO). Homozygous Rnf113a2KO/KO mice were generated by intercrosses between heterozygous mice. Genotyping was performed by PCR using primers Rnf113-F and Rnf113-Chk-wt-R for the wildtype allele and Rnf113-F and Rnf113- Chk-mut-R for the mutant allele (Sup. Table 3). All animals were housed at the Mouse Facility of Instituto Gulbenkian de Ciência (MF-IGC). Wildtype FVB/N mice were generated in the production area of the MF-IGC. Transgenic animals were produced and kept in the barrier experimental area of the MF-IGC.
Characterization of mouse reproductive parameters
Upon sexual maturation, 8-week-old male F2 Rnf113a2KO/KO homozygous mice and corresponding wildtype littermates (Rnf113a2WT/WT) were caged individually and mated with two 8-week-old wildtype FVB/N females. Five males were tested for each genotype to meet the requirements for statistical validity. After the mating period, females with a confirmed copulatory plug were separated. To promote animal welfare and to reduce the number of animals used, we avoided the normal delivery of pups as a means of assessing fertility. For this, separated females were euthanized by cervical dislocation 10 days post copulation and male fertility was evaluated by the presence of developing embryos. Whenever present, the number of embryos was recorded to quantitate the number of progeny sired by the tested males. All males were considered tested after at least two females had confirmed copulatory plug. To complement the in vivo fertility tests, 10-week-old males were euthanized by cervical dislocation for testicular collection and to examine internal anatomy. Testes were weighted before being processed for histology (see Human and mouse testicular imaging). All experiments followed the Portuguese (Decreto-Lei n° 113/2013) and European (Directive 2010/63/EU) legislations, concerning housing, husbandry, and animal welfare.
RNA-Seq in human and Drosophila testes
Total RNA from two snap-frozen testicular tissue samples of the azoospermic individual M1911 (RNF113B variant) as well as from one sample each from three unrelated individuals with normal spermatogenesis [M1544 (obstructive azoospermia), M2224 (anorgasmia) and M2234 (previous vasectomy)]32, were extracted using the Direct-zol RNA Microprep kit (Zymo Research), following the manufacturer’s instructions. RNA quality was estimated by electrophoresis (Agilent Technologies), with all samples having a RNA integrity number (RIN) >4.5 (range: 4.5 - 5.6), except one replicate of M1911 (with RIN= 2.0). rRNA depletion was performed with the NEBNext rRNA Depletion Kit v2 (New England Biolabs), and total RNA libraries were prepared with the NEBNext Ultra II Directional RNA Library Prep Kit for Illumina (New England Biolabs), using 700ng of RNA per sample. Paired-end sequencing with 150bp per read was performed using the NextSeq2000 system (Illumina) at the University of Münster, Germany. An average of 44 million reads was generated per sample. The expression pattern of RNF113B was determined in a control testis by retrieving the corresponding data from our recently published single-cell RNA- Seq dataset of normal human spermatogenesis33.
For fruit flies, RNA was extracted from 40 pairs of adult testes per sample per condition (3 to 7 days post-eclosion) using the PureLink RNA Mini Kit (Thermo Fisher Scientific) following the manufacturer’s instructions. Three conditions were analyzed: dRNF113 RNAi (meiotic arrest), Prp19 RNAi (meiotic arrest) and the mCherry RNAi negative control (normal spermatogenesis), with three independent biological replicates per condition. Extracted RNA was treated with DNAse (Thermo Fisher Scientific), with a RNA quality number of 10 for all samples, as assessed by electrophoresis (Advanced Analytical Technologies). Total RNA libraries were prepared using the Zymo-Seq RiboFree Total RNA-Seq Library Kit (Zymo Research), using 1000ng of RNA per sample. Paired-end sequencing with 150bp per read was performed using DNBSEQ technology (BGI Group), with an average of 46 million reads per sample. Confirmation of knockdown efficiency via quantitative RT-PCR was performed for both the dRNF113 and Prp19 RNAi. For this, cDNA was synthetized using the Transcriptor First Strand cDNA Synthesis Kit (Roche), according to the manufacturer’s instructions. Reactions were performed in a QuantStudio 6 Real- Time PCR System (Thermo Fisher Scientific), using SYBR Green chemistry (Applied Biosystems). All reactions were set up in triplicates and Act57B and RpL32 were used as internal controls. Primers are listed in Sup. Table 3.
RNA-Seq reads were aligned against the human and fruit fly genomes (GRCh38 and BDGP6.32, respectively) using the STAR aligner34. Gene-level counts were obtained using feature-counts, taking into account the strandedness. Counts were normalized with the TMM method35, and differential gene expression analysis was performed using a quasi-likelihood F-test36, as implemented in the edgeR package. DEGs were defined using the previously listed criteria (see Differential gene expression analysis of the meiotic transitions). To determine to what extent the human patient sample with lower RIN could be affecting our results, we repeated this analysis correcting for sample degradation (using DegNorm37) This revealed that 97.8% of all DEGs, as defined by our criteria, were also listed as such after taking into account sample degradation. This fits with a largely similar number of expressed genes (at TPM >1) detected in the control and in the experimental groups (average of 53.9% and 46.8% expressed protein-coding genes, respectively).
Whole exome sequencing and orthoBackbone analysis
Whole exome sequencing (WES) in the 1,021 azoospermic men included in the MERGE study cohort was performed as previously described38. All men provided written informed consent, in agreement with local requirements. The study protocol was approved by the Münster Ethics Committees/Institutional Review Boards (Ref. No. Münster: 2010-578-f-S) in accordance with the Helsinki Declaration of 1975.
Exome sequencing data were filtered for high-confidence loss of function (LoF) variants in the 27 orthoBackbone genes that were similarly downregulated in the testicular transcriptome of M1911 and of the dRNF113 RNAi, and for which there were available functional data supporting a role in spermatogenesis. These genes were: BAZ1A, BOLL, CDC14A, CFAP61, CFAP91, CNOT6, CNOT7, CUL3, DHX36, DNAH7, DNAH8, DNAI1, DNAI2, DRC7, EIF2A, HSPA2, KPNA2, LIPE, MTMR7, PPP1R2, PSMF1, SKP1, SRPK1, TOB1, TPGS2, TUBA4A and UBE2K. From this list, all those previously associated with known causes of human infertility or with congenital syndromic conditions were excluded (CFAP61, CFAP91, DNAH7, DNAH8, DNAI1, DNAI2 and LIPE). The inheritance mode for each candidate gene was predicted using three sources: the DOMINO algorithm39, the gnomAD observed/expected (o/e) constraint score for LoF variants40, and the recent dosage sensitivity map of the human genome41. Filtering criteria were: stop-gain, frameshift and splice site variants with a minor allele frequency (MAF) <0.01 in gnomAD, a maximum occurrence of 10x in our in-house database and a read depth >10. Variants were only considered when detected in accordance with the predicted mode of inheritance.
Supplementary Text, Figures and Tables
Clinical genetics
The homozygous frameshift variant c.556_565del;p.(Thr186GlyfsTer119) in RNF113B in individual M1911 was previously identified by filtering whole exome sequencing data of 74 men with complete bilateral meiotic arrest for rare LoF variants (see Methods for criteria). Two first degree cousins of M1911 were also infertile, but genetic analyses were not possible. It is likely that both are also affected by the homozygous RNF113B frameshift variant, as the parents of both men are consanguineous (therefore, family ascendants of M1911). The heterozygous LoF variants in the HSPA2 and KPNA2 genes (2 cases for each, for a total of four unrelated men) were identified by analyzing data from 1,021 azoospermic men from the MERGE cohort. For HSPA2, the heterozygous stop-gain variant c.175C>T;p.(Gln59Ter) was identified in individual M1678, and the heterozygous frameshift variant c.1806dup;p.(Glu603ArgfsTer81) in M2190. For KPNA2, the heterozygous splice acceptor variant c.667-2A>G;p.? was detected in individual M1645, and the heterozygous frameshift variant c.1158_1161del;p.(Ser387ArgfsTer14) in M2098.
The inheritance for HSPA2 and KPNA2 is predicted to be “very likely autosomal dominant” by DOMINO. The gnomAD o/e scores are 0.32 (0.17 - 0.67) for HSPA2 and 0.21 (0.11 - 0.44) for KPNA2. The pHaplo scores are 0.350 for HSPA2 and 0.747 for KPNA2. The proposed “hard”/conservative thresholds for the “loss of function observed/expected upper bound fraction” or “LOEUF” score (the upper limit of the confidence interval) is <0.35 and for pHaplo >0.86. Thus, both LOEUF and pHaplo hint towards likely autosomal dominant inheritance, more so for KPNA2 than for HSPA2, but the formal cut-offs are not met in both cases.
We tested the possibility that the actual cause of infertility in M1911 and his infertile brother was due to another homozygous variant located in the vicinity of RNF113B. To this end, we queried all genes within a range of 1 million base pairs in both directions from RNF113B. The following
12 genes were analyzed: FARP1, IPO5, RAP2A, MBNL2, MIR3170, STK24, STK-24AS1, SLC15A1, DOCK9-AS1, DOCK9, LINC00456 and DOCK9-DT. Of these, only IPO5 is robustly expressed in the testis (data: GTEx Project, NIH). Variants in this gene have been associated with keratoconus, but not with male infertility. Additionally, no coding variants affecting any of these 12 genes (minor allele frequency <1%) were detected in M1911. These results render extremely unlikely the possibility that M1911’s infertility phenotype is caused by a co-segregating variant.
To rule out possible alternative monogenic causes for the infertility phenotype in all 5 selected individuals, the sequences of 230 genes previously identified with an at least limited level of evidence of being associated with male infertility1, were screened in the exomes of M1911 (RNF113B variant), M1678 (HSPA2 variant), M2190 (HSPA2 variant), M1645 (KPNA2 variant) and M2098 (KPNA2 variant; see Methods for detailed criteria). Of these variants, only those affecting genes with a quantitative impact on spermatogenesis, that were consistent with the reported mode of inheritance of the respective gene, and occurring <10x in our in-house database, were considered for further analysis. These criteria identified the heterozygous missense variant c.2377G>T;p.(Ala793Ser) in DNMT1 as a possible alternative cause for male infertility in M1911, a possibility that was disproven by subsequent familial studies. Briefly, DNA from M1911’s mother, a fertile brother and an azoospermic brother were available for segregation analysis, which was performed by Sanger sequencing. The azoospermic brother of M1911 was also homozygous for the frameshift variant in RNF113B, while the mother and the fertile brother were heterozygous for the variant. Of note, the missense variant c.2377G>T p.(Ala793Ser) in DNMT1 was also identified in M1911’s fertile brother in a similarly heterozygous state (primers are listed in Sup. Table 3). The latter result strongly suggests that the DNMT1 variant has no pathogenic effect on male fertility. Segregation analysis was also performed in M1654’s parents, revealing that the splice-site variant c.667-2A>G;p.? was not present either in his father or in his mother (i.e., indicating that it was a de novo variant). Paternity was confirmed by short-tandem- repeat analysis. No chromosomal aberrations or Y-chromosomal microdeletions in the AZF regions were identified in any of the selected individuals.
Clinical data of individuals M1678 and M2190 (HSPA2 variants)
The testicular phenotype of M1678 is unknown, as this individual did not undergo a testicular biopsy. However, FSH levels of 26.7 U/L (normal range 1-7 U/L) and a bi-testicular volume of 10 mL (normal value >24 mL) are clearly indicative of non-obstructive azoospermia. M2190 underwent a testicular biopsy at the Urology and Andrology department of the University Hospital of Giessen, Germany. He was diagnosed with complete bilateral SCO and sperm retrieval was unsuccessful.
Clinical data of individuals M1645 and M2098 (KPNA2 variants)
Individual M1645 underwent a testicular biopsy in a different clinical unit, with the corresponding pathology report stating complete bilateral SCO without successful sperm retrieval. Histology data was not available for this man. M2098 was also diagnosed with complete bilateral SCO and sperm retrieval was equally unsuccessful.
Supplementary Figures
References
- 1.Construction of a human cell landscape at single-cell levelNature 581:303–309
- 2.The Strategy of the Genes. (RoutledgeLondon https://doi.org/10.4324/9781315765471
- 3.Evolutionary Dynamics of Gene and Isoform Regulation in Mammalian TissuesScience 338:1593–1599
- 4.Panorama of ancient metazoan macromolecular complexesNature 525:339–344
- 5.The evolution of cell types in animals: emerging principles from molecular studiesNat. Rev. Genet 9:868–882
- 6.The origin of Metazoa: a unicellular perspectiveNat. Rev. Genet 18:498–512
- 7.Evolution of the bilaterian germ line: lineage origin and modulation of specification mechanismsIntegr. Comp. Biol 47:770–785
- 8.The rapid evolution of reproductive proteinsNat. Rev. Genet 3:137–144
- 9.Sperm competition and the evolution of spermatogenesisMol. Hum. Reprod 20:1169–1179
- 10.SpermTree, a species-level database of sperm morphology spanning the animal tree of lifeSci. Data 9
- 11.Gene expression across mammalian organ developmentNature 571:505–509
- 12.Cellular Source and Mechanisms of High Transcriptome Complexity in the Mammalian TestisCell Rep 3:2179–2190
- 13.Origins, evolution, and phenotypic impact of new genesGenome Res 20:1313–1326
- 14.New genes often acquire male-specific functions but rarely become essential in DrosophilaGenes Dev 31:1841–1846
- 15.Relaxed Selection and the Rapid Evolution of Reproductive GenesTrends Genet 36:640–649
- 16.Comparative transcriptomic analysis reveals conserved programmes underpinning organogenesis and reproduction in land plantsNat. Plants 7:1143–1159
- 17.The Conservation of the Germline Multipotency Program, from Sponges to Vertebrates: A Stepping Stone to Understanding the Somatic and Germline OriginsGenome Biol. Evol 9:474–488
- 18.Single-Cell RNA Sequencing of the Cynomolgus Macaque Testis Reveals Conserved Transcriptional Profiles during Mammalian SpermatogenesisDev. Cell 54:548–566
- 19.Single-Cell RNA Sequencing of Human, Macaque, and Mouse Testes Uncovers Conserved and Divergent Features of Mammalian SpermatogenesisDev. Cell 54:529–547
- 20.The molecular evolution of spermatogenesis across mammalsNature 613:308–316
- 21.Sex and suicide: The curious case of Toll-like receptorsPLOS Biol 18
- 22.Widespread Transcriptional Scanning in the Testis Modulates Gene Evolution RatesCell 180:248–262
- 23.Higher Germline Mutagenesis of Genes with Stronger Testis Expressions Refutes the Transcriptional Scanning HypothesisMol. Biol. Evol 37:3225–3231
- 24.A subcellular map of the human proteomeScience 356
- 25.A phylostratigraphy approach to uncover the genomic history of major adaptations in metazoan lineagesTrends Genet 23:533–539
- 26.A phylogenetically based transcriptome age index mirrors ontogenetic divergence patternsNature 468:815–818
- 27.Transcriptional regulation during Drosophila spermatogenesisSpermatogenesis 2:158–166
- 28.Using a meiosis detection toolkit to investigate ancient asexual “scandals” and the evolution of sexBioEssays 30:579–589
- 29.Super-enhancer switching drives a burst in gene expression at the mitosis-to-meiosis transitionNat. Struct. Mol. Biol 27:978–988
- 30.Transcriptome and translatome co-evolution in mammalsNature 588:642–647
- 31.Unraveling transcriptome dynamics in human spermatogenesisDevelopment dev 152413https://doi.org/10.1242/dev.152413
- 32.STRING v11: protein–protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasetsNucleic Acids Res 47:D607–D613
- 33.Assessment of network module identification across complex diseasesNat. Methods 16:843–852
- 34.Evolutionary conservation of motif constituents in the yeast protein interaction networkNat. Genet 35:176–179
- 35.Unequal evolutionary conservation of human protein interactions in interologous networksGenome Biol 8
- 36.The effective graph reveals redundancy, canalization, and control pathways in biochemical regulation and signalingProc. Natl. Acad. Sci 118
- 37.The distance backbone of complex networksJ. Complex Netw 9
- 38.Contact networks have small metric backbones that maintain community structure and are primary transmission subgraphsPLOS Comput. Biol 19
- 39.Shortest path counting in probabilistic biological networksBMC Bioinformatics 19
- 40.Protein Interaction Mapping in C. elegans Using Proteins Involved in Vulval DevelopmentScience 287:116–122
- 41.Evolution of biological interaction networks: from models to real dataGenome Biol 12
- 42.Transcription Factor RFX2 Is a Key Regulator of Mouse SpermiogenesisSci. Rep 6
- 43.Chromodomain Protein CDYL Acts as a Crotonyl-CoA Hydratase to Regulate Histone Crotonylation and SpermatogenesisMol. Cell 67:853–866
- 44.Widespread Presence of Human BOULE Homologs among Animals and Conservation of Their Ancient Reproductive FunctionPLoS Genet 6
- 45.Identifying gene expression programs of cell-type identity and cellular activity with single-cell RNA-SeqeLife 8
- 46.The adult human testis transcriptional cell atlasCell Res 28:1141–1157
- 47.Attenuated chromatin compartmentalization in meiosis and its maturation in sperm developmentNat. Struct. Mol. Biol 26:175–184
- 48.Tissue, cell type and stage-specific ectopic gene expression and RNAi induction in the Drosophila testisSpermatogenesis 2:11–22
- 49.Q. midlife crisis encodes a conserved zinc-finger protein required to maintain neuronal differentiation in DrosophilaDevelopment 140:4155–4164
- 50.Primate and Rodent Specific Intron Gains and the Origin of Retrogenes with Splice VariantsMol. Biol. Evol 28:33–37
- 51.Single-cell RNA-seq unravels alterations of the human spermatogonial stem cell compartment in patients with impaired spermatogenesisCell Rep. Med 2
- 52.Human RNF113A participates of pre-mRNA splicing in vitroJ. Cell. Biochem 120:8764–8774
- 53.A systematic review of the validated monogenic causes of human male infertility: 2020 update and a discussion of emerging gene–disease relationshipsHum. Reprod. Update 28:15–29
- 54.Genetic studies in intellectual disability and related disordersNat. Rev. Genet 17:9–18
- 55.A unique view on male infertility around the globeReprod. Biol. Endocrinol 13
- 56.Bi-allelic Mutations in M1AP Are a Frequent Cause of Meiotic Arrest and Severely Impaired Spermatogenesis Leading to Male InfertilityAm. J. Hum. Genet 107:342–351
- 57.Role of heat shock protein HSP70-2 in spermatogenesisRev. Reprod 4:23–30
- 58.Heat shock cognate 70 genes contribute to Drosophila spermatocyte growth progression possibly through the insulin signaling pathwayDev. Growth Differ 63:231–248
- 59.Targeted gene disruption of Hsp70-2 results in failed meiosis, germ cell apoptosis, and male infertilityDev. Biol 93:3264–3268
- 60.Importin Alpha2-Interacting Proteins with Nuclear Roles During Mammalian SpermatogenesisBiol. Reprod 85:1191–1202
- 61.The Meiosis-Specific Hop2 Protein of S. cerevisiae Ensures Synapsis between Homologous ChromosomesCell 94:375–386
- 62.The Hop2 Protein Has a Direct Role in Promoting Interhomolog Interactions during Mouse MeiosisDev. Cell 5:927–936
- 63.Achiasmy: Male Fruit Flies Are Not Ready to MixFront. Cell Dev. Biol 4
- 64.Drosophila melanogaster Importin ɑ1 and ɑ3 Can Replace Importin ɑ2 During Spermatogenesis but Not OogenesisGenetics 161:157–170
- 65.A de novo paradigm for male infertilityNat. Commun 13
- 66.Genetic Architecture of Azoospermia—Time to Advance the Standard of CareEur. Urol 83:452–462
- 67.Molecular mechanisms of sperm motility are conserved in an early-branching metazoanProc. Natl. Acad. Sci 118
- 68.Discovery of archaeal fusexins homologous to eukaryotic HAP2/GCS1 gamete fusion proteinsNat. Commun 13
- 69.A conserved node in the regulation of Vasa between an induced and an inherited program of primordial germ cell specificationDev. Biol 482:28–33
- 70.Preformation and epigenesis converge to specify primordial germ cell fate in the early Drosophila embryoPLOS Genet 18
- 71.Is it possible to overcome issues of external validity in preclinical animal research? Why most animal models are bound to failJ. Transl. Med 16
- 72.The influence of evolutionary history on human health and diseaseNat. Rev. Genet 22:269–283
- 73.Protein complex prediction with AlphaFold-MultimerbioRxiv https://doi.org/10.1101/2021.10.04.463034
- 74.Structural Basis of Egg Coat-Sperm Recognition at FertilizationCell 169:1315–1326
- 1.Salmon provides fast and bias-aware quantification of transcript expressionNat. Methods 14:417–419
- 2.OrthoFinder: phylogenetic orthology inference for comparative genomicsGenome Biol 20
- 3.Fast and sensitive protein alignment using DIAMONDNat. Methods 12:59–60
- 4.Broad phylogenomic sampling improves resolution of the animal tree of lifeNature 452:745–749
- 5.Premetazoan genome evolution and the regulation of cell differentiation in the choanoflagellate Salpingoeca rosettaGenome Biol 14
- 6.Revisiting metazoan phylogeny with genomic sampling of all phylaProc. R. Soc. B Biol. Sci 286
- 7.A phylostratigraphy approach to uncover the genomic history of major adaptations in metazoan lineagesTrends Genet 23:533–539
- 8.SciPy 1.0: fundamental algorithms for scientific computing in PythonNat. Methods 17:261–272
- 9.A phylogenetically based transcriptome age index mirrors ontogenetic divergence patternsNature 468:815–818
- 10.STRING v11: protein–protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasetsNucleic Acids Res 47:D607–D613
- 11.eggNOG 5.0: a hierarchical, functionally and phylogenetically annotated orthology resource based on 5090 organisms and 2502 virusesNucleic Acids Res 47:D309–D314
- 12.Contact networks have small metric backbones that maintain community structure and are primary transmission subgraphsPLOS Comput. Biol 19
- 13.The distance backbone of complex networksJ. Complex Netw 9
- 14.Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotypeNat. Biotechnol 37:907–915
- 15.featureCounts: an efficient general purpose program for assigning sequence reads to genomic featuresBioinforma. Oxf. Engl 30:923–930
- 16.Differential expression analysis of multifactor RNA- Seq experiments with respect to biological variationNucleic Acids Res 40:4288–4297
- 17.edgeR: a Bioconductor package for differential expression analysis of digital gene expression dataBioinforma. Oxf. Engl 26:139–140
- 18.Targeted gene expression as a means of altering cell fates and generating dominant phenotypesDev. Camb. Engl 118:401–415
- 19.Tissue, cell type and stage-specific ectopic gene expression and RNAi induction in the Drosophila testisSpermatogenesis 2:11–22
- 20.The Transgenic RNAi Project at Harvard Medical School: Resources and ValidationGenetics 201:843–852
- 21.An accurate and interpretable model for siRNA efficacy predictionBMC Bioinformatics 7
- 22.A systematic review of the validated monogenic causes of human male infertility: 2020 update and a discussion of emerging gene–disease relationshipsHum. Reprod. Update 28:15–29
- 23.The mouse Gene Expression Database (GXD): 2021 updateNucleic Acids Res 49:D924–D931
- 24.FlyBase: updates to the Drosophila melanogaster knowledge baseNucleic Acids Res 49:D899–D907
- 25.Assessment of Spermatogenesis Through Staging of Seminiferous TubulesSpermatogenesis: Methods and Protocols Humana Press :299–307https://doi.org/10.1007/978-1-62703-038-0_27
- 26.Preparation of live testis squashes in DrosophilaCold Spring Harb. Protoc 2011
- 27.Background illumination correction – Novel context-based segmentation algorithms for intelligent microscopy
- 28.User-friendly semiautomated assembly of accurate image mosaics in microscopyMicrosc. Res. Tech 70:135–146
- 29.Immunostaining of Whole-Mount Drosophila Testes for 3D Confocal Analysis of Large SpermatocytesJ. Vis. Exp. JoVE https://doi.org/10.3791/61061
- 30.The Trithorax group protein dMLL3/4 instructs the assembly of the zygotic genome at fertilizationEMBO Rep 19
- 31.Manipulating the Mouse Embryo: A Laboratory ManualCold Spring Harbor Laboratory Press
- 32.Transcriptome analyses in infertile men reveal germ cell–specific expression and splicing patternsLife Sci. Alliance 6
- 33.Single-cell RNA-seq unravels alterations of the human spermatogonial stem cell compartment in patients with impaired spermatogenesisCell Rep. Med 2
- 34.STAR: ultrafast universal RNA-seq alignerBioinforma. Oxf. Engl 29:15–21
- 35.A scaling normalization method for differential expression analysis of RNA-seq dataGenome Biol 11
- 36.It’s DE-licious: A Recipe for Differential Expression Analyses of RNA-seq Experiments Using Quasi-Likelihood Methods in edgeRMethods Mol. Biol. Clifton NJ 1418:391–416
- 37.DegNorm: normalization of generalized transcript degradation improves accuracy in RNA-seq analysisGenome Biol 20
- 38.Bi-allelic Mutations in M1AP Are a Frequent Cause of Meiotic Arrest and Severely Impaired Spermatogenesis Leading to Male InfertilityAm. J. Hum. Genet 107:342–351
- 39.DOMINO: Using Machine Learning to Predict Genes Associated with Dominant DisordersAm. J. Hum. Genet 101:623–629
- 40.The mutational constraint spectrum quantified from variation in 141,456 humansNature 581:434–443
- 41.A cross-disorder dosage sensitivity map of the human genomeCell 185:3041–3055
- 1.A systematic review of the validated monogenic causes of human male infertility: 2020 update and a discussion of emerging gene–disease relationshipsHum. Reprod. Update 28:15–29
- 2.The Mammalian Spermatogenesis Single-Cell Transcriptome, from Spermatogonial Stem Cells to SpermatidsCell Rep 25:1650–1667
Article and author information
Author information
Version history
- Preprint posted:
- Sent for peer review:
- Reviewed Preprint version 1:
- Reviewed Preprint version 2:
- Version of Record published:
Copyright
© 2024, Correia et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics
- views
- 2,445
- downloads
- 139
- citations
- 0
Views, downloads and citations are aggregated across all versions of this paper published by eLife.