Patients were recruited from Khyber Medical University, Pakistan (n = 24), University of Pittsburgh, United States (n = 110), and Institute of Human Genetics, Polish Academy of Sciences, Poland (n = 39). We analyzed WES data from a previously unpublished cohort of 924 sporadic cases of non-obstructive azoospermia sequenced by the GEMINI Consortium (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC9792524/). Informed consent was obtained from patients and participating family members after briefing about the research project. 1,033 patients within the overall SPGF cohort were diagnosed with azoospermia (no sperm in the ejaculate), 7 patients were diagnosed with cryptozoospermia (<0.1 million sperm/mL), 38 were diagnosed with severe oligozoospermia (0.1–5 million sperm/mL), and 19 patients were diagnosed with oligozoospermia (5–15 million sperm/mL) in accordance with American Urology Association and American Society for Reproductive Medicine guidelines (Schlegel et al., 2021). Patient evaluation also included serum reproductive hormone levels, karyotyping, ultrasound, and testicular biopsy (where appropriate). Individuals with a history of non-genetic causes for SPGF (e.g., trauma, surgery, or medication), obstructive azoospermia (e.g., CBAVD), abnormal sex chromosome evaluation on karyotype or Y chromosome microdeletions were excluded from the study.

DNA was isolated from peripheral blood and whole exome sequencing and analysis were performed as previously described (Hardy et al., 2021). Briefly, DNA was isolated from blood with the Gentra Puregene kit (Qiagen, United States). For samples prepared in Pittsburgh (Khyber Medical University and University of Pittsburgh), whole exome sequencing libraries were constructed with the SOPHiA Whole Exome Solution (SOPHiA GENETICS, Inc., United States) and sequenced on the Illumina NovoSeq 6,000 (Novogene, Sacramento, CA). Paired-end fragments were sequenced with target read length of 150 bp and average depth coverage of ∼110x per target interval. Raw data quality obtained was evaluated with FastQC software (Babraham Bioinformatics, Babraham Institute, UK). FASTQ files were aligned with the GRCh37/hg19 reference genome using SOPHiA Genetics DDM platform (SOPHiA GENETICS, Inc., Boston, MA, United States), utilizing proprietary algorithm with BaseSpace Burroughs-Wheeler Aligner (BWA) with the BWA-MEM algorithm. Among many polymorphisms with minor allele frequencies (MAF) greater than 1%, the following variants with MAF>15% were identified at the expected frequencies across our cohort: p. Leu1724Val (gnomAD MAF: 0.18297; SPGF MAF: 0.16044), p. Asn1698Ser (gnomAD MAF: 0.18309; SPGF MAF: 0.16089), p. Ile1422Val (gnomAD MAF: 0.25998; SPGF MAF: 0.23701), p. Cys491Arg (gnomAD MAF: 0.26515; SPGF MAF: 0.24157), and p. Thr2862 = (gnomAD MAF: 0.99269; SPGF MAF: 0.99407). Since the majority of the patients were recruited without ethnic identification, we are unable to reliably ethnically match our cohort population to the gnomAD reported frequencies, this may contribute to the slight discrepancy between the observed and expected MAFs of our variants in the SPGF cohort and the general population. We have estimated matched frequency of polymorphic variants based on general population structure to validate the power of our variant detection softwares.

Whole-genome sequencing (WGS) was performed on an Illumina HiSeq X (Illumina San Diego, CA) with minimum coverage of 30 × (100–120 Gb per sample) (Malcher et al., 2022). The SPGF patients from the Polish Academy of Sciences (n = 39) were sequenced with WGS methodology. Our WGS analysis did not identify new or significant deep intronic variants.

WES variants were annotated using the SOPHiA Genetics DDM (SOPHiA GENETICS, Inc., Boston, MA, United States) and Fabric Enterprise (Fabric Genomics, Oakland, CA) platforms. Variants were filtered through a list of 336 genes curated from previously reported gene candidate lists (Oud et al., 2019; Alhathal et al., 2020; Houston et al., 2021) with a moderate or above level of evidence. To infer haplotype for sporadic patients without familial segregation data, we utilized gnomAD’s phasing program, which makes use of the Expectation-Maximization algorithm to infer haplotypes of compound heterozygous variants (https://gnomad.broadinstitute.org/variant-cooccurrence) (Niu, 2004). WGS variants were annotated using Ensembl Variant Effect Predictor 30. Databases including 1,000 Genomes (NCBI browser), ESP6500 (National Heart, Lung, and Blood Institute), ExAC and gnomAD (Broad Institute) were used for identification of variant allele frequency in the general population. For evaluating potential clinical significance of variants, the following databases were used: HGMD (Stenson et al., 2014) (Qiagen) (http://www.hgmd.cf.ac.uk/ac/index.php), OMIM (Amberger et al., 2015) (Johns Hopkins University) (https://www.omim.org/), and ClinVar (NCBI) (https://www.ncbi.nlm.nih.gov/clinvar/). The MGI database (Blake et al., 2021) (Jackson Laboratory, United States) (http://www.informatics.jax.org/) was used for evaluation of gene variants from animal models. AceView (National Center for Biotechnology Information, NCBI) (https://www.ncbi.nlm.nih.gov/), GTEx Portal (Lonsdale et al., 2013) (Broad Institute) (https://gtexportal.org/home/), Human Protein Atlas (Uhlen et al., 2010) (https://www.proteinatlas.org/), and BioGPS (Wu et al., 2009) (Scripps Research Institute) (http://biogps.org/#goto=welcome) databases were used for evaluation of gene expression at the tissue level. Clustal Omega (Madeira et al., 2019) (EMBL-EBI) (https://www.ebi.ac.uk/Tools/msa/clustalo/) and PhyloP (Pollard et al., 2010) (Cornell University) (http://compgen.cshl.edu/phast/) were used for assessing conservation of variants. CADD (Kircher et al., 2014) (University of Washington, Hudson-Alpha Institute for Biotechnology and Berlin Institute of Health) (https://cadd.gs.washington.edu/), NNSplice (Berkley Drosophila Genome Project) (https://www.fruitfly.org/seq_tools/splice.html), PolyPhen-2 (Adzhubei et al., 2010) (Harvard University) (http://genetics.bwh.harvard.edu/pph2/), SIFT (Kumar, Henikoff, and Ng, 2009) (J. Craig Venter Institute) (https://sift.bii.a-star.edu.sg/), MaxEntScan (Yeo and Burge, 2004) (Massachusetts Institute of Technology) (http://hollywood.mit.edu/burgelab/maxent/Xmaxentscan_scoreseq.html), MutationTaster (Schwarz et al., 2010) (Charité) (https://www.mutationtaster.org/), and VVP (Flygare et al., 2018) were used for prediction of consequence of amino acid changes on protein function. The R programming environment was used to identify WES variants with causative effect for non-obstructive azoospermia within the GEMINI cohort.

Variants obtained by WES were confirmed by Sanger sequencing using BigDye sequencing kit (ThermoFisher Scientific). Primer3+ software was used for designing primers. Sequencher (GeneCodes) was used for analysis of Sanger sequencing results.

Previously generated, processed merged human digital expression matrix data file from Drop-seq experiments on human testicular cells from 4 adult males was obtained from Gene Expression Omnibus (GEO: GSE142585) as previously described (Shami et al., 2020). The data was generated through global analysis of a total of 35,941 transcripts from 13,597 cells. The following known RNA markers for human testicular cell types were used for expression profiling: spermatogonia (GFRA1, HORMAD1, ID4, ITGA6, LY6K, STRA8, SYCP2, UCHL1, and UTF1), spermatocytes (PIWIL1 and SYCP3), spermatids (ACRV1, PRM1, TNP1, and TSSK6), Leydig cells (IFG1/2 and STAR), endothelial cells (NOSTRIN and VWF), testicular macrophages (CD52, CD163, LYZ, and TYROBP), pericytes (ADIRF, MCAM, PDGFRB, and STEAP4) or myoid cells (ACTA2 and MYH11). The normalized PCA and UMAP images developed from this data collection process were scaled and reported to support the findings in this study (Figure 1A).

The sequence of TEX15 (accession #Q9BXT5, A0A1W2PS94 isoform producing a 3,176 amino acid polypeptide) was obtained from Uniprot (UniProt Consortium, 2021). Three dimensional structural models were generated using AlphaFold (Jumper et al., 2021) and Phyre2 (Kelley et al., 2015). Because the length of TEX15 exceeds the size AlphaFold can query, we folded shorter regions for each prediction. Individual domain boundaries were initially informed by pfam (Mistry et al., 2021) and then manually curated and iteratively remodeled based on sequence alignments and structural modeling from earlier rounds of prediction. Likewise, residue sequences used to predict interdomain interactions were manually dissected and modeled with AlphaFold (Jumper et al., 2021). Structural alignment was performed in COOT (Emsley and Cowtan, 2004) and structural model figures made in PyMOL (Schrödinger, LLC).

A 20-year-old, married male of Pashtun ethnicity from the Khyber Pakhtunkhwa (KPK) region of Pakistan presented with unexplained infertility. The parents were first cousins and had 9 children (Figure 2). The proband’s semen analysis indicated cryptoozoospermia with normal semen volume, and viscosity (Table 1). Secondary factors affecting reproduction like ejaculatory defects, immunological irregularities, congenital disorders or environmental exposures were ruled out. Physical examination including weight, height and secondary sexual characteristics were normal. Hormonal analysis for follicle stimulating hormone, luteinizing hormone, prolactin, and testosterone all showed normal values (Table 1). T-cell cytogenetic analysis showed normal 46, XY male karyotype (Table 1). Y chromosome microdeletions as well as abnormal vas deferens associated mutations were not detected (Table 1).

DNA from the proband, both parents, and one unaffected brother were subjected to WES analysis in order to identify genetic causes of SPGF. Candidates were first screened against a panel of 336 genes with evidence of association with male infertility in humans and/or mice (Oud et al., 2019; Alhathal et al., 2020; Houston et al., 2021). Among identified genes from this list, a homozygous variant in TEX14 was ruled out due to its presence in the unaffected brother. Based on likely modes of inheritance, potential candidate variants for the family found in other relevant genes were a compound heterozygous mutation in AKAP12, homozygous mutation in CHAF1B, homozygous mutation in KIF12, homozygous mutation in RIMBP3, and hemizygous mutation in ERCC6L. These variants were not considered further due to overall lack of sufficient criteria (e.g., MAF>1% polymorphic in one of the ethnic subpopulations) leaving TEX15 as the top candidate. Following thorough analysis, the homozygous TEX15 missense variant c.6835G>A, p. Ala2279Thr identified in the proband (see Case 6 in Table 2), was found to co-segregate with cryptozoospermia as it was not found in the unaffected brother and is heterozygous in both parents, indicating autosomal recessive inheritance (Supplementary Figure S1).

Single cell RNAseq data (GEO: GSE142585) analysis from 4 adult men confirmed the expression of TEX15 RNA in all cell types assessed, except pericytes, in adult human testes (Figure 1B) (Shami et al., 2020). TEX15 RNA shows high expression in male germ cells including both undifferentiated and differentiating spermatogonia, spermatocytes, and spermatids. Loss of function mutations are associated with spermatogenic failure 25 (OMIM#617960 - SPERMATOGENIC FAILURE 25; SPGF25). The homozygous c.6835G>A, p. Ala2279Thr variant found in the proband affects an uncharacterized region of the 3,176 amino acid protein (ENST00000638951.1, ENSP00000492713.1) (Figure 3). Three-dimensional modeling suggests that this residue is located in a predicted alpha-helix whose position is important for interhelical packing within its TEX15 repeat (Figure 4A). Substitution of this residue with threonine results in steric clashes with surrounding amino acids; in one putative threonine rotamer, the methyl group in the side chain would clash with the side chain of lysine 2,314 (Figure 4B). In the other rotamer, the hydroxyl group of the side chain would approach the backbone of phenylalanine 2,276 too closely (Figure 4B). In each of these scenarios, the steric clash would be predicted to be propagated through the tertiary structure, destabilizing this TEX15 repeat and its interaction with the neighboring repeat. Additionally, this position is highly conserved among mice, chimpanzees, rhesus macaques, and humans (Supplementary Figure S2). These predicted steric hindrances and the relevance of position conservation are supported by the unanimously deleterious predictions in all in silico variant assessments performed (Table 2).

Follow-up analysis of WES and WGS results from the remaining 1,096 unrelated male patients (see Materials and Methods) with unexplained SPGF identified 6 cases with a total of 13 unique, potentially causative, likely in trans compound heterozygous/unphased deleterious TEX15 variants (Figure 3; Table 2). A c.700–6T>C splice region variant of uncertain significance (VUS) in the N-terminal “domain of unknown function” along with a likely pathogenic c.6860_6861del, p. Arg2287Asnfs*22 frameshift resulting in a premature stop prior to the second TEX15 domain were identified in azoospermic Case 1; the former of which may affect a predicted ADP ribosylation site as observed with 3D modeling (Supplementary Figure S3).

Case 2, diagnosed with azoospermia due to late meiotic arrest, presented with 3 TEX15 missense VUSes: c.1636A>T, p. Ile546Phe, c.1786A>G, p. Ile596Val, c.8820T>C, p. Ser2940Pro. The first variant, p. Ile546Phe, is moderately conserved (e.g., in some species instead of an isoleucine, it is a methionine or serine); none of the aligned sequences have an aromatic amino acid at this position (Supplementary Figure S2). For the second variant, p. Ile596Val, the conserved position is either isoleucine, valine, or leucine in the alignment (Supplementary Figure S2), so the small difference in the human sequence is a slight size difference, with valine being smaller. The third variant, p. Ser2940Pro, could have a strong effect on the structure of the disordered region, given the nature of proline’s side chain. The additional proline residue could lead to decreased flexibility and/or increased steric clash in the backbone. Given the similarity of the MAFs for the p. Ile546Phe and p. Ile596Val variants (.0001 and .00013 respectively), it is likely that these variants appear as in cis compound heterozygotes (Table 2). The novel p. Ser2940Pro variant is currently not reported in gnomAD, and as a result, we could not use gnomAD’s phasing algorithm (See Material and Methods) to predict haplotype for this variant in relation to the other two variants. However, we believe that it is possible that the p. Ile546Phe and p. Ile596Val variants are in trans compound heterozygotes with the PhyloP/MutationTaster software-predicted deleterious p. Ser2940Pro variant (Table 2), which may have either contributing or pathogenic associations with the azoospermic phenotype observed in the patient. This belief is supported by the notion that within Intrinsically Disordered Regions (IDRs), any position may be important for intra- or inter-molecular interactions. These computational tools are particularly suggestive as they evaluate not only amino acid substitution effect, but more specifically the effect at the relatively conserved position as mentioned above in Materials and Methods. Notably in Case 2, two believed in trans compound heterozygous (assumed on their significantly different MAF frequencies, see Materials and Methods for inference validation) Meiotic Double-Stranded Break Formation Protein 1 (MEI1) variants, consisting of one inframe deletion and one predicted highly deleterious missense variant, were also seen (Supplementary Table S1). MEI1 variants (including frameshifts, non-senses, and missenses) have been previously associated with non-obstructive azoospermia in male patients (Ben Khelifa et al., 2018; Nguyen et al., 2018; Malcher et al., 2022). Given that Case 2 is a sporadic patient without familial segregation, we are unable to confidently rule out the MEI1 variants as potentially causal; however, without further information, they also do not have sufficient evidence for pathogenic classification than the identified TEX15 variants, and so we retained them in Case 2 in analysis.

Case 3, diagnosed with azoospermia, had a missense c.2561G>T, p. Trp854Leu substitution in addition to a c.4131C>A, p. Ser1377Arg, both of uncertain significance. The first variant, p. Trp854Leu, is a large change from a hydrophobic aromatic amino acid to a hydrophobic aliphatic amino acid. While both amino acids are classified as hydrophobic, due to the additional NH group in the tryptophan side chain, it is notably less hydrophobic than leucine. The alignment suggests that the position is not conserved for tryptophan specifically, as seen in Supplementary Figure S2, the other amino acids in that position are cystine, histidine, and asparagine, which are all polar amino acids; this is further supported by the PhyloP benign prediction (Table 2). However, again as suggested previously, due to the intrinsically disordered nature of the protein, specific position conservation is not conclusive of variant effects. The second variant, p. Ser1377Arg, is in a mostly conserved region, with only the Western clawed frog (X. tropicalis) and the Zebrafish (D. rerio) having alternative amino acids at the position: glutamine and valine respectively (Supplementary Figure S2). The change is from the tiny uncharged polar amino acid, serine, to the much larger and positively charged arginine, which could cause problems in the structure of the disordered region.

A likely pathogenic c.3970C>T, p. Arg1324* premature stop upstream of both TEX15 domains and a c.4957G>A, p. Val1653Ile VUS were found in Case 4, diagnosed with azoospermia as a result of early meiotic arrest. The missense variant p. Val1653Ile is in a partly conserved position; however, the position is conserved as either valine or isoleucine, which are both aliphatic amino acids, however, valine is slightly smaller (Supplementary Figure S2). From the alignment, the following amino acid is leucine for the species that have valine at the position of interest, whereas in species with isoleucine at this position, the subsequent leucine is replaced by a valine, except in M. musculus, where it is a histidine (Supplementary Figure S2). Given the fact that the region has unique differences in M. musculus compared to the others, it could mean the structure in most of the aligned species requires a small amino acid near that position.

Azoospermic Case 5 had a c.5602G>C, p. Glu1868Gln VUS as well as a pathogenic c.6886_6889del, p. Ser2296Lysfs*11 premature termination prior to the second TEX15 domain. For the missense variant, p. Glu1868Gln, the alignment shows that the position is highly conserved for glutamic acid (Supplementary Figure S2), suggesting that the glutamic acid is important for the proper function and/or structure of the protein.

Additionally, azoospermic Case 7 had a pathogenic c.7426C>T, p. Arg2476* premature stop in the second TEX15 domain as well as an intronic c.8276-17G>A VUS.

In addition to these compound heterozygous variants, single heterozygous variants were observed in several patients (Supplementary Table S2). Interestingly, the oligozoospermic patient Case 18 had in cis compound heterozygous variants that included a pathogenic frameshift c.6860_6861del, p. Arg2287Asnfs*22 and a likely pathogenic frameshift c.6886_6889del, p. Ser2296Lysfs*11 in the uncharacterized region between the 2 TEX15 domains, both of which resulted in premature stops (Supplementary Table S2; Supplementary Figure S4). While 14 of the remaining 1,090 patients (excluding the 7 described above) observed had single or in cis compound heterozygous variants, these variants were classified as carrier status and overall are not thought to be pathogenic in these patients, considering TEX15-associated disorders are known to exhibit autosomal recessive inheritance.