Identification of 22 susceptibility loci associated with testicular germ cell tumors

Testicular germ cell tumors (TGCT) are the most common tumor in young white men and have a high heritability. In this study, the international Testicular Cancer Consortium assemble 10,156 and 179,683 men with and without TGCT, respectively, for a genome-wide association study. This meta-analysis identifies 22 TGCT susceptibility loci, bringing the total to 78, which account for 44% of disease heritability. Men with a polygenic risk score (PRS) in the 95th percentile have a 6.8-fold increased risk of TGCT compared to men with median scores. Among men with independent TGCT risk factors such as cryptorchidism, the PRS may guide screening decisions with the goal of reducing treatment-related complications causing long-term morbidity in survivors. These findings emphasize the interconnected nature of two known pathways that promote TGCT susceptibility: male germ cell development within its somatic niche and regulation of chromosomal division and structure, and implicate an additional biological pathway, mRNA translation.

TGCTs are the most common cancers in young men of European ancestry, and incidence of TGCT has doubled over the past 20 years1,2. Family history and cryptorchidism are the strongest known risk factors3,4,5, but no robust environmental risk factors have been identified1. Despite the high heritability of TGCT, estimated at 37–49%6,7, CHEK2 is the only moderate penetrance gene in which pathogenic variants have been associated with risk of TGCT8.

In contrast, genome-wide association studies (GWAS) have succeeded in identifying common variation associated with TGCT susceptibility9,10,11,12,13,14,15,16,17,18,19,20,21. Most risk variants map to loci containing genes encoding proteins implicated in critical pathways for male germ cell development, chromosomal segregation, sex determination, and DNA maintenance. Biologically these findings complement the current understanding of disease pathogenesis involving in utero transformation of fetal germ cells into germ cell neoplasia in situ (GCNIS), the common precursor of TGCT22,23.

To gain further insight into the genetic underpinnings of TGCT, the Testicular Cancer Consortium (TECAC) present results from a large meta-analysis of 10,156 men with TGCT and 179,683 men without TGCT that combined summary data from numerous extant TGCT GWAS and de novo genotyping from men with and without TGCT. We identify 22 independent loci for TGCT (P < 5 × 10−8), many of which map to genes that encode proteins in pathways related to male germ cell development, sex determination and chromosomal segregation, as well as mRNA translation. Polygenic risk score (PRS) analysis of all 78 identified risk loci to date reveals a 6.8-fold increase in TGCT risk for men in the top 5% of PRS score compared to those at the median.

Our meta-analysis incorporated estimates from our published TGCT analysis9, genotyping data from deCODE genetics24 and the UK Biobank25, and summary statistics from genotypes collected from 14 studies collaborating as part of the Testicular Cancer Consortium (TECAC) (Supplementary Tables 1, 2; Supplementary Methods). Initial findings were extended by incorporating results from targeted genotyping of 1039 men with TGCT and 1398 men without TGCT (Supplementary Tables 3, 4).

GWAS meta-analysis of TGCT

Our final meta-analysis identified 22 independent susceptibility loci for TGCT (P < 5 × 10−8) (Table 1, Fig. 1, Supplementary Fig. 1), including four independent signals at previously identified genetic regions (Supplementary Table 5; Supplementary Fig. 2) and four loci on the previously disregarded X chromosome (Supplementary Data 1). The Q–Q plot (Supplementary Fig. 3) and estimated genomic inflation factor (λ = 1.03) suggested minimal systematic bias. Only three signals (rs9987332, rs8104804, and rs4898474) showed effect heterogeneity (I2 > 50). Forty-four of the 56 previously identified TGCT susceptibility loci9,10,11,12,13,14,15,16,17,18,19,20,21 replicated at P ≤ 5 × 10−8 (Supplementary Data 2; Supplementary Data 3). Possible reasons for not replicating all known loci include differences in underlying population substructure, prior overestimation of genetic effect size, effect size heterogeneity, and low r2 between the current and previously published loci (Supplementary Data 3). Multiple independent signals were observed at BAK1 (2), TKTL1 (2), TERT (3), DMRT1 (4), and the 19p11-p12 (6) region (Table 1; Supplementary Data 2), a complex region containing multiple KRAB-zinc finger proteins (Supplementary Fig. 4). Minimal overlap is present between the 66 novel and replicated independent loci for TGCT and susceptibility loci identified in GWAS of other cancers (Supplementary Data 4). Only four (6%) loci were associated with risk of another cancer type, each with consistency in direction of effect: BCL2L11 (rs6708784–rs1439287, r2 = 0.93) with chronic lymphocytic leukemia, TERT (rs2735940) with colorectal cancer, HEATR3 (rs2160570–rs10852606, r2 = 0.99) with glioblastoma, and HNF1B (rs11263762–rs12601991, r2 = 1.00) with cancer (pleiotropy).

Novel markers identified in the current meta-analysis are shown as blue squares () with lowercase letters corresponding to column 1 of Table 1. Susceptibility markers identified in previous studies that surpassed genome-wide significance (P ≤ 1 × 10−8) in the current meta-analysis are shown as green circles () with numbers corresponding to column 1 of Supplementary Data 2. Susceptibility markers identified in previous studies that failed to attain genome-wide significance (P > 1 × 10−8) in the current meta-analysis are shown as red diamonds () with numbers corresponding to column 1 of Supplementary Data 2. a Markers are plotted against a full range y-axis that incorporates rs4474514 at KITLG (P = 1.42 × 10−154). b Markers are plotted against a partial range y-axis capped at P= 1.42 × 10−40 to allow for better visualization and discrimination of most associations.

Full size image

Stratified analyses by histology, family history, or cryptorchidism (Supplementary Table 6) did not identify subgroup associations. All 22 susceptibility signals displayed marked differences in minor allele frequency between men of European and African ancestry (Supplementary Table 7) likely explaining some of the observed racial differences in TGCT risk. The 22 identified loci explain 7.0% of father-to-son heritability and 4.7% of heritability among siblings, increasing the overall heritability estimates to 44.0% and 29.1%, respectively.

To generate a polygenic risk score (PRS) for TGCT, we modeled all 78 identified TGCT susceptibility markers, including those that did and did not achieve genome-wide significance in the current study. We found that men in the 95th percentile of PRS had a 6.8-fold increased disease risk (3.4% lifetime risk) compared to men with median scores (Fig. 2). This model identifies men with TGCT with 78.1% accuracy.

Polygenic risk scores (PRS) were calculated for independent samples from n = 5602 men with TGCT and 5006 men without disease from a model incorporating the 22 novel and 56 previously identified markers and effect size estimates from the meta-analysis. Odds ratios are relative to the median risk, composed of subjects in 45–55th percentile of PRS. Men in the top 95th percentile had a 6.8-fold increase (odds ratio (OR) = 6.75, 95% confidence interval (CI) 4.92–9.26; P = 2.84 × 10−32) in risk of developing TGCT compared to men at the 45–55th percentile. Dashed line indicates OR = 1; error bars represent 95% CI.

Full size image

Assessment of credible risk variants (CRV)

We defined a credible risk variant (CRV) as a SNP in strong LD (r2 ≥ 0.8) with any of the 66 novel or replicated signals to determine if among the set of 4755 CRVs there are potential functional variants that influence function or expression of the target gene (Supplementary Table 8 and Supplementary Data 5). A total of 108 unique genes were in regions demarcated by the CRVs on the autosomes and X chromosome. Most GWAS have implicated noncoding variation that work through gene regulation (e.g., enhancers, promoters), but coding variation can also influence target gene function. Seventy-three (1.5%) CRVs were located in coding regions; 34 (0.7%) were synonymous and 39 (0.8%) were missense variants (Supplementary Data 6). None were predicted to be pathogenic using REVEL and VEST426,27. Seven (0.1%) CRVs were annotated at a splice site; but only one, rs1060604 at PMF1 was predicted to influence splicing28. These results align with those from other GWAS and support that most susceptibility functional variants affect the regulation of target genes rather than directly altering gene function.

Inference of autosomal genes associated with TGCT

To identify highly and moderately likely target genes on autosomes, we assessed the gene regions delimited by 4484 CRVs corresponding to 61 top signals (Supplementary Data 5). The total number of target genes evaluated was 108, corresponding to 101 unique genes. As further detailed below, we evaluated (i) the number of genes in the region, (ii) location of the most significantly associated signal, (iii) results from colocalization eQTL analysis29, (iv) gene expression in fetal germ cells30, and (v) results from promoter Capture-C analysis of the TGCT cell line NT2-D1 (NTERA2)31,32 evaluated in conjunction with data from ATAC-seq (Fig. 3). The number of genes in each region ranged from one to eight. For 46 (75%) signals, the gene region included only one or two genes; and for six (10%) signals the gene region encompassed no genes (Supplementary Data 5). Forty-three (70%) of the top signals were in an exon, an intron, or within 10 Kb of a start site (Supplementary Data 5). The colocalization analysis found an eQTL in at least two (non-testis) tissues for 23 (21%) genes, and in testis tissue for 4 (4%) genes (Supplementary Data 5, 7).

Highly and moderately likely target genes were determined by evaluating information derived from GWAS results (blue) and external data sources (orange) including eQTL and promoter Capture-C analyses, and gene expression in fetal testis. To explore potential functional variants, Empirical Bayes modeling in PAINTOR (green) was conducted for all credible risk variants after annotation from multiple publicly available and locally derived data sources.

Full size image

Type 2 TGCTs originate from either fetal primordial germ cells or gonocytes and then develop from the noninvasive precursor GCNIS22. In the absence of available RNA sequencing data on GCNIS, we used single-cell RNA sequencing data from Li et al.30 to evaluate candidate genes for expression in fetal gonads measured at various timepoints. We included male and female germ cells and soma to get a complete picture of the potential expression of genes that may be influencing TGCT development. Transcript levels were categorized as low expression (≤698) for 33 (31%) genes, medium expression (699–2348) for 37 (34%) genes, and high expression (≥2349) for 38 (35%) genes based on tertiles of expression values (Supplementary Fig. 5; Supplementary Data 5). We conducted Assay for Transposase-Accessible Chromatin analysis using sequencing (ATAC-seq) on four TGCT cell lines (Supplementary Data 8, available at https://genome.ucsc.edu/s/jpluta/TECAC2020). The CRVs were significantly enriched in open chromatin regions in all cell lines (EP2102, P = 0.0015; NT2-D1 [NTERA2], P = 2.63 × 10−10; NCCIT, P = 4.37 × 10−8; TCAM2, P = 1.04 × 10−14), consistent with a potential effect on gene regulation. We further evaluated data from ATAC-seq in the context of promoter Capture-C data available on one of the cell lines, NT2-D1 to determine whether the promoter region of a target gene demonstrated a connection with a CRV located in an open chromatin region. Seventeen (16%) genes demonstrated these connections (Supplementary Data 5). Connections appearing in two or more cell lines were scored more highly than a connection found in just one cell line. Based on this evaluation of the potential target genes on the autosomes, we classified 37 (37%) genes as highly likely, 25 (25%) genes as moderately likely, and 39 (39%) as unlikely to be associated with TGCT; genes with multiple classification levels were counted in the highest likelihood group (Table 1, Supplementary Data 2, 5).

Inference of sex chromosome genes associated with TGCT

On the X chromosome, we assessed the gene regions delimited by 271 CRVs corresponding to five top signals (Supplementary Table 8). The total number of unique target genes interrogated was seven. Due to the absence of available eQTL data for X chromosome genes and the lack of expression data in fetal gonads for one target gene, it was not possible to create an equivalent schema to evaluate candidate target genes on the X chromosome. Still, based on our reduced evaluation scheme, one (14%) gene was scored as highly likely and two (29%) genes as moderately likely to be associated with TGCT (Table 1; Supplementary Table 8). However, should eQTL and expression data become available, the four (57%) genes unlikely to be associated with TGCT could be scored as highly or moderately likely (and similarly, the two moderately likely genes could be scored as highly likely); thus, we considered all genes as possible target genes (Supplementary Table 8).

Testis-specific gene enrichment

Genes selected for enrichment analysis included target genes (n = 62) on autosomes that scored moderately or highly likely to be associated with TGCT and all target genes (n = 7) on the X chromosome; two of these genes did not have available expression data. There was enrichment of testis-specific expression (P = 0.00067) in this gene set with three genes having at least 5-fold greater expression in testis compared to all other tissues (Supplementary Fig. 6). The expression of three other genes was enhanced in testis as indicated by five-fold or greater expression in testis compared to the average in all other tissues.

Functional assessment of variants by PAINTOR analysis

We also explored potential functional variants determined by PAINTOR, a Baysian approach that combines genetic association, linkage disequilibrium and enriched genomic features (Fig. 3)33. We annotated all 4755 CRVs with information from 36 datasets relevant to TGCT, including publicly available data and locally generated data from TGCT cell lines (histone marks, open chromatin marks, transcription factor binding sites, methylation), adult testis (histone marks, open chromatin marks, transcription factor binding sites, methylation, transcription start sites), embryo testis (open chromatin marks), and fetal testis (open chromatin marks) (Supplementary Table 9; Supplementary Fig. 7). PAINTOR analysis prioritized 100 variants as potentially functional, the majority of which had high posterior probabilities (≥95%); four (4%) variants had a posterior probability between 90 and 95%, and only one (1%) fell below 90%34 (Supplementary Data 9). Potentially functional variants were found for 57 (86%) of the 66 top signals. Two top signals, rs55873183 in DMRT1 and rs17336718 in TKTL1, contained only one CRV and thus could not be evaluated by PAINTOR. Most variants identified through the PAINTOR analysis were intronic (67%), one was exonic, and most (20%) of the remainder fell within 10 kb of the target gene start site. Of the 102 variants, 83 (81%) disrupted transcription factor binding sites.

Our meta-analysis has increased the number of susceptibility loci for TGCT by one-third. Men in the 95th percentile of the PRS have a 6.8-fold increased disease risk compared to men at the median PRS (Fig. 2); and these men have a 3.4% lifetime risk as compared to 0.4% in the general population2. The PRS for TGCT contains fewer SNPs than those available for most other common cancers, yet with a larger effect. For example, women in the 95th percentile of the PRS for breast cancer (313 SNPs) had a 2.4-fold increased disease risk compared to women at the median PRS35. The performance of the PRS derived from TGCT susceptibility loci suggests that men at highest risk of disease can be identified.

Evaluation of top association signals from our meta-analysis identified 65 target genes that were evaluated as moderately or highly likely to be associated with TGCT. Many of these genes encode proteins that fall into biological pathways relevant to TGCT susceptibility, including those that influence male germ cell specification and migration, sex determination and maturation, and regulation of the mitotic cell (HSA-69618, FDR 8.5 × 10−5; Fig. 4). For several target genes, findings from murine models support their direct role in the development of TGCT or TGCT-related phenotypes.

Deletion variants at the Steel locus (Sl) on the murine 129/Sv background are associated with increased incidence of TGCT; and the etiological gene has been demonstrated to be Kitl36,37. KITLG rs4474514 is the most statistically significant signal in our meta-analysis with a per-allele odds ratio over 2.0. Multiple other target genes implicated by top association signals influence male germ cell development in the mouse. Prdm14 is critical for the specification of primordial germ cells from somatic cells, participating in the reacquisition of potential pluripotency and successful epigenetic reprogramming38. The identified region on 12q13.2 contains two candidate target genes, SP1 and AMHR2 (Supplementary Fig. 1n). SP1 is a transcription factor that regulates cellular processes, including inhibition of mouse embryonic stem cell differentiation39. eQTL analysis suggest that the potential functional variant is associated with SP1 upregulation, thus similarly favoring developmental arrest by maintaining fetal germ cells in a relatively dedifferentiated state. AMHR2 is the receptor for anti-Mullerian hormone (AMH) which, in addition to testosterone (and hence involvement also of AR), results in male sex differentiation, preventing the development of Mullerian ducts into the uterus and fallopian tubes40. In the Japanese rice fish (medaka), knockout of amrh2 is associated with sex reversal and excessive proliferation of germ cells41.

Although we did not define AR as a moderately or highly likely target gene due to the lack of available data to inform eQTL analysis or the evaluation of gene expression in fetal testis, the top marker at Xq12 suggests that AR may be involved in the etiology of TGCT. Disruption of AR leads to androgen insensitivity syndrome and partial sex reversal, depending on the degree of disruption42. Furthermore, high linkage disequilibrium (r2 = 1) exists between the AR locus and variants associated with a decrease in male-pattern baldness43(Supplementary Data 4), a phenotype previously associated with risk of TGCT44. Immunohistochemical investigations also identified AR protein to be present in 40–50% of seminoma and GCNIS samples45. Further evaluation of this gene is warranted, results from which may further support the long-held hypothesis that a relative decrease in androgen compared to the overall population contributes to risk of TGCT46.

We identified a fourth independent susceptibility allele at DMRT1, which plays a critical role in sex determination and maintenance of the male somatic niche47. Expression of DMRT1 is enriched in testis tissue. Loss of Dmrt1 on the murine 129/Sv background leads to an over 90% incidence of testicular teratomas, due to a lack of ability to silence regulators of pluripotency48,49. Knockout of Dazl, a master transcriptional regulator essential for spermatogenesis, causes spontaneous gonadal teratomas, likely due to prolonged expression of pluripotency genes50,51. Expression of DAZL is also enriched in testis tissue.

BAK1 and BCL2L11 are both members of the BCL-2 family, which together tightly regulate the mitochondrial apoptotic response to either facilitate or prevent cell death depending upon intercellular stimuli52. Bak (BAK1) is a pro-apoptotic effector of mitochondrial outer membrane permeabilization, which allows release of cytochrome C and other apoptogenic factors leading to cell death53. Bim (BCL2L11) is a pro-apoptotic BH3-only protein that can activate Bak, but preferentially activates pro-apoptotic effector Bax54,55. Interestingly in mouse models, Bim and Bik cooperate to initiate early germ cell apoptosis in a biological pathway that appears to require Bax, but not Bak56. Bax also controls apoptosis of fetal germ cells during their migration, and in Bax null mice ectopic germ cells with retained primitive markers are observed57,58. Further 60% of NestinCreBaxfl/flBak−/− mice develop high-grade tumors within the testis that have expression profiles consistent with germ cell tumors59. Our eQTL analysis suggests downregulation of BCL2L11, implying improper survival of arrested germ cells and their transformation to pre-GCNIS.

We also identified multiple target genes encoding proteins involved in chromosomal segregation and heterochromatin organization. Inherited alterations in these genes likely contribute to unique hallmarks of TGCT that has one of the highest aneuploidy scores among cancers, characterized by near universal 12p isochromosome or amplification and frequent genome doubling60,61. PPP2R5A, a Ser/Thr phosphatase enriched at kinetochores and regulates chromosome-spindle interactions62, is an implicated target gene. Similar to AR, CENPI could not be defined as a top ranking target gene because of lack of available data for the X chromosome; but the top signal at Xq22.1 suggests that CENPI, a centromere protein and part of the CENPA-NAC (nucleosome-associated) complex responsible for chromosome alignment and segregation and mitotic progression important for gametogenesis63,64, may play a role in TGCT risk. At 9q34.3 the 29.5 kb haplotype block (rs28393706) contains two putative effector genes with overlapping promoter regions, ANAPC2, an E3 ligase enzyme that promotes metaphase-anaphase transition as part of the anaphase-promoting complex (APC), and SSNA1 (SS nuclear autoantigen 1), a centrosomal protein regulating the microtubule-severing activity of spastin65,66. Six implicated genes (PMF1, PPP2R5A, ANAPC2, SSNA1, TEX14, and MCM3AP) have an eQTL associated with downregulation, consistent with a more permissive phenotype for chromosomal mis-segregation; and expression of TEX14 is enriched in testis tissue. Further, multiple TGCT-implicated proteins in the chromosomal segregation pathway interact with TGCT-implicated male germ cell development proteins, demonstrating a biological network underlying TGCT susceptibility (Fig. 4).

After pathway analysis of moderately and highly ranking target genes, several were found to encode proteins that interact in mRNA translation, including one of the ribosomal proteins (RPL4), translation termination protein eRF3A (GSTP1) and translocon-associated protein subunit gamma (TRAP-gamma, encoded by SSR3), which is the general ribosomal interactor participating in the co-translational translocation of proteins into the endoplasmic reticulum67. Finally, multiple DNA-binding transcription factors are implicated in TGCT susceptibility, including HNF1B, PITX1, PKNOX2, PRDM14, SP1, TFCP2L1, ZFPM1, ZNF64, and ZNF217. Several are zinc finger proteins (ZNF) (including KRAB-ZNF) critical for proper germ cell development, such as male primordial germ cells specification and epigenetic reprogramming68.

Results from our investigation provide further understanding of the genetic architecture of TGCT, enhance comprehension of the biology of male germ cell development, and highlight biological pathways important to TGCT that are not noted in other cancers. Our findings implicate potentially important pathways, including regulation of apoptosis beyond the BAK1-BCL2L11 axis (AIFM3, CLPTM1L), enzymatic functions (MPV17L, TKTL1, and UCK2) and several genes involved in actin, cytoskeleton, and microtubule organization (CYTH1, ENOSF1, TNXB, and ARL14EP). The latter may contribute to errors in germ cell migration or chromosomal segregation, likely enhancing the dysregulation of genes directing the germ cell-somatic niche interaction during early development (KITL, DMRT1).

Our meta-analysis has identified 66 validated susceptibility loci for TGCT. Many of these loci have a stronger effect size than those observed in adult epithelial cancers, which results in a high fraction of explained heritability of TGCT. Many TGCT risk alleles have higher frequencies in men of European compared to African genetic ancestry, concordant with the known difference in disease incidence between these groups. Importantly, we have established a PRS that identifies men at highest risk of disease. This TGCT PRS could be potentially applied in men with other risk factors, such as cryptorchidism or infertility, to be targeted for early detection and disease mitigation.