γ-Glutamyl hydrolase modulation significantly influences global and gene-specific DNA methylation and gene expression in human colon and breast cancer cells

Abstract

γ-Glutamyl hydrolase (GGH) plays an important role in folate homeostasis by catalyzing hydrolysis of polyglutamylated folate into monoglutamates. Polyglutamylated folates are better substrates for several enzymes involved in the generation of S-adenosylmethionine, the primary methyl group donor, and hence, GGH modulation may affect DNA methylation. DNA methylation is an important epigenetic determinant in gene expression, in the maintenance of DNA integrity and stability, and in chromatin modifications, and aberrant or dysregulation of DNA methylation has been mechanistically linked to the development of human diseases including cancer. Using a recently developed in vitro model of GGH modulation in HCT116 colon and MDA-MB-435 breast cancer cells, we investigated whether GGH modulation would affect global and gene-specific DNA methylation and whether these alterations were associated with significant gene expression changes. In both cell lines, GGH overexpression decreased global DNA methylation and DNA methyltransferase (DNMT) activity, while GGH inhibition increased global DNA methylation and DNMT activity. Epigenomic and gene expression analyses revealed that GGH modulation influenced CpG promoter DNA methylation and gene expression involved in important biological pathways including cell cycle, cellular development, and cellular growth and proliferation. Some of the observed altered gene expression appeared to be regulated by changes in CpG promoter DNA methylation. Our data suggest that the GGH modulation-induced changes in total intracellular folate concentrations and content of long-chain folylpolyglutamates are associated with functionally significant DNA methylation alterations in several important biological pathways.

Introduction

Folate, a water-soluble B-vitamin, mediates the transfer of one-carbon units involved in thymidylate and purine biosynthesis and biological methylation reactions (Kim 2007; Shane 2010). Intracellular folate homeostasis is maintained by folylpolyglutamate synthase (FPGS) that facilitates intracellular retention of folate by polyglutamylation and by γ-glutamyl hydrolase (GGH) that catalyzes the hydrolysis of polyglutamylated folate into monoglutamates, thereby facilitating the export of folate out of the cell (Shane 2010). Polyglutamylated folates are better retained within cells than monoglutamate counterparts and are better substrates for folate-dependent enzymes (Moran 1999).

Folate, in the form of 5-methyltetrahydrofolate, is involved in remethylation of homocysteine to methionine, which is a precursor of S-adenosylmethionine (SAM), the primary methyl group donor for most biological methylation reactions including DNA methylation (Kim 2004, 2005; Ly et al. 2012). Folate deficiency and excess have been shown to affect DNA methylation in a gene, site, and cell-specific manner (Ly et al. 2012). DNA methylation is an important epigenetic determinant in gene expression (an inverse relation except for few exceptions), in the maintenance of DNA integrity and stability, in chromatin modifications, and in the development of mutations (Kulis and Esteller 2010). DNA methylation of cytosine in the cytosine–guanine dinucleotide sequences (CpG) is a heritable, tissue- and species-specific, post-synthetic epigenetic modification of mammalian DNA (Kulis and Esteller 2010). Seventy to 80 % of all CpG sites in human DNA are normally methylated (Kulis and Esteller 2010). However, this methylation occurs primarily in the bulk of the genome where CpG density is low, including exons, noncoding regions, and repeat DNA sites, and allows correct organization of chromatin in active and inactive states (Kulis and Esteller 2010). By contrast, most CpG-rich areas clustered in small stretches of DNA termed “CpG islands”, which span the 5′ end of approximately half of the human genes including the promoter, untranslated region, and exon 1, are unmethylated in normal cells, thereby allowing transcription (Herman and Baylin 2003; Kulis and Esteller 2010). When methylated, CpG islands can correlate with stable heritable transcriptional silencing (Herman and Baylin 2003; Kulis and Esteller 2010). DNA methylation is a dynamic process between active methylation, mediated by CpG DNA methyltransferases (DNMT) using SAM as the methyl donor, and removal of methyl groups from 5-methylcytosine residues by both passive and active mechanisms (Kulis and Esteller 2010; Li 2000). After transfer of the methyl group, SAM is converted to S-adenosylhomocysteine, a potent inhibitor of most SAM-dependent methyltransferases (Kulis and Esteller 2010; Li 2000).

Polyglutamylation is important in DNA methylation considering that polyglutamylated folates are better substrates for methylenetetrahydrofolate reductase and methionine synthase, both of which are involved in the generation of SAM (Kim 2007; McGuire and Bertino 1981; Moran 1999). Therefore, GGH modulation may affect DNA methylation by changes in total intracellular folate concentrations and by alterations in contents of polyglutamylated folate cofactors with consequent functional ramifications. Aberrant or dysregulation of DNA methylation is mechanistically related to the development of several human diseases including cancer, and hence, GGH overexpression/inhibition may have significant implications in human health and disease via its effects on DNA methylation (Ly et al. 2012). We have recently generated a novel in vitro model of GGH modulation in HCT116 colon and MDA-MB-435 breast cancer cells with predictable functional consequences (Kim et al. 2013). Using this model, we tested the hypothesis that GGH modulation would affect global and gene-specific DNA methylation with consequent functional effects on gene expression.

Materials and methods

In vitro model of GGH overexpression and inhibition

We have previously developed and functionally characterized an in vitro model of GGH overexpression and inhibition in HCT116 colon and MDA-MB-435 breast cancer cells (Kim et al. 2013). Cells overexpressing GGH had significantly higher GGH protein expression and activity, lower total intracellular folate concentrations, lower content of long-chain folylpolyglutamates, slower growth rate, lower thymidylate synthase activity, and lower dihydrofolate reductase protein expression and activity compared with controls expressing endogenous GGH (Kim et al. 2013). In contrast, cells in which GGH is inhibited had significantly lower GGH protein expression and activity, higher concentrations of total intracellular folate, higher content of long-chain folylpolyglutamates, faster growth rate, higher thymidylate synthase activity, and higher dihydrofolate reductase protein expression and activity compared with controls expressing endogenous GGH (Kim et al. 2013). Cells were grown in RPMI-1640 medium (Invitrogen) containing 2.3 µmol/L folic acid supplemented with 10 % fetal bovine serum, 500 µg/mL Geneticin^®, 50 U/mL penicillin with 50 µg/mL streptomycin, and 0.25 µg/mL fungizone amphotericin B. Cell cultures were maintained at 37 °C in 5 % CO₂. Cells were harvested at 80 % confluence and were processed for subsequent analyses.

Global DNA methylation analysis

Total genomic DNA was extracted by a standard technique using proteinase K followed by organic extraction (Laird et al. 1991). The size of DNA estimated by agarose gel electrophoresis was >20 kb in all instances. The final preparations had an A₂₆₀ to A₂₈₀ ratio between 1.8 and 2.0 and were free of RNA and protein contaminations. The concentration of each DNA sample was determined as the mean of three independent spectrophotometric readings. Global CpG DNA methylation was determined by the in vitro methyl acceptance assay using [³H-methyl] SAM (New England Nuclear) as a methyl donor and a prokaryotic CpG DNMT, Sss1 (New England Biolabs), as previously described (Sohn et al. 2004). The manner in which this assay is performed produces an inverse relationship between the endogenous DNA methylation status and exogenous [³H-methyl] incorporation. All analyses were performed in quadruplicate and repeated using two independent cell lysates.

DNMT activity assay

Total cellular CpG DNMT activity was measured by incubating cell lysate containing 10 μg of protein with 0.5 μg of poly[d(I–C)·d(I–C)] template (Sigma-Aldrich), 3 μCi [³H]-SAM (New England Nuclear), and lysis buffer in a total volume of 20 μL for 2 h at 37 °C as previously described (Sohn et al. 2004). Each reaction was performed in triplicate, and the assay was repeated three times.

Gene-specific promoter CpG DNA methylation analyses

The Illumina Infinium HumanMethylation27 (HM27) BeadChip (Illumina) was used to interrogate the DNA methylation status of 27,578 individual CpG sites located at promoter regions of 14,495 genes (Bibikova et al. 2009). Briefly, 1 µg of genomic DNA was bisulfite-converted using the EZ-96 DNA Methylation Kit (Zymo Research) according to the manufacturer’s protocol. Unmethylated cytosines are deaminated to uracil in the presence of bisulfite, while methylated cytosines are refractory to the effects of bisulfite and remain as cytosine. The bisulfite conversion included a thermocycling program with a short denaturation step (16 cycles of 95 °C for 30 s followed by 50 °C for 1 h). The amount of bisulfite-converted DNA and completeness of bisulfite conversion were assessed using a panel of MethyLight-based quality control reactions as previously described (Campan et al. 2009). All of the samples passed the quality control tests and were used for the Infinium DNA methylation assay. A measure of the level of DNA methylation at each CpG site was scored as beta (β)-values. DNA methylation β-values represent the ratio of the intensity of the methylated bead type to the combined locus intensity ranging from 0 to 1. Values close to 0 indicate low levels of DNA methylation, while values close to 1 indicate high levels of DNA methylation (Bibikova et al. 2009). The detection P-values measure the difference of the signal intensities at the interrogated CpG site compared with those from a set of 16 negative control probes embedded in the assay. We identified all data points with a detection P-value >0.05 as not statistically significantly different from background measurements and therefore were not considered trustworthy measures of DNA methylation. These data points were replaced by “NA” values as previously described (Noushmehr et al. 2010). Statistical analysis and data visualization were carried out using the R/Bioconductor software packages (http://www.bioconductor.org). The Illumina Infinium DNA methylation β-values were represented graphically using a heatmap, generated by the R/Bioconductor packages called heatmap.plus and MATLAB. The annotation for Illumina Infinium HM27 probes used in this study is presented in Supplement 1.

DNA hyper- or hypomethylation was calculated for each HM27 probe by subtracting the β-value of the corresponding control from the β-value of cells expressing the sense GGH cDNA (Sense) or cells transfected with the GGH-targeted small-interfering RNA (siRNA). We determined the β-value difference of 0.2 as having 99 % confidence based on intra- and inter-assay variations (Bibikova et al. 2009).

Gene expression analysis

RNA was isolated using the RNeasy Microarray Tissue Mini Kit (Qiagen) according to the manufacturer’s protocol. Total RNA was assessed for the RNA quality verification and microarray hybridization. The Agilent 2100 Bioanalyzer (Agilent Technologies), a microfluidics-based platform, was used for sizing, quantification, and quality of RNA. The RNA Integrity Number score was generated on the Agilent software. For the microarray analysis, the RNA quality for all of the samples had a RNA Integrity Number score ≥7.

The Illumina^® TotalPrepTM-96 RNA Amplification Kit (Ambion) was used for generating biotinylated, amplified cRNA according to the manufacturer’s protocol. The cRNA yield was quantified by NanoDrop (NanoDrop Technologies). A total of 750 ng of purified biotinylated cRNA generated from the samples were randomized in triplicate and used to hybridize onto the Illumina HumanHT-12 v4.0 BeadChip (Illumina). Each array on this BeadChip targets 31,335 annotated genes and includes 47,231 probes designed to cover content from NCBI RefSeq Release 38 (November 7, 2009), as well as legacy UniGene content. Twelve samples were hybridized to one slide for higher throughput and reduced sample-to-sample variability since gaskets separate each array. The BeadChip was incubated at 58 °C, with rotation speed 5 for 18 h for hybridization. After washing and staining, each BeadChip was scanned on the iScan (Illumina), a laser-based imaging system with two lasers (red: Cy5 and green: Cy3) for detecting fluorescence information. The intensity files were quantified in GenomeStudio^® (Illumina, version 2010.2) to generate intensity measurements without normalization algorithms.

For the gene expression analysis, data were checked for overall quality using R version 2.13.2 with the Bioconductor framework and the LUMI package (http://www.bioconductor.org). Data were imported in GeneSpring GX 11.5 (Agilent Technologies) and normalized using a standard quantile-based normalization followed by a “per probe” median-centered normalization. All data were log₂-transformed for analysis and visualization. Data were initially filtered in order to remove any confounding effects on subsequent analyses of probes with no signal. Only probes in the upper 80th percentile of the distribution of intensities in 100 % of samples from any of the groups were retained following filtering. Normalization and data filtering were performed separately for each subset of samples being used in a particular analysis. For comparisons between cells expressing the sense GGH (Sense) and control (Control-S) and between cells transfected with the GGH-targeted siRNA (siRNA) and control (Control-si), an unpaired t-test using a false discovery rate (FDR) Benjamini and Hochberg multiple testing correction with a P-value cutoff of 0.05 was performed.

Integrated analysis of DNA methylation and gene expression data

We merged the DNA methylation and gene expression data sets using Entrez Gene IDs for the integrated analysis to identify genes for which the differential expression might have been influenced by DNA methylation in response to GGH modulation. We used a β-value difference (|Δβ|) of 0.20 as a threshold for differential DNA methylation between Sense and Control-S and between siRNA and Control-si. This threshold of |Δβ| = 0.20 was determined previously as a stringent estimate of Δβ detection sensitivity across the range of β-values (Bibikova et al. 2009). Gene expression data with (1) a fold change greater or less than 1.3 and (2) a one-way ANOVA with a Benjamini and Hochberg-corrected P-value ≤0.05 were used for integrated analysis. We set 1.3 as a fold change not to overlook small changes in response to GGH modulation as we identified a small number of genes differentially expressed especially in the GGH-modulated HCT116 cells.

Functional analysis

The functional analysis was performed using Ingenuity Pathway Analysis (IPA, Ingenuity^® Systems; http://www.ingenuity.com) to identify biological functions and/or diseases that were most significant to genes differentially methylated and/or regulated in each system. Genes with a threshold of |Δβ| = 0.20 and/or fold change of |1.3| and that were associated with biological functions and/or diseases in the Ingenuity Knowledge Base were considered for the analysis. The right-tailed Fisher’s exact test was used to calculate P-values in determining the probability that each biological function and/or disease assigned to that data set is due to chance alone.

Quantitative reverse transcriptase-PCR (qRT-PCR)

qRT-PCR was performed to confirm the data obtained from gene expression analysis using the Illumina HumanHT-12 v4.0 BeadChip. Template RNA was reversely transcribed to cDNA using the QuantiTect Reverse Transcription Kit (Qiagen). Selected primer sequences were synthesized by the Integrated DNA Technologies and are presented in Supplement 2. The reactions were run in triplicate on MicroAmp Optical 384-well plates (Applied Biosystems, Life Technologies), and their amplifications were tracked by SYBR Green fluorescent dye (Applied Biosystems, Life Technologies). Completed plates were spun at 1,200 rpm for 2 min at 4 °C and then placed in the ViiA-7TM Real-Time PCR System (Applied Biosystems, Life Technologies). The reaction conditions for stage one were as follows: 2 min at 50 °C, followed by 10 min at 95 °C to activate the polymerase. This stage was followed by 40 cycles beginning with 15 s at 95 °C to denature the target strand followed by 1 min at 60 °C to allow for the polymerase to anneal and extend the target strand. Finally, a melt curve followed consisting of 15 s at 95 °C then 1 min at 60 °C. Relative gene expression data were analyzed using the comparative threshold (C _t) method as described previously (Livak and Schmittgen 2001).

Statistical analysis

For global DNA methylation, DNMT activity, and qRT-PCR analyses, comparisons between Sense and Control-S (GGH overexpression) and between siRNA and Control-si (GGH inhibition) were determined using the Student’s t-test function of SPSS Statistics 17.0 (IBM SPSS, Chicago, IL). The results were considered statistically significant if two-tailed P-values were <0.05.

Results

Global DNA methylation and DNMT activity

In both HCT116 and MDA-MB-435 cells, GGH overexpression was associated with significantly lower global DNA methylation (HCT116, 16 % lower, P = 0.024, Fig. 1a; MDA-MB-435, 22 % lower, P < 0.001, Fig. 1c) and lower DNMT activity (HCT116, 66 % lower, P < 0.001, Fig. 1e; MDA-MB-435, 26 % lower, P = 0.003, Fig. 1g) compared with controls expressing endogenous GGH.

Fig. 1

Effects of GGH modulation on global DNA methylation and DNMT activity in HCT116 and MDA-MB-435 cells. In both cell lines, GGH overexpression was associated with significantly lower global DNA methylation (HCT116, a; MDA-MB-435, c) and lower DNMT activity (HCT116, e; MDA-MB-435, g) compared with controls expressing endogenous GGH, whereas GGH inhibition showed significantly higher global DNA methylation (HCT116, b; MDA-MB-435, d) and higher DNMT activity (HCT116, f; MDA-MB-435, h) compared with controls expressing endogenous GGH. The in vitro methyl acceptance assay for global DNA methylation measurement produces an inverse relationship between the endogenous DNA methylation status and exogenous [³H-methyl] incorporation into DNA. All analyses were performed in quadruplicate and repeated using two independent cell lysates (a–d). The assay for DNMT activity produces a positive relationship between the endogenous enzyme activity and exogenous [³H-methyl] incorporation into DNA. Each reaction was performed in triplicate, and the assay was repeated three times (e–h). Control-S, cells expressing endogenous GGH; Sense, cells transfected with the sense GGH cDNA; Control-si, cells expressing endogenous GGH; siRNA, cells transfected with the GGH-targeted siRNA. *P < 0.05; **P < 0.01; ***P < 0.001 compared with corresponding control by the Student’s t-test. Values are mean ± SD

In contrast, GGH inhibition showed significantly higher global DNA methylation (HCT116, 15 % higher, P = 0.013, Fig. 1b; MDA-MB-435, 7 % higher, P = 0.013, Fig. 1d) and higher DNMT activity (HCT116, 47 % higher, P < 0.001, Fig. 1f; MDA-MB-435, 27 % higher, P = 0.002, Fig. 1h) compared with controls expressing endogenous GGH in both cell lines.

Effect of GGH modulation on gene-specific CpG promoter DNA methylation

We then determined whether GGH modulation would affect gene-specific DNA methylation, thereby identifying biological pathways that are primarily affected by the GGH modulation-induced DNA methylation alterations, in HCT116 and MDA-MB-435 cell lines. Scatter plots of DNA methylation β-values showed differentially methylated loci between cells overexpressing GGH (Sense) and corresponding controls (Control-S) and between cells with GGH inhibition (siRNA) and corresponding controls (Control-si). MDA-MB-435 breast cancer cells showed more CpG methylation alterations in response to GGH modulation than HCT116 colon cancer cells (Supplement 3).

Genes differentially methylated in the GGH-modulated HCT116 cells

In the HCT116 cell line, we identified 905 genes that were differentially methylated (546 hypermethylated and 359 hypomethylated) in response to GGH overexpression, while 1,869 genes were differentially methylated (998 hypermethylated and 871 hypomethylated) in response to GGH inhibition (Fig. 2a). We performed a functional analysis using IPA to identify biological and disease processes involving the GGH modulation-induced differentially methylated genes in each system, and these results are presented in Supplement 4. In the GGH-overexpressed HCT116 cells, the hypermethylated genes were associated with cell morphology, cellular development, and antigen presentation, while the hypomethylated genes were involved in molecular transport, cellular assembly and organization, and cell-to-cell signaling and interaction (Supplement 4). As for the GGH-inhibited HCT116 cells, major function categories of the hypermethylated genes included cell-to-cell signaling and interaction, amino acid metabolism, and drug metabolism, whereas those of the hypomethylated genes consisted of lipid metabolism, molecular transport, and small molecule biochemistry (Supplement 4).

Fig. 2

Number of genes differentially methylated and expressed in the GGH-modulated HCT116 and MDA-MB-435 cells. Sixty-one hyper- and 54 hypomethylated genes were common between the HCT116 and MDA-MB-435 cell lines in response to GGH overexpression, whereas 117 hyper- and 129 hypomethylated genes were common between these cell lines in response to GGH inhibition (a). DNA hyper- or hypomethylation was calculated by subtracting the β-value of corresponding control from the β-value of Sense or siRNA. We determined the β-value difference of 0.2 as having 99 % confidence based on intra- and inter-assay variations (a). Nine down- and 10 upregulated genes were common between two cell lines after GGH overexpression, whereas 11 down- and 12 upregulated genes were common in these cell lines in response to GGH inhibition (b). We determined the number of genes with a fold change >1.3 or <−1.3 using an unpaired t-test with the FDR corrected P-value ≤0.05 (b)

Genes differentially methylated in the GGH-modulated MDA-MB-435 cells

In the MDA-MB-435 cell line, we identified 2,394 genes that were differentially methylated (1,058 hypermethylated and 1,336 hypomethylated) in response to GGH overexpression, while 2,666 genes were differentially methylated (1,122 hypermethylated and 1,544 hypomethylated) in response to GGH inhibition (Fig. 2a). Genes with functions relating to cell-to-cell signaling and interaction, cellular movement, and molecular transport were differentially methylated in the GGH-overexpressed MDA-MB-435 cells (Supplement 4). In the GGH-inhibited MDA-MB-435 cells, the hypermethylated genes were associated with cellular movement, cell-to-cell signaling and interaction, and cell morphology, while the hypomethylated genes were involved in cell signaling, molecular transport, and vitamin and mineral metabolism (Supplement 4).

Genes differentially methylated in both the GGH-modulated HCT116 and MDA-MB-435 cells

We also identified genes displaying differential DNA methylation changes in response to GGH modulation in both HCT116 and MDA-MB-435 cell lines. In GGH overexpression system, 61 genes involved in cell death, cell cycle, cell morphology, cellular assembly and organization, and cellular compromise were hypermethylated in both cell lines (Fisher’s exact test P = 1.0 × 10⁻³), while 54 genes associated with cellular assembly and organization, cellular development, cellular growth and proliferation, amino acid metabolism, and cell cycle were hypomethylated in both cell lines (P = 3.0 × 10⁻⁴) (Fig. 2a). In the GGH inhibition system, 117 hypermethylated genes with functions relating to cell signaling, small molecule biochemistry, carbohydrate metabolism, cell-to-cell signaling and interaction, and lipid metabolism were common between two cell lines (P = 5.0 × 10⁻⁶), whereas 129 genes associated with cellular movement, antigen presentation, lipid metabolism, molecular transport, and small molecule biochemistry were commonly hypomethylated in both cell lines (P = 8.9 × 10⁻⁵) (Fig. 2a). The list of top networks generated by mapping the focus genes that displayed differential DNA methylation in both the GGH-modulated HCT116 and MDA-MB-435 cells is presented in Supplement 4. The list of genes commonly differentially methylated in both cell lines in response to GGH modulation is presented in Supplement 5.

Effect of GGH modulation on gene expression

We investigated whether GGH modulation would affect gene expression profile in HCT116 and MDA-MB-435 cell lines using the Illumina HT-12 platform. Similar to DNA methylation changes, MDA-MB-435 breast cancer cells showed more number of genes with altered expression in response to GGH modulation than HCT116 colon cancer cells (Fig. 2b).

Genes differentially expressed in the GGH-modulated HCT116 cells

In the HCT116 cell line, we identified 152 genes that were differentially expressed (91 downregulated and 61 upregulated) in response to GGH overexpression, while 321 genes were differentially expressed (139 downregulated and 182 upregulated) in response to GGH inhibition (Fig. 2b). As a result of the classification according to function using IPA, genes involved in cellular movement, cell death, and carbohydrate metabolism were differentially expressed in the GGH-overexpressed HCT116 cells, while genes associated with cell death, cell cycle, and cellular movement were identified in the HCT116 cells in which GGH is inhibited (Supplement 6).

Genes differentially expressed in the GGH-modulated MDA-MB-435 cells

In the MDA-MB-435 cell line, we identified 1,383 genes that were differentially expressed (628 downregulated and 755 upregulated) in response to GGH overexpression, while 859 genes were differentially expressed (402 downregulated and 457 upregulated) in response to GGH inhibition (Fig. 2b). Genes with functions relating to cell-to-cell signaling and interaction, cellular movement, and cell death were identified in response to GGH overexpression, while genes associated with cell death, cellular development, and cellular growth and proliferation were differentially expressed in response to GGH inhibition (Supplement 6).

Genes differentially expressed in both the GGH-modulated HCT116 and MDA-MB-435 cells

We determined genes displaying differential gene expression alterations associated with GGH modulation in both HCT116 and MDA-MB-435 cell lines. We identified nine genes that were downregulated and ten genes that were upregulated in HCT116 and MDA-MB-435 cell lines in response to GGH overexpression (Fig. 2b). We also identified 11 genes that were downregulated and 12 genes that were upregulated in both cell lines in response to GGH inhibition (Fig. 2b). In the GGH-overexpressed HCT116 and MDA-MB-435 cell lines, the commonly downregulated genes were associated with cellular growth and proliferation, cell cycle, gene expression, carbohydrate metabolism, and cellular function and maintenance, whereas the commonly upregulated genes were related to cell cycle, cellular development, molecular transport, small molecule biochemistry, and cellular growth and proliferation (Fig. 2b). In the GGH-inhibited HCT116 and MDA-MB-435 cell lines, the major function categories of the commonly downregulated genes included energy production, lipid metabolism, small molecule biochemistry, cell morphology, and cellular development, while those of genes upregulated in common consisted of cell cycle, cell death, cellular development, cellular growth and proliferation, and cell morphology (Fig. 2b). The list of top networks generated by mapping the focus genes that were commonly differentially expressed in both the GGH-modulated HCT116 and MDA-MB-435 cells is presented in Supplement 6. The list of genes commonly differentially expressed in both cell lines in response to GGH modulation is presented in Fig. 2b.

Integrated analysis of gene expression and DNA methylation changes

We performed an integrated analysis of differentially methylated and expressed genes in response to GGH modulation in order to identify genes, expression of which was regulated by promoter DNA methylation. We detected 21 and 148 genes whose expression was inversely regulated by promoter DNA methylation changes in the GGH-modulated HCT116 and MDA-MB-435 cells, respectively (Fig. 3; Table 1). There were no common genes for which gene expression alterations were correlative with promoter DNA methylation changes in response to GGH modulation between the two cell lines. The list of genes with altered promoter DNA methylation and expression in the GGH-modulated HCT116 and MDA-MB-435 cell lines is presented in Supplement 7.

Fig. 3

Integrated analysis of gene expression and promoter DNA methylation changes in the GGH-modulated HCT116 and MDA-MB-435 cells. In HCT116 colon cancer cells, we identified 5 (a) and 16 genes (b) for which gene expression was inversely regulated by CpG promoter DNA methylation changes in response to GGH overexpression and inhibition, respectively. In MDA-MB-435 breast cancer cells, we identified 101 (c) and 47 genes (d) for which gene expression was inversely regulated by CpG promoter DNA methylation changes in response to GGH overexpression and inhibition, respectively. The β-value difference and log₂-transformed gene expression value difference between Sense and Control-S and between siRNA and Control-si are plotted on x- and y-axes, respectively. Red data points highlight those genes that are hypermethylated with β-value difference >0.2 and show <−1.3 fold change in their expression levels, while blue data points indicate those genes that are hypomethylated with β-value difference <−0.2 and show >1.3 fold change in their expression levels

Table 1 Number of genes and the top molecular and cellular functions associated with genes with altered expression and promoter DNA methylation in the GGH-modulated HCT116 and MDA-MB-435 cells

In HCT116 cells that overexpressed GGH, FGFBP1 (fibroblast growth factor binding protein 1), a gene involved in cellular growth and proliferation, displayed promoter DNA hypomethylation and increased gene expression (Supplement 7). In the GGH-inhibited HCT116 cells, epigenetically silenced genes were involved in cell-to-cell signaling and interaction, cellular development, cellular growth and proliferation, gene expression, and lipid metabolism, while hypomethylated and upregulated genes were associated with lipid metabolism, small molecule biochemistry, carbohydrate metabolism, cellular movement, and drug metabolism (Table 1).

In the MDA-MB-435 cells that overexpressed GGH, downregulated genes with promoter CpG hypermethylation were associated with cell signaling, cellular assembly and organization, cellular movement, drug metabolism, and lipid metabolism, while hypomethylated and upregulated genes were involved in cellular assembly and organization, cell-to-cell signaling and interaction, cellular growth and proliferation, cellular development, and cell death (Table 1). In the MDA-MB-435 cells in which GGH is inhibited, hypermethylated and downregulated genes were involved in cell morphology, cellular development, gene expression, cellular assembly and organization, and cell-to-cell signaling and interaction, while hypomethylated and upregulated genes were associated with cell morphology, cell cycle, cellular growth and proliferation, cellular development, and cellular movement (Table 1). Cellular growth and proliferation was the common function of hypomethylated and upregulated genes in both the GGH-overexpressed HCT116 and MDA-MB-435 cells. In both cell lines with GGH inhibition, cell-to-cell signaling and interaction, cellular development, and gene expression were the common functions associated with hypermethylated and downregulated genes, while cellular movement was the common function associated with hypomethylated and upregulated genes (Table 1).

Validation of gene expression by qRT-PCR

We performed qRT-PCR to validate the gene expression results of selected genes that were inversely regulated by promoter DNA methylation changes. We selected genes that were associated with a great magnitude of fold change in gene expression on microarray analysis and possessed relevant biological function of interest including cancer, cell cycle, proliferation, and apoptosis. As presented in Supplement 8, although the magnitude of change was different, the direction of change in gene expression in response to GGH modulation was consistent between Illumina gene expression arrays and qRT-PCR analyses in both cell lines (P < 0.05), thereby validating the microarray data.

GGH-specific gene expression analysis

To investigate genes whose altered expression might be GGH modulation-specific, we identified genes differentially expressed in the opposite direction between GGH overexpression and inhibition (Supplement 9). In the HCT116 cell line, 28 genes that were downregulated in response to GGH overexpression and upregulated in response to GGH inhibition were associated with cellular movement, cell death, cellular growth and proliferation, cell-to-cell signaling and interaction, and RNA post-transcriptional modification (Fig. 4a). Eleven genes were upregulated in response to GGH overexpression and downregulated in response to GGH inhibition in HCT116 cells, and these genes were involved in carbohydrate metabolism, cellular function and maintenance, small molecule biochemistry, lipid metabolism, and molecular transport (Fig. 4a). The list of the top genes associated with the GGH-specific altered expression in the GGH-modulated HCT116 cells is presented in Table 2. The list of top networks generated by mapping the focus genes associated with the GGH-specific altered expression in the GGH-modulated HCT116 cells is presented in Supplement 10.

Fig. 4

Number of genes and the top molecular and cellular functions associated with the GGH-specific gene expression. We identified genes differentially expressed in the opposite direction between GGH overexpression and inhibition to investigate genes associated with the GGH-specific altered expression in the GGH-modulated HCT116 and MDA-MB-435 cell lines. Twenty-eight (a) and 80 genes (b) that were downregulated in response to GGH overexpression and upregulated in response to GGH inhibition were identified in HCT116 and MDA-MB-435 cell lines, respectively. We identified that 11 (a) and 133 genes (b) were upregulated in response to GGH overexpression and downregulated associated with GGH inhibition in HCT116 and MDA-MB-435 cell lines, respectively

Table 2 List of the top genes associated with the GGH-specific altered expression in the GGH-modulated HCT116 colon cancer cells

In the MDA-MB-435 cell line, 80 genes were downregulated in response to GGH overexpression and upregulated in response to GGH inhibition, and these genes were involved in cellular compromise, cell death, cell morphology, cellular movement, and cellular development (Fig. 4b). One hundred and thirty-three genes that were upregulated in response to GGH overexpression and downregulated in response to GGH inhibition were associated with cellular compromise, cell cycle, cellular assembly and organization, cellular function and maintenance, and cellular movement in MDA-MB-435 cells (Fig. 4b). The list of the top genes associated with the GGH-specific altered expression in the GGH-modulated MDA-MB-435 cells is presented in Table 3. The list of top networks generated by mapping the focus genes associated with the GGH-specific altered expression in the GGH-modulated MDA-MB-435 cells is presented in Supplement 10.

Table 3 List of the top genes associated with the GGH-specific altered expression in the GGH-modulated MDA-MB-435 breast cancer cells

We identified several genes with GGH-specific altered expression changes which were regulated by promoter DNA methylation. TYR was one of the most downregulated (fold change −31.95) and hypermethylated (β-value difference 0.26) genes in the GGH-overexpressed MDA-MB-435 cells, and it was upregulated (fold change 4.07) in the GGH-inhibited MDA-MB-435 cells (Table 3; Supplement 7). TYR encodes tyrosinase, a melanosomal enzyme that catalyzes the rate-limiting steps of melanin biosynthesis (Spritz et al. 2003). BCHE, which encodes butyrylcholinesterase, was associated with upregulation and hypomethylation in the GGH-inhibited MDA-MB-435 cells, and it was downregulated in response to GGH overexpression in the same cell line (Table 3; Supplement 7). Furthermore, in the GGH-overexpressed MDA-MB-435 cells, we identified upregulated and hypomethylated genes that were downregulated in the GGH-inhibited MDA-MB-435 cells including S100A4 and NNMT (Table 3; Supplement 7). In MDA-MB-435 cells, CDC42EP5, CDC42 effector protein (Rho GTPase binding) 5, was hypomethylated (β-value difference −0.58) and upregulated (fold change 3.76) in GGH overexpression, while it was hypermethylated (β-value difference 0.24) and downregulated (fold change −2.53) in GGH inhibition (Table 3; Supplement 7). No common genes were identified in both HCT116 and MDA-MB-435 cell lines.

Discussion

In the present study, we investigated the effects of GGH modulation on global and gene-specific DNA methylation and gene expression using a recently developed in vitro model of GGH overexpression and inhibition in HCT116 and MDA-MB-435 cells with predictable functional consequences (Kim et al. 2013). Generally, most of the observed functional consequences of GGH overexpression and inhibition in this model were consistent with the known biological function of GGH and provided an appropriate in vitro model to test the effect of GGH modulation on DNA methylation (Kim et al. 2013).

We first investigated whether the changes in total intracellular folate concentrations and in the content of polyglutamylated folate cofactors resulting from GGH modulation would affect global DNA methylation and DNMT activity in HCT116 colon and MDA-MB-435 breast cancer cells. In both HCT116 and MDA-MB-435 cells, GGH overexpression was associated with decreased global DNA methylation and DNMT activity, which is likely related to the lower total intracellular folate concentrations and lower content of long-chain folylpolyglutamates in response to GGH overexpression (Kim et al. 2013). In contrast, GGH inhibition demonstrated increased global DNA methylation and DNMT activity. This finding is likely accounted for by the higher total intracellular folate concentrations and higher content of long-chain folylpolyglutamates in response to GGH inhibition (Iacopetta et al. 2008; Kim et al. 2013). Indeed, low GGH expression was associated with CpG island methylator phenotype (CIMP+) colorectal cancer, a subgroup of colorectal cancer with concurrent hypermethylation of a large number of CpG islands mostly of tumor suppressor genes (Toyota et al. 1999), and increased folate intermediates in colorectal cancer (Kawakami et al. 2008). These observations suggest that low GGH expression may be linked to increased promoter methylation in CIMP+ tumors by causing elevation of folate concentrations and that aberrant GGH in the folate metabolic pathway may be involved in CIMP+ colorectal cancer carcinogenesis (Kawakami et al. 2008). Collectively, expression of GGH appears to be one of the important factors in determining DNA methylation status, probably through an effect on folate metabolism (Kawakami et al. 2003, 2008).

Epigenomic and gene expression analyses data suggest that GGH modulation influenced promoter CpG DNA methylation and gene expression involved in important biological pathways. MDA-MB-435 cells revealed more promoter CpG methylation and gene expression alterations in response to GGH modulation compared with HCT116 cells. GGH inhibition demonstrated greater promoter CpG methylation changes compared with GGH overexpression in both cell lines. In both HCT116 and MDA-MB-435 cells, differentially methylated genes associated with cell cycle and cellular assembly and organization were affected by GGH overexpression, while genes involved in small molecule biochemistry and lipid metabolism were differentially methylated in response to GGH inhibition. In both cell lines, expression of genes involved in cellular growth and proliferation and cell cycle were affected by GGH overexpression, whereas genes related to cell morphology and cellular development were differentially expressed in response to GGH inhibition. However, some of the pathways affected by GGH modulation might be indirect consequences of changes in total folate content and polyglutamylated folate cofactors or of changes in cell proliferation or morphology. Overall, only a small number of genes were associated with the inverse relationship between promoter DNA methylation and gene expression in response to GGH modulation. This observation is similar to reports from other studies (Houshdaran et al. 2010; Pike et al. 2008). It appears that the genetic, other epigenetic mechanisms such as histone modifications, chromatin remodeling, and RNA interference, and/or changes in activities of other genes for which expression was directly influenced by promoter DNA methylation are likely responsible for the observed gene expression changes in response to GGH modulation in the present study.

In MDA-MB-435 cells, we identified genes that were upregulated and hypomethylated in response to GGH overexpression and downregulated in response to GGH inhibition including S100A4 and NNMT. S100A4 encodes a member of the S100 family of calcium-binding proteins. S100 family members have a wide range of intracellular functions, including the regulation of homeostasis, protein phosphorylation, cytoskeletal rearrangements, and transcriptional activity, and extracellular functions such as the regulation of cell proliferation and activation, apoptosis, and chemotaxis (Garrett et al. 2006; Santamaria-Kisiel et al. 2006). NNMT encodes nicotinamide N-methyltransferase, an enzyme that catabolizes nicotinamide and other pyridine compounds in a reaction that uses the methyl group generated during the conversion of SAM to S-adenosylhomocysteine, involved in the biotransformation of many drugs and xenobiotic compounds (Aksoy et al. 1994). A recent study found that downregulation of NNMT inhibited proliferation in KB cancer cells, suggesting NNMT might be a target for therapeutics and could alter the efficacy of standard chemotherapeutic drugs (Pozzi et al. 2011). Furthermore, in MDA-MB-435 cells, CDC42EP5 was hypomethylated and upregulated in GGH overexpression, while it was hypermethylated and downregulated in GGH inhibition. CDC42EP5, CDC42 effector protein 5, is involved in the organization of the actin cytoskeleton and acts downstream of CDC42 to induce actin filament assembly leading to cell shape changes (Hirsch et al. 2001).

In MDA-MB-435 cells overexpressing GGH, differential expression of several genes involved in folate biosynthesis and one-carbon metabolism including DPYD, ABCC5, and SLC25A32 is in line with the GGH modulation-induced changes in cellular folate homeostasis and drug efficacy reported in the previous study (Kim et al. 2013). DPYD encodes dihydropyrimidine dehydrogenase that is the rate-limiting enzyme involved in 5-fluorouracil (5FU) metabolism and a major determinant of 5FU efficacy (Oguri et al. 2005; Yoshinare et al. 2003). DPYD overexpression in cancer cell lines is associated with 5FU resistance (Takebe et al. 2001), and high DPYD mRNA expression in colorectal tumors has been shown to correlate with resistance to 5FU (Salonga et al. 2000). In contrast, in CIMP+ colorectal cancers, promoter methylation-induced silencing of DPYD was associated with an enhanced response to 5FU (Iacopetta et al. 2008). Thus, it appears that the observed upregulation of DPYD in the GGH-overexpressed MDA-MB-435 cells might be associated with 5FU resistance. Indeed, we have shown that GGH overexpression decreased 5FU efficacy in MDA-MB-435 cells (Kim et al. 2013).

Multidrug-resistance-associated protein 5 (MRP5) encoded by the ABCC5 gene is associated with resistance to antifolates and 5FU. MRP5 can efflux mono- and diglutamate forms of methotrexate (MTX) and transport 5-fluoro-2-deoxyuridine-5-monophosphate, a metabolite of 5FU (Assaraf 2007; Hooijberg et al. 2006; Pratt et al. 2005; Wielinga et al. 2005). MRPs contribute to drug resistance or increase drug efficacy depending on polyglutamylation of antifolates and intracellular folate concentrations (Assaraf 2006). In addition, the SLC25A32 gene encodes a folate transporter that shuttles folate from the cytoplasm into the mitochondria (Titus and Moran 2000). A reduced, monoglutamylated form of cytoplasmic folate (probably tetrahydrofolate or 5-formyltetrahydrofolate) is transported to the mitochondria by the mitochondrial folate transporter, followed by the mitochondrial FPGS-induced polyglutamylation resulting in mitochondrial folate accumulation (Chen et al. 1996). Taken together, downregulation of ABCC5 and SLC25A32 associated with GGH overexpression supports the role of MRP5 and the mitochondrial folate transporter in the modulation of the intracellular folate levels as well as cellular folate homeostasis since we have previously found that GGH overexpression decreased total intracellular folate concentrations in MDA-MB-435 cells (Kim et al. 2013).

In addition, CDK2 encodes cyclin-dependent kinase 2 that maintains a balance of S-phase regulatory proteins and thereby coordinates subsequent p53-independent G2/M checkpoint activation (Chung and Bunz 2010). CDK2 expression was downregulated in the 5FU resistant cell lines, suggesting that decreased CDK2 activity may delay the transition of resistant cells from G1 into S-phase (Guo et al. 2008; Wang et al. 2004). Delayed S-phase entry and/or reduced S-phase traverse may provide resistant cells with enough time to repair 5FU-induced damage before progressing to G2-M phase. Thus, 5FU resistance may be, at least partially, reversed by specific targeting of the G1-S checkpoint arrest in the resistant cells (Wang et al. 2004). Accordingly, downregulation of CDK2 in the GGH-overexpressed MDA-MB-435 cells is likely to be associated with decreased 5FU efficacy (Kim et al. 2013). Restoration of the G1 checkpoint by targeting CDK2 is currently one of the major strategies for anticancer drug development (Malumbres and Barbacid 2001).

GGH, along with a tightly coupled counter-regulation by FPGS, plays an important role in maintenance of optimal intracellular folate concentrations and polyglutamylated forms for critical folate-dependent one-carbon transfer reactions involved in nucleotide biosynthesis and biological methylation reactions (Shane 2010). Intracellular folate depletion/excess as well as altered folylpolyglutamate distribution leads to perturbations in the nucleotide synthesis and biological methylation pathways (Shane 1989, 2010). Indeed, folate deficiency has been linked to the development of anemia, coronary heart disease, neural tube defects and other congenital disorders, cognitive impairments, and cancers, primarily through aberrant DNA synthesis, stability, integrity, repair, and methylation (Kim 2005, 2007). Furthermore, folate excess has been shown to exert adverse health effects, including tumor promotion and metabolic syndrome, likely via aberrant nucleotide biosynthesis and biological methylation reactions (Kim 2007; Yajnik et al. 2008). Given the critical role of GGH in intracellular folate homeostasis, dysregulation of GGH and consequent perturbations in intracellular folate concentrations and folylpolyglutamate distribution will likely contribute to the development and progression of certain human diseases that are linked to aberrant folate and one-carbon metabolism.

The role of GGH in human health and disease, however, has not yet been clearly demonstrated. An emerging body of studies has reported several single nucleotide polymorphisms (SNPs) in GGH significantly influence GGH expression and activity (Chave et al. 2003; Cheng et al. 2004; Dervieux et al. 2004; DeVos et al. 2008). It has been recently reported that some GGH SNPs are associated with an increased risk of cancer as well as a decreased risk of cardiovascular disease (DeVos et al. 2008; Organista-Nava et al. 2010; Wernimont et al. 2012). In addition to SNPs, CpG methylation in the human GGH promoter region has been shown to significantly alter GGH expression and activity (Cheng et al. 2006). Our data demonstrate that GGH modulation associated with significant changes in intracellular folate concentrations and folylpolyglutamate distribution can lead to altered expression and CpG DNA methylation of genes involved in important biological pathways that might contribute to the development and progression of diseases that are linked to aberrant intracellular folate and one-carbon metabolism. Whether or not recently identified and characterized GGH SNPs or epigenetic variants with functionally significant effects on GGH expression and activity may be associated with altered expression and CpG methylation of genes involved in important biological pathways needs to be confirmed in future clinical and epidemiologic studies. This is an important issue as individuals with certain GGH SNPs may be at risk of developing diseases that are linked to aberrant intracellular folate and one-carbon metabolism and specific forms and amount of folate or folic acid may be required to nullify the risk.

The role of GGH in cancer risk, progression, and treatment has been under intense investigation. Similar to folates, intracellular homeostasis of antifolates such as MTX are regulated by GGH (Cheng et al. 2004; Panetta et al. 2002; Rhee et al. 1999; Yao et al. 1995). Alterations in concentrations and polyglutamylation of both folate and antifolates resulting from GGH modulation can significantly influence cancer risk, progression and treatment response to antifolates (Kim et al. 2013). In addition to the previously mentioned association of GGH expression with CIMP+ colorectal cancer (Kawakami et al. 2008), GGH expression may be a prognostic indicator of cancer. For example, a recent study has reported that high GGH protein level is associated with poor prognosis and unfavorable clinical outcomes in patients with invasive breast cancer (Shubbar et al. 2013). Furthermore, clinical studies investigating the role of GGH in modulating chemosensitivity to 5FU, the cornerstone of colorectal and breast cancer chemotherapy, and antifolates have begun to emerge. Low GGH expression was reported to correlate with an enhanced response to 5FU-based chemotherapy in patients with metastatic colorectal cancer (Nakajima et al. 2008). In patients with advanced pancreatic cancer treated with an oral fluoropyrimidine derivative, S-1, low GGH expression was associated with improved overall survival (Nakamura et al. 2011). In acute myelogenous leukemia, it was found that high GGH activity may play a role in inherent drug resistance to MTX (Rots et al. 1999). High GGH expression was shown to be associated with a higher risk of developing advanced toxicity to pemetrexed, a multi-targeted antifolate, in patients with advanced breast cancer (Llombart-Cussac et al. 2007). Furthermore, several recently identified and characterized functionally significant genetic and epigenetic polymorphisms of GGH have been reported to predict response to and toxicity of antifolate-based treatment in patients with several cancers (Cheng et al. 2004; Kim et al. 2008; Koomdee et al. 2012; Silva et al. 2013; Smit et al. 2012; Wang et al. 2014) and inflammatory arthritis (Dervieux et al. 2004; Hayashi et al. 2009; Jekic et al. 2013; Owen et al. 2012; van der Straaten et al. 2007; Yanagimachi et al. 2011). Our data herein provide evidence that GGH modulation significantly influences expression and CpG DNA methylation of genes involved in important biological pathways that might account for the observed effects of GGH modulation on cancer risk, prognosis, and treatment response. Furthermore, our data provide a framework for future studies aimed at interrogating specific biological pathways associated with GGH modulation and at exploring upstream and downstream targets of GGH modulation.

In conclusion, our results suggest that the GGH modulation-induced changes in total intracellular folate concentrations and content of long-chain folylpolyglutamates significantly influenced global DNA methylation and DNMT activity as well as promoter DNA methylation and gene expression. We demonstrated that GGH overexpression was associated with decreased global DNA methylation and DNMT activity, while GGH inhibition showed increased global DNA methylation and DNMT activity. Furthermore, we showed that GGH modulation was associated with differential promoter DNA methylation and gene expression involved in important biological pathways and some of the observed altered gene expression appeared to be regulated by promoter DNA methylation changes. In the GGH-overexpressed MDA-MB-435 cells, we identified several differentially expressed genes involved in folate biosynthesis, one-carbon pool by folate, and cell cycle, which might in part have contributed to the observed decreased total intracellular folate concentrations and 5FU efficacy in response to GGH overexpression. The potential role of GGH modulation in DNA methylation and its associated downstream functional effects needs further exploration.