Multiple gene-specific DNA methylation in blood leukocytes and colorectal cancer risk: a case-control study in China

The relationship between gene-specific DNA methylation in peripheral blood leukocytes and colorectal cancer (CRC) susceptibility is unclear. In this case-control study, the methylation status of a panel of 10 CRC-related genes in 428 CRC cases and 428 cancer-free controls were detected with methylation-sensitive high-resolution melting analysis. We calculated a weighted methylation risk score (MRS) that comprehensively combined the methylation status of the panel of 10 genes and found that the MRS_10 was significantly associated with CRC risk. Compared with MRS-Low group, MRS-High group and MRS-Medium group exhibited a 6.51-fold (95% CI, 3.77-11.27) and 3.85-fold (95% CI, 2.72-5.45) increased risk of CRC, respectively. Moreover, the CRC risk increased with increasing MRS_10 (Ptrend < 0.0001). Stratified analyses demonstrated that the significant association retained in both men and women, younger and older, and normal weight or underweight and overweight or obese subjects. The area under the receiver operating characteristic curves for the MRS_10 model was 69.04% (95% CI, 65.57-72.66%) and the combined EF and MRS_10 model yielded an AUC of 79.12% (95% CI, 76.22-82.15%). Together, the panel of 10 gene-specific DNA methylation in leukocytes was strongly associated with the risk of CRC and might be a useful marker of susceptibility for CRC.


Sample size
We estimated the sample size according to a logistic regression test assessing whether aberrant DNA methylation in blood leukocytes was associated with the risk of colorectal cancer (CRC). A sample size of 652 participants was needed to achieve 90% power (at the 5% level of statistical significance) to detect odds ratios (ORs) of 1.8 or more with a 20% prevalence in the control group. In addition, taking into consideration incomplete questionnaires and the failure rate for MS-HRM detection, we included about 20% more samples and finally targeted a total sample size of 800 participants.

DNA extraction
After the peripheral blood samples were centrifuged at 1,600 g for 10 minutes to separate plasma, leukocyte-derived DNA was extracted from buffy coats using the QIAamp DNA Blood Mini Kit (Qiagen, Hilden, Germany) according to the manufacturer's protocol. Briefly, buffy coats were resuspended in 400 μl of phosphate-buffered saline. After the addition of 40 μl of proteinase K and 400 μl of buffer AL (Qiagen), the mixture was incubated at 56°C for 10 min, treated with 400 μl of ethanol (99%), and mixed. The mixture was transferred to a QIAamp spin column and centrifuged at 6,000 g at room temperature for 1 minute, followed by washing with buffers AW1 and AW2 and centrifugation. DNA was eluted in 150 μl of buffer AE and was quantified using a NanoDrop 2000c bioanalyzer (Thermo-Fisher, USA). Then, the DNA was divided into aliquots and stored at -80°C.

Bisulfite modification
From each subject, 500 ng of DNA was bisulfitemodified using the EpiTect Plus DNA Bisulfite Kit (Qiagen) according to the manufacturer's protocol. The bisulfite-treated DNA was eluted in 30 μl of elution buffer and quantified using a NanoDrop 2000c bioanalyzer set to RNA as the sample type [1]. Then, the bisulfite-modified DNA sample was diluted to a final concentration of 10 ng/μl and was distributed into three aliquots for storage (-20°C).

Pyrosequencing verification
To verify the results of MS-HRM, we performed pyrosequencing in a subset of the samples. For DAPK1, pyrosequencing was performed for quantitative methylation analyses of 48 HRM-positive samples (36 CRC case and 12 control samples that were selected from 106 cases and 43 controls, respectively) and 48 HRMnegative samples (24 CRC case and 24 control samples that were selected from 322 cases and 385 controls, respectively); for MLH1, 24 HRM-positive samples (16 CRC case and 8 control samples) and 48 HRM-negative samples (24 CRC case and 24 control samples that were selected from 412 cases and 420 controls, respectively).
The regions analyzed via MS-HRM included 7 and 21 CpG dinucleotide sites for DAPK1 and MLH1, respectively. Subsequently, we designed pyrosequencing primer sets to analyze the corresponding regions for DAPK1 and MLH1 using PyroMark Assay Design Software (version 2.0.1.15, Qiagen). PyroPCR was performed using a PyroMark PCR Kit (Qiagen) according to the manufacturer's protocol. The primer sets used are listed in Supplementary Table 11 (the reverse primer in each set was biotinylated). The standard conditions for PyroPCR were as follows: 12.5 μl of PyroMark PCR Master Mix, 2.5 μl of CoralLoad Concentrate, 0.5 μl of each primer (10 pmol/μl) and 2 μl (nearly 20 ng) of bisulfite-modified DNA in a final volume of 25 μl. The PyroPCR thermocycling protocol was as follows: initial PCR activation (95°C for 15 min), 45 3-step cycles (94°C for 30 s, 56°C for 30 s, and 72°C for 30 s) and final extension (72°C for 10 min). For each sample, 10 μl of PyroPCR product was used for pyrosequencing reactions and methylation quantification, which were performed using a PyroMark Q24 Advanced System version 2.0.6 (Qiagen) according to the manufacturer's protocol. The pyrosequencing results were analysed using PyroMark software version 3.0.0 (Qiagen) to obtain a percentage of methylation at each CpG site. In the MS-HRM analysis, we could detect only the mean methylation percentage of the analyzed sequences and could not precisely quantify the methylation percentage of individual CpG sites. Therefore, the mean methylation percentage of all CpG sites based on pyrosequencing was calculated for each analyzed gene, and this value was used for comparison with the results of MS-HRM. The concordance between the methylation levels derived from the two DNA methylation assessment methods was assessed using Spearman correlation coefficients, ROC curves and AUC analyses. Additionally, we used the Bland-Altman plot as a graphical method to compare the two methods.

Missing data analysis and imputation
All questionnaire-derived variables were analyzed via missing value analysis (Little's missing completely at random (MCAR) test) to assess whether the pattern of missing data was related to the observed data (MCAR or missing at random (MAR)). The MAR assumption fit our data. Therefore, missing data were imputed using multiple imputation (five imputations) via the expectationmaximization method, with all questionnaire-derived variables (including outcomes) included in the imputation model. We used the method of multiple imputation to impute the missing dataset rather than single imputation because the former is generally considered to be superior to the latter for solving the problem of missing data [2, 3].

Statistical analysis
The sample size was calculated using PASS version 11.0.7 (NCSS LLC., USA). The Bland-Altman plot analysis was performed using MedCalc software version 15.4 (Ostend, Belgium).

Pyrosequencing verification
For DAPK1 and MLH1, the methylation status obtained via MS-HRM was compared with the mean methylation level based on quantitative pyrosequencing. Using Bland-Altman plots, we found that for MLH1, all examinations of methylation fell within the limits of agreement except for one sample displaying a high methylation level, in which MS-HRM reported a lower methylation percentage than pyrosequencing (Supplementary Figure 1). For DAPK1, all samples fell within the limits of agreement except for three samples with relatively high methylation, in which MS-HRM reported lower levels of methylation than pyrosequencing. In general, only among samples with a higher methylation percentage, MS-HRM tended to indicate lower levels of methylation than pyrosequencing. These resluts indicated that the MS-HRM results were well confirmed by the pyrosequencing results.