Genome-wide DNA methylation assay reveals novel candidate biomarker genes in cervical cancer

Abstract

The oncogenic human papilloma viruses (HPVs) are associated with precancerous cervical lesions and development of cervical cancer. The DNA methylation signatures of the host genome in normal, precancerous and cervical cancer tissue may indicate tissue-specific perturbation in carcinogenesis. The aim of this study was to identify new candidate genes that are differentially methylated in squamous cell carcinoma compared with DNA samples from cervical intraepithelial neoplasia grade 3 (CIN3) and normal cervical scrapes. The Illumina Infinium HumanMethylation450 BeadChip method identifies genome-wide DNA methylation changes in CpG islands, CpG shores and shelves. Our findings showed an extensive differential methylation signature in cervical cancer compared with the CIN3 or normal cervical tissues. The identified candidate biomarker genes for cervical cancer represent several types of mechanisms in the cellular machinery that are epigenetically deregulated by hypermethylation, such as membrane receptors, intracellular signaling and gene transcription. The results also confirm extensive hypomethylation of genes in the immune system in cancer tissues. These insights into the functional role of DNA methylome alterations in cervical cancer could be clinically applicable in diagnostics and prognostics, and may guide the development of new epigenetic therapies.

Keywords: :

• DNA methylation
• Infinium Human Methylation 450 BeadChip
• HPV
• cervical cancer
• immune system

Introduction

Cervical cancer is the third most frequent cancer among women world-wide.¹ The causative agent is a persistent infection with one of 15 oncogenic human papilloma viruses (HPVs) that leads to cervical precancerous lesions which may progress to cancer.² There are more than 100 identified HPV genotypes, 40 of them are considered to infect the genital tract and the types 16, 18, and 45 are found to be the most common in cervical cancer.³

The majority of the viral load is cleared by the cell-mediated immunity, but if not cleared the persistent infection affects the host cell apoptosis and cell cycle control, cell adhesion and DNA repair mechanisms.⁴ It may also activate inflammatory pathways, which was suggested to be critical for tumor development.⁵

The integration of HPV virus in the host genome often occurs in the transcribed genomic region and this was suggested as a mechanism adapted by the virus to improve the expression of some viral products, notably E6 and E7 viral oncogenes.^6,7 In addition, the integration of HPV type 16 was strongly associated with progression of precancerous lesions, i.e., cervical intraepithelial neoplasia grade 1 to 3 (CIN1‒CIN3).⁸

Genomic instability and other DNA alterations including epigenetic ones, causing changed gene expression, are common features in many cancers.⁹ Regions differentially methylated in cancers often co-localize with tissue-specific DNA methylated regions (T-DMR).^10,11 In cervical cancer, DNA methylation can occur on the integrated viral DNA but aberrant DNA methylation is also induced in the host cell genome.¹² Increased methylation of HPV type 16 has been correlated with a more severe cancer progression.¹³ Aberrant methylation status of several host genes related to cell cycle, apoptosis, development, cell adhesion and cellular signaling have been analyzed in clinical studies by Szalmaz et al. and Wentzen et al.^14,15 These meta-studies also showed that there is a large variation in DNA methylation frequency, implying that new biomarkers are needed for early detection and progression of cervical cancer.

In this study, we employed the Illumina Infinium HumanMethylation450 BeadChip array to identify novel biomarker genes that are differentially methylated in cervical squamous cell carcinoma (SCC) DNA samples compared with DNA from normal cervical scrapes and CIN3 tissue. This approach enabled us to acquire a comprehensive epigenetic signature between the different cervical tissues in terms of global methylation, regional methylation in CpG islands (CGI), shores, gene related regions (TSS, 5′UTR, gene-body, 3′UTR), and also methylation of distinct CpG sites that discriminates between these tissues.

Results

DNA methylation profile in cervical cancer or CIN3 compared with normal cervical tissue

The probe call rate was > 99% for all of the samples and 454 215 CpG sites out of 485 577 were included in the analysis. The principal component analysis (PCA) analysis showed that all samples were close to each other. The majority of the cancer samples tended to hold together, followed by the cluster of CIN3 samples and the other cluster of the normal cervical samples (Fig. S1). The frequency distribution of the ∆β-values showed that cervical cancer tissue, compared with CIN3 and normal cervical tissue harbored hypermethylation at specific regions (data not shown). The frequency distribution of differentially methylated CpG sites across chromosomes showed that cancer tissue compared with the CIN3 and normal tissues had most differentially methylated CpG sites on chromosomes 1, 6, and 19 (Fig. S2).

Unsupervised hierarchical clustering analysis distinguished normal cervical tissues from tissues with CIN3 diagnosis and cancer tissues (Fig. 1). In this heatmap, 16 of the 25 statistically significant CpG sites belong to an annotated gene. The first ten CpG sites, located in different genes, are mostly hypomethylated in the cancer tissues as well as in the tissues with CIN3 diagnosis compared with the normal tissues (Fig. 1). In contrast, the remaining 15 CpG sites in specific genes are mostly hypermethylated both in cancer and tissues with CIN3 diagnosis compared with the normal tissues. The mean β-values for all the 25 CpG sites are shown in Table 1.

Figure 1. Unsupervised clustering analysis of the normal cervical tissue (N), CIN3 tissue (CIN3), and the cervical cancer tissue (CC). The color gradient green to red displays the β-value and can range from 0–1.

Table 1. Mean β-values for CpG sites represented in the heat map in Figure 1

		Normal tissue			CIN3 tissue			Cancer tissue
CpG ID	Gene name	Mean	SD	n	Mean	SD	n	Mean	SD	n
cg18803079	EFCAB7, DLEU2L	0.78	0.04	20	0.58	0.10	18	0.53	0.03	6
cg21581312	LOC723972	0.84	0.04	20	0.62	0.11	18	0.50	0.07	6
cg17461271	NA	0.79	0.05	20	0.60	0.11	18	0.54	0.05	6
cg12050434	NA	0.65	0.06	20	0.44	0.12	18	0.28	0.07	6
cg13070650	STK3	0.57	0.12	20	0.79	0.08	18	0.76	0.10	6
cg01262337	TKT	0.39	0.12	20	0.21	0.09	18	0.17	0.03	6
cg00470565	NA	0.62	0.11	20	0.85	0.19	18	0.65	0.09	6
cg02116612	MGLL	0.45	0.10	20	0.26	0.12	18	0.24	0.08	6
cg24662154	EXT2	0.06	0.06	20	0.16	0.08	18	0.13	0.06	6
cg14156405	RGS7	0.03	0.01	20	0.16	0.16	18	0.40	0.28	6
cg16241932	ZDHHC14	0.39	0.31	20	0.77	0.21	18	0.63	0.36	6
cg08272731	LHX8	0.05	0.04	20	0.22	0.18	18	0.47	0.22	6
cg03861379	NA	0.59	0.14	20	0.77	0.12	18	0.78	0.09	6
cg23243867	ST6GALNAC5	0.02	0.02	20	0.13	0.14	18	0.40	0.22	6
cg21632158	NA	0.06	0.05	20	0.19	0.14	18	0.40	0.17	6
cg00532449	NA	0.12	0.11	20	0.30	0.17	18	0.40	0.26	6
cg20784693	HDAC4	0.71	0.15	20	0.50	0.17	18	0.63	0.22	6
cg21183256	NA	0.07	0.04	20	0.22	0.17	18	0.53	0.17	6
cg17985646	TBX20	0.05	0.05	20	0.21	0.17	18	0.45	0.22	6
cg22617898	RAB3A	0.80	0.13	20	0.61	0.17	18	0.56	0.15	6
cg22762215	NA	0.06	0.05	20	0.21	0.17	18	0.48	0.22	6
cg20302133	KCNA3	0.01	0.01	20	0.18	0.23	18	0.42	0.31	6
cg22721334	ZSCAN18	0.01	0.02	20	0.14	0.16	18	0.39	0.21	6
cg23385248	NA	0.34	0.18	20	0.14	0.11	18	0.11	0.12	6
cg18075755	C1QTNF7	0.48	0.28	20	0.17	0.22	18	0.32	0.29	6

NA, not annotated. Bold-face numbers indicate genes selected as candidate biomarkers (see Table 4).

Differentially methylated CpG loci and gene regions

We analyzed the location of the hypo- and hyper-methylated CpG loci in relation to the CpG islands and gene context. Most of the differential methylation occurred in the context of a gene when comparing cancer tissue to precancerous and normal tissue (Fig. 2). A large fraction of the hypermethylated CpG loci were located in the CGI in contrast to the hypomethylated ones that were mostly located in the CGI shore or further away (Fig. 2).

Figure 2. Distribution of hypermethylated and hypomethylated CpG locus in CpG island (CGI), CpG island shore, CpG island shelves and other regions.

In Table S1 an overview of the number of differentially methylated regions and the range of the differences are shown comparing the three diagnostic groups. Most of the differentially methylated regions are found when comparing cancer tissue and normal cervical tissue. The majority of the events occur in the traditional gene promoter region (TSS1500, TSS200, 5′UTR, and 1st exon) and in the CpG islands. The comparison between cancer tissue and CIN3 tissues follow the same pattern. There are no differentially methylated regions reaching the criterion of ∆β of ≥ |0.2| when comparing the CIN3 tissues and normal cervical tissues.

Analysis of the differentially methylated regions in the vicinity of the traditional promoter region (TSS1500, TSS200, 5′UTR, and 1st exon) in cancer tissues compared with the normal tissues showed that many of the hypomethylated genes in cancer were found to be involved in the immune system response (Table 2). The hypermethylated genes are mostly receptors and transcription factors involved in cell development (Table 3).

Table 2. Hypomethylated gene promoter regions in cancer tissues compared with the normal tissues

Gene name	Function	Region	∆β	Adjust P value	Abberrantly methylated in
AIM2	Inflammasome activity	TSS1500 TSS200	-0.34 -0.33	0.004 0.002	Leukemia^31 pancreratic cancer^30
B3GALT6	Transmembrane protein	1st Exon	-0.32	0.001
BST2	Expressed in B-cells	1st Exon	-0.31	0.002	Lupus,^54 glioblastoma^55
BTN3A3	Inhibitor of immune response	5′UTR	-0.28	0.000
C1orf150	Not known	TSS200	-0.36	0.038
CASP10	Cell apoptosis	TSS200	-0.33	0.026
CLEC2B	Cell adhesion, cell-cell signaling, immune response	TSS1500	-0.31	0.033
CST7	Immune regulation	TSS1500	-0.31	0.028
CSTA	Epidermal development	TSS1500	-0.33	0.034
DAPP1	B-cell adhesion	TSS200 5′UTR	-0.32 -0.30	0.024 0.034
EDARADD	Ectodermal dysplasia	1st Exon	-0.28	0.035
FAM26F	Not known	5′UTR	-0.28	0.001
GBP4	Cell growth, anti-viral activity	TSS200	-0.37	0.007
GPR55	Cell proliferation	5′UTR	-0.29	0.042
HAVCR2	Regulates T-cells	5′UTR 1st Exon	-0.30 -0.30	0.030 0.028
HBBP1	hemoglobin, β pseudogene 1	TSS200	-0.35	0.010
HLA-DPB1	Immune response	1st Exon	-0.36	0.001
IGLL1	B-cell regulation	TSS200	-0.35	0.017
IL12RB1	Immune response	TSS200	-0.31	0.006	Atherosclerosis^56
MARCO	Innate antimicrobial immune system	TSS200	-0.30	0.032
MIR15B	Regulation of gene express	TSS1500	-0.35	0.007
NAGS	Urea cycle	1st Exon	-0.27	0.028
OR8B12	G protein-mediated transduction	1st Exon	-0.31	0.036
PATE3	Fertility and reproduction	5′UTR 1st Exon	-0.28 -0.28	0.049 0.046
PIGB	GPI-anchor biosynthesis	TSS1500	-0.32	0.005
PIK3R6	Leukocyte chemotaxis	5′UTR 1st Exon	-0.29 -0.28	0.030 0.016
PRTN3	Regulates immune cells	5′UTR 1st Exon	-0.29 -0.29	0.006 0.005	myeloid tissue^57
S1PR4	Cell-signaling in lymphoid tissue	TSS1500	-0.46	0.017
SDF4	Regulates calcium	TSS1500	-0.31	0.002
SIGLEC12	Immunoglobulin	5′UTR	-0.28	0.030
SPTA1	Cell shape	TSS1500	-0.35	0.002
STAT5A	JAK/STAT pathway	5′UTR	-0.30	0.003	lymphoma^58
VPREB1	B-cell differentiation	TSS200	-0.30	0.013

The list is merged of top-ten differentially methylated (according to Δβ-value) genes in the TSS1500-, TSS200-, 5′UTR- and 1st exon regions.

Table 3. Hypermethylated gene promoter regions in cancer tissues compared with the normal tissues

Gene name	Function	Region	∆β	Adjust p-value	Abberrantly methylated in
ACAN	Part of the extracellular matrix in cartilagenous tissue	TSS1500	0.41	< 0.001
ADCYAP1	Neurotransmitter	TSS1500	0.42	< 0.001	PTSD^59 cervical cancer^60
AJAP1	Adheren junctions	TSS1500	0.41	< 0.001	gastric cancer^61 oligidendrioma^62
AKAP2	Transcription factor	TSS1500	0.47	< 0.001
BARHL2	Homeobox-containing gene	TSS1500 5′UTR	0.46 0.57	< 0.001 < 0.001	Lung cancer^63
BOLL	Cell development	TSS200	0.52	< 0.001	Colorectal cancer,^64 lung cancer^65
C1orf114	ND	1st Exon	0.52	< 0.001
CTNND2	Encodes an adhesive junction associated protein	TSS200	0.56	< 0.001
CXCL12	Stromal cell-derived α chemokine member of the intercrine family	TSS200	0.55	< 0.001	Breast cancer,^66 rheumatoid arthr.^67 lung cancer,^68 glioma,^69 colorectal cancer,^70 pancreatic cancer^71
FAM159B	ND	5′UTR 1st Exon	0.55 0.55	< 0.001 < 0.001
FBXL7	F-box family member and function in phosphorylation-dependent ubiquitination	TSS200 5′UTR 1st Exon	0.52 0.56 0.56	< 0.001 < 0.001 < 0.001
GALR1	Inhibits adenylyl cyclase via a G- protein of the Gi/Go family	5′UTR	0.57	< 0.001	HNSCC^28
GYPC	Integral membrane glycoprotein regulates the mechanical stability of red cells	5′UTR 1st Exon	0.56 0.56	< 0.001 < 0.001
HYDIN	Involved in the ciliary function	1st Exon	0.54	< 0.001
KIF19	ND	TSS200	0.50	< 0.001
MARCH11	Ubiquitin-mediated protein sorting in the trans-Golgi network	TSS1500 TSS200	0.41 0.50	< 0.001 < 0.001
MIR129–2	Targets SOX4	TSS200	0.50	< 0.001	Gastric cancer^72
MIR663	Targets TGFB1	TSS1500	0.44	< 0.001
NAP1L3	NAP family are involved in neuronal differentiation	5′UTR	0.53	< 0.001
NPAS4	Transcriptional regulator	5′UTR	0.53	< 0.001
PTGDR	G-protein-coupled receptor	5′UTR	0.55	< 0.001	bladder cancer^73 neuroblastoma^74 asmathic pat.^75
SCT	Endocrine hormone	5′UTR 1st Exon TSS200	0.52 0.52 0.53	< 0.001 < 0.001 < 0.001
SORCS1	Membrane receptor	TSS200	0.61	< 0.001
SOX8	Transcription factor	TSS1500	0.41	< 0.001
TCERG1L	Transcription factor	1st Exon	0.53	< 0.001	colrectal cancer^76
TRIM58		1st Exon	0.52	< 0.001	Age-associated^77 hepatocellular carcinoma^78
TTYH1		1st Exon	0.57	< 0.001
VSTM2B		TSS1500	0.48	< 0.001
ZFP42		TSS1500	0.42	< 0.001	Renal cancer^79 esophageal cancer^80
ZIK1		5′UTR 1st Exon	0.61 0.61	< 0.001 < 0.001	Esophageal cancer^80
ZNF582		TSS201	0.58	< 0.001	Cervical cancer^27

The list is merged of top-ten differentially methylated (according to Δβ-value) genes in the TSS1500, TSS200, 5′UTR and 1st exon regions. ND, not determined; PTSD, post-traumatic stress disorder; HNSCC, head and neck squamous cell carcinoma.

Identification of functionally related gene groups

A total number of 74 gene ontology term clusters were generated when analyzing the hypomethylated CpG sites in cancer compared with normal cervical tissues. The GO_BP cluster with the highest enrichment score of 6.42 was related to the immune system (Table S2). The GO_FAT terms included several KEGG pathways involving the immune system, such as antigen processing and presentation, autoimmune thyroid disease and asthma.

The gene ontology analysis of the hypermethylated CpG sites created three clusters, and GO_BP had the highest enrichment score of 36.52 (Table S3). These gene ontology terms belong to the categories of cell development and cell differentiation. The GO_FAT and KEGG classification included categories of transcriptional regulation and DNA binding with the highest enrichment score of 24.41.

Candidate biomarker genes

The selection of candidate biomarker genes was based on both differentially methylated CpG sites generated with unsupervised hierarchical clustering and differentially methylated gene regions.

From the list of genes generated with unsupervised hierarchical clustering (Table 1) the selection criteria were 1) Δβ-value > 0.2 between the cancer and normal tissues, and 2) a mean baseline (normal tissue) β-value < 0.2 of a CpG site hypermethylated in cancer or a baseline (normal tissue) β-value > 0.4 for a CpG site hypomethylated in cancer. This approach led to the selection of six genes, the best candidate biomarkers: RGS7, LHX8, STGALNAC5, TBX20, KCNA3, and ZSCAN18 (Table 1 and 4), all with ≥ 9 significant CpG sites.

Table 4. Methylation status and gene expression of candidate biomarker genes

Gene name	Function	Methylation status, CpGs on the array^1	Gene expression data set, FC GSE6791 GSE9750 GSE7803			Aberrantly methylated in cancers
ACAN	Part of the ECM	Hyper, (10/18)	-1.34*	-0.47	-0.38	−
AJAP1	Membrane protein	Hyper, (30/69)	-1.39*	-1.24	-1.31*	gastric cancer^61 oligidendrioma^62
BARHL2	TF	Hyper, (15/15)	−	−	−	lung cancer^63
BOLL	TF	Hyper, (8/17)	-1.55*	−	−	CRC^64 lung cancer^65
C1orf114	NA	Hyper, (9/13)	-1.38*	-1.08	-1.18*	−
FBXL7	ubiquitin	Hyper, (26/28)	-1.26	-1.89*	-1.66*	−
GALR1	galanin receptor	Hyper, (29/34)	-1.16*	-1.42	1.18	HNSCC^28
GYPC	erythrocyte surface protein	Hyper, (4/18)	-1.69*	-1.04	1.17	−
KCNA3	ion channel	Hyper, (9/19)	-1.48*	-1.01	1.46*	−
KIF19	microtubule-depolymerizing activity	Hyper, (12/22)	-1.24*	−	−	−
LHX8	TF	Hyper, (43/46)	-1.09	−	−	−
MIR663	micro RNA	Hyper, (11/11)	−	−	−	−
PTGDR	prostaglandin D(2) receptor	Hyper, (13/19)	-1.25*	-0.16	0.02	bladder cancer^73 neuroblastoma^74
RGS7	regulator of G protein signaling	Hyper, (18/35)	-1.38*	-0.30	-1.06	−
S1PR4	1-phosphate receptor	Hypo, (5/12)	-1.40*	1.26	-1.10	−
SORCS1	VPS10 domain containing receptor	Hyper, (10/20)	-1.42*	−	−	−
ST6GALNAC5	modify proteins on the cell surface	Hyper, (18/29)	1.24	1.08	-1.22	CRC^81
TBX20	TF	Hyper, (19/20)	-1.17*	−	−	−
TRIM58	NA	Hyper, (14/16)	-1.13*	1.16	-1.18	−
TTYH1	chloride anion channels	Hyper, (12/16)	-1.32*	1.05	-1.13	−
VSTM2B	NA	Hyper, (17/20)	−	−	−	−
ZIK1	TF	Hyper, (12/15)	-1.79*	−	−	esophageal cancer^80
ZNF582	TF	Hyper, (10/12)	-1.41*	−	−	cervical cancer^27
ZSCAN18	TF	Hyper, (24/38)	-2.20*	-3.68*	-3.34*	renal carcinoma^82

The gene expression was analyzed in three independent data sets. ¹Number of statistically significant CpG sites / total number CpG sites on the array. The statistically significant (P < 0.05) gene expressions are marked with (*). The negative value means lower expression in cancer compared with the normal tissue. The (−) sign indicates no data available. TF, transcription factor; CRC, colorectal cancer; HNSCC, head and neck squamous cell carcinoma.

From Table 2 and 3, displaying differentially methylated gene regions, we selected genes that were differentially methylated with a Δβ-value > |0.4|. The gene regions were then sorted according to the adjusted p-value and the top 20 regions were selected as candidates corresponding to 18 genes: ACAN, AJAP1, BARHL2, BOLL, C1orf114, FBXL7, GALR1, GYPC, KIF19, MIR663, PTGDR, S1PR4, SORCS1, TRIM58, TTYH1, VSTM2B, ZIK1, and ZNF582. The enlarged list of candidate biomarker genes is presented in Table 4. With exception of S1PR4, the candidate genes were hypermethylated in cancer compared with CIN3 and normal cervical tissues. As illustrative examples, the mean methylated fractions of all interrogated CpG sites in three selected genes, S1PR4, TBX20, and ZNF582 are shown in Figure S3. There is a stepwise increment in methylation (β-value) when comparing normal tissue, CIN3 tissues, and cervical cancer. We further validated the differential methylation of four selected genes, TBX20, RGS7, KCNA3, and S1PR4, with the Pyrosequencing assay technology. Figure S4 shows the percentage methylation of the analyzed CpG sites in cervical cancer and normal tissues. The results are in very good agreement with our array data. The mean methylation of the cervical cancer tissues was higher compared with the normal cervical tissues in the TBX20, RGS7, and KCNA3 genes, while the S1PR4 gene was hypomethylated in cervical cancer tissues. We found a higher variation in the methylated fraction within the cervical cancer group compared with the normal cervical tissues, as detected also by the array technique (Table 1).

External validation

External validation of our candidate biomarker genes was done using the mRNA expression data sets from Peyon et al.,¹⁶ Scotto et al.,¹⁷ and Zhai et al.,¹⁸ which fulfilled quality criteria of tissue samples and employed high quality whole-genome expression arrays. Peyon et al.¹⁶ reported probe sets for 21 of our 24 genes of interest, whereas the two other, reported sets for 14 of the genes. There was a good agreement between our DNA methylation data and the gene expression array: fold changes ranged between −1.09 and −2.20 for all hypermethylated genes, see Table 4. The hypomethylated S1PR4 gene was, as expected, upregulated in cervical cancer in one gene expression data set, but it was downregulated in the two other.

Discussion

The aim of this study was to acquire a comprehensive epigenetic signature of cervical cancer, and to identify new potential biomarker genes that are differentially methylated in cervical SCC compared with CIN3 or the normal cervical tissue. Our results showed that different cellular pathways are hypo- and hyper-methylated in cancer tissue. A short list of 24 novel candidate biomarker genes emerged, which was consistent with the external validation sets and could ultimately be clinically applicable.

Fifteen of the 24 potential biomarker genes (ACAN, Clorf114, FBXL7, GYPC, KCNA3, KIF19, LHX8, MIR663, RGS7, S1PR4, SORCS1, TBX20, TRIM58, TTYH1, and VSTM2B) have not yet been correlated with any cancer type, while eight (AJAP1, BARHL2, BOLL, GALR1, PTGDR, ST6GALNAC5, ZIK1, and ZSCAN18) have been implicated in cancers other than cervical cancer (Table 4). Nevertheless, all of these genes are included in crucial cellular functions, such as cell signaling, gene transcription, immunity, glycosylation, membrane activity and ubiquitination. Therefore, the hypermethylation of genes involved in such processes may contribute to the accumulation of damages that can progress to cervical cancer. By external validation using three independent data sets, the gene expression of 23 potential biomarkers was downregulated. This correlated with their hypermethylated state, confirming a link between DNA methylation and mRNA expression.

In cervical carcinogenesis, both global DNA hypomethylation and hypermethylation of specific genes promoters has been suggested to occur.^19,20 We found here that most of the analyzed CpG loci were equally methylated (94%) in cervical cancer cells and normal cervical cells (Table S4). However, the genome-wide DNA methylation array confirmed that hypermethylation in cervical cancer occurs in fact at specific genes (24 199 CpG sites; Table S4). Even though we included more CIN3 cases precisely to increase the power to detect differences between CIN3 and normal tissues, the differences between the methylome of cervical cancer and normal tissues was greatest in all analyses (p- and β-values) compared with the methylome of CIN3 vs. normal tissues; therefore, our report focuses on these differences.

Distinctly different cellular pathways appear to be hypo- or hyper-methylated. In cervical cancer, the hypomethylated genes were related to the immune system response (Table S2), whereas the hypermethylated genes were involved in the biological processes of cell development, cell differentiation, and transcriptional regulation (Table S3). Shen et al.²¹ interrogated a smaller number of CpG sites (Infinium Human Methylation27 BeadChip) in colon, kidney, stomach, lung and breast cancers as well as the corresponding normal tissue; they reported results consistent with our methylation data performed on cervical cancer.²¹

Aberrant DNA methylation signatures in cervical cancer tissue, compared with the normal cervical tissue, are located in different regions depending on the methylation status. Hypermethylated CpG sites are predominantly located in the CGIs and hypomethylated CpG sites in the CGI shores (Fig. 2). In the gene context, most of the changes occurred in the TSS200 region (Table S1). Aberrant DNA methylation has been mostly studied in promoters with a CpG island. A few years back, CpG island shores were identified as potentially important regions harboring DNA methylation changes in colorectal cancer.¹¹ CpG island shore methylation was statistically significantly correlated with clinical characteristics in breast cancer when we analyzed the SLC25A43 gene.²² Feber et al.²³ analyzed pooled DNA from benign and malignant nerve sheath tumors with methylated DNA immunoprecipitation and next generation sequencing (MeDIP-seq), and suggested that both hypo- and hyper-methylated differentially methylated regions were located in the CGI shore regions.²³ The MeDIP-seq methodology yields a semi-quantitative methylation percentage and differs from the Illumina Infinium HumanMethylation450 BeadChip that detects selected CpG sites at a single-nucleotide resolution in virtually all genes.^24,25

Furthermore, methylation changes that are tumor driving were suggested to be present in premalignant lesions and to occur early in the process of carcinogenesis.²⁶ From a clinical perspective, methylation biomarkers for cervical cancer detection and prognosis are highly needed.¹⁵ Lendvai et al.¹⁹ recently used MeDIP-seq to perform a genome-wide analysis of high-grade CIN3 lesions compared with normal cervical epithelium, and identified 80 differentially methylated regions. They suggested COL25A1 and KATNAL2 genes as potential biomarkers for early detection of high-grade cervical intraepithelial neoplasia.¹⁹ Our study also identified these two genes as differentially methylated in cervical cancer tissue compared with CIN3 or normal cervical tissue, but not when comparing CIN3 tissues and the normal cervical tissues (data not shown), indicating that these two genes might be late biomarkers, however they did not qualify as biomarkers according to our selection criteria. Our data extends the study of Lendivai et al.¹⁹ by comparing the methylome of cervical cancer tissues compared with the CIN3 lesions and the normal tissues. By this approach we generated a list of 24 novel potential biomarker genes (Table 4).

One of our candidate biomarker gene, ZNF582, was previously shown to be hypermethylated in cervical cancer suggesting silencing of the ZNF582 gene in cervical neoplasms (Fig. S3).²⁷ Although the biological function of ZNF582 is not yet well characterized, it could be a tumor suppressor gene as well as potential candidate gene for molecular cervical cancer screening.²⁷ Another highly interesting gene is GALR1, which inhibits a key regulatory enzyme adenylyl cyclase. This gene was suggested to be hypermethylated in a subset of head and neck cancer caused by oncogenic HPV types, notably HPV 16.²⁸ Thus, frequent promoter hypermethylation of GALR1 and cellular growth suppression after re-expression of the gene supports the hypothesis that GALR1 could be a tumor suppressor gene in head and neck cancers,²⁸ as well as in cervical cancer (our study).

Our findings of hypomethylated genes in the immune system provide strong evidence that the immune response is under control of hypermethylation in normal cells, while in cancer cells, these genes become active by losing the methylation mark in the gene promoter region, which represses gene expression. One example is the hypomethylated AIM2 gene, whose functional protein forms an inflammasome, resulting in caspase activation in inflammatory cells. This induction of inflammatory cell response could occur in cervical cancer cells due to the hypomethylation of genes such as the AIM2 gene. Furthermore, AIM2 protein has been found overexpressed in a wide variety of tumor types, including neuroectodermal tumors,²⁹ pancreatic cancer,³⁰ and leukemias.³¹

Another interesting hypomethylated gene is HLA-DPB1, whose gene product belongs to the HLA class II molecules, which plays a central role in the immune system by presenting peptides derived from extracellular proteins. There is evidence that HLA-DP participates in the functional T-cell responses against HPV in cervical cancer.³² The host immune response to HPV is generally at low-level because HPV is shielded in the basal epithelial cells from the circulating immune cells during initial stages of the infection. Furthermore, HPV employs several mechanisms to downregulate innate and cell-mediated immunity.³³ This is confined to the early stages of infection since it is known that functionally active immune responses are generated at later stages of the HPV infection.³⁴ Thus, our findings in cervical cancer are consistent with an immune response against HPV, which is activated through hypomethylation of specific genes included in the immune response and inflammatory reactions.

The immune system, upon activation by hypomethylation, can prevent tumor development by suppressing the viral infection, prevent an inflammatory environment leading to tumorigenesis, and identify tumor cells based on their antigen expression and eliminate them.³⁵ Inflammatory cells are suggested to affect the tumor microenvironment, which in turn can promote or repress tumor growth.³⁶ The theoretical model of aberrant DNA methylation induced by chronic inflammation involves cytokine production TNFα, IL-1β, and oxidative stress in epithelial cells.^36,37 Thus, without adequate immune responses, transformation to cancer might be quicker. This is the case of the adenovirus 12, which has the capability to turn off the immune system genes, which may lead to tumor progression.³⁸ Tumor tissue may contain increased amounts of leukocytes of various types, which could partly account for the observed hypomethylation of immune-system-related genes. The recent introduction into clinical practice of liquid-based cervical cytology may offer opportunities to address this issue in future studies.

In order to determine if DNA methylation changes occur early or late in cervical carcinogenesis, genome-wide DNA methylation studies should extend to lower CIN grades, CIN1 and CIN2 as well as atypical squamous cells of undetermined significance (ASCUS). However, our study clearly indicates that DNA methylation profile of specific genes in samples with high-grade cervical lesions (CIN3) is located between samples with normal diagnosis and samples with invasive cervical cancer (Fig. S1), but closer to cluster of normal samples, as exemplified by specific CpG sites in the ZNF582, TBX20, and S1PR4 genes (Fig. S3).

We have validated our array results with an independent method analyzing the methylated fraction of four genes (Fig. S4). Our results suggest as previously reported by Cachill et al.³⁹ and Deneberg et al.,⁴⁰ that the Illumina HumanMethylation450 array measures differential DNA methylation comparable to the golden standard pyrosequencing assay technology.

Conclusion

In this study, we selected candidate genes with CpG sites highly methylated in cervical cancer tissue and without methylation in normal cervical tissue. This combination provides an optimal condition for a successful diagnostic test. The results highlight the extensive differential methylation signature in HPV-induced cervical cancer tissue compared with CIN3 or normal cervical tissues. The 24 candidate biomarkers that we identified represent several types of mechanisms in the cellular machinery that are epigenetically deregulated, such as membrane receptor function, intracellular signaling, and gene transcription. External validation confirmed a link between DNA methylation and mRNA levels for these genes. Considering the extensive hypomethylation of genes involved in the immune system, epigenetic therapy and immune therapy might need to be combined when treating cervical cancer. Understanding of the functional role of DNA methylation alterations in cancer genomes may prove to be clinically applicable in disease diagnostics and prognostics, and may guide the development of new epigenetic therapies.

Materials and Methods

Study group

The cervical specimens were collected in the Sisters of Mercy Hospital, Zagreb, Croatia during the period from 2004 to 2011. The clinical characteristics of the specimens are presented in Table S5. The cytological diagnoses were obtained for normal cervical smears, high-grade squamous intraepithelial lesion (HSIL)/CIN3, while cervical squamous cell carcinoma (SCC) samples were defined by histopathology.⁴¹ Subsequently, cervical specimens were tested and typed for HPV.⁴²

The cervical sample collection is regulated through laboratory service request forms that have to be signed and approved by the practicing physician. The extracted DNA from cervical specimens was processed without initial knowledge of patient data (age, diagnosis, HPV detection, and typing result). The study was approved by the Ethical Board of the Rudjer Boskovic Institute, as well as the Ethical Board of Sisters of Mercy Hospital, and conforms to the Helsinki declaration (DoH/Oct2008). The Ethical Review Board, Uppsala, approved the Swedish participation in this study.

DNA preparation

DNA from cervical cell samples was isolated on the BioRobot EZ1 (Qiagen, Hilden, Germany) according to the manufacturer’s instruction. After DNA extraction, the purified DNA was dissolved in 50–100 μl of tri-distillate sterile water and stored at –20 °C until further analysis. Each DNA was analyzed by electrophoresis on 1% agarose gels and spectrophotometrically.⁴³

HPV detection and typing

Three sets of consensus primers for HPV detection were used: PGMY09/PGMY11,⁴⁴ L1C1/L1C2-1/L1C2-2⁴⁵ and GP5+/GP6+.⁴⁶ The quality of the isolated DNA was tested by amplification of the 268 bp sequence of the β-globin gene using PC04/GH20 primers⁴⁷ in a multiplex PCR with PGMY primers. Type-specific primers for HPV types 6/11, 16, 18, 31, 33, 45, 52, and 58 were used for HPV typing according to Milutin-Gasperov et al., 2007.⁴² Aliquots of each PCR product (10 μl) were analyzed by electrophoresis on 2% agarose gels stained with ethidium bromide. The amplified products were visualized by UV irradiation of the gels either by ImageMaster VDS (Pharmacia Biotech) or UVItec Cambridge (Alliance 4.7).

Methylation array

The Illumina Infinium HumanMethylation450 BeadChip (Illumina, Sweden) includes 485,577 CpG sites located all over the genome.²⁴ The quality and integrity of the samples were evaluated with Nano Drop and gel electrophoresis; samples were included if the ratio A_260/280 was between 1.7–1.9, the ratio A _260/230 was > 1 and the samples resolved on gel electrophoresis showed a distinct DNA band > 10 000 bp. In this study, 20 normal cervical samples (HPV negative), 18 samples with CIN3 lesions (HPV positive) and 6 cervical cancer tissues (HPV positive) were included. The DNA methylation assay was performed as follows. Briefly, approximately 500 ng of DNA from cancer, CIN3 lesions, and normal cervical tissues were bisulfite treated with Zymo EZ DNA Methylation Gold kit (Zymo Research) according to manufacturer’s procedure. After the bisulfite treatment, 200 ng DNA was subjected to the whole genome amplification (WGA) and enzymatic digestion with reagents provided within the Infinium HumanMethylation450 kit (Illumina). The hybridization of the samples on the BeadChips and washing procedures followed the standard procedures obtained from Illumina. The iScan scanner (Illumina) was used to read the BeadChips and the data collection was performed in the GenomeStudio software (version 1.0).

Data processing

The Infinium Human 450K Methylation Bead chip has several built in sample dependent and independent controls that measure the processing of the chip. They were analyzed in the Genome Studio software with the Dash-Board module, prior to further statistical analyzes. The probe call rate for a passed sample should be > 98%, and the detection p-value for the probes should be < 0.01, otherwise they were excluded. The β-value generated for each CpG locus measures the intensity of methylated (β = 1) and unmethylated probes (β = 0). It is calculated as β = [M/(M + U + 100)], where M = methylated allele and U = unmethylated allele. The study cohort contains only women; therefore the Y-chromosome was excluded from the analysis. The raw β-values were generated with the GenomeBead studio using the settings normalization controls and background subtraction. This data has been deposited to the Gene Expression Omnibus database under the accession number GSE46306.

A CpG locus was considered differentially methylated if the Δβ-value was ≥ ∣0.2∣ and the adjusted P value < 0.05. This cut-off value of ∣0.2∣ represents the 99% confidence interval of the detection limit.²⁵

Unsupervised clustering of CpG loci and DAVID analysis

Heat map showing clustering of the patient samples was performed based on P value < 0.05 and ∆β-value ± 0.5 between CIN3 tissues and normal tissues. They were produced with the Illumina Methylation Analyzer (IMA) package version 2.1.1 (Illumina). The Database for Annotation, Visualization and Integrated Discovery (DAVID)⁴⁸ was used to analyze biological features associated with genes hypo- or hypermethylated in cervical cancer compared with normal cervical tissue. The gene lists were obtained with the criteria P value < 0.05 and a more stringent Δβ-value of > |0.3| (between cancer and normal tissues) due to the input limit of genes to the DAVID tool. The gene ontology (GO) groups were clustered into functional annotation referring to biological process (GO_BP), cellular component (GO_CC) and molecular function (GO_MF). GO_FAT terms were also retrieved to filter the broadest GO categories. In addition, the Kyoto Encyclopedia of Genes and Genomes (KEGG) database was used to identify functional gene groups and ontology terms that are significantly overrepresented.

Differentially methylated CpG loci and gene regions

The CpG loci were categorized as either hypermethylated (positive Δβ-value) or hypomethylated (negative Δβ-value). The CpG loci were then categorized in the context of (I) CpG islands, CpG shores, CpG shelves or others and (II) genomic location of TSS1500 (1.5 kb within transcriptional start site), TSS200 (200 bp within transcriptional start site), 5′UTR (5′ untranslated region) and 1st exon, gene body and 3′UTR (3′ untranslated region) or others. To compare the regions (CpG island, CpG shore, CpG shelves, TSS1500, TSS200, 5′UTR, 1st exon, gene body, and 3′UTR) between the different diagnostic groups, the median β-value was calculated including all the CpG sites in the analyzed region. Each region contains on average 1.53–9.92 probes.²⁵

External validation of the candidate biomarker genes

Three data sets GSE6791,¹⁶ GSE9750,¹⁷ and GSE7803¹⁸ deposited to the Gene Expression Omnibus (GEO) database were used for external validation of mRNA expression of candidate biomarker genes. The data sets measured gene expression by the Affymetrix GeneChip Human genome U133 Plus 2.0 (GSE6791) array and Affymetrix HG-U133A array (GSE9750 and GSE7803). The GSE6791data set was composed of 20 cervical cancer samples (17 HPV+, 3 HPV-) and 8 normal cervical samples (HPV-); the GSE9750 data set of 27 cervical cancers (HPV+) and 24 normal cervical tissues (HPV-); and GSE7803 data set of 24 cervical cancer tissues (HPV+) and 10 normal cervical tissues (HPV-). Using the online tool GEO2R (www.ncbi.nlm.nih.-gov/geo/geo2r), the log fold change (FC) was calculated to estimate the differences in the gene expression in cervical cancer compared with the normal cervical samples. The unadjusted P values < 0.05 were considered as statistically significant, since these analyses were hypothesis driven.

Pyrosequencing validation of differential methylation

We developed assays for the following genes: TBX20, RGS7, KCNA3, and S1PR4. The PCR and sequencing primers were designed to assess the exactly same CpG sites analyzed by the array. All primers were purchased from www.biomers.net (Ulm/Donau). In the Table S6, the sequences, amplicon sizes, and the optimal annealing temperatures are indicated.

The analysis was performed on 5 normal and 5 cancer cervical tissues, which were already tested by Illumina Infinium HumanMethylation450 BeadChip array. Briefly, approximately 500 ng extracted DNA was used for the bisulfite treatment performed with the EZ DNA Methylation Gold kit according to the instructions by the manufacturer and eluted in 20 µl elution buffer (Zymo Research, Orion Diagnostica). The PCR reactions were performed in a total volume of 40 µl with the KAPPA2G Robust HotStarTaq DNA Polymerase Kit (Kapa Biosystems, Inc.), containing 0.10 µmol/L of each primer, 0.5 units of Taq polymerase, 1.5 mM MgCl₂, and 0.1 mM each of dGTP, dATP, dTTP, dCTP and 50 ng of bisulphite treated DNA was added as template. The PCR reaction for the RGS7 gene contained also the 1X Enhancer solution provided in the kit. The PCR program was as follows: initial denaturation step of 1 min at 95 °C, followed by 45 cycles of 30 s denaturation at 95 °C, specific annealing T_a for 30 s and extension for 30 s at 72 °C and one cycle of final extension for 30 s at 72 °C. Sequencing was performed using a Pyromark Gold Q96 Reagent Kit and a PSQ 96ID system (Qiagen) as described previously by Farkas et al.⁴⁹ The nucleotide addition order was optimized by the Pyro Q-CpG software version 1.0.9 (Qiagen) and the results were automatically analyzed using the same software.

Statistical analysis

All statistical analyses were performed using the R software with the IMA package.⁵⁰ The raw β-values were arcsine square root transformed and the empirical Bayes moderated t-statistic was used to generate the P values.⁵¹ The Benjamini-Hochberg method was used to adjust the P values for multiple testing.⁵² For site-level analysis the mean Δβ-values were used and for the regional analysis median Δβ-values value was calculated including all CpG sites in one region. The principal component analysis (PCA) was performed to visualize the difference in DNA methylation between the cancer and normal tissue samples. In the DAVID analysis the statistic used was the modified Fishers exact P value, adjusted for multiple testing with the Bonferroni method.⁵³

Abbreviations:
CGI	=	CpG island
CIN3	=	cervical intraepithelial neoplasia grade
HPV	=	human papilloma virus
PCA	=	principal component analysis
SCC	=	squamous cell carcinoma
T-DMR	=	tissue-specific differentially methylated regions

Acknowledgments

Methylation profiling was performed by the SNP&SEQ Technology Platform in Uppsala. The platform is part of Science for Life Laboratory at Uppsala University and supported as a national infrastructure by the Swedish Research Council. This work was supported by Lions Cancer Foundation, Nyckelfonden, Örebro läns landsting (Farkas SA and Nilsson TK), and partially by the Croatian Ministry of Science, Education and Sports (Grant number: 098-0982464-2510; Grce, M and Milutin-Gašperov, N).

10.4161/epi.26346

Disclosure of Potential Conflicts of Interest

No potential conflicts of interest were disclosed.

Supplemental Materials

Supplemental materials may be found here: www.landesbioscience.com/journals/epigenetics/article/26346