Approximately 7% of the male population worldwide suffer from infertility, and one of the main causes is azoospermia, which is a condition contributing to male infertility with the highest frequency of known genetic factors (about 25%) (Krausz and Riera-Escamilla 2018; Rodrigues, Polisseni et al., 2020). Given that nearly 50% of azoospermia cases are estimated to be associated with genetic defects, the etiology of azoospermia is still very elusive (Rodrigues, Polisseni et al., 2020). In terms of pathobiological mechanism, azoospermia usually is categorized into two groups: non-obstructive azoospermia (NOA) which is defined as no spermatozoa due to failure of spermatogenesis and obstructive azoospermia which occurs due to an obstruction of the seminal ducts and is characterized by the absence of spermatozoa in the ejaculate despite the normal spermatogenesis.

NOA is the most severe form of azoospermia; currently, no treatment can restore spermatogenesis in the majority of patients; however, some patients can benefit from treatment with assisted reproductive technologies. NOA is a heterogenous condition with variable histopathological phenotypes: the Sertoli cell-only syndrome (SCOS) is characterized by a complete absence of germ cells, while spermatogenic arrest is a condition originating from the physiological spermatogenesis defect and is known as the interruption of germinal cell formation of specific cellular type.

Approximately 20% of NOA cases were caused by detectable structural anomalies including chromosomal abnormalities such as the sex-chromosome aneuploidy, the trans-locational aberrations, and microdeletions of azoospermia factors (AZFs) (Krausz and Riera-Escamilla 2018). In addition, monogenic variations have been identified to be a causative for or associated with NOA; the attested genes include X chromosome linked genes TEX11 (Yatsenko, Georgiadis et al., 2015), FAM47C (Chertman, Arora et al., 2019), MAGEB4 (Okutman, Muller et al., 2017), and AR (Mou and Gui 2016); and autosomal genes include PIWIL1 (Gou, Kang et al., 2017), BRCA2 (Zhoucun, Zhang et al., 2006), BOLL and DAZL (Zhang, Xue et al., 2016), TEX15 (Okutman, Muller et al., 2015), SYCP3 (Miyamoto, Hasuike et al., 2003), NR5A1 (Safari, Zare-Abdollahi et al., 2014), SOHLH1 (Choi, Jeon et al., 2010), WT1 (Wang, Li et al., 2013; Xu, Jiang et al., 2017), TAF4B and ZMYND15 (Ayhan, Balkan et al., 2014), FANCM (Kasak, Punab et al., 2018), SPINK2 (Kherraf, Christou-Kent et al., 2017), TEX14 (Gershoni, Hauser et al., 2017), DNAH6 (Gershoni, Hauser et al., 2017), MEIOB (Gershoni, Hauser et al., 2017), RNF212 and STAG3 (Riera-Escamilla, Enguita-Marruedo et al., 2019), XRCC2 (Zhang, Li et al., 2019), TDRD7 and TDRD9 (Arafat, Har-Vardi et al., 2017; Tan, Tu et al., 2019), DMC1 (He, Tu et al., 2018), SYCE1 (Maor-Sagie, Cinnamon et al., 2015), NPAS2 (Ramasamy, Bakircioglu et al., 2015), MCM8 (Tenenbaum-Rakover, Weinberg-Shukron et al., 2015), MEI1 (Ben Khelifa, Ghieh et al., 2018), and STX2 (Nakamura, Kobori et al., 2018). Despite these efforts, the etiology remains unknown in an enormously large part of NOA patients (Krausz and Riera-Escamilla 2018). Considering that a surprisingly large proportion of all genes (roughly 1 in 25 of all mammalian genes) are specifically expressed in the testis, hundreds of genes may be involved in the spermatogenesis process (Schultz, Hamra et al., 2003; Bell, Mello et al., 2020). Therefore, the physiological pathway complexity may contribute to a large number of unresolved patients.

However, current tools for understanding the unidentified genetic etiology of NOA are very limited due to the lack of large pedigrees with infertility; therefore, the discovery of novel genetic factors for unexplained NOA is still the main challenge. In this study we performed whole exome sequencing (WES), one of the most promising and powerful approach for exploring the pathogenic factors of inherent diseases, on 133 unrelated patients and 343 independent controls for rare variant association study (RVAS) (Lee, Abecasis et al., 2014; Artomov, Stratigos et al., 2017; Locke, Steinberg et al., 2019), which was a kind of analysis derived from genome-wide association study (GWAS); meanwhile, 11 trios were analyzed in a maternal inherited or a de novo pattern to delineate the causative genes. Protein-protein interaction (PPI) network functional enrichment analysis was also applied to provide potential evidence to further address the abovementioned questions.

The workflow diagram of this study is shown in Figure 1, a total of 1739 patients were enrolled, amidst which 299 cases (about 17%) with abnormal karyotyping results (such as sex-autosome translocation, and Y chromosome polymorphisms (Yqh–or Yqh+)) were excluded, followed by 281 cases (about 16%) exclusion for harboring the AZF region (AZFa, AZFb, or AZFc) microdeletions; then 1,023 cases who did not consent to accept further sequencing test were not included. Furthermore, a custom panel for detecting the known genes of azoospermia (Supplementary Table S2) was developed, from which three cases with positive results (two with TEX11 mutation, one with AR mutation) were excluded. Eventually, 133 patients and 343 controls were subjected to WES; linkage analysis for 11 trios was also performed in this study.

Among the 133 patients, WES identified 7 LOFs in six known genes (Table 1). Since FAM47C was identified to be highly associated with NOA previously, the X chromosome-linked (XL) FAM47C mutation (p.Val118Leufs*40, frameshift mutation similar with previous literature (Chertman, Arora et al., 2019)) from patient 58 in our study was also identified as the potential genetic etiology accordingly; the BRCA2 mutation in patient 58 (p.Ser 1955*) and patient 82 (p.Thr2197_Glu2198delinsLys) may be associated with NOA. The causative gene SYCP3 mutation in patient 88 (c.-13-2A > C) may be pathogenic (functional evidence needed) according to an autosomal dominant (AD) inheritance pattern. In autosomal recessive (AR) inheritance, the other genes are uncertainly causative concerning the individuals owing to another pathogenic allele not determined. Pathogenic CNVs were not identified among the 133 patients under WES.

Whole exome-wide RVAS results were summarized in Supplementary Table S3, including the p-values and the OR values, the normalized expression (NE) of genes in testis tissue from the Human Protein Atlas database were also provided. NE was calculated according to the transcript profiling based on a combination of three transcriptomics datasets HPA, GTEx, and FANTOM5 (genes with no NE data were ruled out in this study). We quantified genes with the OR value >1 or with OR value not estimable due to the absence of the rare variations among the controls, meanwhile we employed a flexible cutoff of p < 0.05 due to the statistical power limitation of the relatively small sample size (133 patients) and the low proportion of rare variations among the entire samples (a lot of rare variations absent in controls). As a result, a total of 648 genes (65 genes with p < 0.01) were identified to be associated with NOA (three of which were previously reported). The most significant genes reported were DHRS4, WARS1, PICK1, RRBP1, and ENTPD2 (red), whereas the other 179 significant genes were shown in blue; meanwhile, three known genes BRCA2, SYCP3, and TDRD7 (yellow) were also identified by the RVAS test in our study (Figure 2).The 187 genes identified by RVAS (red, blue, and yellow) and 28 previously reported proteins (green) were analyzed for the PPI network as shown in Figure 3. Interestingly, the previously reported proteins, except for only MAGEB4 and FAM47C, tended to be involved in the meiosis-related network, which is centered around PIWIL1, TEX11, TEX15, DAZL, and TEX14, etc. The majority of RVAS identified proteins (117), such as WARS, RRBP1, DHRS4 (red), and PIWIL2, DDX39B, SLX4, HNRNPC (blue), were involved and enriched in the meiosis-related PPI network (FDR value for the male gamete generation pathway was 9.56 × 10^–31), thereby betraying the functional potential to be strongly associated with NOA; whereas 70 proteins, like PICK1 and ENTPD2 (red), were not interacted with the PPI network, attenuating the possibility that these genes play major functional roles in the development of NOA.

For association study, OR values were not estimable for most of the 187 genes (Supplementary Table S3) due to the absence of rare variations among the controls, we summarized a short list of genes with OR values to evaluate the genetic risk (Figure 4). All the genes with p values lower than 0.05, with the testis expression (NE > 16), and with the OR values higher than 1 were included. Consistent with the previous literature, BRCA2 (OR, 7.48; 95% CI, 1.37–40.89; p value, 8.79 × 10^–3) was identified to be associated with NOA in our study. The following genes with high OR values were also involved in the meiosis-related PPI network (Figure 3), including ENTPD2 (OR, 23.34; 95% CI, 2.78–195.61; p value, 8.79 × 10^–5) which had the highest OR; while RRBP1 (OR, 15.32; 95% CI, 3.23–72.57; p value, 1.86 × 10^–5), WARS (OR, 8.83; 95% CI, 3.29–23.71; p value, 2.71 × 10^–7), DHRS4 (OR, 7.54; 95% CI, 3.46–16.42; p value, 4.61 × 10^–8), NOL7 (OR, 15.32; 95% CI, 1.70–138.23; p value, 2.28 × 10^–3), HLA-A (OR, 15.03; 95% CI, 1.67–134.84; p value, 9.13 × 10^–4), EBNA1BP2 (OR, 11.28; 95% CI, 1.17–108.90; p value, 1.06 × 10^–2), CC2D1B (OR, 11.22; 95% CI, 1.17–108.12; p value, 1.07 × 10^–2), DDX39B (OR, 7.64; 95% CI, 2.28–25.58; p value, 4.20 × 10^–4), and COPB1 (OR, 7.64; 95% CI, 1.38–42.19; p value, 8.55 × 10^–3) were shown.

In total, 11 trio families were recruited as shown in Figure 1, and their data were analyzed and summarized in Table 2. For linkage analysis, according to a maternally inherited pattern or the de novo mutation mechanism to which the gender-specific disease NOA condition was supposed to be genetically attributed, the paternal inherited variations were excluded; meanwhile, the probands’ LOFs were preserved only if the LOFs were not found in the 343 male controls. Sanger sequencing was performed for validating the variations. To facilitate the screening for the mutations, the p-value and OR value of RVAS for each gene were also rendered. Eventually, six genes (in bold and double-underlined in Table 2), ACTL8 and TRA2B for patient 62, IDE for patient 117, FIBP for patient 124, NUP37 for patient 130, and PIGT for patient 131 were regarded as the top candidates of associated genes (p < 0.05; OR values not estimable due to the absence of variations among the controls) based on the RVAS test (Figure 2) and their involvement in the meiosis-related PPI network (Figure 3).

These six novel genes potentially accounted for a fraction (3.76%, 5 out of 133 patients) of the heritability of NOA in question. Combined with the three known genes (yellow), the six novel genes potentially accounted for an overall 6.77% (9 out of 133 patients) heritability of NOA in this study. In addition, the 19 following genes (p < 0.05; OR > 1 or values not estimable; NE values not less than 1) were also regarded as expansive candidate genes for each proband, including NDUFAF1 for patient 38, TMEM63A and SEMA3F for patient 52, BICC1, SPRR2F, and GSTM3 for patient 60, ABCB8, PRRG4, RPUSD4, and ZNF221 for patient 62, SELP and GNPTAB for patient 98, DDX11 and POLD2 for patient 117, LRRFIP1 and LMF1 for patient 124, ULK4 for patient 130, ACSF2 and PPP1R3C for patient 131.

Whole exome-based RVAS analysis identified a series of genes associated with NOA in this study. Especially, the most significant proteins, WARS1 (tryptophanyl-tRNA synthetase), DHRS4 (dehydrogenase/reductase 4), and RRBP1 (ribosome binding protein 1), were also involved in the meiosis related pathway, indicating they are associated with NOA. Meanwhile, three previously reported genes, TDRD7, SYCP3, and BRCA2, were also identified by RVAS and involved in the meiosis-related pathway, which partially supported the efficacy of the RVAS test flow in this study.

Combining RVAS analysis and trios linkage analysis showed the overlapped genes, ACTL8 (encodes actin-like 8 protein, restricted expression toward testis), TRA2B (encodes transformer 2 beta homolog protein), IDE (encodes insulin-degrading enzyme), FIBP (encodes FGF1 intracellular binding protein), NUP37 (encodes nucleoporin 37, a constituent of the nuclear pore complexes required for mitosis and probably for meiosis) (Loiodice, Alves et al., 2004), and PIGT (encodes phosphatidylinositol glycan anchor biosynthesis class T protein), were associated with NOA, as well as involved in the meiosis-related PPI network. Especially, the gene NUP37 was previously identified as indispensable for mitosis, we speculate that NUP37 (p.Trp80*, a frameshift mutation found in one out of 133 patients while absent in controls) may also play important roles during meiosis and spermatogenesis. In addition, the expression of ACTL8 gene is restricted in testis, the actin-like 8 protein ACTL8 (p.Val254Aspfs*4, a frameshift mutation found in one out of 133 patients, while being absent in controls) is assumed to execute the functions during spermatogenesis and play a role in NOA.

The other genes identified by RVAS or parental-proband trios analysis, especially those involved in the meiosis-related PPI network like PIWIL2, HNRNPC, and DDX39B may also play a role in NOA. It has not escaped our notice that PIWIL1 (encodes Human P-element-induced wimpy testis 1 protein, a paralog of human PIWIL2) was a previously identified gene causative for NOA; Human P-element-induced wimpy testis (PIWI) proteins act as protectors of germline, and are expressed mainly in the germline cells (Gou, Kang et al., 2017), the phenotype of Piwil2-deficient (knockout) male mice exhibited azoospermia with complete spermatogenic arrest according to the Mouse Genome Database, hence PIWIL2 (p.Gln11Serfs*78, a frameshift mutation found in one, P129, out of 133 patients while absent in controls) may play a similar role during spermatogenesis and in NOA.

Despite significant advances have been achieved in understanding the etiology of NOA, this study has certain limitations. First, prospectively more functional evidence like the cellular or physiological experiments should be investigated to address this question better. Second, the limitation lies in that non-LOFs like the missense mutations and small insertions or deletions which may be pathogenic were not considered. Finally, the false-positive rate of the de novo mutations (Table 2) is relatively high, which is in part attributable to sequencing errors; therefore, we suggest that de novo mutations should be treated with caution in clinical practice.

We demonstrate that RVAS identified a pool of genes, especially the most significant novel genes WARS, DHRS4, and RRBP1 which were also involved in the meiosis-related PPI network, indicating they are associated with NOA. In addition, six novel genes, ACTL8, TRA2B, IDE, FIBP, NUP37, and PIGT were revealed to be associated with NOA in both the RVAS and trios linkage analysis, as well as involved in the meiosis-related pathway; we conclude that these six novel genetic risk factors potentially account for a fraction (3.76%, 5 out of 133 patients) of the heritability of NOA in question, together with the three previously reported genes would potentially account for 6.77% (9 out of 133 patients) heritability of NOA in question in this study.

This study not only sheds light on the underlying pathological mechanism of NOA, but also offers valuable insight into the genomic landscape of NOA to constitute a potential basis for more efficient diagnosis yield in the future clinical application, and might have broad implications on men’s health.