DNA methylation in multiple myeloma is weakly associated with gene transcription

Abstract

Previous studies have now demonstrated that both genic and global hypomethylation characterizes the multiple myeloma (MM) epigenome. Whether these methylation changes are associated with global and corresponding increases (or decreases) in transcriptional activity are poorly understood. The purpose of our current study was to correlate DNA methylation levels in MM to gene expression. We analyzed matching datasets generated by the GoldenGate methylation BeadArray and Affymetrix gene expression platforms in 193 MM samples. We subsequently utilized two independent statistical approaches to identify methylation-expression correlations. In the first approach, we used a linear correlation parameter by computing a Pearson correlation coefficient. In the second approach, we discretized samples into low and high methylation groups and then compared the gene expression differences between the groups. Only methylation of 2.1% and 25.3% of CpG sites on the methylation array correlated to gene expression by Pearson correlation or the discretization method, respectively. Among the genes with methylation-expression correlations were IGF1R, DLC1, p16, and IL17RB. In conclusion, DNA methylation may directly regulate relatively few genes and suggests that additional studies are needed to determine the effects of genome-wide methylation changes in MM.

Figure 1. Methods used to identify genes with methylation-expression correlations.

(A) Pearson correlation was used to measure linear relationships between DNA methylation and gene expression levels for 1505 CpG probes represented on the GoldenGate Methylation BeadArray. The panels represent examples of a gene with high (left) and low (right) Pearson correlation coefficients when analyzing DNA methylation levels (x axis) against gene expression levels (y axis). (B) A discretization approach was used to classify samples into methylated (M) or unmethylated (U) groups based on the mean (μ) methylation value and standard deviation (σ) of a given probe. Statistically significant gene expression differences between M and U groups indicated a methylation-expression correlation for the gene in question.

Figure 2. Confirming expression trends by qRT-PCR for correlated genes identified by discretization approach.

Box plots represent gene expression levels generated by either microarray or qRT-PCR. Data are shown for samples classified as U or M based on the methylation status of p16 (A), DLC1 (B), IGF1R (C), or IL17RB (D). For microarray data, probe intensities are plotted on the y-axis. Relative fold-change differences are plotted for data generated by qRT-PCR. The number of samples in each group is displayed above each plot. The GoldenGate BeadArray probe names are indicated above each pair of box plots.

Figure 3. A diagrammatic representation of sample class enrichments for selected genes.

Samples classified as M or U by the discretization approach are indicated by green and blue bars respectively. The genomic context of samples in each group are shown as white, black or red bars. A representative locus is shown for p16 (A), DLC1 (B), IGF1R (C) and IL17RB (D).

Figure 4. Comparison of Pearson and discretization approaches used to identify methylation-expression correlations.

The Venn Diagram displays the total number of overlapping loci with positive/negative methylation-expression correlations identified by computing either a Pearson correlation or applying the discretization method.

Figure 5. Validation of methylation-expression relationships by Pearson correlation for p16, DLC1, IGF1R and IL17RB using pyrosequencing and qRT-PCR in an independent sample set.

Scatter plots depict the extent of linear correlation of DNA methylation (x axis) to gene expression (y axis) for p16 (A), DLC1 (B), IGF1R (C), IL17RB (D) in 43 MM samples. Methylation data was generated by bisulfite pyrosequencing and the average percent methylation for all CpG loci interrogated is shown. Gene expression relative fold-change was obtained by qRT-PCR. Pearson correlation coefficients are shown and P values are generated by random permutation tests.

Figure 6. Pearson correlation coefficients calculated for methylation of individual CpG loci.

Graphs show the Pearson correlation coefficient (y-axis) of individual CpG loci analyzed by pyrosequencing to the expression level for the corresponding gene. The distance to TSS is shown (x-axis) to demonstrate the proximity of each CpG locus to each other. Each CpG locus analyzed for p16 (A), DLC1 (B), IGF1R (C) and IL17RB (D) is depicted as a point on the graph.

Figure 7. Validation of methylation-expression relationships using a discretization approach for p16, DLC1, IGF1R and IL17RB by pyrosequencing and qRT-PCR in an independent sample set.

Bisulfite pyrosequencing data for 43 independent samples was used to discretize samples into methylated (M) or unmethylated (U) groups based on the percent methylation values obtained for p16 (A), DLC1 (B), IGF1R (C) and IL17RB (D). The average of all CpG loci interrogated by pryosequencing for each gene was used in the analysis. Differential gene expression analysis was conducted by qRT-PCR to compare the expression of each gene between U and M groups. Positive methylation-expression correlations were confirmed for p16 and DLC1. A negative correlation was validated for IGF1R and IL17RB. The number of patients in U and M groups is given above each box plot. Y-axis represents expression levels by plotting relative fold-change (2^−ΔΔCT).