Identification of nine new susceptibility loci for testicular cancer, including variants near DAZL and PRDM14.

Testicular germ cell tumor (TGCT) is the most common cancer in young men and is notable for its high familial risks. So far, six loci associated with TGCT have been reported. From genome-wide association study (GWAS) analysis of 307,291 SNPs in 986 TGCT cases and 4,946 controls, we selected for follow-up 694 SNPs, which we genotyped in a further 1,064 TGCT cases and 10,082 controls from the UK. We identified SNPs at nine new loci (1q22, 1q24.1, 3p24.3, 4q24, 5q31.1, 8q13.3, 16q12.1, 17q22 and 21q22.3) showing association with TGCT (P < 5 × 10(-8)), which together account for an additional 4-6% of the familial risk of TGCT. The loci include genes plausibly related to TGCT development. PRDM14, at 8q13.3, is essential for early germ cell specification, and DAZL, at 3p24.3, is required for the regulation of germ cell development. Furthermore, PITX1, at 5q31.1, regulates TERT expression and is the third TGCT-associated locus implicated in telomerase regulation.

Testicular germ cell tumor (TGCT) is the most common malignancy in men aged 15–45 years[1]. Studies in families have estimated the risk to brothers of TGCT cases to be increased 8- to 10- fold and to sons of cases 4- to 6-fold, and these risks are substantially higher than the equivalent familial relatives risk of ~2 typical for common cancers such as breast, colorectal and prostate[2]. Recent genome-wide association studies (GWAS) of TGCT have identified eight associated SNPs at six loci, which together account for >11% of the genetic risk of TGCT[3-7]. These loci lie on chromosome 12q21, encompassing KITLG, 5q31 (SPRY4), 6p21 (BAK1), 5p15 (TERT and CLPTM1L), 12p13 (ATF7IP) and 9p24 (DMRT1)[11].

To identify further susceptibility loci for TGCT, we identified the most strongly associated 1,050 SNPs from our previous GWAS, which comprised 307,291 SNPs genotyped in 986 cases of TGCT and 4,946 controls from the UK[3,4]. These SNPs were included on a custom Illumina iSelect genotyping array (the iCOGS array). The iCOGS array includes 211,155 cancer-related SNPs selected for further evaluation in different cancers and has already yielded 80 new loci for breast, ovarian and prostate cancer[12,13]. Using the iCOGS array, we genotyped 1,094 cases of TGCT from the UK, which had not been included in the original GWAS. We compared these genotypes to those of three cancer-free UK control series, also genotyped on the iCOGS array. Following quality control exclusions, we obtained data for 694 of the 1,050 SNPs in 1,064 cases and 10,082 controls (the ‘iCOGS replication series’).

We tested association between each SNP and TGCT risk using a 1df trend test, adjusted for six principle components. Inflation in the test statistic, quantified using 36,879 uncorrelated non-TGCT-related SNPs, revealed evidence of only modest inflation of the test statistics (Supplementary Figure 1, λ=1.14, λ1000=1.07; following adjustment for six principle components, λ=1.05, λ1000=1.02). There was clear evidence of association of the six previously reported loci in the iCOGS replication series (Ptrend= 3.88×10−9 to 9.31×10−28), with odds ratios consistent with those previously published (Supplementary Table 1).

In a combined fixed-effects meta-analysis of the GWAS and the adjusted iCOGS replication for the remaining 640 SNPs outside of these six loci, we identified association of P<5×10−8 for a further 18 SNPs at nine loci: 1q22, 1q24.1, 3p24.3, 4q24, 5q31.1, 8q13.3, 16q12.1, 17q22 and 21q22.3. Following multiple logistic regression to test for independence of effects, there remained one SNP at each of the nine loci showing association with TGCT (Table 1). The associations at 1q24.1 and 17q22 have been replicated concurrently in other independent TGCT series[14,15].

Table 1

Summary results for nine SNPs showing genome-wide association (P<5×10−8) in the combined analysis.

SNP1	Location2	Alleles3	RAF4	GWAS		iCOGS replication		Combined
				OR5(95% CI)	Ptrend6	adjusted OR7(95% CI)	Ptrend6	Pmeta8
rs2072499	154,436,2341q22	G/A	0.35	1.24(1.12-1.37)	1.84×10−5	1.19(1.08-1.3)	2.98×10−4	2.57×10−8
rs3790672	164,140,0161q24.1	C/T	0.28	1.26(1.14-1.4)	1.10×10−5	1.20(1.09-1.33)	2.92×10−4	1.59×10−8
rs10510452	16,600,0523p24.3	A/G	0.69	1.24(1.12-1.39)	7.06×10−5	1.24(1.12-1.37)	4.69×10−5	1.28×10−8
rs2720460	104,274,1354q24	A/G	0.61	1.28(1.16-1.42)	2.25×10−6	1.24(1.12-1.36)	1.97×10−5	2.18×10−10
rs3805663	134,394,0995q31.1	T/C	0.63	1.20(1.08-1.33)	4.86×10−4	1.25(1.13-1.38)	7.14×10−6	1.56×10−8
rs7010162	71,139,0598q13.3	G/A	0.61	1.21(1.09-1.34)	2.03×10−4	1.22(1.11-1.34)	6.05×10−5	4.60×10−8
rs8046148	48,700,44516q12.1	G/A	0.79	1.29(1.13-1.46)	1.03×10−4	1.32(1.17-1.48)	8.75×10−6	3.68×10−9
rs9905704	53,987,54217q22	T/G	0.67	1.29(1.16-1.43)	3.74×10−6	1.21(1.1-1.34)	1.50×10−4	3.06×10−9
rs2839186	46,514,49621q22.3	T/C	0.46	1.20(1.09-1.32)	3.01×10−4	1.26(1.15-1.38)	7.43×10−7	1.18×10−9

dbSNP rs number

Build 36 position, chromosomal region

Risk/non-risk associated alleles.

RAF:Frequency of the risk allele

OR: per allele odds ratio

Ptrend: p-value for trend, via logistic regression

Adjusted OR: per allele odds ratio, adjusted for six principle components

Pmeta: P-value for fixed effects meta-analysis

For each SNP, we examined for evidence of departure from a log-additive (multiplicative) model: only rs2072499 at 1q22 showed evidence for departure from the log-additive model, indicating a recessive effect (heterozygote OR (ORhet): 1.00 (0.87-1.16), homozygote OR (ORhom): 1.53 (1.27-1.84), P2df=5.18×10−6, Supplementary Table 2). This observation was supported by genotypic (2df) analysis of the GWAS data for rs2072499 (ORhet 1.09 (0.94-1.26) ORhom 1.64 (1.24-2.01), P2df 1.25 ×10−6). We found no evidence for statistical interaction between any pair of these new nine SNPs or with the previous eight SNPs; for every combination of SNPs, the combined risk was consistent with the product of the individual risks.

We investigated whether the newly identified loci were associated with different risks in subgroups of TGCT cases characterized by specific phenotypic characteristics (Supplementary Tables 3,4,5). Consistent with our findings at the previous six loci, none of the newly identified loci showed a significant difference in effect when comparing cases of the two TGCT histological subtypes, seminoma and nonseminoma. Likewise, there was no clear evidence of differential association by age of diagnosis, presence of testicular maldescent, family history of TGCT or unilateral versus bilateral disease, although power to detect differences in effects was limited on account of the case distribution.

Encyclopedia of DNA Elements (ENCODE) data indicated that rs2072499, rs3790672, 3805663, rs7010162, rs9905704 and rs2839186 may have effects on transcription on account of evidence of location within a DNase I hypersensitivity site, local modification of histones indicating promoter/enhancer activity, influence on binding of transcription factors and/or effect on a regulatory motif. Functional annotations from the ENCODE data for all correlated variants within each of the nine LD block (r2>0.8), are shown in Supplementary Table 6. We then investigated whether the genotype at these nine SNPs was associated with differential expression of genes lying within 500kB of the SNP: there were no compelling patterns evident from expression data available from fibroblast, skin, adipose and lymphoblastoid cells (Supplementary Table 7). Further examination in testicular germ cells both for evidence of transcriptional regulation and of gene expression would be of interest.

Evidence supports TGCT developing from primordial germ cell or gonocytes, with tumor initiation described in utero through a pre-invasive stage termed intra germ cell neoplasia unclassified (IGCNU) or carcinoma in situ (CIS)[16]. It is therefore noteworthy that a number of the genes identified at loci identified in this study relate to pathways of early differentiation and development of the testicular germ cell[11] (Table 1).

rs10510452 lies in an LD block of 82kb at 3p24.3 which contains only one gene, DAZL (deleted in azoospermia-like). DAZL encodes an RNA-binding protein which has been shown from expression studies in human primordial germ cells to have a central role in early differentiation of primordial germ cells[17]. Dazl knockout mice are infertile, with differentiation of the germ cells halted at the A spermatogonia level [9] and expression of DAZL is testis-specific [18,19]. Sequence analysis indicates that the DAZ cluster on the Y-chromosome arose during primate evolution via transposition of the DAZL gene[20]. Of note, DAZ has been implicated as the critical gene within the Y-chromosome AZF (azospermia factor) region, so-named as deletions of this region cause non-obstructive spermatogenetic failure (azoospermia) in humans[21].

rs7010162 lies in an LD block of 61kb at 8q13.3 in which there is one gene, PRDM14, which is required for both the reacquisition of potential pluripotency and the genomewide epigenetic reprogramming necessary for primordial germ cell specification[8,22,23]. PRDM14 encodes a transcriptional regulator which controls expression of key pluripotency genes such as POU5F1 (OCT4), NANOG and SOX2, which are highly expressed in IGCNU and TGCT[24,25]. Whilst normal PRDM14 expression is restricted to germ cell and stem cell lineages, amplification and/or over-expression of PRDM14 have been reported in a number of different cancer types, including TGCT[26-29].

rs3805663 is located in a 51kb LD block at 5q31.1 which contains two genes: CATSPER3 and PITX1. CATSPER3 (Cation-Channel, Sperm-associated 3) is one of a family of genes which together form a functional hetero-tetrameric cation channel, exclusively expressed in the testis and essential for hyperactivated motility of the differentiated germ cell (spermatozoa); accordingly Catsper3−/−male mice are completely infertile [30-32]. An equally plausible candidate gene at this locus is PITX1, (paired-like homeodomain transcription factor 1), which binds to specific PITX1-binding sites in the promoter region of TERT, thus regulating expression of telomerase[10]. Thus, along with the previously identified TGCT-associated loci at 5p15 (TERT) and 12p13 (ATF7IP, which regulates expression of TERT), the signal at 5q31.1 may also be mediated via effects on telomerase regulation. Whilst the 5p15 (TERT) region is notable for its associations with multiple cancers[33-39], to date there have been no associations reported with other cancers for the loci at 5q31.1 (PITX1) or 12p13 (ATF7IP).

rs9905704 lies at 17q22 in an LD block of 813kB in which there are several genes including RAD51C and TEX14 (testis expressed 14). TEX14 is a protein kinase highly expressed in the human male germ cell [40,41]. Male Tex14 −/− mice are infertile whilst the female knockouts are normally fertile[42]. Cellular studies have demonstrated that TEX14 regulates kinetochoremicrotubule assembly and spindle assembly checkpoint in the germ cells of the testis and requires recruitment of key proteins such CENPE, a kinesin-like motor protein which associates with the centromere [41-43]. It is thus noteworthy that CENPE (Centromere-associated protein E) is located within the 214kb TGCT-associated LD block at 4q24. Furthermore, PMF1, one of the two genes located in the TGCT-associated region at 1q22, is part of the MIS12 complex which is required for normal chromosome alignment, segregation and kinetochore formation during mitosis[44]. Therefore, regulation of microtubule assembly may be the mechanism linking the signals at 17q22, 4q24 and 1q22.

In summary, we have identified nine new loci for TGCT, which brings the total number of TGCT-associated loci identified to date to fifteen. These fifteen loci have provided considerable new insights into TGC tumourigenesis, implicating genes involved in germ cell specification and/or differentiation (DAZL, PRDM14) including the KIT-KITL signalling pathway (KITL, SPRY4, BAK1) and genes involved in sex-determination (DMRT1), microtubule assembly (TEX14, CENPE, PMF1) and telomerase regulation (TERT, ATF7IP, PITX1). The nine new susceptibility alleles account for ~4% of the excess familial risk to brothers and ~6% to sons of men with TGCT: this brings the cumulative totals to ~15% and ~22% respectively. The power to detect these loci in this combined analysis was modest (29-84%): accordingly there are likely to be several additional loci for TGCT of equivalent or lesser effects, which may be identifiable via still larger follow-up studies of these data and meta-analyses with other GWAS data. Moreover, imputation and fine-scale association mapping of the loci may uncover more of the missing heritability, either through identifying more strongly associated variants or by identifying additional signals, including rare variants, at these loci. Integration of these genetic variants together with non-genetic TGCT risk factors into risk models may enable clinically useful risk profiling for TGCT in unaffected men; which in turn may provide a rational basis for application of strategies for screening and targeted prevention of TGCT.

Methods

Samples

See Supplementary Note

GWAS analysis

Cases were genotyped on the Illumina HumanCNV370-Duo bead array and controls were genotyped on the Illumina Infinium 1.2M array at the Wellcome Trust Sanger Institute (see [3] and [4] for full methods). The genotypes were re-called since the analysis reported in [4], such that SNPs on both sex chromosomes were included in the current analysis. We used data on 314,861 SNPs that were successfully genotyped on both arrays. We excluded individuals: with low call rate (<95%), with abnormal autosomal heterozygosity or with >10% non-European ancestry (based on multi-dimensional scaling). We filtered out all SNPs with (i) minor allele frequency <1%, (ii) a call rate of <95% in cases or controls or (iii) minor allele frequency of 1–5% and a call rate of <99% or (iv) deviation from Hardy-Weinberg equilibrium (10−12 in controls and 10−5 in cases). 307,291 SNPs passed the QC filters and we put forward for replication the 1050 SNPs with the lowest P-value (P<1.3 × 10−3), based on the 1-d.f. Cochran-Armitage trend test.

iCOGS Genotyping

Genotyping was conducted using a custom Illumina Infinium array (iCOGS array) comprising 211,155 SNPs selected across multiple consortia within the COGS (Collaborative Oncological Gene-environment Study, see [12] for details of SNP inclusion). Of the 1,050 SNPs submitted from our analysis, 740 attained an Illumina design score of ≥0.8, and were included on the array. Genotyping was performed at Genome Quebec, one of the seven genotyping centres providing genotyping for the COGS collaboration. Genotypes were called using Illumina’s proprietary GenCall algorithm. Pan-COGS initiatives were employed across experiments to optimize calling algorithms for the iCOGS array: firstly, initial calling employed a cluster file generated using samples from Hapmap2; secondly, a cluster file based on 3,018 individuals selected across genotyping centers, across participating consortia and across ethnicities, combined with 380 samples of European, Asian or African ancestry genotyped as part of the Hapmap and 1000 genomes project was applied to call the genotypes for the subsequent samples. Thirdly, the performance of GenCall was compared to two other calling algorithms: Illuminus54 and GenoSNP55: all three algorithms were >99% concordant in their calling for 91% of the SNPs on the array. Manual inspection of a sample of the discrepant SNPs indicated that the GenCall calls were almost invariably superior (because Illuminus and GenoSNP frequently attempted to call SNPs that clustered poorly) [12].

Quality Control

Using the full SNP set of 211,155 SNPs on the iCOGS array, we applied quality control exclusions as follows to subjects: i) subjects with overall call rate <95% or low or high heterozygosity (P<10−6) (5 cases), ii) using identity-by-state estimates based on 37,046 uncorrelated SNPs, we identified “cryptic” duplicates and related samples and the sample with the lower call rate was excluded (7 cases), iii) we identified ethnic outliers by multi-dimensional scaling by combining the iCOGS data with the three Hapmap2 populations using 37,046 uncorrelated markers and removed individuals with >10% non-Western European ancestry (18 cases). We included 1,064 cases and 10,082 controls in the final analysis.

We applied quality control exclusions as follows to SNPs: i) discrepant calls in more than 2% of duplicate samples across COGS consortia, ii) call rate <95%, MAF<1%, call rate <99% if MAF=1-5%, iii) deviation from Hardy-Weinberg (P<10−5 in controls, P<10−12 in cases ). Following quality control exclusions applied individually to case and control data, we included genotypes from 694 SNPs in subsequent analyses. We checked genotype intensity cluster plots manually for SNPs in each new region for which a genome-wide significant association was achieved.

Statistical Analysis

For the 694 SNPs in the TGCT iCOGS replication experiment, we estimated the per-allele odds ratio (OR), 1df, using logistic regression, adjusting for principal components as covariates. We generated eight principle components and examined reduction in the inflation factor using iterative combinations of principle components. Inclusion of principle components 1-6 achieved optimal reduction in inflation with no further reduction achieved by incorporating additional principle components into the regression model. We found no significant differences in the point estimates of effect size on performing sensitivity analysis between the control groups used in the TGCT iCOGS replication experiment.

The inflation factor (λ) was derived by dividing the median of the lowest 90% of the 1-d.f. statistics by the 45% percentile of a 1-d.f. χ2 distribution (0.357), utilizing 36,879 uncorrelated SNPs selected to not include TGCT-associated SNPs. The inflation was converted to an equivalent inflation for a study with 1000 cases and 1000 controls (λ1,000): We excluded from further analysis the 54 SNPs lying within the six previously identified TGCT-associated loci, defining these LD blocks via the Oxford recombination rates. We included the 640 remaining SNPs in the replication analysis.

Because the inflation factor in the GWAS data had been modest (λ = 1.078, λ1000 = 1.045), we did not adjust the analysis of the GWAS data. We obtained overall significance levels by combining the estimates from the GWAS (unadjusted) and the iCOGS replication (adjusted for principle components 1-6) using a fixed effects meta-analysis, to derive a 1df test, calculating an I2 statistic to evaluate heterogeneity between the GWAS and replication. We used a threshold of 5×10−8 to denote genome-wide significance. Of note we had previously examined the loci at 1q24.1 and 4q24 in our first GWAS replication experiment[3]: whilst there had been evidence of replication of these signals via Taqman genotyping of 565 cases and 1,758 controls, these associations had failed to achieve genome-wide significance.

We also computed genotype-specific ORs (2 df) for the iCOGS replication, adjusted for principle components 1-6 (Supplementary Table 2). We assessed each SNP for dose response (Phetergeneity) in the replication series by comparing 1-d.f. and 2-d.f. logistic regression models,both adjusted for principle components 1-6, using a likelihood ratio test (evaluated using a significance threshold of Phetergeneity < 0.006 to account for 9 tests). We examined for statistical interaction between TGCT predisposition SNPs (nine new SNPs and the previously published eight SNPs) by evaluating the effect of adding an interaction term to the regression model, adjusted for stage, using a likelihood ratio test (using a significance threshold of P < 0.0003 to account for 144 tests).

We assessed for modification of the odds ratios by covariate phenotype and risk factors using a case-only analysis incorporating both GWAS and iCOGS replication data (unadjusted). We evaluated the effect of age on SNP genotype via polytomous regression: we divided age at diagnosis into 5 ordinal categories and fitted maximum-likelihood multinomial logistic models, executed using the mlogit command in Stata12 (using a significance threshold of P < 0.006 to account for 9 tests, Supplementary Tables 3, 4, 5).

We identified variants from within each LD block reported within the 1000 genomes project (r2>0.8 and <200kB from sentinel SNP) and used HaploReg and ENCODE data to apply functional annotations relevant to the regulation of transcription: (i) whether the variant lies in a region in which modification of histone proteins is suggestive of enhancer and other regulatory activity (H3K4Me1 and H3K27A histone modification) or promoter activity (H3K4Me3 histone modification), (ii) whether the variant lies in a region where the chromatin is hypersensitive to cutting by the DNase enzyme (suggestive of regulatory region), (iii) whether the variant lies in a region of binding of transcription factor proteins (as assayed by chromatin immunoprecipitation with antibodies specific to the transcription factor followed by sequencing of the precipitated DNA (ChIP-seq)), (iv) whether the variant affects a specific regulatory motif, as evaluated from position weighted matrices assembled from TRANSFAC, JASPAR and protein-binding microarray experiments (Supplementary Table 6)[45,46].

We investigated associations of the nine TGCT-associated SNPs with gene expression using GENEVAR, which includes gene expression profiling and genotypic data from four datasets: (i) lymphoblastoid cell lines (LCL) from 726 HapMap individuals adipose (HapMap3), (ii) LCL and skin collected from 856 healthy female twins of the MuTHER resource (Muther pilot Twin 1), (iii) adipose, LCL and skin derived from a subset of 160 MuTHER healthy female twins (Muther pilotTwin2) and (iv) fibroblast, LCL and T-cells derived from umbilical cords of 75 Geneva GenCord individuals (Gencord). Spearman Rank Correlation between normalized gene expression levels and the count of one of the alleles of the SNP (0, 1 or 2) was analysed with significance assessed by permutation (10,000 permutations). Expression of genes within 500kB of the sentinel SNP was examined (using a significance threshold of P < 0.00001 to account for 440 tests, Supplementary Table 7).

LD matrices between SNPs reported in HapMap were based on Data Release 27/phaseII + III (Feb 2009) on NCBI B36 assembly, dbSNP b126, viewed using Haploview software (v4.2) and plotted using SNAP. LD blocks were evaluated using the HapMap recombination rates (cM/Mb) and defined using the Oxford recombination hotspots[47]. GWAS associations for other diseases were evaluated using the NHGRI (National Human Genome Research Institute) Catalog of Published Genome-Wide Association Studies. All genomic references are based on NCBI Build 36. Analyses were performed using R (v2.6), Stata12 (State College) and PLINK (v1.07) software.