A meta-analysis of 87,040 individuals identifies 23 new susceptibility loci for prostate cancer

Genome-wide association studies (GWAS) have identified 76 variants associated with prostate cancer risk predominantly in populations of European ancestry. To identify additional susceptibility loci for this common cancer, we conducted a meta-analysis of >10 million SNPs in 43,303prostate cancer cases and 43,737 controls from studies in populations of European, African, Japanese and Latino ancestry. Twenty-three novel susceptibility loci were revealed at P<5×10-8; 15 variants were identified among men of European ancestry, 7 from multiethnic analyses and one was associated with early-onset prostate cancer. These 23 variants, in combination with the known prostate cancer risk variants, explain 33% of the familial risk of the disease in European ancestry populations. These findings provide new regions for investigation into the pathogenesis of prostate cancer and demonstrate the utility of combining ancestrally diverse populations to discover risk loci for disease.

Prostate cancer is the most common non-skin cancer in men in the Western world and epidemiological studies have shown strong evidence for genetic predisposition to prostate cancer, based on two of the most important factors, ancestry and family history. Genomewide association studies (GWAS) have identified 76 common risk loci (reviewed in ref 1 ); however, over 1,000 additional common SNPs are estimated to contribute prostate cancer risk. 2,3 Previous prostate cancer GWAS have been conducted primarily in populations of European ancestry 2,[4][5][6][7] , with the majority of risk loci discovered also found to be associated with prostate cancer risk in other racial/ethnic populations. 8,9 The generalizability of risk associations for a large fraction of loci suggests that combining GWAS across ancestral populations could increase power to detect risk loci that are shared among diverse populations.
To search for additional genetic risk factors for prostate cancer, we combined data from studies with existing high-density SNP genotyping in prostate cancer GWAS discovery or replication efforts in the following populations: European ancestry[34,379 cases and 33,164 controls from UK/Australia 4 , Cancer of the Prostate in Sweden (CAPS) 10 , Breast and Prostate Cancer Cohort Consortium (BPC3) 6 , PEGASUS, and iCOGS/PRACTICAL 2 ]; African ancestry[5,327 cases and 5,136 controls from the African Ancestry Prostate Cancer GWAS Consortium (AAPC) 11 and the Ghana Prostate Study 12 ]; Japanese ancestry[2,563 cases and 4,391 controls from a GWAS in Japanese in the Multiethnic Cohort (MEC) 8 , and Biobank Japan 13,14 ]; and, Latino ancestry[1,034 cases and 1,046 controls from the MEC 8 ]. Imputation was performed in each study using a cosmopolitan reference panel from the 1000 Genomes Project (1KGP; March, 2012). Across the various studies, 5.8-16.8M genotyped and imputed SNPs, as well as insertion/deletion variants ≥1% frequency were examined in association with prostate cancer risk (Online Methods, Supplementary Tables 1-3, Supplementary Information).
We first conducted ethnic-specific meta-analyses, with the large European ancestry sample providing the strongest statistical power for discovery of novel loci, followed by a multiethnic meta-analysis of all populations to identify additional loci with pan-ethnic effects. For these primary analyses we employed a P-value threshold of 5×10 -8 to define genome-wide significance. Secondary meta-analyses focused on a) aggressive disease in the large European ancestry sample; b) aggressive disease in the combined multiethnic sample; and c) prostate cancer diagnosed at ≤55 years of age in the European ancestry sample only. Aggressive prostate cancer was defined as a Gleason score ≥8, disease stage as 'distant', a prostate-specific antigen (PSA) level >100 ng/ml, or death from prostate cancer. For these two secondary phenotypes, we utilized a more stringent P-value threshold of 5×10 -8 /2=2.5×10 -8 for genome-wide significance. In each study, we tested for gene dosage effects via a 1-d.f. test for trend from logistic regression models adjusted for genetic ancestry (principal components). We observed little evidence of inflation in the test statistics in any single study or population (λ/λ 1000 : European, 1.14/1.00; African, 1.03/1.01; Japanese, 1.06/1.02; Hispanic, 1.00/1.00) or in the multiethnic analysis (λ=1.08,/λ 1000 =1.00; Online Methods, Supplementary Table 4, Supplementary Figure 1).
In the meta-analysis of the European ancestry studies,20 novel signals in18 regions ±500 kb outside of previously associated loci were observed to be associated with prostate cancer risk at P<5×10 -8 ( Figure 1; Supplementary Figure 2 Supplementary Figure 3). The most significant associations in each region were observed with imputed variants and we were able to confirm the imputed genotypes for 15 variants which had high imputation information scores (r 2 range, 0.76-1) through direct genotyping or sequencing across multiple studies (Table 1;Online Methods, Supplementary Tables 5-8). Two of the variants were located within 370kb of each other on chromosome Xq13and are independent signals based on conditional analyses (rs6625711, P=6.1×10 -10 and rs4844289, P=2.0×10 -8 ; r 2 <0.01 in EUR 1KGP; Supplementary Table 9). All 15 variants were common, with minor allele frequencies (MAFs) ≥0.09, in the European ancestry population, and all but three (rs80130819/12q13, rs76939039/10q11 and rs17694493/9p21) were also common (MAF≥0.05) in African, Japanese and Latino populations. Evidence of heterogeneity in the per-allele OR was noted with 4 variants (P het =0.01-8.4×10 -6 ; rs17599629/1q21, rs115306967/6p21, rs17694493/9p21 and rs6625711/Xq13). Four of the 15 variants (rs10009409/4q13, rs4713266/6p24, rs80130819/12q13 and rs2807031/Xp11) had directional effects that were consistent with men of European ancestry and were nominally statistically significant (P<0.05) in at least one other population (Table 1) and for 3 SNPs, combining data across populations strengthened the statistical significance of the association ( Table 1). In this large European ancestry sample we also confirmed the reported signal at 22q13 with variant rs58133635 (P=5.8×10 -9 ; r 2 =0.74 with rs9623117 in 1KGP European ancestry populations (EUR); Supplemental Figure 2; Supplementary Figure 3). 15 No novel risk loci were revealed in ethnic-specific analyses within the African, Japanese or Latino ancestry populations possibly due to lack of power (Supplemental Figure 2). However, in combining results across populations in a multiethnic meta-analysis (43,303 cases, 43,737 controls), 11additional variants were identified in association with prostate cancer risk in novel risk regions at P<5×10 -8 (Table 1; Supplemental Table 5; Figure 2). We confirmed the imputed genotypes for 7 variants which had high imputation information scores (r 2 range, 0.81-1) through additional genotyping and sequencing (Online Methods, Supplementary Tables 6-8). All 7 variants were nominally associated with risk (P<0.05) in at least one of the non-European ancestry populations and the per-allele effects were directionally consistent across all 4 populations for 6 of the 7 variants. All variants had MAFs≥0.05 in all four populations, and no significant evidence of population heterogeneity was noted with any of these 7 variants (Table 1).
In secondary GWAS analyses, we detected an association with variant rs636291 at 1p36 (risk allele frequency, 0.16; OR=1.18; P=2.1×10 -8 ; Table 1) and early-onset disease among men of European ancestry (4,147 cases ≤55 years of age and all controls, n=27,212). The association with this variant was weaker for cases diagnosed >55 years of age (23,564 cases versus all controls, n=27,212: OR=1.04; p=0.004; P het =2.2×10 -4 ; Supplementary Table 10). We did not detect any genome-wide significant associations with aggressive disease in the European population (n=7,903 cases) or in the combined multiethnic sample (n=10,209 cases; Supplemental Figure 4).
For the 23 novel risk variants (15 in European, 7 in multiethnic and 1 in the early onset analysis), the per-allele effects ranged from 1.06-1.14 and were consistent with log-additive effects (Supplemental Table 11). The association of each variant was noted for both aggressive and non-aggressive prostate cancer (Supplemental Table 12); for only one variant, rs7153648 at 14q23, there was suggestive evidence of a difference by disease severity (OR=1.17 for aggressive and OR=1.09 for non-aggressive disease; P het =0.03).
These results confirm what has been observed in prostate cancer GWAS to date; risk loci appear to confer risk for prostate cancer overall and not discriminate between the aggressive and indolent disease. In analyses stratified by age, 17 of the 23 variants demonstrated larger effects at younger ages (≤55 versus >55 years), although only 6 had evidence of a significant difference (p<0.05) (Supplemental Table 9). Only two of the 23 variants was modestly associated with PSA levels among controls (rs9287719 at 2p25, P=0.03 and rs115306967 at 6p21, P=0.05; Supplemental Table 13).
Of the 23 novel risk variants, 13 are located in intronic regions of genes and 2 are correlated with non-synonymous variants in adjacent genes (rs12051443/16q22, r 2 =0.98 with rs4788821/E60Kin MARVELD3; rs2238776/22q11, r 2 =0.67 with rs72646967/ N397HinTBX1). Based on functional annotations of transcription factor (TF) occupancy, response element disruption, histone marks and DNaseI sensitive regions in prostate cancer cell lines (Online Methods), 12 of the risk variants are either directly located within putative functional elements or are correlated (at r 2 >0.9 in 1KGP EUR) with such variants (Supplementary Table 14). Using gene expression data for 145 prostate cancer tumor samples from The Cancer Genome Atlas (TCGA) (Online Methods) we also examined the cis-associations between the index SNP and expression of gene transcripts within a 1Mbregion. Among the 23 loci, 5 cis-associations were observed, albeit the associations were modest (Supplemental Table 14;Online Methods).
A number of the novel susceptibility regions are located in close proximity to genes which have either an established role, or have been directly implicated, in cancer ( Table 1). The most notable is rs1041449 on chromosome 21q22, which is situated 20kb 5′ of theTMPRSS2 gene which encodes a member of a serine protease family. 16 Expression of TMPRSS2 is highly specific to prostate tissue and chromosomal translocation resulting in fusion of the TMPRSS2 promoter/enhancer region with the ETS transcription factors ERG and ETV1 are frequently observed in prostate cancer. 17 In analyzing data of 552 tumors characterized for the TMPRSS2-ERG fusion (46% positive) (Online Methods), we found no evidence of an association between the risk allele and fusion status (p=0.53; Supplementary Table 15). The variant risk rs1041449 is located within a number of histone marks and TF occupancy sites in the predicted enhancer region of TMPRSS2 ( Figure 3) however we found little evidence that this variant influences TMPRSS2 expression in prostate tumors (n=244, P=0.60), or in normal prostate tissue (n=87, P=0.62) (Online Methods).
Another region of notable importance is on chromosome 9p21. The risk variant, rs17694493, is intronic in CDKN2B-AS1, which encodes a long non-coding RNA -ANRIL, and is part of the CDKN2B-CDKN2A gene cluster ( Figure 3). The region contains highly penetrant alleles for familial melanoma and common susceptibility alleles for melanoma, breast cancer, basal cell carcinoma, lung cancer and glioma. [18][19][20][21][22][23][24] The index SNP, rs17694493, falls within chromatin bio features and is predicted to disrupt two TF motifs (STAT1 and RUNX1) suggesting that it may have a functional effect on the regulation of the CDKN2B-AS1 or CDNK genes ( Figure 3, Supplementary Table 14), however, the variant was not found to be strongly associated with expression of either CDKN2A (P=0. 19) or CDKN2B (P=0.40) in the 145 TCGA prostate tumors.
Variant rs4713266 at chromosome 6p25, is located in intron 1 of NEDD9, a gene that participates in cell adhesion, motility, the cell cycle and apoptosis, and has been implicated in progression and metastasis of several cancer types. 25 Variant rs9443189 on chromosome 6q14 is intronic in MYO6, a modulator of androgen-dependent gene expression which has been found to be overexpressed in prostate cancer tumors and enhance prostate tumor growth and metastasis. [26][27][28] Variant rs636291 on chromosome 1p36, which we found in association with early-onset prostate cancer, is located in intron 2 of PEX14 and is correlated with rs616488 (r 2 =0.66 in 1000 Genomes Project, EUR population), a variant reported in a GWAS of breast cancer. 29 The identification of novel risk loci for prostate cancer through a multiethnic analysis demonstrates the value of combining genetic data across populations to increase statistical power for discovery. As further support for conducting multiethnic analyses, we examined the genome-wide evidence for consistency in the direction of the allelic associations between populations. Excluding SNPs ± 500kb of index signals at known loci (n=77), we defined independent signals (r 2 <0.2) for the European ancestry population of nominal significance at various P-value thresholds between <10 -2 -10 -5 . For the sets of SNPs defined for men of European ancestry, 53-64% had ORs that are directionally concordant for African (p=0.04-0.003, dependent on the p-value threshold bin), Asian (p=0.31-0.02) or Hispanic men (p=0.04-0.002) with the ORs in Europeans. This same observation remained once we removed the 23 risk loci identified by the current study (Supplementary Figure 5). The excess of directionally consistent associations between populations implies that additional common risk loci for prostate may be revealed through discovery efforts in multiethnic studies. These 23 novel loci (includingrs58133635 at 22q13) 15 bring the total number of susceptibility variants for prostate cancer to 100 (Supplementary Table 16). In total, we estimate these 100 risk loci account for ∼33% of the familial risk of prostate cancer in populations of European ancestry, with these additional 23 loci, with effect sizes ranging from 1.06 to 1.14, explaining ∼3.1% of the familial risk (Online Methods). Based on a polygenic risk score comprising these 100 variants for men of European ancestry (Online Methods), the top 10% of men in the highest risk stratum have a 2.9 fold (95% CI 2.8-3.1) relative risk of prostate cancer and the top 1% of men have a 5.7 fold (95% CI 4.8-6.6) relative risk compared with the population average (Supplemental Table 17). The top 10% is at a RR compared with the average of the population where it will be important to examine whether targeted screening based on family history genetic risk may reduce the overdiagnosis of indolent disease, which is a main limitation of PSA screening. Our findings demonstrate the importance of conducting large-scale genetic studies in diverse populations for the discovery of novel risk loci which continue to provide novel insights into disease mechanisms for complex traits.

Online Methods
Primary genotype data were used from four prostate cancer GWAS in men of European ancestry (UK/Australia Stages 1 and 2; CAPS 1 and 2; BPC3 and Pegasus), and a ∼200K custom replication array (iCOGS), two GWAS in men of African ancestry (AAPC and Ghana Prostate Study), two GWAS in Japanese men (JAPC and BBJ) and a single scan in Latinos (LABC). 2,[4][5][6][7][8][10][11][12][13][14] (Supplementary Tables 1-3; Supplementary Information). Genotypes in all scans were imputed for ∼17 M SNPs/indels using the 1000 Genome Project (March 2012 release) as a reference panel. UK/Australia stages 1 and 2, CAPS 1 and 2, Pegasus, iCOGS, AAPC, Ghana Prostate Study, LABC and JAPC were imputed using IMPUTE V2. 30 BPC3, BBJ and Pegasus were imputed using Minimac. Betas and standard errors for each SNP were estimated stratified by study adjusting for principal components. In addition to analyses of overall prostate cancer risk, we performed secondary analyses of aggressive and early onset disease (age at diagnosis ≤55). Aggressive prostate cancer was defined as a Gleason score ≥8, disease stage as 'distant', a prostate-specific antigen (PSA) level >100 ng/ml, or death from prostate cancer. We included imputed data for SNPs with quality information scores >0.3 (IMPUTE V2) or with estimated correlation between the genotype scores and the true genotypes (r 2 )>0.3 (Minimac). We limited the analysis to SNPs/indels on chromosomes 1-22 as well as the X with minor allele frequency greater than 1%, except in iCOGS and Pegasus, which utilized arrays with coverage of less common alleles, where the MAF threshold was reduced to 0.5%.
Tests of homogeneity of the ORs across populations and study were assessed using likelihood ratio tests. Risk heterogeneity by disease aggressiveness and age was assessed using a case-only analysis. The associations between SNP genotypes and PSA level were assessed using linear regression, after log-transformation of PSA level to correct for skewness. Analyses were performed using SNPTEST, ProbABEL 31 , PLINK, Stata and an in-house C++ program (Supplementary Table 2). METAL was used to perform fixed effect ethnic-specific and multi-ethnic meta-analyses for overall prostate cancer, as well as secondary meta-analyses of aggressive and early-onset disease. 32

Inflation
We excluded SNPs with ± 500kb distance of any previously known prostate cancer risk locus and estimated the inflation for each study based on the 45 th percentile of the test statistic. The inflation was estimated to be 1.00 in the Latino, 1.03 in the African, 1.06 in the Japanese and 1.14 in the European ancestry studies, and, 1.07 in the European ancestry studies when SNPs at known risk loci and the iCOGS and UK2 studies were removed (see Supplementary Table 4). The inflation was converted to an equivalent inflation for a study with 1000 cases and 1000 controls (λ 1000 ) by adjusting by effective study size, namely where n k and m k were the number of cases and controls, respectively, for study k. Following the conversion the study-specific lambdas ranged from 0.995-1.083.

Genotyping and Concordance
The most significant associations in the meta-analyses were observed with imputed SNPs. To validate the accuracy of the imputed genotypes we genotyped each variant in ≥1847 samples (except rs9443189 and rs12051443 which were sequenced in 183 and 265 samples, respectively) that were included in the meta-analysis, and estimated the correlation between imputed and genotyped alleles. A correlation of ≥0.75 was used as the confidence threshold for imputation quality (Supplemental Table 6).

Functional Annotation
We used a number of publicly available prostate epithelia and prostate cancer ENCODE datasets of chromatin features to identify putative enhancer/regulatory regions at each risk locus. 33,34 The integration of chromatin bio feature annotations with the index SNPs and correlated markers (r 2 >0.9) from 1KGP EUR populations was performed using FunciSNP. 35 These datasets included LNCaP and RWPEI DnaseI HS sites (GSE32970) ENCODE; PrEC DNaseI HS sites (GSE29692) ENCODE; LNCaP CTCF ChIP-seq peaks (GSE33213) ENCODE; LNCaP H3K27ac and TCF7L2 (GSE51621) 33 , H3K4me3 and H3K4me1 histone modification ChIP-seq peaks GSE27823 36 ; FoxA1 ChIP-seq peaks (GSE28264) 37 ; Androgen Receptor (AR) ChIP-seq peaks 38 and AR binding sites (GSE28219) 39 ; NKX3-1 ChIP-seq peaks (GSE28264). 37 We also used the highly conserved set of predicted targets of microRNA targeting (miRcode 11, June 2012 release) 40 . To determine whether any of the putative functional SNPs potentially affect the binding of known transcription factors, position-specific frequency matrices were employed from Factorbook. 33,41

cis-eQTL analysis
Each risk locus is represented by an index SNP. For each index SNP, we retrieved all the correlated (r 2 ≥0.9) variants EUR populations from 1KGP. The genotypes of the correlated variants in 145 prostate tumor samples and 33 normal tissue samples were downloaded from TCGA database (Feb 2013). If a variant was not represented in the TCGA data, the genotypes were imputed using IMPUTE2. 30 A cis-eQTL analysis was performed for these variants and any transcript within a 1 Mb interval (500 kb on either side). Gene expression values were adjusted for somatic copy number and CpG methylation as previously described (ref. 42 ). Each risk variant was corrected for the number of transcripts in the interval. Significant associations were defined as a nominal p-value <0.05 and a false discovery rate <0.05 based on Benjamini-Hochberg method.
For the TMPRSS2 locus, we also used gene expression data generated from formalin-fixed paraffin embedded (FFPE) tissue in the Physicians' Health Study cohort. 43 RNA was extracted with the Agencourt Form a Pure FFPE kit (Beckman Coulter, Indianapolis, IN) and amplified using the WT-Ovation FFPE System V2 (NuGEN, San Carlos, CA). cDNA was hybridization on the GeneChip Human Exon 1.0 ST microarray (Affymetrics, Santa Clara, CA). The residuals were shifted to have the original mean expression values and normalized using the RMA method. 44,45 The SNP (rs1041449) was available in the BPC3 GWAS samples 6 ; 99 participants had both tumor expression and genotype data; 54 had both normal prostate expression and genotype data.

Determination of TMPRSS2-ERG fusion status
The TMPRSS2-ERG fusion was assessed in a subset of 552 cases from study samples of FHCRC, UKGPCS, TAMPERE, ULM and IPO-PORTO. The majority of cases were typed for TMPRSS2-ERG rearrangements on FFPE tumor materials using FISH techniques according to Summersgill, et al. 46 (for UKGPCS and FHCRC), Perner, et al. 47 (for ULM), or Saramaki, et al. 48 (for TAMPERE). The IPO-PORTO group applied qRT-PCR on RNA from fresh-frozen tumor tissues using a TaqMan gene expression assay (Hs03063375_ft, Life Technologies, Carlsbad, CA) for the fusion transcript T1G4, which is present in approximately 90% of all TMPRSS2-ERG positive prostate cancer.

Comparison of Number of Associated Loci among populations
We used the meta-analysis results from each population to evaluate the excess fraction of directionally consistent effect estimates (ORs) across populations, as evidence for additional shared susceptibility loci. We excluded the previously known prostate cancer risk regions as well as those identified in the current study (±500kb of index SNP) and compared the direction of association of SNPs defined in the European ancestry population with the other populations for several p-value thresholds. The p-values provided are based on a Chi-square binomial test for comparing proportions versus 50% chance to be in the same direction for each p-value cut-off.

Contribution to Familial Risk and Risk Stratification
The contribution of the known SNPs to the familial risk of prostate cancer, under a multiplicative model, was computed using the formula where λ 0 is the observed familial risk to first degree relatives of prostate cancercases, assumed to be 2, and λ k is the familial relative risk due to locus k, given by: where p k is the frequency of the risk allele for locus k, q k =1 − xp k and r k is the estimated per-allele odds ratio. 2 Based on the assumption of a log-additive model, we constructed a polygenic risk score (PRS) from the summed genotypes weighted by the per-allele log-odds ratios. 3 Thus for each individual j we derived:  Manhattan Plot of results from the multiethnic meta-analysis of overall prostate cancer risk. All SNPs within 500kb of known GWAS SNPs are omitted. The green line represents P=5×10 -8 . This figure shows all new variants with P<5×10 -8 , regardless of the confirmation results, as well as signals that were reported in the European meta-analysis that also reached 5×10 -8 in the multiethnic meta-analysis (see Table 1 and Supplementary  Regional plots of two novel genome-wide significant loci associated with prostate cancer risk. rs1041449/21q22 (TMPRSS2 region, left) and rs17694493/9p21 (CDKN2B-AS1region, right). Top: SNPs are plotted by their position 500kb on either side of the index SNP (purple diamond) on the chromosome against their association (-log 10 P) with prostate cancer from the multiethnic meta-analysis (rs1041449) and European meta-analysis (rs17694493). SNPs surrounding the index SNP are colored to indicate the local LD structure using pairwise r 2 data from the EUR panel of the 1000 Genomes (March 2012). MIDDLE: Significant peaks from TF and histone modification ChIP-seq experiments in the same genomic window (see Online Methods). All ChIP-seq in LNCaP unless otherwise indicated. BOTTOM: Genomic sequence (enclosed in black box) surrounding the SNP (red box) aligned to a LOGO graphic representing the proposed motif disruption. f Analysis limited to European ancestry populations as only small numbers of early onset cases (≤55 years) were available in the other populations.