Abstract
Colorectal cancer has a strong epigenetic component that is accompanied by frequent DNA methylation (DNAm) alterations in addition to heritable genetic risk. It is of interest to understand the interrelationship of germline genetics, DNAm, and colorectal cancer risk.
We performed a genome-wide methylation quantitative trait locus (meQTL) analysis in 1,355 people, assessing the pairwise associations between genetic variants and lymphocytes methylation data. In addition, we used penalized regression with cis-genetic variants ± 1 Mb of methylation to identify genome-wide heritable DNAm. We evaluated the association of genetically predicted methylation with colorectal cancer risk based on genome-wide association studies (GWAS) of over 125,000 cases and controls using the multivariate sMiST as well as univariately via examination of marginal association with colorectal cancer risk.
Of the 142 known colorectal cancer GWAS loci, 47 were identified as meQTLs. We identified four novel colorectal cancer–associated loci (NID2, ATXN10, KLHDC10, and CEP41) that reside over 1 Mb outside of known colorectal cancer loci and 10 secondary signals within 1 Mb of known loci.
Leveraging information of DNAm regulation into genetic association of colorectal cancer risk reveals novel pathways in colorectal cancer tumorigenesis. Our summary statistics-based framework sMiST provides a powerful approach by combining information from the effect through methylation and residual direct effects of the meQTLs on disease risk. Further validation and functional follow-up of these novel pathways are needed.
Using genotype, DNAm, and GWAS, we identified four new colorectal cancer risk loci. We studied the landscape of genetic regulation of DNAm via single-SNP and multi-SNP meQTL analyses.
Introduction
Large genome-wide association studies (GWAS) have identified 142 genetic variants associated with colorectal cancer risk to date (1, 2). Although large GWAS have been essential for variant discoveries, novel information can be gained by combining GWAS with functional studies of genomic characteristics. Studying association of genetic variants with molecular features such as gene expression and DNA methylation (DNAm) may reveal underlying disease mechanisms and discover additional novel loci. With the availability of genomic studies, results from large-scale GWAS studies, and summary statistics-based methods, it is now feasible to study these relationships. In this paper, we examine the impact of genetic variants on DNAm and utilize this information to study the association of methylation quantitative trait loci (meQTL) with colorectal cancer risk.
DNAm refers to the addition of a methyl group to a cytosine phosphate guanine (CpG) site of a DNA strand. It plays a critical role in colorectal carcinogenesis (3–6) and is associated with survival time in colorectal cancer (7, 8). In addition, a well-studied subtype of colorectal cancer is the CpG Island Methylator Phenotype (CIMP), which is defined by the degree of methylation at specific CpG islands in the tumor genome (9–11). Other measures of DNAm and epigenetic modifications have also been used as screening and early detection, such as ColoGuard (NDRG2 and BMP3) and Epi ProColon (SEPT9; ref. 12).
The interrelationships between DNAm, germline genetics, and complex phenotypes have been examined previously. Freytag and colleagues (13) developed a genetic predictive model for DNAm at each CpG site using elastic net penalization and assessed the association of predicted DNAm with complex phenotypes, including schizophrenia, Alzheimer's disease, and irritable bowel disease. Yang and colleagues developed genetic predictive weights of DNAm within the Framingham Heart Study and used these weights to assess associations of predicted DNAm with epithelial ovarian (14) and breast cancers (15). However, to our knowledge, there has been no similar study relating genetically regulated methylation with colorectal cancer risk.
Many methods have been developed to evaluate how genetic variants that impact a genomic characteristic are associated with an outcome. Typically, these methods assess the associations between SNPs and a genomic feature and then use the resulting effect sizes as weights to impute the genomic feature into a much larger genetic data set for association with the outcome (16–19). This has been expanded to use only GWAS summary statistics, alleviating both computational and data accessibility burdens. These approaches serve as a vital tool in leveraging genomic information for novel discoveries. They are closely related to two-stage least-squares Mendelian randomization approaches, when the genetic variants affect the outcome only through the genomic feature. Mixed-effects score test (MiST; refs. 19, 20) and the summary statistics-based sMiST (21) allow for SNP effects through not only the predicted genomics features such as in Mendelian randomization but also alternative pathways by modeling the residual genetic effects. By this, MiST or sMiST tests the total effect of SNPs on the outcome.
In this paper, we evaluated the relationship between germline SNPs and DNAm based on the GWAS and methylation data of 1,355 study participants by both single-SNP and multi-SNP meQTL analyses. We then assessed the association of meQTLs with colorectal cancer risk by utilizing summary statistics from over 125,000 colorectal cancer cases and controls.
Materials and Methods
Description of study with DNAm and GWAS data
The DNAm and GWAS data were from the Ontario Familial Colon Cancer Registry (22). Cases were identified from the Provincial Cancer Registry between 1997 and 2000, and controls were recruited randomly via household telephone interviews (22, 23). Informed consent was obtained from participants. All study protocols were approved by the relevant Institutional Review Boards, and informed consent was obtained from study participants in accordance with the Helsinki Accord. After quality control (QC) and filtering, a total of 1,355 cases and controls had both genome-wide genotyping and methylation data available. All participants were of European descent as confirmed by GWAS data. DNAm was measured in lymphocytes. Detailed information on DNAm preprocessing and QC is provided in Supplementary Methods. After QC, there were 298,248 CpG sites. The majority of cases and controls were plated separately, except for 267 study participants (115 cases and 152 controls).
Study participants were genotyped on the Affymetrix or the Human1M/1Mduo platforms (22, 23) and imputed to the Haplotype Reference Consortium panel (1). Variants were excluded if the imputed minor allele frequency (MAF) was < 1%, missing rate > 10%, or imputation R2 < 0.3. In total, there were 7,672,848 imputed genetic variants, which were coded as imputed dosage ranging from 0 to 2.
GWAS summary statistics and linkage disequilibrium (LD) reference panel
The GWAS summary statistics of over 125,000 colorectal cancer cases and controls were obtained from Huyghe and colleagues (1), including data from the Genetic and Epidemiology Colorectal Cancer Consortium (GECCO), the Colon Cancer Family Registry (CCFR), and the Colorectal Transdisciplinary Study (CORECT). A full list of studies that contribute to the summary statistics is given in Supplementary Table S1 with more information published previously (1, 24–26). We removed variants with MAF < 0.01. The LD reference panel is based on 29,303 individuals of European ancestry from GECCO (1, 24).
Single-SNP MeQTL analysis
We assessed the pairwise association of SNPs and CpGs to identify meQTLs. We undertook a two-stage discovery and replication design by splitting the participants genotyped on Affy Chip in the discovery (n = 856) and those genotyped on Human1M/1Mduo (n = 499) in the replication. We adjusted for sex, case–control status, sample plate, position on plate, age, two principal components (PC) for technical effects, and 6 PCs from a subset of CpG sites to adjust for unobserved cell types. More details on the PC are provided in Supplementary Methods. DNAm was logit transformed so that the range would not be restricted (27). The analysis was done using the matrixEQTL R package (28).
We classified SNP–CpG pairs into cis-meQTL when the SNP was within 1 Mb of the CpG site, and trans-meQTL otherwise. A total of 1,710,673,549 cis-associations and 2,286,700,896,755 trans-associations were tested. We adjusted for multiple testing using the Bonferroni correction with type I error 0.05 for cis- and trans-meQTL associations, separately. For discovery, we set the genome-wide significance threshold 0.05/1,710,673,549 ≈ 2.92e–11 for cis-associations, and 0.05/2,286,700,896,755 ≈ 2.19e–14 for trans-associations. For replication, the threshold was based on the number of cis (ncis) or trans (ntrans) that passed the discovery stage, i.e., 0.05/(ncis or ntrans). The direction of the association in the replication stage also needed to be the same as in the discovery stage. Top meQTL associations were defined as having the lowest P value in the replication set.
We performed LD pruning on the meQTLs to account for the high SNP LD. For each CpG site, we selected the lead meQTL (based on the replication P value) and removed any meQTL with LD R2 > 0.2 with the lead meQTL. We then moved to the next most significant remaining meQTL and repeated until no meQTLs left. For computational ease, we assumed no cross-chromosomal LD.
Association of single-SNP meQTL with colorectal cancer
We examined associations of meQTLs identified from the single-SNP meQTL analysis with colorectal cancer risk based on GWAS summary statistics (1). For each CpG site, we picked its lead meQTL with the largest test statistic and examined their associations with colorectal cancer risk. We adjusted for multiple testing using Bonferroni correction based on the number of SNPs being tested. An association was considered novel if there were no known colorectal cancer loci within 1 MB. The known loci are those that reached the genome-wide significance level 5e–08 as previously reported (Supplementary Table S2).
Multi-SNP meQTL analysis
To assess the joint effects of SNPs on DNAm, we fit a penalized regression model on each CpG site. We did not include trans-SNPs as the majority of meQTLs are cis and substantial computational burden if we include trans-SNPs. To improve power, we combined both discovery and replication methylation data and adjusted for the genotyping platform in addition to the covariates in the single-SNP meQTL analysis. These covariates were not penalized. We fit four different elastic net models with α values of 0.25, 0.5, 0.75 and 1.0 (lasso), following previous literature that showed elastic net models tend to provide better prediction when some sparsity was induced (16, 18, 29, 30). We performed 10-fold cross-validation and determined the optimal penalty λ value based on cross-validation MSE. We estimated the weights at the optimal penalty value using the entire data set. Analyses were done using the glmnet package in R (31). We evaluated performance by the mean cross-validated partial R2 of DNAm variation explained by SNPs across 10 folds.
Association of multi-SNP meQTLs with colorectal cancer
We tested the joint association of meQTLs identified from the penalized regression with colorectal cancer risk for each CpG site using sMiST (21). sMiST combines both the association of the genetically predicted DNAm at that CpG site using the weights obtained from the penalized regression and the association of residual effects of individual meQTLs (random effects; ref. 21; Supplementary Methods). We used the Fisher combination procedure to combine the tests for these two components. We restricted our analysis to CpG sites with a mean cross-validated partial R2 above 1%. A P < 0.05/# of CpG sites tested was considered genome-wide statistically significant.
We performed conditional analysis of CpG sites that passed the genome-wide significance, adjusting for all known colorectal cancer GWAS loci. An association with P < 0.05/# of genome-wide significant CpG sites was considered novel if there were no known loci within 1 MB or as a novel secondary signal otherwise. As the participants with methylation data were part of GECCO, which contributed to the GWAS colorectal cancer summary statistics, we conducted a sensitivity analysis of significant associations by excluding the GECCO data.
Results
Table 1 provides descriptive statistics of the people in the methylation analysis data set. There were 1,355 participants, with 856 genotyped using the Affy Chip assay and 499 using the Human1M/1Mduo genotyping array. Roughly 50% were men, 47.5% were colorectal cancer cases, and the mean age was 62 years old.
. | Affy Chips . | Human1M/1Mduo . | Overall . |
---|---|---|---|
. | (n = 856) . | (n = 499) . | (n = 1,355) . |
Sex | |||
Female | 454 (53.0%) | 225 (45.1%) | 679 (50.1%) |
Male | 402 (47.0%) | 274 (54.9%) | 676 (49.9%) |
Case/control | |||
Case | 456 (53.3%) | 187 (37.5%) | 643 (47.5%) |
Control | 400 (46.7%) | 312 (62.5%) | 712 (52.5%) |
Age reference | |||
Mean (SD) | 62.3 (7.62) | 60.3 (9.22) | 61.6 (8.30) |
Median (min–max) | 63.0 (29.0–77.0) | 62.0 (27.0–76.0) | 63.0 (27.0–77.0) |
. | Affy Chips . | Human1M/1Mduo . | Overall . |
---|---|---|---|
. | (n = 856) . | (n = 499) . | (n = 1,355) . |
Sex | |||
Female | 454 (53.0%) | 225 (45.1%) | 679 (50.1%) |
Male | 402 (47.0%) | 274 (54.9%) | 676 (49.9%) |
Case/control | |||
Case | 456 (53.3%) | 187 (37.5%) | 643 (47.5%) |
Control | 400 (46.7%) | 312 (62.5%) | 712 (52.5%) |
Age reference | |||
Mean (SD) | 62.3 (7.62) | 60.3 (9.22) | 61.6 (8.30) |
Median (min–max) | 63.0 (29.0–77.0) | 62.0 (27.0–76.0) | 63.0 (27.0–77.0) |
Landscape of MeQTL
A total of 5,228,463 cis SNP–DNAm associations passed the genome-wide threshold in the discovery stage and 3,444,113 were replicated. Of these 3,444,113 cis-associations, there were 32,686 unique CpG sites and 1,280,770 unique SNPs, with many SNPs associated with multiple CpG sites. For trans-associations, a total of 559,695 SNP–DNAm associations passed the genome-wide threshold in the discovery stage and 292,919 trans-associations were replicated. Of these 292,919 associations, 160,752 were CpG–SNP pairs on the same chromosome but over 1 Mb apart. Across the 2,288,411,570,304 tests performed (every CpG–SNP pair) in the discovery, the median P value was between 0.499 and 0.501, suggesting no overall inflation.
Of the 142 colorectal cancer known variants (Supplementary Table S2), all but two variants (rs1391441 and rs11087784) were either meQTLs (n = 47) or within 1 Mb of a meQTL (n = 93). For the 47 known loci that are meQTLs, the majority (n = 43) were associated with CpGs on the same chromosome (Supplementary Tables S3 and S4). SNP rs3131043 (6:30758466_A/G) was a meQTL for 28 CpG sites including 18 cis, 8 trans on the same chromosome, and 2 trans on the different chromosomes. Among the 94 that were not meQTLs but within 1 Mb of meQTLs (one variant was not in our data set), only three were in high LD with a meQTL (R2 above 0.7): rs7121958, rs12594720, and rs6983267.
We examined the number of CpGs associated with each SNP stratified by cis-associations and trans-associations that on the same and different chromosomes (Fig. 1). Not surprisingly, there was a peak of cis-associations on chromosome 6 that map to the MHC region (Fig. 2). In total, there were 584,770 unique CpG–SNP associations within 5 Mb on either side of HLA-A (the peak), accounting for 15.6% of all CpG–SNP associations and 42.8% of all trans-associations (459,341 cis, 97,930 trans on the same chromosome, and 27,499 trans on the different chromosomes). This high amount of trans on the same chromosome is likely due to the extensive LD in the MHC region.
As many cis- and trans-associations are in LD, we performed LD pruning for each CpG site. This resulted in 70,106 associations with 67,620 cis-associations, 529 trans on the same chromosome, and 1,957 trans on different chromosomes. Due to the uniqueness of the HLA region, the associations where CpGs or SNPs were within 5 Mb of HLA-A on either side were excluded in tallying the number of cis- and trans-associations. The top cis-association was between rs7083920 (10:45071890_C/T) and cg02113055 (10: 45072520; replication P < 1e−100), which were 630 bp apart. The top trans on the same chromosome association was rs3764855 (17:66244578_T/C) and cg27123423 (≈ 3.27 Mb, replication P < 1e−100). The top trans-association on different chromosomes was rs12222783 (11:108531761_T/C) and cg23247274 (MRP63/SKA3, 13:21750649; replication P < 1e−100). These associations were in excellent agreement between the replication and discovery data sets (Fig. 3). Figure 3 also displays the next three strongest associations in all three categories (list of associations in Fig. 3 is given in Supplementary Table S5).
Novel association of single-SNP meQTLs with colorectal cancer
We assessed the association of single-SNP meQTLs with colorectal cancer risk. Across all CpG–SNP associations, there were 33,845 unique CpG sites. Taking the top meQTL for each of these CpG sites resulted in 27,287 unique SNPs associated with DNAm, of which 3,917 SNPs were the top meQTL for multiple CpGs. Using the colorectal cancer summary statistics, 70 SNPs passed the Bonferroni significance threshold 0.05/27,287 = 1.83e−06 (Supplementary Table S6). Of these 70 SNPs, 68 were within 1 Mb of 142 known GWAS variants. The two SNPs outside 1 Mb were in high LD (R2 = 0.99), and they were associated with increased colorectal cancer risk (SNP rs2749870, 14:52493068_reference allele/effect allele [G/A], OR = 1.04, 95% CI = 1.03–1.06, P = 1.2e−06; SNP rs2516600, 14:52487766_T/C, OR = 1.04, 95% CI = 1.03–1.06, P = 1.3e−06). Conditioning on known loci did not change the associations greatly (rs2749870 P = 1.37e−05, rs2516600 P = 1.37e−05).
SNP rs2749870 upregulated methylation at cg26923084 (14: 52536893, replication effect = 0.39, 95% CI = 0.35, 0.41), whereas SNP rs2516600 downregulated methylation at cg20550154 (14: 52487779, replication effect = −0.71, 95% CI = −0.76, −0.69) with DNAm at these two CpG sites negatively correlated (Pearson correlation = −0.5; Supplementary Fig. S1). We performed colocalization for the region from 52.4 to 52.5 Mb using the coloc package (32). There was strong evidence of a common signal for colorectal cancer risk and methylation at cg20550154 (posterior probability = 0.96), but not at cg26923084 (posterior probability = 0.24).
Multi-SNP meQTLs and predictive weights
We built genetic prediction models for DNAm at each CpG site using cis SNPs. Of the 298,248 CpG sites, 179,087 CpG sites had at least one SNP selected by an elastic net model (178,390, 175,709, 174,647, and 174,015 for α = 0.25, 0.50, 0.75, and 1, respectively). All models had similar partial R2 (Supplementary Fig. S2), but lasso had the smallest model size with mean number [standard deviation (SD), range] of selected SNPs 15.5 (19.2, 1–755). For the remaining analyses, we focused on lasso-predictive weights, as it provided the most parsimonious models. A total of 170,340 CpG sites had a mean cross-validated partial R2 for predictive DNAm above 1% and 45,965 above 10%. Beyond that the numbers decreased rapidly (Fig. 4), with only 4,422 having a partial R2 above 50%. The mean and median partial R2 were 9.8% and 5.4%, respectively. All genetic weights are available upon request from the authors.
Novel multi-SNP meQTLs with colorectal cancer
We ran sMiST on 168,963 CpG sites with partial R2 ≥ 0.01, after removing 1,377 CpG sites where all selected SNPs had MAF < 1% in the GWAS summary statistics. A total of 853 CpG sites passed the significance threshold (0.05/168,963 ≈ 2.96e−07; Supplementary Fig. S3). After adjusting for the known loci, 13 CpG sites were still significant. Of these 13, 3 CpG sites cg10226744 (22: 46068542, ATXN10), cg05743278 (7:129709803, KLHDC10), and cg19037357 (7: 130081024, CEP41) had no colorectal cancer known loci within 1 Mb (Table 2). In our sensitivity analysis, removing GECCO from the GWAS summary statistics did not change the result (Supplementary Table S7).
CpG . | Gene . | NSNPa . | R2b . | Unadj. Pc . | Cond. Pd . | Burden Pe . | Random Pf . | Nearest known CRC SNPg . |
---|---|---|---|---|---|---|---|---|
cg18541609 (5:430,359) | AHRR | 26 | 0.082 | 1.50E−07 | 2.50E−05 | 3.90E−03 | 4.40E−04 | rs78368589 |
cg23330450 (6:31,548,285) | 56 | 0.122 | 2.00E−17 | 3.30E−05 | 5.10E−04 | 4.70E−03 | rs2516420 | |
cg14426347 (6:31,632,827) | BAT4/CSNK2B | 20 | 0.042 | 7.60E−13 | 1.70E−05 | 3.60E−05 | 3.10E−02 | rs2516420 |
cg13653456 (6:31,763,819) | VARS | 20 | 0.074 | 3.80E−10 | 3.70E−06 | 4.70E−05 | 4.80E−03 | rs2516420 |
cg23945481 (6:31,938,984) | STK19/DOM3Z | 14 | 0.053 | 4.20E−11 | 5.50E−05 | 1.30E−02 | 3.00E−04 | rs3830041 |
cg22767321 (6:32,143,681) | AGPAT1 | 12 | 0.043 | 1.20E−08 | 2.40E−06 | 1.90E−06 | 7.60E−02 | rs3830041 |
cg07588430 (6:41,561,444) | FOXP4 | 4 | 0.014 | 5.80E−10 | 3.80E−08 | 2.20E−07 | 8.30E−03 | rs62396735 |
cg05743278 (7:129,709,803) | KLHDC10 | 5 | 0.031 | 1.40E−07 | 8.50E−07 | 3.30E−06 | 1.40E−02 | — |
cg19037357 (7:130,081,024) | CEP41 | 17 | 0.067 | 1.60E−07 | 2.30E−06 | 6.10E−01 | 2.30E−07 | — |
cg07005513 (11:61,595,956) | FADS2 | 10 | 0.076 | 3.70E−11 | 1.00E−06 | 5.50E−02 | 1.10E−06 | rs174533 |
cg18213653 (12:112,605,734) | C12orf51 | 7 | 0.016 | 2.10E−08 | 2.20E−09 | 2.40E−03 | 3.80E−08 | rs597808 |
cg19193595 (15:67,396,487) | SMAD3 | 5 | 0.022 | 4.40E−14 | 2.30E−05 | 6.50E−01 | 2.50E−06 | rs56324967 |
cg10226744 (22:46,068,542) | ATXN10 | 7 | 0.093 | 1.10E−07 | 2.80E−09 | 1.60E−04 | 7.30E−07 | — |
CpG . | Gene . | NSNPa . | R2b . | Unadj. Pc . | Cond. Pd . | Burden Pe . | Random Pf . | Nearest known CRC SNPg . |
---|---|---|---|---|---|---|---|---|
cg18541609 (5:430,359) | AHRR | 26 | 0.082 | 1.50E−07 | 2.50E−05 | 3.90E−03 | 4.40E−04 | rs78368589 |
cg23330450 (6:31,548,285) | 56 | 0.122 | 2.00E−17 | 3.30E−05 | 5.10E−04 | 4.70E−03 | rs2516420 | |
cg14426347 (6:31,632,827) | BAT4/CSNK2B | 20 | 0.042 | 7.60E−13 | 1.70E−05 | 3.60E−05 | 3.10E−02 | rs2516420 |
cg13653456 (6:31,763,819) | VARS | 20 | 0.074 | 3.80E−10 | 3.70E−06 | 4.70E−05 | 4.80E−03 | rs2516420 |
cg23945481 (6:31,938,984) | STK19/DOM3Z | 14 | 0.053 | 4.20E−11 | 5.50E−05 | 1.30E−02 | 3.00E−04 | rs3830041 |
cg22767321 (6:32,143,681) | AGPAT1 | 12 | 0.043 | 1.20E−08 | 2.40E−06 | 1.90E−06 | 7.60E−02 | rs3830041 |
cg07588430 (6:41,561,444) | FOXP4 | 4 | 0.014 | 5.80E−10 | 3.80E−08 | 2.20E−07 | 8.30E−03 | rs62396735 |
cg05743278 (7:129,709,803) | KLHDC10 | 5 | 0.031 | 1.40E−07 | 8.50E−07 | 3.30E−06 | 1.40E−02 | — |
cg19037357 (7:130,081,024) | CEP41 | 17 | 0.067 | 1.60E−07 | 2.30E−06 | 6.10E−01 | 2.30E−07 | — |
cg07005513 (11:61,595,956) | FADS2 | 10 | 0.076 | 3.70E−11 | 1.00E−06 | 5.50E−02 | 1.10E−06 | rs174533 |
cg18213653 (12:112,605,734) | C12orf51 | 7 | 0.016 | 2.10E−08 | 2.20E−09 | 2.40E−03 | 3.80E−08 | rs597808 |
cg19193595 (15:67,396,487) | SMAD3 | 5 | 0.022 | 4.40E−14 | 2.30E−05 | 6.50E−01 | 2.50E−06 | rs56324967 |
cg10226744 (22:46,068,542) | ATXN10 | 7 | 0.093 | 1.10E−07 | 2.80E−09 | 1.60E−04 | 7.30E−07 | — |
aNumber of SNPs selected by the lasso model.
bPartial R2 of DNAm explained by genetic variants adjusting for covariates.
csMiST P value not adjusting for known GWAS loci.
dsMiST P value adjusting for all known GWAS loci.
eP value for the conditional burden component (predicted DNAm) adjusting for all known GWAS loci.
fP value for the conditional random effects component adjusting for all known GWAS loci.
gNearest known GWAS SNP (based on Table S1).
Upon closer examination, for most CpG sites, both the genetically predicted CpG component and the random effects component contributed to the overall association, highlighting the power of sMiST for testing the total effects of meQTLs and the idea that the SNPs may contribute to colorectal cancer risk beyond the CpG site considered. The most significant conditional association was cg18213653 (C12orf51) on chromosome 12 (overall P = 2.23e−09) with burden (genetically predicted CpG) P = 2.42e−03 and random effects P = 3.82e−08. The association signal of the novel CpG site cg05743278 (KLHDC10) was driven mainly by the predicted methylation (P 3.3e−06). A total of five meQTLs in this CpG explained 3% of the variation. On the other hand, cg19037357, cg07005513, and cg19193595 were driven mainly by the random effects, suggesting potential individual variant effects, in addition to the effect of the predicted methylation effect. Interestingly, we did not observe much correlation among the individual DNAm of these 13 sites (Supplementary Fig. S4).
As a follow-up analysis, we examined sMiST results of the two CpG sites from the meQTL colorectal cancer. The meQTL rs2516600 was selected as one of the 11 SNPs in the lasso regression, which together strongly predicted methylation of cg20550154 (partial R2 = 0.72). The sMiST P value for this CpG site is 1.04e−05 with the P value for predicted methylation 6.9e−07 and the random effects P = 0.99, suggesting the colorectal cancer association was likely driven by methylation. The top meQTL of the other novel CpG site, cg26923084, was also selected as one of cg26923084’s 28 SNPs by lasso. The sMiST for this CpG site had an overall P value of 4.2e−07, with the P value for predicted methylation at 2.29e−06 and for random effects at 0.01, indicating the association was also likely driven by methylation.
Discussion
Herein, we provided a systematic examination of the relationship between genetic variants and DNAm and their subsequent association with colorectal cancer risk. We used two complementary approaches: (i) examining pairwise association between SNPs and DNAm and then testing whether the significant meQTLs were associated with colorectal cancer risk; and (ii) examining joint association of cis-variants with DNAm and using these weights to test the overall association of predicted methylation and residual variant effects with colorectal cancer risk via sMiST. The single-SNP meQTL approach revealed a novel colorectal cancer risk locus and that about one third of known loci were meQTLs. The multi-SNP approach revealed three novel loci, as well as several novel secondary signals of regions within 1 Mb of known loci.
In our single-SNP meQTL approach, the two novel SNPs, rs2516600 and rs2749870, are in high LD and likely represent one signal. However, these two SNPs are not the most associated SNPs in the cis-region of CpG sites and rs1051069 is (Supplementary Fig. S5; colorectal cancer P = 6.59e−08). The CpG sites for which both rs2516600 and rs2749870 are top meQTLs are mapped to Nidogen 2 (NID2). Interestingly, these two CpG sites are located within two different regions of the NID2 gene: cg20550154 is within the gene body, whereas cg26923084 is within 1,500 bp of the transcription start site. NID2 may be involved in maintaining the structure of the basement membrane (33). The encoded cell-adhesion protein NID2 binds to laminins (33). Members of the laminin family, including LAMA5 and LAMC1, have been implicated in colorectal cancer risk by GWAS (34). The NID2 promotor is aberrantly methylated in human gastrointestinal cancers, including colon cancer, and methylation has been shown to silence gene expression (35). NID2 gene methylation is also robustly associated with bladder cancer and urine-based methylation biomarker panels containing NID2 have been proposed for the early detection of primary bladder cancer (36, 37).
In our multi-SNP approach, we found an association with cg10226744, mapped to the body of the ATXN10 gene. The Ataxin 10 (ATXN10) gene is expressed in most tissues and identified as a downstream effector of the p53 and pRB tumor suppressor pathways involved in cellular senescence (38). Circumventing senescence has been proposed to represent an essential step in tumor progression (39). Four of the seven lasso-selected variants for cg10226744 are within ATXN10, two near FAM118A, and one near WNT7B. WNT7B is a member of the WNT pathway, which plays a major role in colorectal cancer signaling (40).
Both CpG sites cg05743278 and cg19037357 from the multi-SNP analysis are mapped to the 7q32.2 loci. This is a known imprinted region, with associations with basal-cell carcinoma and type 2 diabetes (41). cg05743278 is within 1,500 bp of the transcription start site of KLHDC10 and cg19037357 is within 1,500 bp of the transcription start site of CEP41. One SNP (rs10277874) was selected for both CpG sites by lasso. KLHDC10 has been found to be a suppressor of PP5 (42). PP5 is an alias for TFPI2, which is a methylation biomarker for colorectal cancer (43, 44). CEP41, though not yet implicated in colorectal cancer, has been found to express in gastrointestinal tissues (45). Besides KLHDC10 and CEP41, the gene KLF14, where nearby SNPs were selected for cg05743278 and cg19037357, may also be a candidate, as it has been found to be associated with colorectal cancer (46).
The 70,106 unique CpG–SNP associations after LD pruning were made up of 55,209 unique SNPs. Of these, 46,254 were associated with only one CpG. Of the remaining 8,955 that were associated with more than one CpG sites, only 8,576 were associated with cis CpG sites, making it difficult to ascertain whether the associations were hotspots or due to correlation among the nearby CpG sites. Among the remaining SNPs, 12 were meQTLs for at least five CpG sites on different chromosome. The SNP rs11673023 was also reported by Huan and colleagues as being in hotspot 22 (H22) of their trans-meQTL hotspots (47). In addition, three variants were within the region of the previously reported hotspot 14 (47, 48).
Although we have found multiple potentially novel loci, there are several limitations to our study. First, to increase the sample size for determining joint associations of SNPs with DNAm, we used the entire methylation data set to infer meQTLs. The weights from the penalized regression model are likely biased (49). Ideally, the weights would be reestimated using an independent validation data set. However, the power loss for identifying novel loci with colorectal cancer risk due to biased weights is likely minimal as a scale shift for genetically predicted methylation does not affect testing. Second, our participants are all of the European ancestry. It is well recognized that many of these SNPs tag causal variants, and due to LD and minor allele differences between populations, the effects of these SNPs may not be generalizable to other populations. It is important to validate these results in other racial and ethnic populations. Third, the DNAm was collected in whole blood as opposed to the relevant colon and rectal tissues. However, it has been observed that cross-tissue meQTL effects are prevalent (50, 51), which, to some extent, are similar to gene-expression quantitative trait loci (52). We also observed that gene-regulatory elements of immune cell types from blood are enriched for variants with colorectal cancer association (53). Interestingly, when we examined the gene expression of the three novel associations (NID2, KLHDC10, and CEP41), there was higher gene expression in the colon tissue compared with whole blood (Supplementary Figs. S6–S9). Regardless, future studies for identifying meQTLs in the relevant colon and rectal tissues and their associations with colorectal cancer risk are needed.
In our meQTL analyses, we adjusted for case–control status, which may introduce collider bias. We examined the association of novel colorectal cancer with methylation at cg20550154 and cg26923084; ignoring case–control status or restricting the analysis to just controls, the associations did not alter however (Supplementary Table S8). Further, there may be differential meQTL effects between cases and controls. Toward this end, we assessed the 3,737,032 unique associations from our discovery for case/control–SNP interaction in a subsample of cases–controls where cases and controls were on the same plate. Despite diagnosis/treatment potentially impacting DNAm, only 615 associations had a significant interaction effect (0.05/3,737,032), and none overlapped the 70 associations detected in the single-SNP meQTL colorectal cancer GWAS. The limited number of interaction effects identified is likely due to low power. We further performed separate cis-meQTL analysis for cases and controls separately. The P values of SNP–DNAm associations were highly correlated between cases and controls, suggesting meQTLs are generally shared between cases and controls (Supplementary Fig. S10). Thus, taking a combined analysis of cases and controls will have the most power in identifying these shared meQTLs. There were only 979 SNP–DNAm (after LD pruning) that reached our previous discovery cutoff (alpha = 2.9e−11) in cases but showed no association in controls (P < 2.9e−11 in cases but >0.5 in controls), and 1,393 SNP–DNAm that showed associations in controls but not in cases, 15 of 979 and 103 of 1,393 of which were associated in our combined analysis. None of the SNPs were known colorectal cancer loci. To maximize the utility of the data, we studied common effects across cases and controls; studies with larger sample size of cases and controls will be needed to study potentially differential effects.
Our study has many strengths. These include large sample sizes for GWAS colorectal cancer and methylation; use of novel powerful statistical methods to account for the total effect of both genetically predicted methylation and individual variant effects via alternative pathways; and the first large-scale examination of methylation association with colorectal cancer risk.
In summary, we have used a novel framework for utilizing DNAm data to discover novel loci as well as secondary signals associated with colorectal cancer risk. The findings could elucidate novel pathways in colorectal cancer tumorigenesis. These findings, in particular, the novel loci NID2, ATXN10, KLHDC10, and CEP41, warrant follow-up in future functional studies. Finally, we applied the summary sMiST analysis, a powerful statistical method to combine both information from the effect through methylation and residual direct effects of the meQTLs. We hope this paper will motivate researchers to utilize functional genomic data such as DNAm or gene expression and large-scale GWAS summary statistics with these statistical tools to discover new genetic associations.
Authors' Disclosures
M. Giannakis reports grants from Bristol-Myers Squibb, Merck, Servier, and Janssen outside the submitted work. H. Hampel reports other support from Invitae, Genome Medical, Promega, and GI OnDemand outside the submitted work. T.J. Hudson reports being currently employed by AbbVie. T.J. Hudson reports that contributions to this manuscript are related to former position at the Ontario Institute for Cancer Research; AbbVie has no connection to the manuscript. V. Moreno reports grants from the Agency for Management of University and Research Grants (AGAUR) of the Catalan Government, Instituto de Salud Carlos III, cofunded by FEDER Funds—a Way to Build Europe, grants from Spanish Association Against Cancer (AECC) Scientific Foundation, and Consortium for Biomedical Research in Epidemiology and Public Health (CIBERESP) during the conduct of the study. R.K. Pai reports personal fees from Alimentiv, Inc., Allergan, Eli Lilly, AbbVie, and PathAI outside the submitted work. P.D. Pharoah reports grants from Cancer Research UK during the conduct of the study. No disclosures were reported by the other authors.
Authors' Contributions
R. Barfield: Conceptualization, formal analysis, investigation, visualization, methodology, writing–original draft, writing–review and editing. J.R. Huyghe: Resources, data curation, investigation, writing–original draft. M. Lemire: Resources, data curation. X. Dong: Data curation, software, methodology. Y. Su: Data curation, software, methodology. S. Brezina: Resources, investigation. D.D. Buchanan: Resources, investigation. J.C. Figueiredo: Resources, investigation. S. Gallinger: Resources, data curation. M. Giannakis: Resources, investigation. A. Gsur: Resources. M.J. Gunter: Resources, investigation. H. Hampel: Resources, investigation. T.A. Harrison: Resources, investigation. J.L. Hopper: Resources, investigation. T.J. Hudson: Data curation. C.I. Li: Resources, investigation. V. Moreno: Resources, investigation. P.A. Newcomb: Resources, data curation, investigation. R.K. Pai: Resources, investigation. P.D. Pharoah: Resources, investigation. A.I. Phipps: Investigation. C. Qu: Resources, data curation. R.S. Steinfelder: Resources, data curation. W. Sun: Resources, investigation. A.K. Win: Investigation. S.H. Zaidi: Data curation, investigation. P.T. Campbell: Resources, investigation. U. Peters: Conceptualization, resources, data curation, funding acquisition, investigation, writing–original draft, writing–review and editing. L. Hsu: Conceptualization, supervision, funding acquisition, methodology, writing–original draft, writing–review and editing.
Acknowledgments
The authors thank the participants across all studies. In addition, we thank the reviewers for valuable feedback. This project was a part of Richard Barfield's postdoctoral work at Fred Hutch. L. Hsu was partially funded by the NIH/NCI grant R01 CA189532. U. Peters was supported by the Genetics and Epidemiology of Colorectal Cancer Consortium (GECCO): NCI, NIH, U.S. Department of Health and Human Services (U01 CA164930, U01 CA137088, R01 CA059045, and R01 CA201407). Analyses were run on the Fred Hutch cluster, which is part of the Scientific Computing Infrastructure funded by ORIP grant S10OD028685. A full list of study-specific acknowledgments is provided in the Supplementary Data.
The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.