Article https://doi.org/10.1038/s41467-022-35661-z Diverse monogenic subforms of human spermatogenic failure Received: 8 March 2022 A list of authors and their affiliations appears at the end of the paper Accepted: 16 December 2022 Non-obstructive azoospermia (NOA) is the most severe form of male infertility and typically incurable. Defining the genetic basis of NOA has proven chal- Check for updates lenging, and the most advanced classification of NOA subforms is not based on genetics, but simple description of testis histology. In this study, we exome- 1234567890():,; 1234567890():,; sequenced over 1000 clinically diagnosed NOA cases and identified a plausible recessive Mendelian cause in 20%. We find further support for 21 genes in a 2-stage burden test with 2072 cases and 11,587 fertile controls. The disrupted genes are primarily on the autosomes, enriched for undescribed human “knockouts”, and, for the most part, have yet to be linked to a Mendelian trait. Integration with single-cell RNA sequencing data shows that azoospermia genes can be grouped into molecular subforms with synchronized expression patterns, and analogs of these subforms exist in mice. This analysis framework identifies groups of genes with known roles in spermatogenesis but also reveals unrecognized subforms, such as a set of genes expressed across mitotic divisions of differentiating spermatogonia. Our findings highlight NOA as an understudied Mendelian disorder and provide a conceptual structure for organizing the complex genetics of male infertility, which may provide a rational basis for disease classification. Non-obstructive azoospermia, or lack of sperm in the ejaculate due considering that the number of genes linked to the trait in mice is to disruption of spermatogenesis, is a multifactorial trait with a over 5009, half of which are specifically implicated in azoospermia10. prevalence of 0.4–2% in the male population1–3 and has an incidence The few genomewide studies published have shown promise of using of around 10% in cohorts of infertile men4. Mounting evidence sug- whole-exome and whole-genome sequencing to identify the missing gests that male infertility is broadly associated with late-onset genetic causes of azoospermia11–15. What is needed now is a con- somatic comorbidities, including cancer, cardiovascular disease, ceptual framework to organize what appears to be a very large col- and other chronic diseases5. The mechanisms by which these risks are lection of monogenic disorders that we collectively call “male linked are generally unknown. From the perspective of patient infertility”. We hypothesized that distinct subforms of NOA exist, counseling, disease prevention and management, it is thus becoming currently invisible by standard histological classification, but increasingly relevant to establish the genetic origins of male inferti- detectable at a molecular level. lity and spermatogenic impairment. During pioneering genetic In this work, we test this hypothesis through analysis of the investigations over 40 years ago, large Y chromosome deletions largest collection of azoospermia cases to date, combined with a were identified as a highly penetrant genetic cause of human new single-cell RNA-sequencing (scRNA-seq) atlas of human testis azoospermia6. The advent of next-generation sequencing has gene expression (Fig. 1). The results presented provide a broad renewed the hunt for Mendelian forms of azoospermia and the and detailed survey of rare genetic changes causing azoospermia number of confirmed monogenic causes has since grown to in humans, and a framework for organizing these changes into encompass at least 40 loci7,8. Nevertheless, our understanding of subgroups, for improved characterization and personalized genetic causes of male infertility phenotypes is still lacking management of infertility patients. e-mail: conradon@ohsu.edu Nature Communications | (2022)13:7953 1 Article https://doi.org/10.1038/s41467-022-35661-z COMBINED GENETIC ANALYSIS GEMINI Genetics of Male Infertility Initiative INTEGRATED SINGLE CELL ANALYSIS Two independent male infertility cohorts Whole exome sequencing scRNAseq atlas of adult human testes (exomes of 1148 cases) (924 NOA cases) (12 donors) tSNE1 Clustering of Recurrence of genetic causes molecular NOA causes tSNE2 MERGE GEMINI genes genes Histology Molecular ‘subforms’ MA Newcastle/ SCO Radboud genes HSG Fig. 1 | Overview of the study. DNA samples from NOA cases recruited in 11 centers analysis with scRNAseq atlas of adult human testis enabled to cluster the molecular across the world underwent WES to detect rare recessive variants underlying the NOA causes into ‘subforms’ that go beyond the visual classification of the disease disease. The identified Mendelian causes of NOA were additionally screened for in based on histological findings. Map image created by Serhii Brovko and used under independent male infertility cohorts to map recurrent causes of NOA. An integrated license from iStock. Variation prioritization identifies potential monogenic causes in as expected due to the lower accuracy of detecting structural rear- 20% of NOA cases rangements from exome data17. Through the GEnetics of Male INfertility Initiative (GEMINI), we per- Using PSAP prioritization, we identified a potential molecular formed whole-exome sequencing (WES) on 1,011 unrelated men diag- cause of the disease for 19.3% (n = 178/924) of the NOA cases, nosed with spermatogenic failure, the vast majority with unexplained when considering the SNVs and INDELs only. The PSAP pipeline is NOA. Carriers of known genetic causes of azoospermia or large a filtering tool; some variants that are completely unrelated to structural variants were excluded using the processed sequencing infertility may pass this filter. In order to estimate what fraction of dataset (Supplementary Fig. 1, Supplementary Discussion). The final prioritized variants are false positives, we also applied PSAP to an GEMINI cohort analyzed here consisted of 924 unrelated cases institutional control cohort of 2665 individuals with non- recruited in eleven centers across seven countries (Fig. 1, Supple- reproductive disorders (Supplementary Table 2). Only 6% of mentary Figs. 1b, 2, Supplementary Discussion, Supplementary controls carried prioritized variants, indicating that the discovery Table 1). Testis biopsies were performed on 42% of the subjects. These rate in GEMINI is much higher than background (P = 1.7 × 10−27), biopsies were used to classify patients into three distinct histological but also indicating that as many as 1/3 of prioritized variants in subtypes- Sertoli Cell Only (SCO), complete lack of germ cells and the GEMINI may be false positives. most prevalent diagnosis (n = 248 men), followed by maturation arrest The NOA cohort included 72 unrelated consanguineous cases that (MA, arrest of spermatogenesis, n = 101) and hypospermatogenesis were expected to carry damaging recessive variation leading to the (HSG, n = 37), where the spermatogenic output is severely reduced due disease8. Consistently, a possible genetic cause of NOA was identified to significant variation in inter-tubule content, with or without com- in 76.4% of the consanguineous men, significantly higher than the plete spermatogenesis (Supplementary Fig. 2). remainder of the cohort 14.9%; Fisher’s exact test P = 8.3 × 10−28; We prioritized rare (minor allele frequency <1% in gnomAD) likely (Supplementary Discussion). damaging lesions compatible with recessive Mendelian disease models Additionally, two unrelated NOA cases from Utah were by using an analysis pipeline based on Population Sampling Probability identified with uniparental isodisomy of chromosome 2 (UPD2) (PSAP, “Methods”, Fig. 2a, Supplementary Fig. 3, Supplementary and chromosome 4 (UPD4), and with prioritized variation in Data 1)16. The prioritized variants were further filtered by taking into INHBB18 (regulation of FSH production) and POLN (involved in account the mutation type and the testicular expression patterns of recombination) on affected chromosomes, respectively (Supple- the underlying genes (Methods) yielding a total of 312 variants most mentary Discussion; Supplementary Fig. 4a). The overall detec- likely to have an impact on spermatogenesis. Experimental validation tion rate of a candidate NOA lesion with all variation types confirmed 98% (n = 82/84) of SNVs and 93% (n = 38/41) of INDELs combined reached 20%, an estimate similar to previously repor- tested, whereas a decreased rate was detected for CNVs (43%, n = 3/7) ted rates in genomic studies of male infertility14,15. Nature Communications | (2022)13:7953 2 Article https://doi.org/10.1038/s41467-022-35661-z a b c Genotype 1.00 150 Novel KO KO Non-KO 3 CNV Proportion of cases 0.75 Number of cases INDEL SNV Variation count 100 2 0.50 1 50 0.25 0 DCAF12L1 FAM104A C1orf146 NOTCH1 PCDHB3 SPATA20 CCDC89 CCDC33 PNLDC1 CELSR3 FBXO15 SH3TC2 PLA2G6 SLC4A2 NFKBIB ZNF610 CSMD1 PACRG TDRD9 TSSK4 U2AF1 TTLL9 HCN4 M1AP PLD6 0 0.00 AR Chet XL YL A X Y Autosomes X d e f GEMINI-295 M P C 7 normalized read depth 5 6 STRA8 0 Z-score of 5 Observed −log10P -5 4 ~5kb 3 -10 2 -15 Deletion 1 E1 E2 E3E4 E5 E6 E7 E8 E9 0 0 1 2 3 4 5 PCR_F PCR_R Expected −log10P NLS NES g h Control GEMINI-295 2.5 deleterious genotypes observed in 11,587 fertile controls −log10(p value) from burden test 2.0 1.5 * * 1.0 DCAF12L1 FBXO15 MROH9 SPANXN1 TGIF2LY YY2 C1orf146 CCDC89 DCAF4 FNDC7 HCN4 INHBB NFKBIB PNLDC1 POTEJ * RNF212 SLC4A2 SPIDR SPOCD1 * TDRD9 TMEM63C BAX BRCA2 CLCN2 CSMD1 MAGEB18 PCDHB3 PMP22 RNF17 TDRD12 * NR0B1 * AKAP4 * NOTCH1 USP9Y Fig. 2 | Variant prioritization and burden testing in NOA cases. a Distribution of ~5 kb product in GEMINI-295 (P) validating the homozygous STRA8 deletion. The the prioritized variation types across inheritance modes. AR autosomal recessive, experiment repeated twice with the same result. C, control; M, ladder. g In control Chet Compound heterozygous, XL X-linked, YL Y-linked. b Predicted hemizygous testis, clear staining of STRA8 was observed in spermatogonia (arrows), together loss-of-function genotypes appear to be enriched on sex chromosomes compared with background staining in peritubular and interstitial cells (asterisk). In GEMINI- to biallelic loss-of-function genotypes on autosomes (A) (Fisher’s exact p = 0.018 for 295, STRA8 protein staining in spermatogonia is much lower (arrows), presumably chrX and p = 0.032 for chrY vs A). *p < 0.05. c Summary of all genes with multiple due to nonsense-mediated decay of the transcript. Three different sections of the case findings in the GEMINI cohort. d PSAP p values of all variants detected in controls (n = 8) and the case (n = 1) showed consistent staining patterns. The bar patient GEMINI-295, the carrier of the biallelic STRA8 deletion. STRA8 CNV (the represents 10 µm. h Burden test results from comparison of a combined cohort of smallest p value) was prioritized as the most likely cause of NOA in GEMINI-295. The 2072 NOA cases with 11,587 fertile controls. Thirty-four genes with prioritized gray area represents the 95% confidence interval. e Z-scores of normalized read variation in the GEMINI cohort were selected for burden testing; of these 21 were depth of exome sequencing data around the STRA8 locus, plotted against the nominally associated with NOA (“Methods”). * indicates genes with at least “limited” STRA8 gene model. NLS nuclear localization; NES, nuclear export signal. Green clinical validity defined by ref. 7. arrows, PCR primers. f PCR spanning the predicted deletion region yielded a short Population-based statistical tests associate 21 genes with not a replacement for association testing, which remains a gold stan- azoospermia dard for definitive identification of disease genes. To perform asso- The PSAP analysis described above is useful as a prioritization filter ciation testing, we assembled two independent cohorts of infertile that is sensitive to monogenic disease mutations in n=1 cases, but it is men from MERGE (n = 817 cases from Germany) and Newcastle/ Nature Communications | (2022)13:7953 3 Article https://doi.org/10.1038/s41467-022-35661-z Radboud studies (total n = 331; n = 286 from Netherlands, n = 45 UK), as and plays a vital role in meiotic initiation upon induction by retinoic well as a large control cohort of 11,587 fertile parents (Fig. 1, “Meth- acid in both female and male germline24. The NOA patient with the ods”). Of the 221 putative disease genes prioritized in the GEMINI STRA8 deletion displayed a maturation arrest phenotype where pre- cohort, 18.9% were disrupted in at least one other NOA cohort (Sup- pachytene spermatocytes were rarely observed and germ cells beyond plementary Table 3). Three genes were affected in all NOA cohorts: the pre-pachytene stage were not detected, largely consistent with the M1AP, recently published in a cross-center report as a novel cause of Stra8−/− mouse model (Fig. 2g; Supplementary Fig. 6)24. These data autosomal recessive meiotic arrest in men, and the uncharacterized indicate that undescribed human KOs contribute to the loss of ferti- genes KCTD19 and YY2 (Supplementary Table 3)19,20. lity in men. Using the combined NOA cohorts (n = 2072 cases) and 11,587 fertile controls, we performed burden testing on 34 genes with Clinical validity assessment and the role of recurrent gene prioritized mutations in at least 2 cases (Fig. 2h). Due to data use mutations in establishing pathogenicity limitations we were unable to obtain data on compound hetero- Although the aim of this study was to explore and compile an zygous changes from controls, thus limiting some of the genes extensive catalog of recessive variation possibly leading to NOA that could be tested (Methods). Burden testing identified 21 genes among unrelated cases, we evaluated the clinical validity of the with nominally significant p-values, providing further evidence variants to identify the genes with the best prospects for clinical that variation in these genes is likely deleterious for human male testing. Potentially diagnostic variants classified as “Pathogenic” infertility (Supplementary Table 4). For all 21 of these genes, we or “Likely Pathogenic” were identified in 29 genes disrupted in 37 identified multiple highly deleterious genotypes exclusively in cases from the GEMINI cohort (Supplementary Data 2) based on cases, and none in 11,587 fertile controls (Fig. 2h). Of the nom- the guidelines from the American College of Medical Genetics inally associated genes, only TDRD9 and RNF212 were identified as and Genomics (ACMG)25. These genes consisted of eight with having at least “limited” clinical validity in a recent review of existing human genetic evidence, whereas three appear to be azoospermia research7. We note that the approach to association implicated in human or mouse male infertility for the first time testing that we use here is considered a “two-stage analysis”, in based on the evidence compiled here: DCAF12L1, DCAF4, and YY2 which genes selected from a discovery stage (GEMINI) are (Supplementary Data 2). brought forward to test in a second stage (combined case cohort Recurrent observations are a key for classifying novel disease and fertile controls). In such a framework, if one combines data variation as ‘(Likely) Pathogenic’ according to the ACMG guidelines. from both stages in the final analysis, the p-values should be We observed only 25/221 (11.3%) of prioritized genes in multiple sub- corrected for the number of genes tested in the first stage21. None jects, with a maximum of three recurrent hits seen for four genes of our genes pass the threshold of exome-wide significance for an (Fig. 2c; Supplementary Data 3). This is in accordance with the biolo- association test when using both stage one and stage two data gical complexity of the testis where 77% of all human protein-coding combined (“Methods”). genes are expressed26. False positives identified by PSAP prioritization should only lightly influence this observation. Depending on assump- Undescribed human “knockouts” in azoospermia tions about how false positives are distributed among genes and cases, We next identified rare human biallelic loss-of-function variants, or the expected rate of recurrence ranges from 8.5%–16.1%. “knock-outs” (KO), enrichment of which has been observed for Men- For the majority of the NOA cases with a detected possible cause delian disease, developmental disorders, and autism22. A high- (77%), a single gene was identified, consistent with a monogenic form confidence complete KO was predicted for 50/221 NOA genes identi- of the disease, whereas the remaining subjects carried recessive var- fied in this study with an enrichment on X and Y chromosomes (p = iation in up to five loci depending on the degree of consanguinity 0.018 and p = 0.032 vs autosomes, respectively; hypergeometric test, (Supplementary Fig. 4b). Based on these observations, and a simple Fig. 2b) that have historically been considered the main source of model of NOA genetic architecture, we estimate the total number of azoospermia defects. On chromosome X, 9/19 (47%) genotypes were monogenic azoospermia genes is around 625 (Fig. 2a; “Methods”). KO, including in genes MAGEA3, MAGEB18 and MAGEC3, the members Incorporating the estimated FDR of PSAP prioritization in the calcu- of the cancer/testis-antigen family modulating reproductive success23. lations, the total number drops 12% to 550. The cohort sizes required On chromosome Y, prioritized LoF mutations were found in USP9Y, for detecting multiple occurrences of at least half of all azoospermia DDX3Y, TGIF2LY, and ZFY, the latter three yet to be linked to male genes are thus an order of magnitude higher among outbred cases infertility in humans. Notably, 38.0% (n = 19) of the KO genes repre- (n~6700) than has been feasible until now (Fig. 3b), while the number sented the first instances of human KO observed across more than of consanguineous cases required for equivalent power is much lower 224,000 individuals in five large human variation datasets, including (n~1300) (Fig. 3b). These results highlight the necessity for large cohort gnomAD and the Human Gene Mutation Database (“Methods”; Sup- studies for identifying NOA disease variation with a diagnostic value as plementary Table 5). For comparison, significantly fewer cases of novel per the ACMG guidelines. predicted KOs were observed among the 2265 institutional controls considering only autosomal genes (n = 11; one-tailed binomial test, P = NOA as an understudied Mendelian disorder 5.4 × 10−7). Three genes were recurrently affected by biallelic or hemi- In spite of the growing number of studies on the genetics of NOA, the zygous loss-of-function variants, including poorly studied X-linked total number of the genes linked to the disease has remained min- DCAF12L1 (Fig. 2c). uscule. Consistently, we observed a low fraction (7.7%) of genes pre- In addition to LoF variants, we prioritized five CNVs (1.8 kb viously reported in association with male infertility in humans duplication and five homozygous/hemizygous deletion events of 0.48 (Supplementary Table 6). The remainder represent novel NOA candi- kb-1.1 Mb) that are expected to abolish the function of the underlying date genes in men, of which at least 29.9% (61/204) were known to genes (Supplementary Fig. 5a; Supplementary Discussion). The rear- cause sub/infertility upon disruption in mice. A protein-protein inter- ranged genes were novel in the context of male infertility, with an action (PPI) network analysis of the prioritized genes identified exception of STRA8 where a homozygous splice variant was recently numerous well defined subgroups of interacting proteins with clear found in an NOA/severe oligozoospermia cohort15. We identified and roles in testis function (Fig. 3c). validated a homozygous ~6 kb deletion in the gene, which was pre- To expand the context to all known Mendelian trait loci, we dicted as the top candidate cause of NOA in the patient (Fig. 2d–f, intersected the prioritized gene list with known Mendelian disease Supplementary Fig. 5b–d). STRA8 is considered a meiotic gatekeeper genes in OMIM, an online database of Mendelian traits. For the fifty- Nature Communications | (2022)13:7953 4 Article https://doi.org/10.1038/s41467-022-35661-z a b d GEMINI Mendelian disease 1.00 1.00 1.00 Consanguinity max expression % genes with GEMINI 0.75 no yes Fraction of recurrent genes 0.50 Fraction of recurrent genes 0.75 0.25 0.75 0.00 0.50 across genes (log10) 0.50 2 Mean expression 0 0.25 0.25 −2 GEMINI 9.4% 0.00 0 250 500 750 1000 102 103 104 105 Other Testis Other Testis Size of azoospermia gene pool Number of samples c Ciliogenesis Ubiquitination Endocrine piRNA biogenesis Overlap with GEMINI Disorders of sex development GEMINI Endocrine gland dysfunction Germ cell loss Anatomical malformation Spermiogenesis defect Fertilization defect DNA recombination and repair Centrosome Fig. 3 | Characteristics of NOA genes in the context of male infertility and the for non-consanguineous cases. c High confidence protein-protein interactions broader Mendelian disease landscape. a Estimation of the total number of between GEMINI genes and known male infertility disease genes from Oud et al. azoospermia genes. To detect 9.4% of genes recurrently in non-consanguineous 20197 (Supplementary Data 5). d Expression patterns of NOA genes are distinct cases with predominantly monogenic NOA, as seen in GEMINI, a gene pool of from that of known recessive Mendelian disease genes (n = 615) in the GTEx atlas of roughly 625 azoospermia genes would be expected. Blue area represents standard 52 human tissues. NOA genes are more likely to have maximum expression in testis deviation around the mean of 1000 iterations. b Projected number of samples (top panel) and greater testis-specific expression (lower panel). Center line, med- required for observing recurrent gene hits in NOA cohorts considering the esti- ian; box limits, upper and lower quartiles; whiskers, 1.5× interquartile range; points, mated pool of 625 azoospermia genes and discovery rates of 76% for con- outliers. sanguineous cases, 20% for GEMINI overall (8% of consanguineous cases) and 15% one (24.9%) NOA genes in OMIM, infertility was the most frequent among genes linked to recessive Mendelian disorders (Fisher’s exact reported disease type (n = 13), followed by developmental disorders test p = 1.0 × 10−40; Fig. 3d). (n = 11) and cancer (n = 7) (Supplementary Data 4). Collectively these findings highlight primary testicular failure as a poorly studied Integrative analysis with testis scRNAseq reveals ‘molecular Mendelian disorder. Cases of idiopathic male infertility can provide subforms’ of NOA knowledge about Mendelian disease that is under-sampled in We next aimed to characterize the testicular cell types and pathways other study populations, likely due to tissue-specific functions of affected by prioritized NOA variants by characterizing scRNAseq data many NOA genes. Indeed, by analysis of expression patterns from twelve human testis samples (“Methods”)28. Testis-expressed across 52 human tissues in the GTEx database27, we find that genes genes were grouped into 70 expression modules (“components”) by with high expression in testis are significantly under-represented soft clustering with Sparse Decomposition of Arrays (SDA), a Bayesian Nature Communications | (2022)13:7953 5 Article https://doi.org/10.1038/s41467-022-35661-z a Genes Components Genes Components Cells Cells Sparse gene loadings matrix Gene expression data Cell scores matrix b 1-Undifferentiated spermatogonia I 9-Sertoli 8 2-Undifferentiated spermatogonia II 10-Endothelial 9 3-Differentiating spermatogonia (early) 11-Leydig 7 4-Differentiating spermatogonia (late) 12-Lymphoid B cell 1 5-Pre-pach spermatocytes 13-Lymphoid T cell 10 2 6-Pach-dip-2nd spermatocytes 14-Macrophage M1 3 14 15 17 7-Round spermatid I 15-Macrophane M2 4 12 8-Round Spermatid 2 16-Myoid 17-Neuro 6 16 5 11 13 Meiosis SPG SPC SPD c (−Inf,−0.457] 35 35 (−0.457,−0.248] (−0.248,−0.0219] (−0.0219,0.175] (0.175,0.314] (0.314, Inf] (−Inf,−0.436] 59 59 (−0.436,−0.275] (−0.275,−0.0794] (−0.0794,0.0195] (0.0195,0.0848] (0.0848, Inf] (−Inf,0.0331] 92 92 (0.0331,0.0674] (0.0674,0.669] (0.669,1.36] (1.36,1.5] (1.5, Inf] Fig. 4 | Decomposition of testis gene expression patterns with sparse decom- assigned to all cells. c By plotting the cell scores for three representative germ cell position of arrays (SDA). a We applied sparse decomposition analysis (SDA) to components on the tSNE, and as a function of pseudotime, it is apparent that testis scRNA-seq data from 12 human donors to identify latent factors (‘compo- transcription during spermatogenesis can be modeled as series of components nents’) representing gene modules. These components are defined by two vectors overlapping in time, coming on and off gradually on different timescales29. Shown – one that indicates the loading of each cell on the component, and one that are components with activity that start in spermatogonia (35), spermatocytes (59), indicates the loading of each gene on the component. b The same scRNA-seq data and spermatids (92). was summarized using a conventional tSNE analysis, and testicular cell types labels method for tensor decomposition (Fig. 4a, “Methods”). Each compo- https://conradlab.shinyapps.io/HISTA/#. These components represent nent is described by two matrices: a cell score matrix and a gene sets of genes with similar expression dynamics independent of cell loadings matrix. Together they provide a quantitative summary of type boundaries and are often driven by shared biology29. The activity which genes are co-expressed across which cells. We provide an of each component can be summarized by plotting the cell scores for interactive browser to view the SDA components in detail online at each component on a map of labeled types (Fig. 4b), and, in the case of Nature Communications | (2022)13:7953 6 Article https://doi.org/10.1038/s41467-022-35661-z germ cell components, plotting cells scores as a function of develop- species31, raising the fascinating possibility that there may be evolvable mental stage (pseudotime, Fig. 4c). In such a way components can be machinery that specifically controls germline, but not somatic, mito- identified that are specifically expressed in spermatogonia (compo- ses. Concordantly, we found that roughly half of the top-loading genes nent 35), spermatocytes (component 59) and spermatids (component on component 46 exhibited testis-enhanced expression (Supplemen- 92) (Fig. 4c). tary Fig. 7d). By comparing the profile of SDA genes loadings for GEMINI NOA Component 98 included the prioritized genes C2orf78, POTEJ, genes to the profile observed for (1) genes essential for spermato- POU4F2, and PROK2, and has a strong, specific expression pattern in genesis in mice based on MGI database and literature searches spermatocytes (Fig. 5d, Supplementary Fig. 7a–c). Loading in the top (denoted ‘MGI genes’; n =197) and (2) genes identified among the 60 genes of this component are 39 lincRNAs, and, in addition to controls of the institutional cohort (n = 131), we observed that the C2orf78, three other uncharacterized protein-coding genes (C17orf96, mouse MGI genes were more similar to NOA rather than to the control C9orf163, C9orf57). Three GO annotations are significant for this mys- genes (R = 0.54 vs R = 0.087; p = 3 × 10−4 Pearson-Filon test, Fig. 5a). terious component, all indicating functions in plasma membrane Unlike controls, NOA and MGI genes showed a very strong, specific cell–cell adhesion (GO:0007156, GO:0098742, GO:0016339, all p < enrichment on SDA components with peak expression in earlier phases 10−10). Component 122 contains the NOA genes INHBB, MAMLD1, and of germ cell development, differentiating spermatogonia and pre- SMIM1, the first two of which have previously characterized roles in pachytene spermatocytes (Fig. 5a). Specifically, component 51 was an gonadal function32,33. This component is specifically expressed in Ser- outlier in the absolute number of human and murine infertility genes toli cells, raising the important possibility that some forms of human (Fig. 5b, c) and showed strong enrichment of gene ontology (GO) spermatogenic impairment may be ascribed to somatic cell defects categories related to meiosis. At least half of the top 50 genes in this (Supplementary Fig. 7a, c). component have already been identified as associated with male infertility in humans, mice, or both, whereas several of the genes with a piRNA processing factors in NOA high specificity towards this component have yet to be characterized We prioritized potential disease variation in six genes that are essential in mammalian spermatogenesis, including RAD51AP2, PRAP1, C18orf63, for the biogenesis of PIWI-interacting RNAs (piRNA). piRNAs are small C5orf47, and the SYCP3 paralogs FAM9B and FAM9C (Fig. 5b). non-coding RNA molecules highly enriched within and critical for the We next compared the SDA profiles in the context of the histo- survival of the germ cell pool through silencing transposable elements logical diagnosis of the NOA cases and observed a different distribu- in fetal germ cells and for transcript storage and degradation in adult tion of gene loadings across the components for patients with MA meiotic and haploid germ cells34,35 (Fig. 6a). In male mice, disruption of compared to patients with SCO or “unknown” histology, the latter of any of the components of the piRNA biogenesis leads to detectable which is likely a mixture of all three histologic subtypes (p = 3.0 × 10-3, changes in the expression of mature piRNAs and causes spermato- hypergeometric test, Fig. 5a). Furthermore, the human SDA profile for genic arrest36–42. Concordantly, spermatogenic arrest was mostly MA genes was negatively correlated with the profile of SCO-associated characteristic to the eleven NOA cases who were affected by rare genes and genes from patients with unknown histology (R = −0.23) recessive variation in six piRNA biogenesis genes (PLD6, PNLDC1, indicating distinct molecular causes for different histological sub- RNF17, TDRD9, TDRD12, TDRKH; Supplementary Fig. 8a). TDRD943 and types. This observation of different loading patterns for genes linked RNF1744 have previously been found to be disrupted in men with to MA compared to SCO (p = 2.8 × 10−3, hypergeometric test) was azoospermia. replicated with data from 197 mouse models of male infertility with Remarkably, when queried against the SDA components, these 6 known histological observations (“Methods”). genes loaded highly and specifically on component 59 expressed in SDA components were classified into “MA-enriched” and “SCO- pre-pachytene spermatocytes (p < 0.01, Fig. 6b) together with other enriched” based on the corresponding case histology of component known piRNA processing genes PIWIL2, TEX19 and GPAT2, and a genes. Among the MA components with multiple gene loadings, potentially novel factor in piRNA processing, TDRD10 (Fig. 6c). Strik- SDA149 and SDA35 were expressed in undifferentiated and differ- ingly, 11 of the 15 top loadings on component 59 were long intergenic entiating spermatogonia and involved multiple known early germ cell non-coding RNAs annotated as pre-pre-piRNAs45, indicating that the markers, such as NANOS3, ID4, and DMRT1, but also DCAF12L1, a gene precursors of piRNA molecules are strongly co-regulated with the recurrently disrupted in NOA cases (Fig. 5d). The SDA94 component transcripts of the piRNA biogenesis genes themselves. represents an alternative molecular time-point affected in MA cases, We set out to validate the functional consequences in testicular meiotic spermatocytes, where cases were identified with mutations in tissue samples available for six NOA cases with prioritized variants in HIST1H1T and C1orf146 (a human ortholog of yeast SPO16). Similarly, PNLDC1, TDRD9 and TDRD12. PNLDC1, a trimmer of pre-piRNAs, was five SCO ‘molecular subforms’ were mapped to the spermatogonial cell recurrently disrupted in three men from different GEMINI centers and population, including SDA67, enriched for DNA replication genes, and one case from the Newcastle/Radboud cohort46, whereas two potentially SDA46, a component with unknown function (see below). Collectively, deleterious variants in TDRD9 and TDRD12, essential factors in second- these findings indicate that within the same histological endpoint, ary piRNA biogenesis, were found in singleton NOA cases from the different ‘molecular subforms’ can be found reflecting the known and GEMINI and MERGE cohorts, respectively (Supplementary Fig. 8a). Dis- novel molecular origins of azoospermia. ruption of all tested piRNA genes predominantly led to spermatogenic arrest seen in testicular biopsy (Fig. 6d, Supplementary Fig. 8a–c, e, f)46 SDA components with uncharacterized genes and the altered piRNA biogenesis was clearly detectable upon small RNA We noted several SDA components that included multiple prioritized sequencing with a shift towards longer immature piRNAs and sig- genes not previously linked to mammalian sperm production. For nificantly reduced fraction of mature piRNA molecules <32 nt compared instance, component 46, which included NOA case genes TUBA3C, to respective controls with normal spermatogenesis (1.1–1.3 fold CDC45, FOXM1, and GTSE1, cycles strongly two times in differentiating decrease, p < 0.01, hypergeometric test; Fig. 6e, f; Supplementary spermatogonia, corresponding to the number of mitotic divisions Fig. 8d, g; Supplementary Discussion)46. Collectively, these findings these cells undergo in humans and which are thought to be the demonstrate that components of piRNAs biogenesis are indispensable equivalent of transit amplification of stem cells in somatic tissues for human spermatogenesis and highlight the benefit of intersecting (Supplementary Fig. 7a–c)30. The number of mitotic divisions that putative disease genes with scRNAseq data for identification of diverse occur during spermatogonial amplification appears to vary among molecular subforms of NOA. Nature Communications | (2022)13:7953 7 Article https://doi.org/10.1038/s41467-022-35661-z a b Mutation previously found in: mouse human mouse+human neither Genes mutated in GEMINI patients TOP2A SMC1B Somatic cells Somat C1orf146 TEX101 Undifferentiated RAD51AP2 TEX12 MEIOB DMC1 spermatogonia SYCP1 TEX19 SPO11 spermatogonia DMRTC2 spermatocytes matocytes spermatocytes ZCWPW1 PRAP1 Differentiating h-dip-2nd Pach-dip-2nd Lymphocyte DPEP3 rmatids Spermatids Endotherial BRCA1 Pre-pach TOPAZ1 PRDM9 SYCP3 CCDC79 Myeloid SYCP2 dig Leydig Sertoli oli STAG3 SYCE1 C14orf39 HORMAD2 CEP85 MCM8 CDC45 Component 51 loading Cell type 0.5 BUB1 VRK1 C9orf84 BRCA2 SYCE3 CEP72 MEIOC 0 5 10 15 SSX3 Gene sets CCDC15 TDRD9 PSMC3IP MAGEA3 TMEM206 M1AP RNF212 MDC1 RTBDN CCDC155 Institutional controls 0.0 NOA cases Mouse −0.5 PMP22 51 AKR1C3 Germ cells Somatic cells 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 20 22 X Y M Histology - human Genomic coordinate 0 MA 1 d Meiosis 2 SPG SPC SPD SCO 3 Primary histology 4 149 HSG 5 MA 6 UNK LUZP4,TFDP3 149 35 51 67 98 52 35 SDA components MA c SSX3, ATP6V0E2, DCAF12L1 Component 51 cell score 67 SCO Sertoli DCAF4, DPEP3, YY2 51 SPG EC SPD Mixed SYCP1, STAG3, DMC1 98 -Inf Unknown -1 PROK2, C2orf78, POTEJ, POU4F2 -0.5 SPC 0 Leydig 52 0.5 Unknown 1 OLAH, SIRT1 Inf Pseudotime Fig. 5 | Using SDA components to define molecular subforms of genetic spermatocytes, SPD spermatids, UNK undetermined histology. b Expression of SDA infertility. a Heatmaps summarizing the count of genes loading on each SDA component 51 visualized on a t-SNE representation of the human testis scRNA-seq component. Top: The distribution of NOA candidate genes across components is dataset. This component is expressed in spermatogonia, and more highly in pre- more similar to the distribution of mouse infertility genes, compared to genes with leptotene spermatocytes. c Gene loadings for SDA component 51. Gene loadings rare damaging genotypes in population controls (p = 3 × 10−4, Pearson and Filon’s Z, reflect which genes are active in a component, with stronger positive or negative two-sided). Bottom: Distribution of NOA candidate genes by testicular histology loadings indicating greater expression changes in the component compared to the across SDA components. Genes found in patients with MA show a different pattern total dataset. d Germ cell expression of components with multiple gene loadings of expression compared to genes in SCO patients. SPG spermatogonia, SPC ordered by pseudotime, shown with representative genes and patient histologies. Discussion implicated over 200 genes in primary spermatogenic failure, including In this study, we performed whole-exome sequencing on an interna- 21 genes with nominally significant association with NOA when com- tional cohort of unrelated NOA cases and reported extensive genetic paring 2072 cases and 11,587 fertile controls. Thus, we provide two lists heterogeneity attributable to the recessive form of the disease. We of genes with two levels of statistical signal: 221 genes prioritized in Nature Communications | (2022)13:7953 8 Article https://doi.org/10.1038/s41467-022-35661-z a b Primary piRNA processing Biogenesis of secondary piRNAs Component 59 cell score NUCLEUS NUCLEUS -Inf -1 -0.5 0 0.5 1 Inf Germ cell Germ cell Sertoli TDRD9 X PIWIL4 SPG EC Primary Transcriptional piRNA PLD6 TDRD12 processing TDRKH silencing precursor 3’ PIWIL2 PNLDC1 PIWIL4 SPD 3’ 5’ piRNA 5’ 3’ intermediate Ping pong cycle PIWIL2 PIWIL4 RNF17 TDRKH 5’ 3’ 3’ 5’ SPC PIWIL1/ PIWIL2 PNLDC1 Leydig 5’ TDRKH PIWIL4 PIWIL2 PNLDC1 Not mutated 5’ 3’ Mature PIWIL1/ PIWIL2 Mutated piRNA c Non-codind RNAs Protein coding components Genes mutated in GEMINI patients d 0.5 TEX19 Component 59 gene loading GPAT2 0.0 RP11−620J15.3 TDRD9 RP11−386M24.9 FAM225B BUB1 TDRD12 NEU3 RNF17 RP11−451H23.3 PNLDC1 TMEM99 TDRD10 SPATA8 PLD6 PIWIL2 STRBP LINC00930 AP000350.5 RP11−844P9.2 ANKRD18A -0.5 LINC01206 LINC01015 RP11−403N16.3 RP13−735L24.1 CTC−523E23.6 RP3−468B3.2 LINC00865 CTD−2036P10.3 LINC00668 TDRG1 LINC01120 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17181920 22 X Y M Genomic coordinate e f 1.11x10-105 3.49x10-40 6.56x10-30 4.71x10-3 6.75x10-8 Normalized read count (x102) Control 15 TDRD12 case 0.8 Fraction of mature piRNAs Control Case 0.6 (<32 bases) 10 0.4 5 0.2 0 0.0 25 30 35 40 45 Portugal Denmark Netherlands Spain Germany piRNA length (nt) PNLDC1 TDRD9 TDRD12 Fig. 6 | Disruption of piRNA biogenesis causes spermatogenic failure in men. d H&E stain of TDRD9 patient biopsy showing spermatogenic arrest. The most a A schematic of piRNA biogenesis with the components affected among NOA cases mature germ cell observed were early round spermatids, which often appeared in this study highlighted in red. piRNAs are produced by two biogenesis pathways. multinucleated (arrows). In addition, many pyknotic cells were observed (arrow- The primary pathway involves transcription of long precursor-transcripts from heads). Scale bar = 50 µm. e Size distribution of piRNAs detected in the NOA case genomic clusters, which are then processed in the cytoplasm. b Cell scores for with a biallelic missense variant in TDRD12 and matching control, derived from component 59, which encodes both piRNA processing genes and target pre-pre- small RNA-sequencing of testis tissue. f Significantly decreased fraction of mature piRNAs, indicate that it is expressed primarily in pre-leptotene spermatocytes. (<32 bases) piRNAs was detected in testis of patients with piRNA pathway muta- c Gene loadings for component 59. Note TDRD10 (black box), a protein-coding tions, indicative of faulty processing of the immature piRNA transcripts. gene that has yet to be characterized for a role in mammalian piRNA processing. Nature Communications | (2022)13:7953 9 Article https://doi.org/10.1038/s41467-022-35661-z GEMINI cases, the majority of which were observed in n = 1 cases, with production across the full spectrum of sperm count, such has already an estimated 30% false discovery rate based on a control cohort been shown with some Y chromosome variants59. Furthermore, it is comparison (Supplementary Data 1). Second, we provide a higher imperative that future studies of NOA should consider oligogenic confidence list of 21 genes which have all been observed in multiple inheritance during case assessment, and population-based tests of NOA cases and have nominally significant association with NOA in a association that include heterozygous genotypes should be used as comparison with 11,587 fertile controls (Supplementary Table 4). The discovery and validation tools. observed detection rate of possible causes of NOA in GEMINI (20%) By intersecting the identified NOA disease genes with an scRNA- remains a conservative estimate, as our approach did not consider seq atlas of human adult testis, we were able to organize the underlying genes exclusively involved in fetal gonad development, excluded molecular pathologies beyond what is currently possible with histo- genes with broad expression across the body, and did not address logical classification. We identified gene expression components with dominant or additive defects likely to contribute to the manifestation clear sets of independent genes, involving distinct cellular processes of NOA as well47,48. A more detailed analysis of X-linked genes is (e.g. mitosis, meiosis, piRNA biogenesis) and distinct cell types (e.g. reported in a companion paper focusing exclusively on the X chro- Sertoli cells, and various germ cell types). In total we found over a mosome, using a different method of prioritization than described dozen clearly distinct components with robust GEMINI gene loadings, here49. and have statistically defined 70 components from total testis One remarkable conclusion of our study is that the Mendelian expression data. Thus the scope is clearly there to use gene expression causes of NOA are well distributed across the vast number of genes components as a rational basis for NOA patient classification with involved in testis function, and not clustered into a small number of higher resolution than our current systems based on histology. We are genes. The highest level of recurrence recognized in the GEMINI not proposing that gene expression signatures like SDA components cohort was n = 3 cases affected by mutations in the same gene, and will provide a comprehensive and definitive classification system for all most cases were singletons, affecting in total 200 genes. This places genetic forms of male infertility, just that such an approach can pro- NOA in a similar space of genetic architecture as congenital hearing vide a useful starting point for organizing the rapidly growing, and, in loss (151 genes with at least “limited” clinical validity)50, inherited some cases, not obviously related genetic causes that are being retinal degenerations (over 260 genes)51, and ciliopathies (over 190 identified. genes)52, but not yet approaching developmental disorders (over The use of components to categorize patients can have benefits to 1500 reportable genes)53. This observation has numerous implica- research and patient management. Most importantly, we showed that tions for the clinical management of genetic forms of male infertility. defects in multiple genes from the same component, 59, produced a First, it seems unlikely that a targeted panel of genes for clinical shared molecular signature of defective piRNA processing that could genetic testing will, in the short term, be an effective replacement for be detected with the same assay of small RNA sequencing. With further whole-exome or genome assays in this disease space. Second, most research, it may be possible to design similar companion diagnostics clinicians will never see multiple unrelated cases of infertility caused for each component, or molecular subform, that can be used in con- by mutations in the same gene. Looking across clinics to share junction with genetic analysis to confirm the pathogenicity of candi- experiences for diagnosis and treatment, and to construct case ser- date variants, even for variants and genes that have yet to be ies, will be essential to optimize patient care. From a research per- characterized. Examples of other molecular markers that could be spective, this scarcity of recurrent cases will challenge variant measured from patient tissue would be DNA methylation (e.g. for interpretation: this will need to be addressed by data sharing con- SPOCD1 mutations found in our study60), metabolite levels, or markers sortia, development of robust assays to evaluate the function of a of immune activity, although the potential number is large and many wide range of variants, and more creative, integrative use of model have been considered to date. The SDA approach simply provides a organisms for functional validation. framework to pair genes of uncertain function with other genes, many It is now clear that what is clinically referred to as “spermatogenic of which may already have well-characterized roles. Finally, we showed impairment” is a collection of numerous rare diseases. However, these that components can be used to map the cell types involved in a rare diseases share a related genetic basis; many of the affected genes defect. This would be critical knowledge to nominate patients with encode proteins that either interact physically, such as in a protein primary defects in somatic cells, who could benefit from powerful complex, or physiologically, in regulating a particular cellular process. in vivo therapies targeted to those somatic cells, without directly In such a system it is expected that oligogenic inheritance will occur, as modifying germ cells. While we are confident that broad germ cell well as variable expression of disease genotypes, where the precise groups (e.g. spermatogonia, spermatocytes, spermatids) can be con- phenotypic effect of a mutation will depend on genetic variation at fidently identified and discriminated with current scRNA-seq analysis other interacting loci, as well as the environment (e.g. patient age, methods, it is still unclear how precisely we can map individual cells to lifestyle)54. We noted that eleven genes prioritized in our study have fine-scale subtypes (e.g. Adark versus Apale spermatogonia) or even previously been linked to male infertility disorders other than NOA, timepoints in germ cell development. Our reservations are that tran- including hypogonadotropic hypogonadism, primary ciliary dyskine- scription and translation are decoupled in mammalian germ cell sia, and multiple morphological abnormalities of the sperm flagella development, and some subpopulations have only been well char- (MMAF; Fig. 3c, Supplementary Table 6, Supplementary Discussion)7,55. acterized with protein biomarkers and morphology. This is especially The latter suggests that problems in sperm assembly may also result in true with human germ cells, which in general are less characterized defects in sperm concentration, consistent with the observation of than mouse germ cells. Thus our description of cell types expressing severe oligozoospermia in some patients with MMAF;56,57,58. Clinical certain SDA components should be considered as preliminary and in data available for the affected GEMINI patients did not indicate mis- need of further investigation. diagnoses, and most of these genes are robustly expressed in adult Infertility treatment is the leading frontier of some genetic tech- testicular cell types relevant to NOA, raising the possibility of an nology, including diagnostic single-cell sequencing61 and therapeutic expanded phenotype for these loci. In total, these observations sug- mitochondrial donation62. The principles outlined in this discovery gest that variant and gene interpretation for NOA should typically be study indicate that genetic forms of male infertility have diverse performed jointly with genetic knowledge of other male infertility mechanisms, and that genetic diagnosis through whole-exome/gen- phenotypes. It should be expected that, while some subtypes of NOA ome sequencing will be an essential tool for developing and targeting may have a unique genetic basis not shared with other traits, many new forms of assisted reproductive therapies but also for risk evalua- variants that have a causal role in NOA will also modulate sperm tion of future comorbidities. Nature Communications | (2022)13:7953 10 Article https://doi.org/10.1038/s41467-022-35661-z Methods the YCMD or 47, XXY testing was not available for all subjects. To Ethics statements & IRBs from centers retrospectively identify and remove cases where YCMD or 47, XXY The study complies with all relevant research regulations and was karyotype was not originally assessed or was missed in the patient approved by the Ethics Committee of all collaborative centers: workup, we analyzed the genotype call set of the cases generated protocols #201107177 and #201109261 approved by the institu- based on the WES data (Supplementary Fig. 1; further details in Sup- tional review board (IRB) of Washington University in St. Louis, plementary Discussion). Structural abnormalities of sex chromosomes USA; IRB_00063950 approved by IRB of University of Utah, USA; (including YCMD) and autosomes were detected from the copy num- PTDC/SAU-GMG/101229/2008 approved by the Ethics Committee ber variation data, whereas cases of Klinefelter syndrome were iden- and Hospital Authority, University of Porto, Portugal; 16030459 tified by a combination of normalized coverage of chromosome Y and and 0102004794 approved by the IRB of Weill Cornell Medical the inbreeding coefficient F of chromosome X (Supplementary Fig. 1c) College, New York, USA; Ref. No.: 2014/04c approved by the IRB and confirmed by the presence of large copy number variants on of Fundació Puigvert, Barcelona, Spain; NL50495.091.14 version 4 chromosome X. approved by the Research Ethics Committee of the Radboud To identify potentially missed cases of obstructive azoospermia, University, Nijmegen, The Netherlands; Ref: 18/NE/0089, Bursa: the genotype call set was additionally screened for pathogenic mono- 05.01.2015/04 approved by the University Research Ethics Com- or biallelic ClinVar variants (curated Expert Panel list, July 2018)65 in the mittee of University of Newcastle, UK; and ethical approvals from CFTR gene, a known cause of obstructive azoospermia. We elected to human ethics committees of Monash Surgical Day Hospital, remove cases who were heterozygous for a single well-acknowledged Monash Medical Centre and Monash University, Australia; the pathogenic CFTR variant as a conservative measure. Due to the lim- ethics committee for the Capital Copenhagen Region (Ref. Nr. H- itations of WES to detect variants in non-exonic regions, we are not 2-2014-103) and the Danish Personal Data Protection Agency able to exclude the presence of other non-coding pathogenic variants (Datatilsynet 2012-58-0004, local Nr. 30-1482, I-Suite 03696) gave in the CFTR gene in trans with the coding variants66,67. Additionally, ethical approval for this work in accordance to the European identity by descent (IBD) analysis was performed using SNPRelate R Commission Directive for the transfer of personal data (MTA/I- package68 to identify and remove one counterpart of each accidental 4728.A1); the Ethics Committee of National Institute of Health Dr sample duplicate or twin pairs and cases with cryptic relatedness to Ricardo Jorge, Lisboa, Portugal gave ethical approval for this avoid ascertainment bias of rare variation. work; and approval 74/54 (last amendment 288/M-13) released by Research Ethics Committee of the University of Tartu, Estonia. Institutional control cohort for development of variation The MERGE study protocol was given ethical approval by the prioritization procedures Ethics Committee of the Ärztekammer Westfalen-Lippe and the A collection of individuals (total n = 2265, 851 of whom are men) University of Münster (Ref. No. 2010-578-f-S). Written informed sequenced with the identical exome platform as the GEMINI cases was consent was obtained from all men. used as institutional controls to refine the variant prioritization pipe- line and to test its ability to specifically aggregate potentially patho- Study subjects genic variation related to testicular function. This comparative dataset Anonymized DNA samples were collected for a cohort of 1011 unre- consisted of 79 patients from Alzheimer study as well as 2186 partici- lated men diagnosed with idiopathic NOA and recruited in four centers pants from a study of skeletal or brain malformations: cases with in the USA, two centers in Portugal (irreversibly anonymized) and a adolescent idiopathic scoliosis (AIS, n = 1245), Chiari malformation (n = total of five centers in Australia, Canada, Denmark, Estonia and Spain. 105), hypermobility (n = 38) and/or limb abnormalities (n = 798), and Thus, determination of NOA and infertility workup varied by practice 43 controls with no malformations. The AIS dataset has previously pattern across centers, but in general followed published clinical been partially reported69 and is fully available through the database of guidelines of American Urological Association/American Society for Genotypes and Phenotypes (dbGaP; https://www.ncbi.nlm.nih.gov/ Reproductive Medicine (AUA/ASRM) and European Academy of gap/; accession ID phs001677.v1.p1). The institutional cohort was Andrology (EAA)63,64. Men were confirmed to have azoospermia (no subjected to exome capture and sequencing methods identical to the spermatozoa detected in the whole ejaculates) or severe oligozoos- GEMINI cases (McDonnell Genome Institute, St. Louis, MO, USA) and permia/cryptozoospermia (n=9 subjects; extremely low concentration similar bioinformatic pipeline was applied to read alignment and joint of spermatozoa, <~1 million/mL; jointly termed as ‘NOA’ in this study) genotyping procedures. The prioritization of variants and genes likely according to the AUA/ASRM guidelines and based on physical exam- disrupting testicular function was performed in parallel with the ination (testis volume), endocrine measures (FSH, LH, and T) and his- GEMINI cases and was aiming to refine the prioritization strategy in an tological findings if available. In the case of mixed findings on inheritance mode-aware manner to maximize patient-specific calls. histology, i.e. the observation of MA tubules and SCO tubules in the Genes with prioritized variants in both cohorts were removed from same patient, we classified the patient as MA if any MA tubules both as a conservative step. were noted. Testis volume and presence and grade of varicocele were eval- Additional cohorts uated by palpation or ultrasound by a trained fertility specialist. The Two independent “replication” cohorts were screened for variation in exclusion criteria included the obstruction or absence of vas deferens, genes linked to NOA in the GEMINI study. Whole-exome sequencing of varicocele of bilateral grade 2–3 or unilateral grade 3, radical pelvic 900 men has been performed in the framework of the Male Repro- surgery, anejaculation, spinal cord injury, radiation treatments or ductive Genomics (MERGE) study and includes cases with NOA or chemotherapy exposure and environmental factors. CBAVD and other extremely low sperm counts (n = 817). Patients were recruited at the forms of obstructive azoospermia were ruled out based on a variety of Centre of Reproductive Medicine and Andrology (CeRA), Münster, factors including semen parameters (semen volume and qualitative Germany and the Clinic for Urology, Pediatric Urology and Andrology, fructose test), endocrine evaluation (normal FSH), normal testis Gießen, Germany. The MERGE cohort was queried for the 221 prior- volume and the presence of vas deferens detected by palpation and/or itized GEMINI genes to identify potential disease variation and retain ultrasound. In addition, men with positive findings of Y chromosome recessive changes that are rare (MAF < 0.01), non-synonymous variants microdeletions (YCMD) or karyotype abnormalities (including Kline- seen in the MERGE cohort less than 20 times, had alternate read count felter syndrome, 47, XXY) were excluded from the study. As a large <5, frequency >25 in all reads and CADD score >10 if the consequence is fraction of the GEMINI cohort are historical cases, the information on missense. The dataset used to compare to the GEMINI study was Nature Communications | (2022)13:7953 11 Article https://doi.org/10.1038/s41467-022-35661-z further filtered to match the GEMINI filtering procedures by excluding and defined as overlapping more than 50% of their respective targets. variation with PSAP P-value >10−3, the positions that were observed in Only CNVs with quality scores greater than or equal to 60 were more than three individuals and the compound heterozygous pairs included in the downstream analysis. As previously recommended81, that were found within a 10 bp proximity or occurred in more than one individuals having a CNV count greater than 3 SD above the mean (27 individual. DELs or 30 DUPs, in this study) were removed from the analysis. Whole-exome sequencing in the Newcastle/Radboud cohort was To evaluate the sensitivity of CNV detection from the WES dataset, performed on a total of 331 patients who presented with idiopathic a comparative CNV call set for the reference sample NOA (n = 164) or severe to extreme oligozoospermia (with or without NA12878 sequenced in parallel with the GEMINI study subjects, was asthenozoospermia; n = 167) at the Radboudumc outpatient clinic available through the whole-genome sequencing (WGS) effort of the between July 2007 to October 2017 (n = 286) and at the Newcastle 1000 Genomes phase 3 project82. The CNV calls of the 1000G dataset Fertility Clinic between January 2018 to January 2020 (n = 45). The were required to overlap with WES CNVs by at least 50% of length and presence of chromosomal anomalies, AZF deletions or pathogenic span the exonic regions targeted in this study. The sensitivity estimate CFTR variants were exclusion criteria, all patients remained idiopathic reached 13% (n = 11/85 of reference CNVs), which is within the expected following thorough clinical evaluation. Variants were filtered for range considering the fragmented coverage of the genome in WES recessive changes that were rare (MAF < 0.01), non-synonymous, seen approach and comparable to previously reported sensitivity estimates in less than 5 patients, with an alternate read count >5, and an alternate (7.67%) for CNV calling from exome sequencing data with XHMM when allele frequency >15%, predicted to be pathogenic by 3 or more of the compared to WGS17. In total, approximately 55% of the autosomal following variant prediction tools: SIFT70, PolyPhen-271, duplications (n = 5/9) and 60% of autosomal deletions (14/23) identi- MutationTaster , FATHMM , Mutation Assessor and CADD75. 72 73 74 fied in NA12878 in this study are supported by previous studies82 and All variants with read depth less than 50 were validated by Sanger are likely true positives. sequencing. Validation of CNVs by PCR Whole-exome sequencing and analysis A subset of likely homozygous/hemizygous deletions was selected for Whole-exome sequencing of the genomic DNA extracted from blood experimental validation by ±PCR, performed with at least one primer was performed at McDonnell Genome Institute of the Washington pair located between the predicted deletion breakpoints (internal) and University in St. Louis, MO, USA (genome.wustl.edu) using an in-house a second primer pair outside, either upstream or downstream from the exome targeting reagent capturing 39.1 Mb of exome and 2 × 150 bp deletion breakpoints in a region unaffected by the predicted deletion paired-end sequencing on Illumina HiSeq 4000. A subset of cases (5%) (external, positive control). CNV was confirmed when the internal set were processed with either Nextera Rapid Capture (Illumina, San of primers did not result in amplification and the external set of pri- Diego, USA) or Nextera Rapid Capture Expanded Exome kit according mers outside the deletion successfully amplified from patient DNA. to the manufacturer’s protocol and sequenced on HiSeq 2500 2 × 101 Both sets were required to amplify from control DNA. Polymerase bp or HiSeq 3000 in 2 × 150 paired-end mode. On average, 99.9% of reaction for ±PCR was performed using 50ng of DNA and Qiagen reads mapped to the target regions yielding an average exome cov- HotStarTaq (Qiagen, Hilden, Germany). Deletions were considered as erage of 80x across sequenced individuals and platforms. true positives when only amplification for the external primer pair(s) Raw sequencing reads were processed and aligned to GRCh38 was obtained. The primer sequences and annealing temperatures are assembly in an alternate contig aware manner using bwa-mem available upon request. v0.7.1776, Genome Analysis Toolkit v3.6.077 (GATK) and Picard tools v2.10.0 (http://broadinstitute.github.io/picard/). All sequenced cases Fine-mapping of the STRA8 deletion were tested for sample contamination with verifyBamID v1.1.378. Joint A homozygous ~6 kb deletion on chromosome 7 genotyping of all samples was performed with GATK tools and fol- (chr7:134,925,271–134,931,454, hg19) residing within the STRA8 gene lowing the best practices of GATK. In order to achieve a high-quality was identified in an NOA case GEMINI-295 by calling CNVs from WES call set, only samples with contamination FREEMIX <5%, average cov- data. We first validated the presence of this deletion in the patient DNA erage >30x, low excess of chimeric reads (<5%) and call rate >85% were using ±PCR with primers anchoring inside the deletion in exons 2–4 as included. Positions with high exome capture kit-specific missingness well as outside the deletion in exon 8 (Supplementary Table 7). Only (≥15%), low InbreedingCoeff (<‘−0.2’) and multi-allelic insertions and the second set of primers outside the deletion amplified from the deletions (INDELs) were filtered out, as well as individual genotypes patient DNA as expected (Fig. S6). We then tested four additional pairs with DP < 10, GQ < 30 or heterozygous variation in non-PAR regions of of primers in the intronic regions flanking the deletion to further chromosome Y. Precision and recall rate of detected variation was delimit the breakpoint. Finally, we designed primers for amplification calculated in relation to the CEPH individual NA12878 run in parallel across the expected breakpoint (gap-PCR) with a fragment of max- with the study cohort and for which a reference call set of high-quality imum expected size of 13 kb without the deletion. This primer set only variants is publicly available (Illumina Platinum Genomes79). The gen- resulted in a single product of ~5 kb in the patient indicative of a otyping achieved precision of 99% for SNVs and 91% for INDELs and a homozygous deletion (Supplementary Fig. 5). The acquired PCR pro- recall rate of 89% and 47%, respectively. duct was sequenced and allowed us to determine the precise locali- zation of the deletion breakpoints at chr7:134925172–134933449 Detection of copy number variation (hg19). The 5′ breakpoint is nearly equidistant (~400 bp) from an Copy number variants were detected from WES data using XHMM as L1ME3 repeat of the L1 family and an MLT1AD repeat of the ERVL-MaLR previously described80. CNVs were called using the Viterbi HMM family. algorithm with default XHMM parameters, and XHMM CNV quality PCR was performed using 50 ng of DNA and Qiagen HotStarTaq scores were calculated as previously described using the (Qiagen, Hilden, Germany). The long PCR spanning the deletion forward–backward HMM algorithm. In addition, all called CNVs were breakpoint was performed with NZYLong DNA polymerase (Nyztech, statistically genotyped across all samples using the same XHMM Lisbon, Portugal) at 62 °C and 6 minutes of annealing (35 cycles). The quality scores and output as a single uniformly called VCF file. ~5 kb PCR product was cloned into TOPO-TA vector (MilliporeSigma, The output of CNV calls was first separated into deletions (DEL) St. Louis, MO, USA) and a clean sequence was obtained, which con- and duplications (DUP) and each subset was then filtered by frequency firmed the breakpoints. Finally, we attempted to find heterozygous to remove common CNVs, which were present in >1% of individuals positions in the regions adjacent to the breakpoints to allow allele- Nature Communications | (2022)13:7953 12 Article https://doi.org/10.1038/s41467-022-35661-z specific analysis and confirm the existence of two deletions with and additionally meet either the prioritization criteria (1) or (2) men- exactly the same breakpoints but we did not find any heterozygous tioned above, whereas Y-linked variants were included if defined as positions in either of the introns sequenced (1, 7 and 8). Primers are LoF. Any NOA gene, which did not display testis-enhanced expression provided upon request. based on the Human Protein Atlas dataset, was required to at least show an appreciable expression specifically in premeiotic germ cells or Variant prioritization somatic cells of the testis (normalized TPM > 0.5 in the testis To prioritize likely pathogenic SNVs, INDELs and CNVs detected from scRNA-seq)28. the WES data, a modified version of the Population Sampling Prob- ability (PSAP) software was utilized, which identifies genotypes that are Burden testing unusual in the context of known human variation (https://github.com/ We performed gene-based burden testing using the combined data conradlab/PSAP/)16. PSAP evaluates the probability of sampling a from GEMINI, MERGE, and Newcastle/Radboud case cohorts, com- genotype or a set of genotypes based on the pathogenicity scores and pared to a control cohort of 11,587 fertile parents sequenced at Rad- frequencies of variants observed in the unaffected population. Unlike boud University. These fertile parents were ascertained from the classical workflows of case-control studies, PSAP enables the identifi- genetic diagnostic lab at Radboud, where they were referred for cation of potential causal variants from a single genome without the genetic testing as part of clinical care for their children. The cohort need for matching control samples. PSAP takes into account the sex consisted of approximately equal numbers of men (5803) and women and the ethnicity (European, African, and ‘other’) of the cases and tests (5,784), with an ethnic makeup reflecting the local Dutch population. for autosomal Mendelian disease models (autosomal dominant, single- Due to data use limitations, it was not possible to analyze compound variant autosomal recessive and compound-heterozygote autosomal heterozygous variants in the control cohort; therefore, burden testing recessive) but also X- and Y-linked inheritance patterns. The impact of was only performed using single-variant genotype calls. Case and fer- nucleotide changes positioned within multiple overlapping genes is tile control genotypes were processed using the same filtering and assessed independently for each gene. The genomic positions of prioritization pipeline described above, so the only variants included detected variation were lifted over to human genome assembly in the analysis were the short list of prioritized thought to represent hg19 supported by PSAP. potential Mendelian causes of NOA. The proportion of individuals with The prioritized variants were subsequently filtered to maximize filtered variants was compared between cases and controls using the true positive calls, very pathogenic rare changes and genes most Fisher Exact test. Testing was performed on the 34 genes with more likely to impair spermatogenesis. Variation was excluded if it had than one prioritized homozygous (autosomal) or hemizygous (X and Y PSAP P-values >10−3, minor allele frequency >0.01 across all popula- chromosomes) variant among all cases (primary and replication). This tions in the gnomAD database v2.1.1 (popmax, https://gnomad. type of selective testing of genes, which is predicated on features of broadinstitute.org/) or was common in the study cohort. Addition- the primary dataset, can be considered a two-stage study design21. We ally, genes enriched for rare pathogenic variation (PSAP P < 10−4, have analyzed both primary and replication datasets jointly in the popmax MAF < 0.01) and identical homozygous or hemizygous burden test. In this context, the resulting p-values should be corrected recessive changes seen in >3 individuals were removed. Compound for the total number of genes screened in the first stage, which is the heterozygous variants were subject to exclusion if identical pair of GEMINI-only analysis. The number of genes truly tested in the first- genotypes was observed in more than one case or were found in cis stage analysis of GEMINI-only data is not easy to define. Here we configuration as determined by pHASER v1.1.183, in-house parsing of conservatively used a value of p < 2.5 × 10−6 as our threshold for sig- variant annotation in the VCF file or manual check of read alignments nificance following multiple-test correction; this corresponds to a using IGV tools84. Bonferroni correction assuming 20,000 single-gene tests in the first In order to further aggregate variation most likely increasing stage (e.g. this is a threshold for ‘exome-wide’ significance). Mendelian disease risk and affecting spermatogenesis pathways, the following information was considered: (1) loss-of-function (LoF; non- Bioinformatic annotation sense, splice site, and frameshift variants) changes, (2) genes known to To determine the ancestry of GEMINI cases, EthSeq software87 was cause male infertility in mice (Mouse Genome Informatics, MGI, http:// implemented using Phase 3 genotype and population data from the www.informatics.jax.org/) or implicated in male infertility in humans 1000 Genomes Project as a reference82. Only variants with MAF > 0.2 (at least ‘limited’ supporting evidence, n = 164)7, and (3) a list of genes and located within the exonic regions targeted in the GEMINI WES with elevated expression in testis retrieved from the Human Protein study were considered. Long runs of homozygosity (ROH) were Atlas v18.185 (n = 2237). This dataset is a compilation of 172 individual detected using the H3M2 tool specifically designed for analyzing WES samples, corresponding to 37 different human tissues analyzed by bulk data88. ROH detection was performed separately for European, African, RNAseq86. The Human Protein Atlas project compared the gene East Asian and South Asian populations with default H3M2 parameters expression of ten normal testis samples to 162 other tissue samples (36 (DNorm = 100.000, p1 = 0.1, p2 = 0.1, F = 5). To identify the longest tissue types) to identify the list of genes with elevated expression in class of autozygous stretches most likely related to disease and testis. The genes were further grouped based on the level of tissue reflecting recent inbreeding89, unsupervised five-component cluster- specificity: tissue enriched (expression in one tissue at least 5-fold ing was performed in a population-specific manner using Mclust higher than all other tissues), group enriched (fivefold higher average (mclust package v.5.4.3 in R). Only the longest class (class 5) of the ROH TPM in a group of two to seven tissues compared to all other tissues) regions reflecting recent consanguinity and most likely contributing to and tissue enhanced (fivefold higher average TPM in one or more disease89 were considered when calculating the fraction of the auto- tissues compared to the mean TPM of all tissues)86. All three groups some being homozygous (FROH). The long ROH regions were inter- were considered for prioritization of genes with a likely impact on sected with CNV calls to exclude hemizygous regions. testicular function. To extract high quality LoF calls predicted to cause gene knock- The criteria for aggregating most plausible NOA variation were outs, homozygous, hemizygous and compound heterozygous LoF evaluated individually for each Mendelian disease model in parallel variation were annotated and filtered using the LOFTEE tool (v1.0.3; with the cohort of institutional controls aiming to maximize patient- https://github.com/konradjk/loftee) implemented in Variant Effect specific calls. All autosomal variants were required to meet at least one Predictor (v99)90 as per Karczewski et al. 201991. Comparative popu- of the three aggregation criteria mentioned above. X-linked variants lation frequency of human “knock-outs” for the respective genes was were expected to be located in genes with testis-enhanced expression acquired by screening large variation databases for LoF variation and Nature Communications | (2022)13:7953 13 Article https://doi.org/10.1038/s41467-022-35661-z performing LOFTEE filtering if not already previously applied intern- of arrays (SDA) framework29. SDA soft clusters genes and cells into ally by the databases. The screened databases included gnomAD latent “components”, which are then manually interpreted and sepa- (v2.1.1. and v3; accessed 02/14/2020), deCODE (Supplemental Table S4 rated into “technical noise” and “signal”. Components corresponding in Sulem et al. 2015)92 and Human Gene Mutation Database (HGMD; 12/ to technical noise are removed from the data, thus providing batch 20/2014 data freeze)93. Additionally, Born in Bradford, Birmingham correction. For the analyses described, here, SDA was run with 150 project, and East London Genes and Health (http://www. components, half of which were removed as technical noise. Gene genesandhealth.org/)94 exome data was aggregated. ontology enrichment analysis was performed on the top 250 genes The functional protein association network with the prioritized from each component (from each side) using the enrichGO function list of disease genes was built using protein-protein interaction data of from the clusterProfiler R package in which p-values are calculated the STRING database v11.0 (stringApp v1.5.0) integrated in the Cytos- based on the hypergeometric distribution and corrected for testing of cape platform (v3.7.2)95. Only high-confidence interactions (>0.7) were multiple biological process GO terms using the Benjamini-Hochberg considered. Genes located on chromosome Y are not mapped by the procedure. STRING database and are not included into the analysis. The resulting network was visualized using a prefuse force directed layout based on Small RNA sequencing and analysis the protein-protein interaction scores. RNA was extracted from the fixed testicular tissue of the NOA cases The median gene expression (TPM) for genes of interest across 52 with prioritized variants in genes involved in piRNA biogenesis and human tissues, including testis, was extracted from the The Genotype- controls with complete spermatogenesis matched to each case by Tissue Expression (GTEx) V8 RNA-Seq dataset27. To explore the fixative. Small RNA sequencing libraries were prepared using CATS expression patterns of GEMINI disease genes across GTEx tissues in Small RNA-seq kit (Diagenode, Cat. #C05010040) and sequenced on comparison to other recessive Mendelian disorders, a list of genes (n = the MiSeq platform (Illumina, San Diego, CA). Raw sequencing reads 2588) linked to Mendelian diseases was retrieved from the Centers for were processed according to the CATS Small RNA-seq protocol and Mendelian Genomics dataset (http://mendelian.org/phenotypes- reads with 25-45 bases in length were mapped to the human assembly genes; accessed Oct 01, 2020). The gene list was then intersected hg19 using bowtie v1.0.199 allowing one mismatch. The mapped reads with an OMIM dataset (retrieved May 15, 2020) to only extract genes were filtered to remove all reads mapping to small non-coding RNAs that follow a recessive inheritance pattern, which includes autosomal other than piRNAs (DASHR v2.0100) and the remaining reads were as well as X-linked and Y-linked recessive modes. GTEx expression intersected with 205 piRNA loci previously identified in the adult values were available for 604 genes out of the 615 defined as recessive human testis45. Individual piRNA counts were normalized to the spike- Mendelian disease genes and the difference in the expression levels, in cel-miR-39-3p (5′-UCACCGGGUGUAAAUCAGCUUG-3′) added to the when compared to the GEMINI disease genes, was evaluated using samples in the first step of the sequencing library preparation. Spike-in Mann-Whitney U test (p-value below 0.05 was considered significant). reads were mapped to the human reference using bowtie v1.0.1 and To explore the functional role of prioritized disease genes in mice, allowing no mismatches. information on the reproductive phenotypes of knock-out mouse models was extracted from the MGI database, International Mouse Sanger sequencing Phenotyping Consortium (www.mousephenotype.org)96 and literature. The regions of interest were amplified using Hot Start Taq 2x To estimate p, the number of monogenic azoospermia genes in Master Mix (New England Biolabs, Ipswich, MA, USA), 0.4 µM non-consanguineous men, we used monte carlo simulations to evalu- primers and a touch-down PCR program. The amplified PCR ate a model that relates p to the expected rate of recurrent gene hits in products were purified with MinElute PCR Purification kit (Qia- a sequencing study, which we denote as r. We assume that for a given gen, Hilden, Germany) or GFX PCR DNA and Gel Band Purification set of azoospermia genes, each gene in the set is equally likely to be a Kit (GE Healthcare, Marlborough, MA) and sequenced at Inte- monogenic cause of disease (e.g. the amount and effect size of disease grated DNA Technologies (Coralville, IA, USA). The primer variation in each gene is the same). Next, we evaluate the likelihood of sequences are available upon request. a particular value of p using the difference in the expected rate of recurrent hits, r, and the observed rate ro (9.4% among non- Immunohistochemistry consanguineous NOA cases) as an objective function. We searched Testicular tissue was fixed in Bouin’s fixative, processed and over a grid of 1000 values of pi, i ranging from 1 to 1000, with p1 = 1, p2 = stained with an antibody directed against STRA8 (Abcam 2, and so forth. For each value of pi, we simulate a set of “gene hits” by ab49405, diluted 1:250) using the ImmPRESS detection kit (Vector random sampling with replacement 221 genes, the actual number of laboratories, CA, USA) as described before101. Antigen retrieval prioritized genes observed in GEMINI, from the pool of pi genes. This is was performed in a microwave in TEG buffer and the staining done 1000 times for each value of pi, each time recording the corre- developed 20 min with AEC. sponding value of ri, and then averaging across all 1000 replicates to arrive at a final estimate for ri. The objective function, ri − rO, was Reporting summary minimized at i = 625. The sensitivity of detecting multiple gene hits as a Further information on research design is available in the Nature function of the cohort size was then calculated considering the esti- Portfolio Reporting Summary linked to this article. mated size of the azoospermia gene pool (n = 625), discovery rates of 15% for outbred samples, 76% for consanguineous samples and 20% for Data availability the GEMINI cohort overall (8% of consanguineous cases) and assuming Sequencing data generated from patients that have been consented the absence of natural selection. for data sharing in this study have been deposited in the dbGaP data- base under accession code phs003115.v1.p1 These data are available to scRNAseq analyses qualified researchers via a standard controlled access procedure We created an aggregated dataset of human testis scRNA-seq from managed by dbGaP. The scRNA-seq data from human testis used in this multiple sources, comprising cells from a total of 12 human donors. study are available in the GEO database under accession codes: These consisted of 7 adults28, 3 adults97 (GEO accession #GSE109037), GSE109037, GSE169062, GSE120506. These data are publicly available and two juvenile samples98 (GEO accession #GSE120506 [https://www. without any access restrictions. The SDA decomposition of human ncbi.nlm.nih.gov/geo/query/acc.cgi]). As described elsewhere, these testis data are publicly available without restriction through the HISTA datasets were integrated and analyzed using the sparse decomposition browser: https://conradlab.shinyapps.io/HISTA/#. Nature Communications | (2022)13:7953 14 Article https://doi.org/10.1038/s41467-022-35661-z Code availability 21. Barrett, J. C. et al. New mutations, old statistical challenges. Pre- The PSAP pipeline used in this manuscript is available from https:// print at bioRxiv https://doi.org/10.1101/115964 (2017). github.com/conradlab/PSAP. 22. Lim, E. T. et al. Rare complete knockouts in humans: population distribution and significant role in autism spectrum disorders. References Neuron 77, 235–242, https://doi.org/10.1016/j.neuron.2012.12. 1. Jarow, J. P., Espeland, M. A. & Lipshultz, L. I. Evaluation of the 029 (2013). Azoospermic patient. J. Urol. 142, 62–65 (1989). 23. Fon Tacer, K. et al. MAGE cancer-testis antigens protect the 2. Willott, G. M. Frequency of azoospermia. Forensic Sci. Int. 20, mammalian germline under environmental stress. Sci. Adv. 5, 9–10 (1982). eaav4832 (2019). 3. Itoh, N., Kayama, F., Tatsuki, T. J. & Tsukamoto, T. Have sperm 24. Anderson, E. L. et al. Stra8 and its inducer, retinoic acid, counts deteriorated over the past 20 years in healthy, young regulate meiotic initiation in both spermatogenesis and Japanese men? Results from the Sapporo area. J. Androl. 22, oogenesis in mice. Proc. Natl Acad. Sci. USA 105, 14976–14980 40–44 (2001). (2008). 4. Tüttelmann, F. et al. Clinical experience with azoospermia: 25. Richards, S. et al. Standards and guidelines for the interpretation aetiology and chances for spermatozoa detection upon biopsy. of sequence variants: a joint consensus recommendation of the Int. J. Androl. 34, 291–298 (2011). American College of Medical Genetics and Genomics and the 5. Giudice, F. D. et al. Clinical correlation among male infertility and Association for Molecular Pathology. Genet. Med. 17, overall male health: a systematic review of the literature. Investig. 405–424 (2015). Clin. Urol. 61, 355 (2020). 26. Djureinovic, D. et al. The human testis-specific proteome defined 6. Tiepolo, L. & Zuffardi, O. Localization of factors controlling sper- by transcriptomics and antibody-based profiling. Mol. Hum. matogenesis in the nonfluorescent portion of the human Y chro- Reprod. 20, 476–488 (2014). mosome long arm. Hum. Genet. 34, 119–124 (1976). 27. Consortium, T. G., The GTEx Consortium. The GTEx Con- 7. Oud, M. S. et al. A systematic review and standardized clinical sortium atlas of genetic regulatory effects across human validity assessment of male infertility genes. Hum. Reprod. 34, tissues. Science 369, 1318–1330, https://doi.org/10.1126/ 932–941 (2019). science.aaz1776 (2020). 8. Kasak, L. & Laan, M. Monogenic causes of non-obstructive 28. Mahyari, E. et al. Comparative single-cell analysis of biopsies azoospermia: challenges, established knowledge, limitations and clarifies pathogenic mechanisms in Klinefelter syndrome. Am. J. perspectives. Hum. Genet. https://doi.org/10.1007/s00439-020- Hum. Genet. 108, 1924–1945 (2021). 02112-y (2020). 29. Jung, M. et al. Unified single-cell analysis of testis gene regulation 9. Nagirnaja, L., Aston, K. I. & Conrad, D. F. Genetic intersection of and pathology in five mouse strains. Elife 8, e43966 (2019). male infertility and cancer. Fertil. Steril. 109, 20–26 (2018). 30. Ehmcke, J. & Schlatt, S. A revised model for spermatogonial 10. Riera-Escamilla, A. et al. Sequencing of a ‘mouse azoospermia’ expansion in man: lessons from non-human primates. Reproduc- gene panel in azoospermic men: identification of RNF212 and tion 132, 673–680 (2006). STAG3 mutations as novel genetic causes of meiotic arrest. Hum. 31. Hess, R. A. & de Franca, Renato L. Spermatogenesis and cycle Reprod. 34, 978–988 (2019). of the seminiferous epithelium. Adv. Exp. Med. Biol. 636, 11. Fakhro, K. A. et al. Point-of-care whole-exome sequencing of 1–15 (2008). idiopathic male infertility. Genet. Med. 20, 1365–1373 (2018). 32. Anderson, R. A. & Sharpe, R. M. Regulation of inhibin production in 12. Krausz, C. et al. Genetic dissection of spermatogenic arrest the human male and its clinical applications. Int. J. Androl. 23, through exome analysis: clinical implications for the management 136–144 (2000). of azoospermic men. Genet. Med. 22, 1956–1966 (2020). 33. Ogata, T., Sano, S., Nagata, E., Kato, F. & Fukami, M. MAMLD1 and 13. Chen, S. et al. Whole-exome sequencing of a large Chinese 46,XY disorders of sex development. Semin. Reprod. Med. 30, azoospermia and severe oligospermia cohort identifies novel 410–416 (2012). infertility causative variants and genes. Hum. Mol. Genet. 29, 34. Ernst, C., Odom, D. T. & Kutter, C. The emergence of piRNAs 2451–2459 (2020). against transposon invasion to preserve mammalian genome 14. Zhang, H. et al. Whole exome sequencing identifies genes asso- integrity. Nat. Commun. 8, 1411 (2017). ciated with non-obstructive azoospermia. Front. Genet. 13, 35. Fu, Q. & Wang, P. J. Mammalian piRNAs: Biogenesis, function, and 872179 (2022). mysteries. Spermatogenesis 4, e27889 (2014). 15. Alhathal, N. et al. A genomics approach to male infertility. Genet. 36. Ding, D. et al. PNLDC1 is essential for piRNA 3’ end trimming and Med. 22, 1967–1975 (2020). transposon silencing during spermatogenesis in mice. Nat. Com- 16. Wilfert, A. B. et al. Genome-wide significance testing of variation mun. 8, 819 (2017). from single case exomes. Nat. Genet. 48, 1455–1461 (2016). 37. Huang, H. et al. piRNA-associated germline nuage formation 17. Tan, R. et al. An evaluation of copy number variation detection and spermatogenesis require MitoPLD profusogenic tools from whole-exome sequencing data. Hum. Mutat. 35, mitochondrial-surface lipid signaling. Dev. Cell 20, 376–387 899–907 (2014). (2011). 18. Lopes, A. M. et al. Human spermatogenic failure purges deleter- 38. Watanabe, T. et al. MITOPLD is a mitochondrial protein essential ious mutation load from the autosomes and both sex chromo- for nuage formation and piRNA biogenesis in the mouse germline. somes, including the gene DMRT1. PLoS Genet. 9, Dev. Cell 20, 364–375 (2011). e1003349 (2013). 39. Wasik, K. A. et al. RNF17 blocks promiscuous activity of PIWI pro- 19. Wyrwoll, M. J. et al. Bi-allelic Mutations in M1AP are a frequent teins in mouse testes. Genes Dev. 29, 1403–1415 (2015). cause of meiotic arrest and severely impaired spermatogenesis 40. Saxe, J. P., Chen, M., Zhao, H. & Lin, H. Tdrkh is essential for leading to male infertility. Am. J. Hum. Genet. https://doi.org/10. spermatogenesis and participates in primary piRNA biogenesis in 1016/j.ajhg.2020.06.010 (2020). the germline. EMBO J. 32, 1869–1885 (2013). 20. Fang, K. et al. Prediction and validation of mouse meiosis-essential 41. Shoji, M. et al. The TDRD9-MIWI2 complex is essential for piRNA- genes based on spermatogenesis proteome dynamics. Mol. Cell. mediated retrotransposon silencing in the mouse male germline. Proteom. 20, 100014 (2021). Dev. Cell 17, 775–787 (2009). Nature Communications | (2022)13:7953 15 Article https://doi.org/10.1038/s41467-022-35661-z 42. Pandey, R. R. et al. Tudor domain containing 12 (TDRD12) is Society for Reproductive Medicine Report on optimal evaluation essential for secondary PIWI interacting RNA biogenesis in mice. of the infertile male. Fertil. Steril. 86, S202–S209 (2006). Proc. Natl Acad. Sci. USA 110, 16492–16497 (2013). 64. Dohle, G. R. et al. EAU guidelines on male infertility. Eur. Urol. 48, 43. Arafat, M. et al. Mutation in TDRD9 causes non-obstructive 703–711 (2005). azoospermia in infertile men. J. Med. Genet. 54, 65. Landrum, M. J. et al. ClinVar: improving access to variant inter- 633–639 (2017). pretations and supporting evidence. Nucleic Acids Res. 46, 44. Li, Z. et al. Excess of rare variants in genes that are key epigenetic D1062–D1067 (2018). regulators of spermatogenesis in the patients with non- 66. de Souza, D. A. S., Faucz, F. R., Pereira-Ferrari, L., Sotomaior, V. S. & obstructive azoospermia. Sci. Rep. 5, 8785 (2015). Raskin, S. Congenital bilateral absence of the vas deferens as an 45. Girard, A., Sachidanandam, R., Hannon, G. J. & Carmell, M. A. A atypical form of cystic fibrosis: reproductive implications and germline-specific class of small RNAs binds mammalian Piwi genetic counseling. Andrology 6, 127–135 (2018). proteins. Nature 442, 199–202 (2006). 67. Ge, B. et al. A rare frameshift variant in trans with the IVS9-5T allele 46. Nagirnaja, L. et al. Variant, defective piRNA processing, and of CFTR in a Chinese pedigree with congenital aplasia of vas azoospermia. N. Engl. J. Med. 385, 707–719 (2021). deferens. J. Assist. Reprod. Genet. 36, 2541–2545 (2019). 47. Gunes, S., Arslan, M. A., Hekim, G. N. T. & Asci, R. The role of 68. Zheng, X. et al. A high-performance computing toolset for relat- epigenetics in idiopathic male infertility. J. Assist. Reprod. Genet. edness and principal component analysis of SNP data. Bioinfor- 33, 553–569 (2016). matics 28, 3326–3328 (2012). 48. Oud, M. S. et al. A de novo paradigm for male infertility. Nat. 69. Haller, G. et al. A missense variant in SLC39A8 is associated with Commun. 13, 154 (2022). severe idiopathic scoliosis. Nat. Commun. 9, 4171 (2018). 49. Riera-Escamilla, A. et al. Large-scale analyses of the X chromo- 70. Kumar, P., Henikoff, S. & Ng, P. C. Predicting the effects of coding some in 2,354 infertile men discover recurrently affected genes non-synonymous variants on protein function using the SIFT associated with spermatogenic failure. Am. J. Hum. Genet. https:// algorithm. Nat. Protoc. 4, 1073–1081 (2009). doi.org/10.1016/j.ajhg.2022.06.007 (2022). 71. Adzhubei, I. A. et al. A method and server for predicting damaging 50. DiStefano, M. T. et al. ClinGen expert clinical validity curation of missense mutations. Nat. Methods 7, 248–249 (2010). 164 hearing loss gene-disease pairs. Genet. Med. 21, 72. Schwarz, J. M., Rödelsperger, C., Schuelke, M. & Seelow, D. 2239–2247 (2019). MutationTaster evaluates disease-causing potential of sequence 51. Duncan, J. L. et al. Inherited retinal degenerations: current alterations. Nat. Methods 7, 575–576 (2010). landscape and knowledge gaps. Transl. Vis. Sci. Technol. 7, 73. Shihab, H. A. et al. Predicting the functional, molecular, and 6 (2018). phenotypic consequences of amino acid substitutions using hid- 52. Wheway, G., Mitchison, H. M. & Genomics England Research den Markov models. Hum. Mutat. 34, 57–65 (2013). Consortium. Opportunities and challenges for molecular under- 74. Reva, B., Antipin, Y. & Sander, C. Predicting the functional impact standing of ciliopathies–the 100,000 Genomes Project. Front. of protein mutations: application to cancer genomics. Nucleic Genet. 10, 127 (2019). Acids Res. 39, e118 (2011). 53. Thormann, A. et al. Flexible and scalable diagnostic filtering of 75. Kircher, M. et al. A general framework for estimating the relative genomic variants using G2P with Ensembl VEP. Nat. Commun. 10, pathogenicity of human genetic variants. Nat. Genet. 46, 2373 (2019). 310–315 (2014). 54. Kousi, M. & Katsanis, N. Genetic modifiers and oligogenic inheri- 76. Li, H. & Durbin, R. Fast and accurate short read alignment tance. Cold Spring Harb. Perspect. Med. 5, a017145 (2015). with Burrows-Wheeler transform. Bioinformatics 25, 55. Chivukula, R. R. et al. A human ciliopathy reveals essential func- 1754–1760 (2009). tions for NEK10 in airway mucociliary clearance. Nat. Med. 26, 77. McKenna, A. et al. The Genome Analysis Toolkit: a MapReduce 244–251 (2020). framework for analyzing next-generation DNA sequencing data. 56. Dong, F. N. et al. Absence of CFAP69 causes male infertility due to Genome Res. 20, 1297–1303 (2010). multiple morphological abnormalities of the flagella in human and 78. Jun, G. et al. Detecting and estimating contamination of human mouse. Am. J. Hum. Genet. 102, 636–648 (2018). DNA samples in sequencing and array-based genotype data. Am. 57. Tang, S. et al. Biallelic mutations in CFAP43 and CFAP44 cause J. Hum. Genet. 91, 839–848 (2012). male infertility with multiple morphological abnormalities of the 79. Eberle, M. A. et al. A reference data set of 5.4 million phased sperm flagella. Am. J. Hum. Genet. 100, 854–864 (2017). human variants validated by genetic inheritance from sequencing 58. Martinez, G. et al. Whole-exome sequencing identifies muta- a three-generation 17-member pedigree. Genome Res. 27, tions in FSIP2 as a recurrent cause of multiple morphological 157–164 (2017). abnormalities of the sperm flagella. Hum. Reprod. 33, 80. Fromer, M. et al. Discovery and statistical genotyping of copy- 1973–1984 (2018). number variation from whole-exome sequencing depth. Am. J. 59. Rozen, S. G. et al. AZFc deletions and spermatogenic failure: a Hum. Genet. 91, 597–607 (2012). population-based survey of 20,000 Y chromosomes. Am. J. Hum. 81. Ruderfer, D. M. et al. Patterns of genic intolerance of rare copy Genet. 91, 890–896 (2012). number variation in 59,898 human exomes. Nat. Genet. 48, 60. Zoch, A. et al. SPOCD1 is an essential executor of piRNA-directed 1107–1111 (2016). de novo DNA methylation. Nature 584, 635–639 (2020). 82. Consortium, T. 1000 G. P. & The 1000 Genomes Project Con- 61. Huang, L., Ma, F., Chapman, A., Lu, S. & Xie, X. S. Single-cell sortium A global reference for human genetic variation. Nature whole-genome amplification and sequencing: methodology 526, 68–74, https://doi.org/10.1038/nature15393 (2015). and applications. Annu. Rev. Genomics Hum. Genet. 16, 83. Castel, S. E., Mohammadi, P., Chung, W. K., Shen, Y. & Lappalai- 79–102 (2015). nen, T. Rare variant phasing and haplotypic expression from RNA 62. Kang, E. et al. Mitochondrial replacement in human oocytes car- sequencing with phASER. Nat. Commun. 7, 12817 (2016). rying pathogenic mitochondrial DNA mutations. Nature 540, 84. Robinson, J. T. et al. Integrative genomics viewer. Nat. Biotechnol. 270–275 (2016). 29, 24–26 (2011). 63. Male Infertility Best Practice Policy Committee of the American 85. Uhlén, M. et al. Proteomics. Tissue-based map of the human pro- Urological Association & Practice Committee of the American teome. Science 347, 1260419 (2015). Nature Communications | (2022)13:7953 16 Article https://doi.org/10.1038/s41467-022-35661-z 86. Pineau, C. et al. Cell type-specific expression of testis elevated Research Foundation Clinical Research Unit ‘Male Germ Cells’ grant genes based on transcriptomics and antibody-based proteomics. DFG CRU326 (C.F., F.T.) The Netherlands Organization for Scientific J. Proteome Res. 18, 4215–4230 (2019). Research VICI grant 918-15-667 (J.A.V.) Investigator Award in Science 87. Romanel, A., Zhang, T., Elemento, O. & Demichelis, F. EthSEQ: from the Wellcome Trust grant 209451 (J.A.V.). FCT/MCTES, through ethnicity annotation from whole exome sequencing data. Bioin- national funds attributed to Centre for Toxicogenomics and Human formatics 33, 2402–2404 (2017). Health—ToxOmics, grant UID/BIM/00009/2016 (J.Go.) National Insti- 88. Magi, A. et al. H3M2: detection of runs of homozygosity from tute of Mental Health of the National Institutes of Health grant T32- whole-exome sequencing data. Bioinformatics 30, MH014677 (W.-L.C). 2852–2859 (2014). 89. Pemberton, T. J. & Szpiech, Z. A. Relationship between deleterious Author contributions variation, genomic autozygosity, and disease risk: insights from L.N., A.M.L., K.I.A., and D.F.C. designed the study. L.N., A.M.L., and the 1000 Genomes Project. Am. J. Hum. Genet. 102, D.F.C. performed computational and statistical analysis. W.-L.C., 658–675 (2018). B.M., T.L., E.M., S.R.C., L.K., and E.F. contributed to computational 90. McLaren, W. et al. The ensembl variant effect predictor. Genome analysis. L.N., R.S., I.G., A.S., L.K. and K.A. performed validation Biol. 17, 122 (2016). experiments. W.-L.C. and C.A.G. provided and analyzed whole- 91. Karczewski, K. J. et al. Variation across 141,456 human exomes and exome data from institutional controls. C.F., M.O., M.J.X., A.P., S.K., genomes reveals the spectrum of loss-of-function intolerance J.A.V. and F.T. participated in the replication studies. L.N., A.M.L. across human protein-coding genes. Nature 581, 434–443 and D.F.C wrote the manuscript. L.N., A.M.L., K.V.-C., K.A. and (2020). D.F.C. generated manuscript figures. S.S., F.C., S.F., A.B., J.Go., 92. Sulem, P. et al. Identification of a large set of rare complete human R.M., K.R.M., E.S.J., N.J., E.R.-D.M., K.A.J., M.P., J.I.S., M.L.E., K.O., knockouts. Nat. Genet. 47, 448–452 (2015). D.J., J.S., J.Gu., E.R.J., T.G.J., J.M.H., D.T.C., P.N.S., M.O’B., M.L., 93. Stenson, P. D. et al. The Human Gene Mutation Database: towards C.K., A.R.-E., S.K., K.A., and K.I.A. provided patient samples and a comprehensive repository of inherited mutation data for medical data. K.V.-C., N.J., E.R.-D.M., K.A.J., M.P., J.I.S., M.L.E., T.G.J., J.M.H., research, genetic diagnosis and next-generation sequencing stu- D.T.C., P.N.S., M.O’B., M.L., and C.K. contributed to the study dies. Hum. Genet. 136, 665–677 (2017). design, discussion and manuscript preparation as the GEMINI 94. Finer, S. et al. Cohort profile: East London genes & health (ELGH), a consortium members. All authors reviewed the manuscript. community-based population genomics and health study of Brit- ish Bangladeshi and British Pakistani people. Int. J. Epidemiol. Competing interests https://doi.org/10.1093/ije/dyz174 (2019). The authors declare no competing interests. 95. Shannon, P. et al. Cytoscape: a software environment for inte- grated models of biomolecular interaction networks. Genome Additional information Res. 13, 2498–2504 (2003). Supplementary information The online version contains 96. Dickinson, M. E. et al. High-throughput discovery of novel devel- supplementary material available at opmental phenotypes. Nature 537, 508–514 (2016). https://doi.org/10.1038/s41467-022-35661-z. 97. Guo, J. et al. The adult human testis transcriptional cell atlas. Cell Res. 28, 1141–1157 (2018). Correspondence and requests for materials should be addressed to 98. Guo, J. et al. The dynamic transcriptional cell atlas of testis Donald F. Conrad. development during human puberty. Cell Stem Cell 26, 262–276.e4 (2020). Peer review information Nature Communications thanks Jon Oatley and 99. Langmead, B., Trapnell, C., Pop, M. & Salzberg, S. L. Ultrafast and the other, anonymous, reviewer(s) for their contribution to the peer memory-efficient alignment of short DNA sequences to the review of this work. human genome. Genome Biol. 10, R25 (2009). 100. Leung, Y. Y. et al. DASHR: database of small human noncoding Reprints and permissions information is available at RNAs. Nucleic Acids Res. 44, D216–D222 (2016). http://www.nature.com/reprints 101. Nielsen, J. E. et al. Characterisation and localisation of the endo- cannabinoid system components in the adult human testis. Sci. Publisher’s note Springer Nature remains neutral with regard to jur- Rep. 9, 12866 (2019). isdictional claims in published maps and institutional affiliations. Acknowledgements Open Access This article is licensed under a Creative Commons We thank all of the study participants and the numerous medical staff Attribution 4.0 International License, which permits use, sharing, that enabled this study. We also acknowledge the contribution of the adaptation, distribution and reproduction in any medium or format, as Lisbon clinical team - Carlos Calhaz-Jorge, Ana Aguiar, Joaquim Nunes, long as you give appropriate credit to the original author(s) and the and Sandra Sousa (Unidade de Medicina da Reproducão, Hospital de source, provide a link to the Creative Commons license, and indicate if Santa Maria, Centro Hospitalar Lisboa Norte, Lisboa, Portugal), Maria changes were made. The images or other third party material in this Graça Pinto and Sónia Correia (Centro de Medicina Reprodutiva), for article are included in the article’s Creative Commons license, unless patients' clinical evaluation.The study has been funded by the follow- indicated otherwise in a credit line to the material. If material is not ing resources. National Institutes of Health of the United States of included in the article’s Creative Commons license and your intended America grant R01HD078641 (D.F.C, K.I.A) National Institutes of Health use is not permitted by statutory regulation or exceeds the permitted of the United States of America grant P50HD096723 (D.F.C.) National use, you will need to obtain permission directly from the copyright Health and Medical Research Council of Australia grant APP1120356 holder. To view a copy of this license, visit http://creativecommons.org/ (M.O., D.C., K.I.A., R.M., and J.A.V.) Spanish Ministry of Health Instituto licenses/by/4.0/. Carlos III-FIS grant FIS/FEDER- PI20/01562 (C.K., A.R.-E.). Estonian Research Council grants IUT34-12 and PRG1021 (M.L., M.P.) ReproUnion © The Author(s) 2022 and the Innovation Fund Denmark grant 14-2013-4 (K.A.) German Nature Communications | (2022)13:7953 17 Article https://doi.org/10.1038/s41467-022-35661-z Liina Nagirnaja 1, Alexandra M. Lopes2,3, Wu-Lin Charng4, Brian Miller1, Rytis Stakaitis 5,6,7, Ieva Golubickaite 5,6,8, Alexandra Stendahl 1, Tianpengcheng Luan9, Corinna Friedrich10, Eisa Mahyari1, Eloise Fadial 1, Laura Kasak 11, Katinka Vigh-Conrad1, Manon S. Oud 12, Miguel J. Xavier 13, Samuel R. Cheers9, Emma R. James14,15, Jingtao Guo 14, Timothy G. Jenkins14, Antoni Riera-Escamilla16,17, Alberto Barros2,18, Filipa Carvalho2,18, Susana Fernandes2,18, João Gonçalves19,20, Christina A. Gurnett4, Niels Jørgensen5,6, Davor Jezek 21, Emily S. Jungheim22, Sabine Kliesch23, Robert I. McLachlan24, Kenan R. Omurtag25, Adrian Pilatz26, Jay I. Sandlow27, James Smith28, Michael L. Eisenberg29, James M. Hotaling14, Keith A. Jarvi 30, Margus Punab31,32, Ewa Rajpert-De Meyts 5,6, Douglas T. Carrell14, Csilla Krausz33, Maris Laan 11, Moira K. O’Bryan 9,34, Peter N. Schlegel35, Frank Tüttelmann 10, Joris A. Veltman 13, Kristian Almstrup 5,6, Kenneth I. Aston14 & Donald F. Conrad 1 1 Division of Genetics, Oregon National Primate Research Center, Oregon Health & Science University, Beaverton, OR, USA. 2i3S - Instituto de Investigação e Inovação em Saúde, University of Porto, Porto, Portugal. 3IPATIMUP - Instituto de Patologia e Imunologia Molecular da Universidade do Porto, Porto, Portugal. 4 Department of Neurology, Washington University, St. Louis, MO, USA. 5Department of Growth and Reproduction, Copenhagen University Hospital - Rigshospitalet, Copenhagen, Denmark. 6International Center for Research and Research Training in Endocrine Disruption of Male Reproduction and Child Health (EDMaRC), Copenhagen University Hospital - Rigshospitalet, Copenhagen, Denmark. 7Laboratory of Molecular Neurooncology, Neuroscience Institute, Lithuanian University of Health Sciences, Kaunas, Lithuania. 8Department of Genetics and Molecular Medicine, Lithuanian University of Health Sciences, Kaunas, Lithuania. 9School of BioSciences, Faculty of Science, The University of Melbourne, Parkville, VIC, Australia. 10Institute of Reproductive Genetics, University of Münster, Münster, Germany. 11Institute of Biomedicine and Translational Medicine, University of Tartu, Tartu, Estonia. 12Department of Human Genetics, Radboud University Medical Centre, Nijmegen, Netherlands. 13Biosciences Institute, Faculty of Medical Sciences, Newcastle University, Newcastle- upon-Tyne, UK. 14Andrology and IVF Laboratory, Department of Surgery (Urology), University of Utah School of Medicine, Salt Lake City, UT, USA. 15 Department of Human Genetics, University of Utah School of Medicine, Salt Lake City, UT, USA. 16Andrology Department, Fundació Puigvert, Universitat Autònoma de Barcelona, Instituto de Investigaciones Biomédicas Sant Pau (IIB-Sant Pau), Barcelona, Catalonia, Spain. 17Molecular Biology Laboratory, Fundació Puigvert, Instituto de Investigaciones Biomédicas Sant Pau (IIB Sant Pau), Universitat Autònoma de Barcelona, Barcelona, Catalonia 08025, Spain. 18 Serviço de Genética, Departamento de Patologia, Faculdade de Medicina da Universidade do Porto, Porto, Portugal. 19Departamento de Genética Humana, Instituto Nacional de Saúde Dr. Ricardo Jorge, Lisboa, Portugal. 20Centre for Toxicogenomics and Human Health, Nova Medical School, Lisbon, Portugal. 21 Department of Histology and Embryology, University of Zagreb School of Medicine, Zagreb, Croatia. 22Department of Obstetrics and Gynecology at Northwestern University, Division of Reproductive Endocrinology, Chicago, IL, USA. 23Department of Clinical and Surgical Andrology, Centre of Reproductive Medicine and Andrology, University Hospital Münster, Münster, Germany. 24Hudson Institute of Medical Research and the Department of Obstetrics and Gynecology, Monash University, Clayton, VIC, Australia. 25Department of Obstetrics and Gynecology at Washington University, Division of Reproductive Endocrinology, St. Louis, MO, USA. 26Clinic for Urology, Pediatric Urology and Andrology, Justus Liebig University, Giessen, Germany. 27Department of Urology, Medical College of Wisconsin, Milwaukee, WI, USA. 28Department of Urology, University California San Francisco, San Francisco, CA, USA. 29 Department of Urology, Stanford University School of Medicine, Stanford, CA, USA. 30Division of Urology, Department of Surgery, Mount Sinai Hospital, University of Toronto, Toronto, ON, Canada. 31Andrology Center, Tartu University Hospital, Tartu, Estonia. 32Institute of Clinical Medicine, University of Tartu, Tartu, Estonia. 33Department of Experimental and Clinical Biomedical Sciences, University of Florence, Florence, Italy. 34School of Biological Sciences, Monash University, Clayton, VIC, Australia. 35Department of Urology, Weill Cornell Medicine, New York, NY, USA. e-mail: conradon@ohsu.edu Nature Communications | (2022)13:7953 18