Methylation of leukocyte DNA and ovarian cancer: relationships with disease status and outcome

Study participants

DNA methylation assays and arrays

Peripheral blood (leukocytes) was used as the source of DNA. DNA was extracted from four milliliters of fresh peripheral blood using the Autogen Flexstar instrument utilizing Flexigene chemistry (salting out methodology). Blood is aliquoted for DNA extraction using an automated liquid handler with barcoding to ensure proper sample placement. Post-DNA extraction, DNA is aliquoted into a permanent storage tube utilizing an automated liquid handler with barcoding, again, to ensure proper sample placement. DNA samples are assessed for quality and concentration using a Trinean DropSense 96 spectrophotometer and DNA is then stored longterm at -80°C. The leukocyte-derived DNA (1 ug) was bisulfite modified (BSM) using the Zymo EZ96 DNA Methylation Kit (Zymo Research, Orange, CA) according to the manufacturer’s protocol. BSM DNA (250 ng) was then assayed on 96 well plates in three batches at the Mayo Clinic Molecular Genome Facility (Rochester, MN): Batch 1 used the Infinium HumanMethylation27 BeadChip on 84 cases and 91 controls, Batch 2 used this array on 172 cases and 176 controls and Batch 3 used Illumina Infinium HumanMethylation450 BeadChip on 156 cases and 157 controls. Methylation status at the target CpG sites was determined by comparing the ratio of fluorescent signal from the methylated allele to the sum from the fluorescent signal from both methylated and unmethylated alleles (i.e., the beta value).

To assess the quality of the DNAm data produced from the Illumina arrays, Centre d'Etudes du Polymorphisme Humain (CEPH) DNA, positive BSM controls (placental DNA) and negative BSM controls (whole genome amplified [WGA] DNA) were assayed within each batch. For the HumanMethylation27 BeadChips (Batch 1 and Batch 2), 9 CEPH DNA, 12 positive control DNA samples and 8 negative control DNA samples were also assayed, in addition to 12 replicate/duplicate samples. Similarly, for the HumanMethylation450 BeadChip batch (Batch 3), 6 CEPH samples, 11 positive control samples, 6 negative control samples and 6 replicate samples were assayed. Lastly, twenty duplicate samples were assayed using both Illumina Infinium HumanMethylation27 and HumanMethylation450 BeadChip in order to compare the methylation levels between the two arrays.

Quality control and normalization

Using Illumina GenomeStudio software, DNAm values from the HumanMethylation27 BeadChip assays were scored as beta values, ranging from 0 (unmethylated) to 1 (methylated), resulting in methylation beta values for 27,578 probes. Quality control was done for Batch 1 and Batch 2 combined and then separately for Batch 3. Probes were then excluded if they were on the Y chromosome, positioned at a single nucleotide polymorphism (dbSNP build 137), had high beta values in BSM negative controls (beyond four standard deviations of mean), or were detected in less than 70% of samples. Quality control was also conducted at the sample level, based on the bisulfite conversion ratio and call rate rates (based on a detection p-value of 0.05). Histograms and scatterplots of these statistics were used to determine which samples to exclude (i.e., “outliers”). Similar quality control steps were completed for the samples assayed using the HumanMethylation450 BeadChips, which contained 485,577 CpG site-specific probes.

For the HumanMethylation27 BeadChip arrays, 25,922 (94%) methylation probes passed quality control; for the HumanMethylation450, 441,716 (91%) methylation probes passed quality control. The pairwise correlations for beta values among CEPH replicates were excellent (≥0.97 for Batches 1 and 2, and >0.99 for Batch 3), as were the intra-class correlations of beta values among CEPH replicates (>0.98 for Batches 1 and 2, and >0.99 for Batch 3) and among duplicated study participant samples (>0.93 for Batches 1 and 2, and >0.81 for Batch 3). For 20 samples assessed across batches, the intra-class correlation for beta values of the 24,520 overlapping probes in the HumanMethylation27 and HumanMethylation450 BeadChips was > 0.88. Of samples in Batches 1 and 2, 6 were excluded based on call rates, and one failed bisulfite conversion; in Batch 3, 10 samples were removed following quality control (9 samples failed the bisulfite conversion, one sample with low mean methylation beta value across probes). Following exclusions, we included 69 cases and 87 controls in Batch 1, 146 cases and 176 controls in Batch 2, and 121 cases and 135 controls in Batch 3.

We assessed possible differences by plate and chips within plates (8 BeadChips per plate were assessed with 12 DNAs each) through principal component analyses. Based on the assessment of technical artifacts using principal component analyses, a plate effect was observed within each of the three batches and a chip within batch effect for the HumanMethylation27 data (Additional file 1: Figure S1). Therefore, we adjusted for a plate effect for batch 3 and for chip within plate effect for batches 1 and 2 using a linear model of the logit-transformed beta value for each CpG site, with the unstandardized residuals saved. The logit-transformed locus mean was added back onto the residuals followed by the transformation of the residual to the 0 to 1 scale, producing an “adjusted beta” value for all CpG sites.

Finally, we restricted analyses to probes in common between the DNAm arrays following quality control, excluding 9,341 CpG probes on the Illumina Infinium HumanMethylation27 shown to associate with cell type distribution at q-value < 0.05 [15, 16], as well as 1,363 CpG probes found by Chen et al. to be non-specific (i.e., mapped to multiple places along the genome) [17]. Thus, analyses focused on the remaining 13,816 CpG probes (i.e., 24,520 probes in common between the two panels following quality control minus 9,341 probes associated with cell type distribution minus 1,363 non-specific probes).

Statistical association analysis

We analyzed each batch separately using Van der Waerden rank, or rank-based inverse Gaussian, transformed beta values and combined results across batch using meta-analysis techniques. This allowed us to examine similarity of effects across batches and to estimate the combined effect. Meta-analysis was completed using a random effect meta-analysis. A Woolf’s test of homogeneity of regression coefficients across batches was performed, i.e. the distribution of regression estimates across batches for each probe is compatible with that expected given a common regression estimate. All statistical tests were 2-sided, and analyses of individual batches were carried out using SAS (version 9.3; SAS Institute Inc., Cary, NC) and R (version 2.14.0). Meta-analyses were carried out using the R package rmeta (http://CRAN.R-project.org/package=rmeta). To control for multiple testing, associations with p < 5 × 10^-6 were considered statistically significant (e.g., Bonferroni adjustment based on number of independent tests). Pathway analysis used Ingenuity Pathway Analysis (IPA) (Ingenuity® Systems, http://www.ingenuity.com) for genes closest to CpG probes associated with disease status or outcome at p < 0.0001.

The following linear model was used to determine if DNAm levels differ between EOC cases and matched controls for each CpG site. Let, ${{\mathit{Y}}_{\mathit{i}}}_{\mathit{j}}={\mathit{\alpha }}_{\mathit{j}}+{\mathit{\beta }}_{\mathit{j}}{\mathit{X}}_{\mathit{i}}+{\mathit{\gamma }}_{\mathit{j}}^{\mathit{T}}{\mathit{Z}}_{\mathit{i}}+{{\mathit{e}}_{\mathit{i}}}_{\mathit{j}},$ where Y _ij represents the adjusted methylation beta value for subject i and CpG probe j (j = 1…, 13816), X _i represents disease status for subject i (1 if case and 0 if control), Z _i represents covariates for subject i and ${{\mathit{\epsilon }}_{\mathit{i}}}_{\mathit{j}}~\mathit{N}\left(0,{\mathit{\sigma }}_{\mathit{j}}^2\right).$ To identify covariates that differ between EOC cases and controls to include in the model (i.e., potential confounders), potential covariates were examined for association with disease status within a stepwise logistic regression model, resulting in the inclusion of parity/age at first live birth combination (nulliparous, 1-2 and age < = 20 years, 1-2 and age > 20 years, 3+ and age < = 20 years, 3+ and age > =20 years, missing), current alcohol use (never, former, current, missing), current smoking status (never or former, current, missing), enrollment year, and recruitment state (MN vs. non-MN). For each CpG probe j, the disease status parameter ($\widehat{{\mathit{\beta }}_{\mathit{j}}}$) was estimated using the rank-transformed adjusted beta methylation values, along with a 95% confidence interval (CI).

We assessed associations of methylation beta values with overall survival (OS) using Cox proportional hazards regression analyses, adjusted for age at diagnosis, tumor stage (III/IV, I/ II), presence of ascites (yes, no, missing) and volume of residual tumor following debulking surgery (<1 cm, >1 cm, missing) based on stepwise Cox regression analysis. The proportionality assumption was assessed by the analysis of scaled Schoenfeld residuals for all covariates included in the statistical analysis and found to be upheld [18]. We accounted for left truncation using start-stop counting process style of input and estimated hazards ratios (HR) and 95% CIs [19].

Disease status and DNA methylation

In a meta-analysis across the three batches (two sets of experiments involving the Illumina Infinium HumanMethylation27 beadchip and one experiment involving the Illumina Infinium HumanMethylation450 beadchip) evaluating association between each of the 13,816 CpG probes and ovarian cancer case-control status (336 cases, 398 controls), 30 CpGs showed p-value ≤ 5×10^-7 (Table 2), where none of the tests for heterogeneity of effects across batches were significant (p > 0.05). We confirmed that these 30 CpGs were also included in the Koestler et al. (2012) analysis, and thus determined not to be associated with cell type distribution. Of these CpGs, the following were also replicated in an independent study (p < 0.001) conducted by Teshendorff et al. [14]: cg04834572 near DUSP13, cg10414058 near HDAC3, cg19280776 near PAG1, and cg24959428 near GBP6. In addition to the replication of specific CpG sites, C19orf18 and MARCH1 contained CpG sites found to be replicated for association with EOC risk [20]. All CpG sites, with the exception of a CpG near PAG1, had negative parameter estimates indicating lower methylation in the cases as compared to controls (e.g., cases were hypo-methylated). Plots of the entire set of results for the 13,816 CpG sites (i.e., sites contained in both the 27K and 450K arrays, specific and not associated with cell type distribution) are presented in Figure 1A. The top association between methylation and disease status, which as also replicated, was observed for the CpG probe cg04834572 located approximately 315 kb upstream of DUSP13 on chromosome 10 (Figure 2A) with a meta-analysis p-value of 1.6 × 10^-14 and individual batch p-values ranging from 2.1×10^-4 to 1.1 × 10^-6. DUSP13 is a member of the protein-tyrosine phosphatase superfamily and interacts with protein kinases involved in the regulation of cell proliferation and differentiation. Other significantly associated CpG sites were near biologically interesting/relevant genes, such as SRC (cg05498681; p = 4.8×10^-7) (Figure 2B), HHIP (cg14580567; p =5.6×10^-11) (Figure 2C), and replicated CpG near HDAC3 (cg10414058; p = 6.3×10^-12) (Figure 2D).

Probe ID	Ch	Position (bp)	Nearest genes	Location of nearest gene (bp)#	Location to Island^	Meta-Analysis§		Batch 1		Batch 2		Batch 3
						$\hat{\mathit{\beta }}$	P	$\hat{\mathit{\beta }}$	P	$\hat{\mathit{\beta }}$	P	$\hat{\mathit{\beta }}$	P
cg04834572*	10	76868766	DUSP13	76854190-76868970	Shelf	-0.65	1.6E-14	-0.82	2.1E-4	-0.55	6.0E-6	-0.73	1.1E-6
cg11722531	19	1449857	APC2	1450148-1473243	Shore	-0.62	2.7E-13	-0.64	3.5E-3	-0.54	6.4E-6	-0.74	1.3E-6
cg10414058*	5	141017903	RELL2	141016517-	Shore	-0.94	6.3E-12	-1.13	3.5E-8	-0.71	7.1E-10	-1.05	9.0E-14
			HDAC3	141020631
				141000443-141016423
cg14580567	4	145567271	HHIP	145567148-145659881	Island	-0.56	5.6E-11	-0.49	0.024	-0.55	3.8E-6	-0.61	9.3E-5
cg08245789	22	40289538	ENTHD1	40139049-40289794		-0.56	3.6E-10	-0.66	4.0E-3	-0.43	3.9E-4	-0.69	6.5E-6
cg27623214	19	58485726	C19orf18 †	58469805-58485902		-0.52	1.3E-9	-0.66	2.9E-3	-0.52	1.4E-5	-0.45	3.4E-3
cg23877385	15	59908652	GCNT3	59903982-59912210		-0.52	2.2E-9	-0.46	0.043	-0.53	1.3E-5	-0.52	6.3E-4
cg26150490	X	47863595	SPACA5	47863734-47869130	Island	-0.50	3.6E-9	-0.34	0.10	-0.54	6.8E-6	-0.52	6.6E-4
			ZNF182	47834250-47863394
cg22336401	9	140336227	ENTPD8	140328816-140335901	Island	-0.50	5.3E-9	-0.27	0.24	-0.51	2.5E-5	-0.60	9.1E-5
cg20775254	2	95940705	PROM2	95940201-95957056		-0.53	6.1E-9	-0.58	5.3E-3	-0.41	8.2E-4	-0.69	5.6E-6
cg19280776*	8	82024586	PAG1	81880045-82024303	Shore	0.49	1.0E-8	0.60	6.2E-3	0.52	1.4E-5	0.38	0.013
cg07634191	8	27850178	SCARA5	27727399-27850369		-0.48	1.8E-8	-0.50	0.025	-0.39	1.2E-3	-0.61	4.5E-5
cg21244955	22	21192955	PI4KA	21061979-21213100		-0.49	2.2E-8	-0.47	0.039	-0.48	6.6E-5	-0.50	1.3E-3
cg18159180	6	43022213	CUL7	43005355-43021683	Shore	-0.48	2.2E-8	-0.35	0.12	-0.57	2.8E-6	-0.41	8.0E-3
			MRPL2	43021767-43027242
cg04439215	13	45768901	KCTD4	45766988-45775175		-0.49	2.4E-8	-0.65	4.4E-3	-0.37	2.2E-3	-0.60	1.0E-4
			GTF2F2	45694631-45858240
cg26787239	5	132008525	IL4	132009678-132018370		-0.48	2.9E-8	-0.50	0.025	-0.50	3.4E-5	-0.43	4.8E-3
cg07259382	4	164536228	MARCH1 †	164445450-165304407		-0.48	3.0E-8	-0.58	8.5E-3	-0.43	3.1E-4	-0.50	1.2E-3
cg14808739	17	17741098	SREBF1	17714663-17740325	Shore	-0.48	3.7E-8	-0.41	0.063	-0.52	2.3E-5	-0.44	3.8E-3
cg00065385	9	111623395	ACTL7A	111624603-111626035		-0.47	6.1E-8	-0.50	0.026	-0.49	5.0E-5	-0.43	6.9E-3
cg04721883	X	103499577	ESX1	103494719-103499599	Island	-0.46	7.6E-8	-0.40	0.068	-0.38	1.7E-3	-0.62	4.6E-5
cg21400640	X	12992967	TMSB4X	12993226-12995346	Shore	-0.46	1.1E-7	-0.50	0.029	-0.46	1.4E-4	-0.44	4.0E-3
cg11871280	12	60082038	SLC16A7	59989821-60183636		-0.45	1.7E-7	-0.24	0.28	-0.51	2.4E-5	-0.46	3.1E-3
cg09261015	X	103499647	ESX1	103494719-103499599	Island	-0.45	2.1E-7	-0.46	0.046	-0.39	1.1E-3	-0.54	4.7E-4
cg18731813	X	100805683	ARMCX1	100805514-100809683	Shelf	-0.44	2.4E-7	-0.23	0.30	-0.46	1.3E-4	-0.51	7.3E-4
cg23279136	X	74375966	ABCB7	74273105-74376132		-0.44	3.2E-7	-0.62	6.3E-3	-0.33	5.7E-3	-0.54	4.9E-4
cg02254461	3	39195904	CSRNP1	39183342-39195102	Shore	-0.44	3.6E-7	-0.37	0.10	-0.48	9.1E-5	-0.42	6.1E-3
cg26246138	X	18372612	SCML2	18257433-18372844	Island	-0.44	3.8E-7	-0.37	0.096	-0.37	2.3E-3	-0.57	1.6E-4
cg24959428*	1	89829951	GBP6	89829436-89853719		-0.44	4.1E-7	-0.50	0.025	-0.42	5.2E-4	-0.44	4.4E-3
cg05498681	20	35973318	SRC	35973088-36033821	Shore	-0.44	4.8E-7	-0.17	0.45	-0.53	1.3E-5	-0.41	8.1E-3
cg01377911	19	49568036	NTF4	49564397-49567124		-0.43	5.0E-7	-0.38	0.086	-0.37	2.3E-3	-0.55	2.5E-4

^# Locations based on NCBI (http://www.ncbi.nlm.nih.gov), build 37.

^Shore CpG sites defined to be within +/- 2 kb from CpG island; Shelf CpG sites defined to be within +/- 2kb of CpG Shore.

^§ All tests for heterogeneity of effects across the three batches were non-significant (p > 0.05); analysis adjusted for age at first live birth, alcohol use, smoking status, enrollment year, and recruitment state; a negative parameter estimate indicates lower methylation in the cases as compared to the control (e.g., cases where hypomethylated).

*Same CpG site found to be associated with EOC risk (p ≤ 0.001) in a prior report [14].

^†CpG sites near this gene found to be associated with EOC risk in a prior report [14].

To identify any commonality of highlighted genes within biological pathways, pathway analysis using Ingenuity Pathway Analysis (IPA) was completed for the 155 genes closest to the CpG probes (based on Illumina provided annotation) that were associated with disease status based on a liberal threshold of p < 0.0001. The top pathways enriched for these 155 genes were the telomerase signaling (five genes in our top 155 were in the list of 99 genes within the telomerase signaling pathway; p = 1.24×10^-3 for enrichment of pathway) and the paxillin signaling (five genes in our top 155 were in the list of 110 genes within the paxillin signaling pathway; p = 1.42×10^-3). The five genes in the telomerase signaling pathway with methylation associated with disease status at p < 0.0001 were HDAC3 (p = 6.33×10^-12), IL2RG (p = 4.33×10^-6), PIK3C2B (p = 1.97×10^-5), PIK3R1 (p = 5.19×10^-5), and POT1 (p = 1.38×10^-6). PIK3C2B has been implicated in development of glioblastoma multiforme, while mutations in PIK3R1 have been seen in ovarian tumors and cancer cell lines and endometrial cancer [21–23]. POT1 has been found to be associated with tumor stage and telomere length in gastric cancer [24–26]. For the paxillin signaling pathway, the five differentially methylated CpGs were near ARFIP2 (p = 4.60×10^-5), ITGB6 (p = 3.95×10^-5), PIK3C2B, PIK3R1 and SRC, with some overlap between the top two pathways (PIK3R1 and PIK3C2B).

Survival following EOC and DNA methylation

Many fewer CpGs were associated with OS among the 366 cases than with case-control status, as illustrated in Figure 1B. None of the associations were statistically significant at the 5×10^-6 level; the top eight CpG probes with meta-analysis p-value < 10^-3 for association with OS are presented in Table 2. The top CpG sites associated with OS were cg10276549 within the promoter region of GABRE (p = 5.8×10^-5) (Figure 3A) and CpG site (cg06171242) within the promoter region of TTRAP/TDP2 (p = 4.4×10^-4). GABRE is a target for many benzodiazepine drugs used in the treatment of pain, insomnia, epilepsy, anxiety and panic related disorders [27–29]. However, little information can be found implicating a role of GABRE in response to chemotherapies (http://www.cancer.gov/clinicaltrials/). In addition to the modest level of association for CpGs near GABRE, there was a trend for association of CpG sites near the following biologically relevant genes: MT1X (p =7.4×10^-4) (Figure 3B), ADORA2B (p = 7.4×10^-4) (Figure 3C), and ABLM3 (p = 9.3×10^-4). These three CpG sites moderately associated with OS were all within CpG islands or shores and within the promoter region of the corresponding gene.

Similar to the analysis of the disease-associated genes, an exploratory pathway analysis using IPA was completed for the 61 genes closest to the CpG probes most associated with OS (meta-analysis p < 0.01). The top canonical pathways enriched for these 61 genes were relaxin signaling (five genes out of 147; GNA12, GNB1, PIK3R4, RAP1A, TDP2; p = 7.09×10^-5 for enrichment of pathway) and CXCR4 signaling (five genes out of 160; GNA12, GNB1, ITPR1, PIK3R4, ROCK1; p = 1.25×10^-4) and IL-8 signaling (five genes out of 192; ARRB2, GNA12, GNB1, PIK3R4, ROCK1; p = 3.05×10^-4). Three genes (GNB1 (p = 0.006), GNA12 (p = 0.009), and PIK3R4 (p = 0.002)) are part of all three of these canonical pathways.