The Impact of Inherited Genetic Variation on DNA Methylation in Prostate Cancer and Benign Tissues of African American and European American Men

Prostate cancer is the second most common cancer and cause of cancer death among men in the United States (1). African American (AA) men are disproportionately affected by prostate cancer, with an incidence rate that is 1.7 times higher (2–4) and a mortality rate that is two to four times higher than European ancestry (EA) men (5, 6). The causes of these disparities are likely complex, with social, environmental, and genetic factors contributing (3).

Genome-wide association studies (GWAS) have identified 269 common risk alleles (SNP) associated with prostate cancer susceptibility, and these account for ∼43% of the familial relative risk for prostate cancer (7, 8). More than 90% of prostate cancer-risk SNPs are located in noncoding regions, suggesting that the causal variants underlying these associations impact gene regulation.

Despite the progress in identifying susceptibility loci, the biological mechanisms by which SNPs impact prostate cancer risk are largely unknown. A common approach for understanding the regulatory mechanisms of disease-associated variants is to assess their association with local gene expression and/or epigenetic features (9, 10). Regions where SNPs affect gene expression and/or DNA methylation (DNAm) are known as expression quantitative trait loci (eQTL) and methylation QTLs (meQTL), respectively. Few studies have identified meQTLs in prostate tissue, and those studies lack adequate representation of individuals of African ancestry. One study reported 7,590 genome-wide cis-meQTLs in prostate cancer tumor samples (11); another focused on 147 prostate cancer-risk SNPs and identified 93 prostate cancer-risk SNPs that were associated with DNAm at nearby CpG sites in prostate cancer tumor tissue (12).

In this study, we attempt to improve our understanding of the mechanisms by which prostate cancer susceptibility SNPs influence prostate cancer biology, by examining these SNPs’ association with DNAm at nearby CpG sites across the genome, in both benign and cancerous prostate tissue. To identify mechanisms that may be relevant to disparities, we analyze data from both AA and EA patients with prostate cancer and conduct analyses stratified by ancestry. We leverage existing GWAS summary statistics to determine whether GWAS and cis-meQTL association signals (identified in ancestry-specific analyses) are likely to share a common causal variant (using colocalization methods). The identification of colocalizing meQTLs and GWAS loci can provide insights on the epigenetic mechanisms by which SNPs influence prostate cancer risk.

Subjects included in this work were male patients with prostate cancer who underwent robotic-assisted laparoscopic radical prostatectomy at University of Chicago Medical Center (UCMC) between 2011 and 2017. Trained interviewers from the Epidemiology Research Recruitment Core consented 74 AA and 74 EA eligible men for the collection of questionnaire data, prostate tissue, and access to medical records. All eligible participants were diagnosed with Gleason scores of at least seven. All enrolled patients provided written informed consent. This study was conducted in accordance with recognized ethical guidelines (U.S. Common Rule) and was approved by the Institutional Review Board of The University of Chicago.

After surgery, prostate specimens were sent to the Human Tissue Resource Center (HTRC) at the University of Chicago. Each prostate specimen underwent histologic examination and Gleason grading by a genitourinary pathologist (G.P. Paner) at the University of Chicago. The presence of adenocarcinoma was confirmed by the overexpression of alpha-methylacyl-coenzyme-A racemase (AMACR), and areas to sample for DNA extraction were marked. Benign tissue and tumor tissue were collected.

The following criteria were used for the selection of benign tissue from FFPE tissue of the resected prostate: (i) samples were selected from blocks that were free of tumor; (ii) tissue from the peripheral zone was prioritized; (iii) if no tumor-free areas were available in the peripheral zone, tissue was collected from the central zone. In 10 cases, neither the peripheral nor the central zone were suitable for sampling, so BPH tissue was collected. Cancer tissue was selected from the index lesion in multifocal tumors. Tissue collection was performed using either a 1 mm biopunch or by laser capture microdissection of ∼100 μm² of tissue (8-μm-thick sections using a Leica LMD 6500 system). In cases where two consecutive diagnostic blocks showed acceptable areas to sample (continuous benign or tumor tissue running through two blocks from base to apex), tissue was collected by punching through the base-most diagnostic block.

The Gentra Puregene Tissue Kit (QIAgen) was used to extract 1 to 2 μg of DNA from benign and cancerous prostate samples. We assessed DNA concentration and quality with the NanoDrop and Agilent BioAnalyser. We excluded DNA samples with a concentration <40 ng/μL or 260/280 ratio outside the range of <1.6 to ≥2.1 and/or fragmented DNA <2 kb. The Illumina Infinium HD FFPE Restoration Kit was used to restore FFPE DNA (according the manufacturer's protocol, including qPCR-based quality check of the DNA samples prior to restoration).

Genome-wide SNP data were generated (from benign tissue DNA) using the Illumina Infinium Multi-Ethnic Global-8 v1.0 array at the University of Chicago Genomics Core Facility. Genotypes were called using a GenCall Threshold 0.15. Genotype data consisted of 148 individuals (74 AA and 74 EA) with 1,707,345 autosomal SNPs measured. We excluded 224,774 SNPs with low call rates (<99%). Minor allele frequency (MAF) and Hardy–Weinberg equilibrium (HWE) thresholds were applied separately for AA and EA samples. For AA samples, we excluded 855,033 SNPs with MAF <0.05 and 53 SNPs with HWE P-values <10⁻⁵, resulting in 627,485 high-quality SNPs. For EA samples, we excluded 932,143 SNPs with MAF <0.05 and 23 SNPs with HWE P-value <10⁻⁵, resulting in 550,405 high-quality SNPs.

We performed imputation using the Haplotype Reference Consortium (HRC, Version r1.1 2016) panel, which includes all samples from the 1,000 Genomes Project (phase III), using the Michigan Imputation Server. Of the 627,485 and 550,405 post-quality check SNPs in AA and EA, respectively, 563,473 SNPs in AA and 498,015 SNPs in EA matched the HRC panel and met the quality check thresholds of the Michigan Imputation Server. A total of ∼39 million SNPs were imputed in each group. In AA, we excluded ∼23M SNPs with imputation accuracy (r²) ≤0.3 and ∼9.5M SNPs with MAF ≤0.07, resulting in 6,463,658 SNPs. In EA, we removed ∼29.3M SNPs with r² ≤0.3 and ∼4.9M SNPs with MAF ≤0.07, resulting in 4,900,500 SNPs. We conducted downstream analyses on the resulting 6.4M SNPs for AAs and 4.9M SNPs for EAs.

The Illumina Infinium MethylationEPIC array was used to interrogate >850,000 CpG sites at the University of Chicago Genomics Core Facility. This array provides dense coverage across CpG islands (>95%), shores (>80%), and shelves (>90%). Methylation data were normalized using the BMIQ function in the ChAMP software, and methylation at each CpG was expressed as β values ranging from 0 (completed unmethylated) to 1 (completed methylated). We processed the methylation data for EA and AA groups combined, but stratified by tissue type. We removed probes with a detection P-value >0.01 (95,481), beadcount <3 in at least 5% of samples (106), non-CpG probes (2,228), nonspecific probes that align to multiple locations (47), and underperforming probes (e.g., low mappability to hg38, unrecognized color channel switch for type I probes, and contain SNPs close to the 3′ end of probe; 69,244; ref. 13). This quality check resulted in 698,812 CpG sites for benign tissue and 682,694 in tumor tissue, and these CpGs were included in all downstream analyses.

We characterized differences in DNAm between tumor and benign samples by applying principal components analysis (PCA) to all CpGs for all 296 samples. PCA demonstrated clear separation of most tumor samples from most benign samples (Supplementary Fig. S1). Within tumor samples, PC1 was also associated with Gleason score (Supplementary Fig. S1). These results demonstrate that tumor/benign status is the largest source of variation in our DNAm data, and our tumor samples vary with respect to the composition/abundance of tumor versus benign cells.

We used FastQTL to conduct genome-wide cis-meQTL analyses. All analyses were conducted separately for AA benign, EA benign, AA tumor, and EA tumor tissue. We tested the local (cis) association of SNPs and CpG sites <500 kb apart. Methylation β values were rank normalized to satisfy linear model assumptions. To identify CpGs affected by a meQTL, we computed CpG-level empirical P values (i.e., the smallest P-value for each CpG) using an approximation of the β-distribution and adaptive permutations as implemented in FastQTL (–permute 1000 10000; ref. 14). To account for multiple testing, we used the Storey/Tibshirani FDR procedure applied at the CpG-level to identify mCpGs (implemented in the R/qvalue package). A FDR of 0.05 was applied in each of the four analyses.

All regression models were adjusted for age, five genotyping principal components (PC), and 10 methylation surrogate variables (SV), to capture variability in cell-type composition (including tumor purity) as well as potential technical variation.

Methylation SVs were estimated using the SVA package (15). To determine the number of SVs that maximized power for meQTL detection, we conducted cis-meQTL analysis of chromosome 1 for each ancestry and tissue type using 5, 10, 15, and 20 SVs.

PC analyses were conducted (separately for AA and EA) using a genome-wide set of independent SNPs in PLINK (–indep-pairwise 50 5 0.2). An additional PC analysis was conducted for AA and EA patients combined to demonstrate that and AA and EA individuals clustered separately, as expected (Supplementary Fig. S2).

An overview of our workflow (and results) is shown in Fig. 1. We first identified meQTLs for benign (6,298 for AA and 6,960 for EA) and tumor tissue (2,641 for AA and 1,700 for EA) and identified the lead SNP for each meQTL. Next, we restricted to the 20,646 SNPs from prostate cancer GWAS that met a threshold of P < 5 × 10⁻⁸ (from Schumacher and colleagues) and identified SNPs that were also a lead meQTL SNP, resulting in 37 AA and 40 EA colocalization candidates for benign tissue and 38 AA and 33 EA candidates for tumor tissue. We conducted colocalization tests (as described in the sections below) at loci where the lead GWAS SNP was in linkage disequilibrium (LD, r² > 0.5) with the lead meQTL SNP in at least one ancestry (8 for AA and 17 for EA in benign, 3 for AA and 14 for EA in tumor), according to LDlink (https://ldlink.nci.nih.gov/). For AA men we used the Americans of African Ancestry in the USA (ASW) reference population, and for EA men we used the European reference populations (EUR). We conducted all GWAS-meQTL colocalization analyses stratified by ancestry and tissue type, and identified one (AA) and seven (EA) colocalizations for benign meQTLs and three (AA) and five (EA) colocalization for tumor meQTLs (as described in the Results section). For the meQTLs identified, mCpG enrichment in genomic features was assessed using X² tests.

We used summary statistics from a prior prostate cancer GWAS (7) and our cis-meQTL results to conduct Bayesian co-localization analysis (using coloc; ref. 16). This approach restricts to SNPs that are present in both sets of summary statistics within a given start and end positions (base pairs) for the region under analysis.

Coloc requires three prior probabilities: the probability a variant is causal for prostate cancer only (p₁), a QTL only (p₂), and both prostate cancer and a QTL (p₁₂). The software uses these priors to calculate the posterior probability of a common causal variant [P(CCV)]. To select a prior for prostate cancer, we used a recent estimate of the number of independent common prostate cancer susceptibility variants (4, 530; ref. 17). Given the ∼20M SNPs tested in recent GWAS, the probability a SNP is causal for prostate cancer is 4,530/20M, approximately 10⁻⁴. This probability is equal to p₁ + p₁₂. To set meQTL priors, we used the previously detected 7,590 cis-meQTLs among 4,894,225 SNPs tested (11), indicating the probability a SNP is a causal meSNP in prostate tissue is approximately 5 × 10⁻³. This probability corresponds to p₂ + p₁₂ (for GWAS-meQTL colocalization). Selection of these priors was informed by the literature but not intended to be exact. Because the true value of p₁₂ is unknown, we varied the value of p₁₂ to correspond to probabilities of a causal prostate cancer SNP being a causal meSNP of 10%, 25%, 50%, and 75%, similar to prior meQTL studies (18, 19). The resulting values for p₁, p₂, and p₁₂ are shown in Supplementary Table S1.

We searched for eQTLs among the GWAS-meQTLs (meSNPs) that reached the colocalization threshold of P(CCV) >80% using the Genotype-Tissue Expression Project (GTEx v8; ref. 20) prostate tissue eQTL results (n = 221). To prepare the eQTL and GWAS summary statistics for coloc, we restricted to the common SNPs and the start/end genomic positions defined by the meQTL analysis. Genomic coordinates for GTEx SNPs were converted from GRCh38/hg38 to GRCh37/h19 using the NCBI Genome Remapping Service.

GWAS-eQTL colocalization was performed only for loci with strong evidence of GWAS-meQTL colocalization [P(CCV) of >80%]. We conducted analyses for GWAS-eQTL pairs using coloc. GTEx identified 7,356 eQTLs in prostate tissue among 11.5M SNPs, suggesting the probability that a SNP is a causal eSNP in prostate tissue is ∼6 × 10⁻⁴. This probability corresponds to p₂ + p₁₂ (for GWAS-eQTL colocalization). We also varied the value of p₁₂ to correspond to probabilities of a causal GWAS SNP being a causal eSNP of 10%, 25%, 50%, and 75% (Supplementary Table S1).

The raw data for this study were generated at University of Chicago, and these data are publicly available through dbGaP (dbGaP Study Accession No.: phs003516.v1.p1). All meQTL summary statistics are available for download from Zenodo (DOIs: 10.5281/zenodo.10304061 and 10.5281/zenodo.10358156).

Characteristics of the 74 AA and 74 EA patients included in our analyses are described in Supplementary Table S2. AA and EA patients were on average 65 years old at diagnosis. Gleason score and tumor volume (defined as the percent of prostate that is tumor) were slightly higher in EAs compared with AAs, but PSA was similar. The slightly less favorable clinical characteristics of EA compared with AA patients may be due to EA men with more advanced prostate cancer residing outside the UCMC catchment area seeking care at the University of Chicago, with our AA patients being more likely to reside in the catchment area.

We identified a cis-meQTL for 6,298 and 6,960 CpGs (FDR 0.05) in benign prostate tissue of AA and EA men, respectively (Table 1). The meQTLs detected were represented by 5,855 unique lead SNPs in AA and 6,496 unique lead SNPs in EA, as some mCpGs had the same lead SNP. The CpGs impacted by meQTLs (i.e., mCpGs) were enriched in non-CpG islands (open sea) and depleted in islands compared with all measured CpGs (Supplementary Fig. S3A). 64% (4,060) and 63% (4,409) of mCpGs in AA and EA were assigned to a gene (based on Illumina annotations) as compared with 73% of all measured CpGs (P < 0.001; Supplementary Fig. S4A). mCpGs were enriched in enhancer regions (P < 0.001) for both AAs (4%) and EAs (4%) compared with all measured CpGs (3%) and depleted in promoters (Supplementary Fig. S4A). In AAs, mCpGs were enriched in DNase hypersensitivity sites (DHS) regions (67%) compared with all CpGs measured (60%; P < 0.001; Supplementary Fig. S4A).

Of the 6,298 cis-meQTLs detected in AA benign tissue, 4,269 (68%) had P < 0.01 and 3,865 (61%) had P < 0.001 in EA benign tissue. The correlation (r) between the β coefficients from AA and EA for the 4,269 cis-meQTLs was 0.97, and 4,258 (68%) were directionally consistent and considered replicated (Supplementary Fig. S5A). Similarly, among the 6,960 cis-meQTLs detected in EA benign tissue, 4,372 (63%) had P < 0.01 and 3,576 (51%) had P < 0.001 in AA benign tissue. The correlation between the β coefficients for the 4,372 cis-meQTLs was 0.95, and 4,343 (62%) were directionally consistent and considered replicated (Supplementary Fig. S5B).

In tumor tissue we identified 2,641 (AA) and 1,700 (EA) cis-meQTLs (FDR < 0.05; Table 1). Tumor mCpGs showed somewhat weaker depletion in islands compared with benign (Supplementary Fig. S3B). In tumor, 69% (AA) and 67% (EA) of mCpGs were assigned to genes compared 73% of all measured CpGs (P < 0.001; Supplementary Fig. S4B). Tumor mCpGs were depleted in promoters (13%) and enriched in DHS regions (67%) compared with all measured CpGs (20% and 61%, respectively, both P < 0.001).

Among the 2,641 cis-meQTLs (FDR <0.05) detected in AA tumor tissue, 1,667 (63%) had P < 0.01 and 1,457 (55%) had P < 0.001 in EA tumor tissue. The correlation among the β coefficients from AA and for the 1,667 cis-meQTLs was 0.96, and 1,659 (63%) were directionally consistent and considered replicated (Supplementary Fig. S5C). Among the 1,700 cis-meQTLs (FDR <0.05) detected in EA tumor tissue, we found 1,131 (67%) at P < 0.01 and 996 (59%) at P < 0.001 in AA tumor tissue. The correlation between the β coefficients in EA versus AA was 0.96 among the 1,131 cis-meQTLs and 1,126 (43%) were replicated with directional consistency (Supplementary Fig. S5D).

Among the 6,298 benign-tissue cis-meQTLs detected in AAs, we replicated 62% (P < 0.01) in AA tumor tissue (Supplementary Table S3). Similarly, for EA, we replicated 55% benign-tissue cis-meQTLs in tumor tissue (P < 0.01). We observed a larger percent of replication of tumor-tissue cis-meQTLs in benign tissue (79% at P < 0.01 for both AA and EA; Supplementary Table S3).

We attempted replication of the 7,590 cis-meQTLs identified previously among 589 localized prostate tumors of EA men (Supplementary Table S4; ref. 11). In AA benign and tumor tissue, 48% and 44% of the 5,231 meQTLs for which we had data were replicated (P < 0.01) with directional consistency. In EA benign and tumor tissue, we observed more replication of specific CpG-SNP pairs compared with AAs: 59% and 49% of the 5,590 previously reported meQTLs were replicated (P < 0.01) with directional consistency (Supplementary Table S4).

In addition to genome-wide meQTLs, Houlahan and colleagues (11) report 75 prostate cancer-risk cis-meQTLs (52 validated in an independent cohort), which describe the association between 27 unique prostate cancer-risk loci and 73 CpG sites. We attempted replication for 41 and 45 of the 52 validated cis-meQTLs in AA and EA tissues, respectively (Supplementary Table S5). In AA benign and tumor tissue samples, we replicated 24% and 41% of the 41 meQTLs (P < 0.01 with directional consistency), respectively. In EA benign and tumor tissue samples, we replicated 29% and 53% of the 45 meQTLs (P < 0.01 with directional consistency), respectively.

We also sought to replicate 110 cis-meQTLs identified in the prostate tumor tissues of 355 EA men (Supplementary Table S6; ref. 12). Of the 110 cis-meQTLs, we had summary results for 60 cis-meQTLs in AAs and 67 cis-meQTLs in EAs. We replicated 15% and 42% of the meQTLs (P < 0.01 with directional consistency) in AA benign and tumor tissues, and 28% and 49% of the meQTLs in EA benign and tumor tissues.

Among our 6,298 (AA) and 6,960 (EA) meQTLs identified in benign tissue, we identified cis-meQTLs residing in the same location as prostate cancer-risk loci (Fig. 1; Supplementary Tables S7–S8). Restricting to lead SNPs with P < 5 × 10⁻⁸ in the Schumacher and colleagues prostate cancer GWAS, we searched for meQTL lead SNPs with an LD (r² > 0.5) with a GWAS lead SNP in AA and/or EA. We identified 20 such SNP–CpG pairs (Supplementary Table S9).

We found evidence of colocalization [P(CCV) >80%] (based on 50% prior probability that a GWAS SNP is an meSNP, see Supplementary Table S1) for two GWAS–meQTL pairs in AA men and four pairs in EA men (Table 2; Supplementary Table S10 for full results). The colocalized signals detected in AA men were in an intergenic region near MLPH and in a promoter of HAUS6 (Fig. 2A and B). The signal in HAUS6 was also suggestive in the EA cohort. The four colocalized signals detected in EA men were located in the gene body of TNS3, promoter near MSMB, promoter near MRPL52/MMP14, and intergenic region of COPRS/UTP6 (Fig. 2C–F). The signals in TNS3 and COPRS are also suggestive in AA and highlight the potential impact of cross-population LD differences on colocalization results. The posterior probabilities that the GWAS and meQTL signal are caused by independent causal variants are shown in Supplementary Table S11.

In tumor tissue, we identified cis-meQTLs located at previously reported prostate cancer-risk loci (P < 5 × 10⁻⁸) in 38 regions (for AA) and 33 regions (for EA), respectively (Fig. 1; Supplementary Tables S12–S13). Restricting to loci with LD (r²) > 0.5 between the lead GWAS and lead meQTL SNP in at least one ancestry, we identified 16 loci to test for colocalization (Supplementary Table S14).

There was evidence of colocalization (based on 50% prior probability that a GWAS SNP is a meSNP, Supplementary Table S1) for three GWAS–meQTL pairs in both AAs and EAs (Table 3; Supplementary Table S15 for full results). Two colocalization signals were shared between AAs and EAs; one near the gene body of IRX4 and the other in the 5′UTR of MMP7 (Fig. 3A and B). In EA, prostate cancer-risk SNP rs12653946 was associated with nine CpGs in the IRX4 region (Supplementary Fig. S6), most showing strong colocalization evidence, including several promoter CpGs. In AA men, evidence of colocalization was observed for a CpG near the gene body of MYO9B (Fig. 3C). This region showed suggestive evidence of colocalization in EA (Table 3). One additional colocalized signal identified in EA was located in the gene body of TNS3 (Fig. 3D), a signal also observed in EA benign tissue and suggestive AA benign and tumor tissue. The posterior probabilities that the GWAS and meQTL signals are caused by independent causal variants are shown in Supplementary Table S16.

Among the nine colocalized GWAS–meQTL pairs across both tissue types and ancestries, we identified eQTLs for four prostate cancer-risk SNPs (rs12653946, rs11568818, rs11666569, and rs10993994) in GTEx normal prostate tissue. We conducted colocalization analyses for the four prostate cancer-risk SNPs and the corresponding eGenes (IRX4, MMP7, MYO9B, and MSMB, respectively). Using the prior sets listed in Supplementary Table S1, we found strong evidence of shared common causal variants affecting both GWAS and eQTL traits for all four SNP-eGenes tested (Supplementary Table S17; Supplementary Fig. S7), suggesting these prostate cancer-risk variants regulate both local methylation and expression.

We performed a genome-wide search for cis-meQTLs in both benign and cancerous prostate tissue of AA and EA men. We identified 6,298 and 6,960 cis-meQTLs in the benign tissue of AA and EA men, respectively, and 2,641 and 1,700 cis-meQTLs in the tumor tissue of AA and EA men, respectively. To determine if known prostate cancer susceptibility loci (7) impact local DNA methylation in the prostate, we used Bayesian colocalization methods to identify nine regions in which prostate cancer risk loci likely share a causal variant with a meQTL in either cancer and/or benign tissue. Four of the nine regions showed evidence of a colocalizing eQTL (IRX4, MMP7, MYO9B, and MSMB), all of which have roles in tumorigenesis (7, 12, 21–23). Our findings highlight potential regulatory mechanisms by which prostate cancer-risk variants impact gene regulation in prostate tissue.

We identified six regions with strong evidence of GWAS-meQTL colocalization in benign tissue (two in AA, four in EA) and four in tumor tissue (two in both AA and EA, one in AA, and one in EA). According to QTLbase (24) and GTEx, only three of our nine colocalized meSNPs are associated with methylation and/or expression in whole blood (rs2292884-cg14458575: eGene = MLPH, rs10993994-cg17030820: eGene = MSMB, rs11666569-cg19418318: eGene = MYO9B; Supplementary Table S18), indicating that some prostate cancer susceptibility SNPs may have prostate-specific regulatory effects.

Interestingly, we observed colocalization for the MLPH region (for AA benign tissue) despite distinct differences in LD between AA (meQTL) and EA (GWAS) in this region (Fig. 2A). The lead meQTL SNP for AA (rs72620822) was among the top GWAS SNPs (Supplementary Fig. S8), reflecting LD between these SNPs in populations of European ancestry (r² = 0.65 in EUR). However, the lead GWAS SNP (rs2292884) is not among the top AA meQTL SNPs (Fig. 2A), reflecting lack of LD between these SNPs in populations of African ancestry (r² = 0.002 in AFR). Evidence for colocalization for the EA meQTL in the MLPH region is much weaker. This this locus requires further investigation in future studies.

Our results suggest that the causal prostate cancer-risk variant represented by rs12653946 affects methylation in tumor tissue at (at least) nine CpGs near IRX4 in EA men (Supplementary Table S14; Supplementary Fig. S6). Six of these CpGs are in a CpG island near the IRX4 start site, and one is in an enhancer (Supplementary Table S14). The rs12653946 risk allele (T) was associated with increased methylation at all seven CpGs within islands (Supplementary Table S14). The meQTL at IRX4 appears substantially weaker (or absent) in benign tissue (P > 0.001 for all nine CpGs) compared with tumor (Supplementary Table S19; Supplementary Fig. S9). In GTEx normal prostate tissue, the rs12653946 risk allele (T) is associated with decreased gene expression of IRX4 (Supplementary Table S17), a previously reported tumor suppressor gene (25). Additional prostate cancer-risk variants showing clear association with both DNAm and gene expression (in GTEx prostate tissue) include rs11568818 (MMP7), rs11666569 (MYO9B), and rs10993994 (MSMB). For MMP7 and MYO9B, the risk allele decreases DNAm (in 5‘ UTR) and increases gene expression. For MSMB, the risk allele increases DNAm (near the TSS) and decreases MSMB expression.

A primary challenge of conducting colocalization using meQTL results from AA men is the LD mismatch with the largely European ancestry GWAS of prostate cancer, which decreases power to detect colocalization using meQTL results from AAs. For example, in AA benign and tumor tissue, we identified a meQTL (cg23694490-rs834603) near TNS3 but the lead meQTL SNP was in low LD (r² = 0.15) with the lead GWAS SNP (rs56232506), whereas in EA men these SNPs were in strong LD (r² = 0.86) and colocalization was observed. Thus, evidence of colocalization was lower for AA likely due LD differences in this region (Fig. 2C). This LD discordance also likely affects the generalizability of associations reported in prior prostate cancer GWAS, as less than half of prostate cancer-risk SNPs (identified in studies of largely European ancestry men) have been replicated in men of African ancestry (26). These examples highlight the need for large prostate cancer GWAS focusing on AAs. Such results would improve our ability to describe the functional impacts of prostate cancer-risk alleles using AA QTL results.

We detected nearly triple the number of cis-meQTLs in benign tissue compared with tumor tissue (in both AA and EA), a finding likely attributable to differences in cell-type composition. Benign tissue samples may be more homogeneous, with patient samples consisting largely of epithelial and stromal cells. In contrast, tumor tissue samples will contain cancer cells, which can differ in various cellular phenotypes across individuals (and within individuals), as well as adjacent normal (epithelial and stromal cells). Assuming meQTL patterns differ to some extent across cell types, power to detect cell type-specific QTLs will be higher in more homogeneous cell-type mixtures. Despite these differences in cell type, we found that 79% of tumor-tissue meQTLs were likely present in benign tissue (P < 0.01) but observed less replication of benign-tissue meQTLs in tumor tissue (62% in AAs, 55% in EAs).

Several of the genes/regions identified in this work have established roles in prostate cancer. For example, studies in MMP7 knockout mice show decreased proliferation, increased apoptosis, and inhibition of angiogenesis and epithelial-to-mesenchymal transition, leading to decreased tumor burden (22). In addition, increased levels of serum MMP7 were recently reported to associate with reduced overall survival in castration-resistant prostate cancer (27). IRX4 is a tumor suppressor gene (25), and has been identified as a target of an eQTL associated with prostate cancer risk (21), with alternatively spliced transcripts potentially playing an important role in prostate cancer prognosis (28). MSMB has been previously identified as a target of a prostate cancer-risk eQTL, with downstream/trans effects on SNGH11 (23). MSMB expression is reduced in tumors and adjacent benign prostate tissue, potentially explaining why serum MSMB are decreased in patients with prostate cancer (29).

Our study has several strengths compared with prior work. First, we employed colocalization tests, which are critical for determining if causal variants for prostate cancer-risk also affect DNAm (and are not simply associated with DNAm due to LD with a nearby meSNP). Although our work appears to suggest that the number of prostate cancer-risk loci potentially explained by meQTLs is smaller than previously reported, we are underpowered to detect weaker meQTL effects, and sample size likely prevented detection of colocalization for some meQTLs, such as in the EA tumor tissue meQTL (cg02493740-rs10187424) near GGCX (Supplementary Table S15; Supplementary Fig. S10). Second, we analyzed equal numbers of AA and EA samples, emphasizing detection of meQTLs in patients of African ancestry. Third, we conducted meQTL analyses using paired tumor and benign samples; however, it is possible that these benign samples may contain cancer-related molecular characteristics and are not representative of prostate tissue from prostate cancer-free men. Finally, although we did not generate gene expression data, we leveraged eQTL results from GTEx normal prostate tissue.

Our study utilized FFPE tissues as a DNA source, and we used a qPCR-based quality check procedure to ensure the quality of DNA samples were appropriate for restoration using the Illumina Infinium HD FFPE Restoration Kit. Prior studies show high concordance of EPIC array DNA methylation measures obtained from FFPE samples (repaired using the Illumina Kit) compared with those from fresh–frozen tissue from the same individuals, reporting correlations >90% (30–32). However, the potential for decreased quality of FFPE samples, and increased noise in the DNAm data generated, may reduce power for meQTL detection compared with frozen tissue.

We adjusted for global ancestry in our meQTL analyses, but not local ancestry. It has recently been shown that local ancestry adjustment can improve power for QTL detection in admixed samples (33), but power gains are likely to be modest (34).

In summary, we conducted a comprehensive search for meQTLs in tumor and paired benign prostate tissues of AA and EA men. We used summary statistics from prior GWAS to identify prostate cancer susceptibility loci whose biological mechanism likely involves alteration of local DNAm. These results are a resource to explore differences in prostate meQTL profiles of AA and EA men. To better understand the biological mechanisms by which prostate cancer susceptibility SNPs influence prostate cancer biology, larger and more diverse studies and meta-analyses of prostate meQTLs and eQTLs are needed, as are larger prostate cancer GWAS of AA men.

K.J. Gleason reports other support from AbbVie outside the submitted work. B.L. Pierce reports grants from Department of Defense during the conduct of the study. No disclosures were reported by the other authors.

D. Delgado: Formal analysis, writing–original draft. M. Gillard: Resources, data curation, formal analysis. L. Tong: Writing–review and editing. K. Demanelis: Methodology. M. Oliva: Methodology. K.J. Gleason: Methodology. M. Chernoff: Data curation, formal analysis, methodology, writing–review and editing. L. Chen: Supervision, methodology. G.P. Paner: Formal analysis, investigation, methodology. D. Vander Griend: Conceptualization, resources, supervision, funding acquisition, investigation, methodology, writing–review and editing. B.L. Pierce: Conceptualization, resources, supervision, funding acquisition, investigation, methodology, writing–review and editing.

U.S Department of Defense CDMRP Health Disparity Research Award (W81XWH-14–1-0529 to B.L. Pierce), National Institute of Environmental Health Sciences award (R35ES028379 to B.L. Pierce and F30ES031858 to M. Chernoff), National Institute of General Medicine (T32GM150375 and T32GM007281 to Marcus Ramsey Clark), NCI (P30CA14599 to Kunle Odunsi), Susan G. Komen Research Training Grant (GTDR16376189 to Eileen Dolan), and the National Institute of Aging (T32AG51146 to David O. Meltzer).

This work was supported by the Epidemiology Research and Recruitment Core, the Human Tissue Research Core, and the Genomics Core Facility at the University of Chicago. This publication was neither originated nor managed by AbbVie, and it does not communicate results of AbbVie-sponsored Scientific Research. Thus, it is not in scope of the AbbVie Publication Procedure (PUB-100).