Characterization of DNA hydroxymethylation profile in cervical cancer

Abstract

Cervical cancer is one of the most common gynecological tumors in females. DNA methylation alteration is a type of epigenetic decoration that controls the gene transcriptional regulation and is essential for the pathological progression of cervical cancer to reflect the prognosis and therapeutic sensitivity in clinical practice. Beyond DNA methylation, DNA hydroxymethylation considered as a more stable biomarker draws the outline of the reversible cycle from DNA methylation and demethylation. However, the landscape of 5-hydroxymethylcytosine (5hmC) distributed in the genome is never characterized in cervical cancer. In this study, we presented the whole 5-methylcytosine (5mC) and 5hmC profile in cervical cancer of I-IIa and IIb to IV stages and cervicitis tissues as control by dot plot assay and immunohistochemistry. We observed that the total 5mC was up-regulated while 5hmC was down-regulated in cervical cancer group compared to the control group. Furthermore, we investigated the distribution of 5mC and 5hmC on genomic DNA by MeDIP- and hMeDIP-Seq. 53 differential methylation/hydromethylation regions (DMRs/DHMRs) displayed a continuously increasing or decreasing trend of 5mC or 5hmC from cervicitis to I–IIa and from I–IIa to IIb–IV stages of cervical cancer. Thirty-seven DMRs and DHMRs have a similar variation trend while the other 8 have the opposite trend compared between CSCC and cervicitis. Moreover, the DMR/DHMR associated genes were closely related to Wnt, MAPK, Rap1 and other important signaling pathways. Finally, 5hmC beyond 5mC at the genes such as ACTG1, SALL3, DNAJA3, SERPINB6, CDC14B and CALN1 were considered as the putative novel hallmarks for cervical cancer diagnosis and prognosis. Altogether, this study first describes the DNA hydroxymethylation atlas of cervical cancer and shows a list of novel genes transcriptionally regulated by DNA methylation and hydroxymethylation.

Keywords:

• Cervical cancer
• hydroxymethylation
• 5hmC
• methylation
• 5mC

Introduction

Cervical cancer is the second most common gynecological tumors among women worldwide with an estimated 520,000 newly diagnosed cases and 260,000 deaths globally every year, which clearly affected youth trends [1]. Almost all of the cervical cancer cases were HPV-associated, including 70% squamous cell carcinomas and 25% adenocarcinomas [2]. So far, vaccines still face the problems of function limitation and unachieved full-type coverage upon HPV infections. Moreover, only about 15% 5-year survival rate among advanced patients indicates the poor prognosis in the late stages of cancers [3]. Therefore, it cries out for elucidating the underlying molecular mechanisms on to understand the development of cervical cancer.

DNA methylation, the most widely studied epigenetic regulation mainly occurs at CpG loci or islands via transferring a methyl group to cytosines to generate 5-methylcytosine (5mC) by methyltransferase DNMT family. The presence of aberrant DNA methylation including hypo- and hypermethylation in cervical cancer and high-grade cervical intraepithelial neoplasia suggests that most genes are hypermethylated [4]. However, direct removal of the methyl from 5mC is infeasible in terms of the energy consumption for DNA demethylation. Recent studies revealed that 5mC could be further oxidized into 5-hydroxymethylcytosine (5hmC), 5-formylcytosine (5foC), and 5-carboxylcytosine (5caC) by ten-eleven translocation (TET) proteins and finally return to regular cytosine for one manner of DNA demethylation process [5–7]. Thus, 5hmC, a more stable epigenetic mark compared to 5mC, plays a crucial role in methylation recycle maintenance as an intermediate [8].

To date, many studies have attempted to bridge 5mC to the diagnostic and prognostic biomarker of cervical cancer [9–12]. Nevertheless, the 5hmC pattern in cervical cancer is less investigated. Only both studies reported that the overall 5hmC level is decreased in cervical squamous cell carcinoma tissues [13,14]. However, current research progress only reflects the global level of 5hmC by immunostaining or high-performance liquid chromatography, but never figures out the 5hmC pattern on the genetic loci from the genome-wide level. Herein, we first describe the DNA hydroxymethylation atlas in different stages of human cervical cancer tissues using genomic high-throughput sequencing. Our data determines that 5hmC in some genes displays a more sensitive and specific characteristic for evaluating the malignancy of cervical cancer compared to the levels of transcription or 5mC.

Materials and methods

Patients and tissues

The tissue specimens from three cervicitis and twelve cervical squamous cell carcinoma (CSCC) patients were selected from the medical records of the Department of Gynecology of Yuhuangding Hospital from January to July 2017. Six tumor specimens were surgically excised from the patients with I–IIa stage, while four cervicitis and the other six tumor specimens with IIb–IV were harvested by cervical biopsy. Tumor stage was evaluated according to the FIGO staging standards [15]. The age of the patients was 38.4 ± 4.82 years for cervicitis and 42.6 ± 7.15 years for CSCC. None of the patients in this study received chemotherapy or radiotherapy before surgery. All the patients understood and signed the informed consent. Signed informed consent and ethics committee documents of Ethics Committee of The Affiliated Yantai Yuhuangding Hospital of Qingdao University were all provided to approve this study.

Dot blot assay

Dot blot assay was performed as previously described [16]. In brief, the patients’ cervicitis or CSCC tissues were extracted by genomic DNA using QIAquick Gel Extraction Kit (Qiagen, Dusseldorf, Germany). After quantifying and diluting the equal DNA concentration, each sample (2 μl) was dropped on nitrocellulose (NC) membranes with double proportion dilution (0, 5, 10, 20 and 40 ng) for dot blot assay. NC membranes were naturally dried at room temperature (RT) and incubated with 5hmC (1:500 dilution, Abcam, ab106918, Cambridge, UK) or 5mC antibodies (1:500 dilution, Abcam, ab10805) in TBST 4 °C overnight, followed by washing by TBST three times for 10 min in room temperature. followed by appropriate secondary antibodies (1:5000 dilution, Beyotime, Shanghai, China) in 10 ml of TBST for 1 h at room temperature with gentle shaking and washed again for three times and washed again for three times. Finally, the spots on NC membranes were developed with 3 ml of DAB substrate/chromogen for 5 min in darkness and analyzed the gray values.

Immunohistochemical (IHC) analysis

Tissues were embedded in paraffin and cut into sections with 5 μm. After deparaffinization, antigen retrieval was conducted within 10 mmol/l citrate buffer (pH 6.0). After blocking with 5% horse serum albumin (Beyotime, China) at 37 °C for 30 min, the sections were washed by PBS and incubated with a rabbit anti-human 5-hmC monoclonal antibody (1:2000 dilution, Abcam) or a mouse anti-human 5-mC monoclonal antibody (1:2000 dilution, Abcam) at 4 °C overnight. After washing, samples were incubated with secondary anti-rabbit or anti-mouse antibodies (1:5000 dilution, Beyotime) in room temperature for 30 min. Sections were washed and incubated with freshly prepared DAB substrate/chromogen (Beyotime) mixture for 5 min and washed again. The sections were counterstained, dehydrated and mount prepared. The positive staining of each specimen was the mean value of the scores. Specimens were scored according to the intensity of the dye color and the number of positive cells. The intensity of the dye color was graded as 0 (no color), 1 (light yellow), 2 (light brown) or 3 (brown), as well as the number of positive cells, was graded as 0 (<5%), 1 (5–25%), 2 (26–50%), 3 (51–75%) or 4 (>75%). Both of these two indices were assessed together to assign the specimens as four categories: 0–1 (-), 2 (+), 3–4 (++) or 5 (+++). 0–2 and 3–5 were considered as low and high expression respectively.

Methylated and hydroxymethylated DNA immunoprecipitation (MeDIP/hMeDIP) assay

Briefly, genomic DNA from the tissues of two cervicitis and four CSCC (two I–IIa and two IIb–IV) were extracted and sonicated to generate DNA fragments less than 500 bp. Illumina barcode adapters (AmpliSeq Library PLUS, Illumina, San Diego, CA, USA) were ligated before MeDIP and hMeDIP. Approximately, 4 mg of adaptor-ligated gDNA from the CSCC and control samples were pooled together in one tube. DNA libraries were denatured and diluted (5% used as input). The denatured DNA was gently shaken with the antibodies of 1 μg 5mC or 5hmC antibodies at 4 °C overnight and captured by protein A/G beads (Thermo, Waltham, MA, USA). The hMeDIP and MeDIP products were purified and sequenced followed by standard Illumina protocols.

DNA libraries were sequenced on an Illumina Nextseq500 platform. The raw sequence reads of input, MeDIP and hMeDIP were trimmed adaptors and filter out low-quality reads using Cutadapt (v1.9.1) and Trimmomatic (v0.35) [17], and checked the quality of clean reads using Fastqc [18]. Next, clean reads were mapped to the human genome (assembly hg38) using the Bowtie 2 (v2.2.6) algorithm [19]. The process of 5hmC and 5mC peak calling (p<.01) were performed by MACS 2 (v2.1.1) [20] and analyzed the different binding domains based on FDR value less than 0.05 and annotated by DiffBind [21]. The peaks on certain genomic loci were visualized by Integrative Genomics Viewer (IGV).

MeDIP and hMeDIP sequencing data were deposited to ArrayExpress assigned with the accession number E-MTAB-7810.

Statistical analysis

The experimental data were processed with SPSS 20 software (IBM, Armonk, NY, USA). One-way ANOVA was used for comparison of the difference among groups. Spearman rank correlation test was used to evaluate the relationship between 5mC/5hmC and tumor malignancy. p Values less than .05 was considered as statistical significance.

Results

The overall pattern of 5mC and 5hmC in cervical cancer tissues

Initially, the contents of 5mC and 5hmC were investigated in twelve human CSCC (six I–IIa and six IIb–IV) and three normal tissues of cervicitis as a control group. Dot blot and IHC were employed to detect the whole DNA methylation and hydroxymethylation levels. We observed that the total 5mC was up-regulated (Figure 1(A–D)) while 5hmC was down-regulated (Figure 1(E–H)) with the increased malignancy of cervical cancer group compared to control group. Furthermore, we also found that 5mC (r = 0.648, p < .05) and 5hmC (r = −0.717, p < .05) levels were significantly correlated with the tumor stages analyzed by Spearman rank test. Altogether, our results presented the aberrant overall DNA methylation and hydroxymethylation levels in cervical cancers.

Figure 1. The overall levels of 5mC/5hmC in CSCC. The dot plot (A) and IHC staining (C) with the corresponding statistical analysis (B, D) of 5mC in CSCC. The dot plot (E) and IHC staining (G) with the corresponding statistical analysis (F, H) of 5hmC in CSCC. Experiments were performed in triplicates. “*” represents the comparison between I-IIa group of CSCC and cervicitis with p values less than 0.05. “#” represents the comparison between IIb and IV group of cervicitis and control with p values less than .05.

The global genomic distribution of 5mC and 5hmC in cervical cancer tissues

Next, genomic DNA from four CSCC tissues of two I–IIa and two IIb–IV stage as well as two cervicitis were isolated and conducted 5mC and 5hmC immunoprecipitation approach combined with deep sequencing to map genome-wide 5mC and 5hmC profiles in cervical cancers. Approximately, 29.6 M reads of 5mC and 28.8 M reads of 5hmC sequencing data were collected, respectively (Table 1). To compare the individuals in cervical cancer and control groups, 5mC and 5hmC densities were normalized. Also, MA plot for a genome-wide comparison of 5mC/5hmC levels between IIb and IV stage CSCC and cervicitis groups displayed that most of the scatter plots are symmetrically distributed at the two sides of the 0 axis, which means that 5mC/5hmC sites were similar between the groups. Thus, no significant genome-wide differences in the 5mC/5hmC densities between the CSCC and cervicitis groups were observed (Figure 2(A,B)). After calling peaks by MACS software (p < .01), we mapped the (h)MeDIP-seq signals of 5mC and 5hmC peaks according to their genomic location and we did observe that 5mC was highly enriched at gene body and TES more than promoter regions (Figure 2(C)), while 5hmC was majorly enriched at the promoter and TSS regions (Figure 2(D)). The variation trend of the 5mC/5hmC levels had no significant difference among the three groups.

Figure 2. The genomic distribution of 5mC/5hmC in CSCC. The MA plot of IIb–IV group of CSCC and cervicitis contrast of 5mC (A) and 5hmC (B) peaks normalized with tag density. The X axis indicates the normalized mean and the Y axis indicates the log2-fold change. The red pixels mean the differential peaks with FDR less than 0.05. The blue pixels mean the peaks with no significant difference. The genome-wide mapping of 5mC (C) and 5hmC (D) in cervicitis and different stages of CSCC. Normalized DMRs and DHMRs tag density distribution across the gene body. Each gene body was normalized to 0–100%. Normalized Tag density is plotted from 20% of upstream of transcription start sites (TSS) and 20% of downstream of transcription terminal sites (TTS). -3 kb to TSS and +3 kb to TTS is shown.

Table 1. The summary of MeDIP and hMeDIP-seq data.

Sample	Total raw reads	Total raw bases	Total clean reads	Total clean bases	Mapped reads	Mapped ratio (%)	Unique mapped reads	Unique mappedratio (%)
Cervicitis_1	4658933	3.5E + 08	4527698	3.4E + 08	4267998	94.26	3077184	67.96
Cervicitis_2	4866549	3.7E + 08	4729759	3.6E + 08	4457715	94.25	3215524	67.98
I–IIa_1	4789317	3.6E + 08	4654346	3.5E + 08	4387511	94.27	3163059	67.96
5hmC
I–IIa_2	4977856	3.8E + 08	4837617	3.7E + 08	4559156	94.24	3285198	67.91
IIb–IV_1	5107214	3.9E + 08	4963388	3.8E + 08	4679174	94.27	3374284	67.98
IIb–IV_2	5166223	3.9E + 08	5020926	3.8E + 08	4733070	94.27	3413340	67.98
Cervicitis_1	4601442	3.5E + 08	4590342	3.5E + 08	4483939	97.68	2973858	64.79
Cervicitis_2	4658933	3.5E + 08	4647606	3.5E + 08	4540046	97.69	3011401	64.79
I–IIa_1	4789317	3.6E + 08	4777779	3.6E + 08	4667329	97.69	3097092	64.82
5mC
I–IIa_2	4866549	3.7E + 08	4854865	3.7E + 08	4742431	97.68	3147006	64.82
IIb–IV_1	4877129	3.7E + 08	4865326	3.7E + 08	4752950	97.69	3151706	64.78
IIb–IV_2	4977856	3.8E + 08	4965935	3.7E + 08	4851746	97.70	3218583	64.81
Cervicitis_1	4658933	3.5E + 08	4637594	3.5E + 08	4311095	92.96	2923315	63.04
Cervicitis_2	4977856	3.8E + 08	4955132	3.7E + 08	4604933	92.93	3123040	63.03
I–IIa_1	5021006	3.8E + 08	4998307	3.8E + 08	4646036	92.95	3149488	63.01
Input
I–IIa_2	4789317	3.6E + 08	4767682	3.6E + 08	4432209	92.96	3007027	63.07
IIb–IV_1	4866549	3.7E + 08	4844124	3.6E + 08	4502630	92.95	3052616	63.02
IIb–IV_2	4089748	3.1E + 08	4071043	3.1E + 08	3784021	92.95	2566918	63.05

Furthermore, the screening condition of significantly different enrichment of 5mC or 5hmC was set to meet the following requirements: cervicitis > I–IIa > IIb–IV or cervicitis  <  I–IIa  <  IIb–IV stages of CSCC in order with FDR  <  0.05. We identified 53 differential peaks of 5mC and 5hmC among these samples, total of 9 down-regulated and 42 up-regulated differential methylation regions (DMRs) as well as 5 down-regulated and 42 up-regulated differential hydroxymethylation regions (DHMRs) (log₂FC >1 or  < -1 and FDR < 0.05) were observed in cervical cancer tissues compared to controls (Table 2). Thirty-seven DMRs and DHMRs have a similar variation trend while 8 have the opposite trend compared between CSCC and cervicitis (Figure 3(A)). The genomic distribution of these DMRs and DHMRs showed that 80% around domains were enrolled in 48 RefSeq genes, in which 31 genes only contained differential enrichment of 5mC at promoter regions (-2k to +2k of TSSs), 9 genes only included showed differential enrichment of 5mC at gene bodies and UTR regions as well as 8 genes had differential 5mC at both of these regions. Also, the rest 20% DMRs and DHMRs were associated with intergenic regions. (Figure 3(B)). In addition, gene ontology (GO) analysis showed to be closely related to Wnt signaling pathway, as well as hormone secretion and cell junction in the differential 5mC associated genes, while it is closely related to ERK/MAPK pathway in the differential 5hmC associated genes (Figure 3(C)). KEGG pathway enrichment analysis revealed that multiple important pathways were enriched, including Wnt pathway, ubiquitin-mediated proteolysis, necroptosis and cervical cancer in 5mC differential genes, while Rap1 pathway, cell junction and cytoskeleton in 5hmC differentia genes (Figure 3(D)). Altogether, we comprehensively analyzed the (h)MeDIP-seq data and determined the specific differential methylated and hydroxymethylated loci of cervical cancer.

Figure 3. The associated functions of DMR/DHMRs in CSCC. The count of selected DMRs/DHMRs with the similar or opposite variation trend of 5mC/5hmC compared among cervicitis, I–IIa and IIb–IV of CSCC (A). Upper and down arrow mean the up- and down-regulated 5mC compared among cervicitis, I–IIa and IIb–IV of CSCC. Left and right arrow mean the up- and down-regulated 5hmC compared among cervicitis, I–IIa and IIb–IV of CSCC. “-” means no variation of 5mC/5hmC among these three groups. The distribution of the selected DMRs/DHMRs in genomic regions (B). GO (C) and KEGG (D) analysis of the genes associated with DMRs/DHMRs.

Table 2. The annotation of RefSeq genes associated with DMRs/DHMRs.

chr	Start	End	Width	Annotation	Genechr	GeneStart	Geneend	Genelength	Genestrand	Distanceto TSS	Ensemble	Symbol	Gene name	5hmC（cancer vs.cervici tis）	5mC（cance r vs.cerviciti s）
chr10	1.04E + 08	1.04E + 08	2182	Promoter (< =1kb)	10	1.04E + 08	1.04E + 08	54717	2	0	ENSG000 00065618	COL17A 1	Collagen type XVII alpha 1 chain	No	Down
chr10	1.15E + 08	1.15E + 08	1131	Intron (uc001lcc. 4/57700, intron 16of 16)	10	1.15E + 08	1.15E + 08	3124	1	14472	ENSG000 00151553	FAM160 B1	Family with sequence similarity 160 member B1	Up	Up
chr10	1209185	1209837	652	Intron (uc009xh q.4/105, intron 7 of 9)	10	1183151	1200385	17235	2	−8801	ENSG000 00185736	ADARB2	Adenosine deaminase, RNA specific B2(inactive)	DOWN	DOWN
chr10	1.32E + 08	1.32E + 08	761	Intron (uc001lkx. 5/282973, intron 15 of 23)	10	1.32E + 08	1.32E + 08	41094	1	6889	ENSG000 00188385	JAKMIP 3	Janus kinase and microtubule interacting protein 3	Up	Up
chr10	17226768	17228949	2181	Promoter (< =1kb)	10	17214239	17226898	12660	2	0	ENSG000 00229124	VIM- AS1	VIM antisense RNA 1	No	Up
chr10	41882611	41884809	2198	Distal intergenic	10	42475543	42495336	19794	1	#####	ENSG000 00185904	LINC008 39	Long intergenic non-protein coding RNA 839	Up	Up
chr10	42104070	42104409	339	Distal intergenic	10	42475543	42495336	19794	1	#####	ENSG000 00185904	LINC008 39	Long intergenic non- protein coding RNA 839	Up	Up
chr10	43408670	43410867	2197	Promoter (< =1kb)	10	43385617	43408882	23266	2	0	ENSG000 00169813	HNRNP F	Heterogeneous nuclear ribonucleoprotein F	Up	Down
chr10	79067549	79069735	2186	Promoter (< =1kb)	10	78943328	79067895	124568	2	0	ENSG000 00224596	ZMIZ1- AS1	ZMIZ1 antisense RNA 1	Up	Up
chr11	1.16E + 08	1.16E + 08	2170	Promoter (< =1 kb)	11	1.15E + 08	1.16E + 08	335288	2	0	ENSG000 00182985	CADM1	Cell adhesion molecule 1	Up	Up
chr11	31817033	31820872	3839	Promoter (< =1 kb)	11	31784792	31817961	33170	2	0	ENSG00000007372	PAX6	Paired box 6	Up	Down
chr11	65558439	65560872	2433	Promoter (< =1 kb)	11	65539411	65558930	19520	2	0	ENSG000 00168056	LTBP3	Latent transforming growth factor beta binding protein 3	Up	Down
chr1	2146821	2147888	1067	Promoter (1–2 kb)	1	2145387	2185116	39730	1	1435	ENSG000 00067606	PRKCZ	Protein kinase C zeta	Down	Down
chr12	91719517	91720391	874	Distal intergenic	12	91986626	91993403	6778	2	273012	ENSG000 00257242	LINC016 19	Long intergenic non-protein coding RNA 1619	Up	Up
chr13	77975030	77977208	2178	Intron (uc001vks.3/100874 222, intron 1 of 2)	13	77989372	77990646	1275	2	13438	NA	LINC010 69	long intergenic non-protein coding RNA 1069	Up	Up
chr14	1.01E + 08	1.01E + 08	2518	Promoter (< =1kb)	14	1.01E + 08	1.01E + 08	49776	1	0	ENSG000 00214548	MEG3	Maternally expressed 3 (non-protein coding)	Up	Up
chr1	4653242	4655419	2177	Promoter (< =1 kb)	1	4654732	4777794	123063	1	0	ENSG000 00196581	AJAP1	Adherens junctions associate d protein 1	Up	Up
chr16	4425172	4427411	2239	Promoter (< =1 kb)	16	4425805	4456772	30968	1	0	ENSG000 00103423	DNAJA3	DnaJ heat shock protein family (Hsp40)member A3	Up	Up
chr16	48931526	48932225	699	Distal intergenic	16	48546097	48620148	74052	2	#####	ENSG000 00102921	N4BP1	NEDD4 binding protein 1	Up	Up
chr16	90031782	90032741	959	Promoter (2–3 kb)	16	90029541	90044961	15421	1	2242	ENSG000 00141013	GAS8	Growth arrest specific 8	UP	UP
chr17	181132	183414	2282	Promoter (< =1 kb)	17	142789	181636	38848	2	0	ENSG000 00272636	DOC2B	DoubleC2 domain beta	No	Up
chr17	31390129	31392485	2356	Promoter (< =1 kb)	17	31391624	31538217	146594	1	0	ENSG000 00131242	RAB11FI P4	RAB11 family interacting protein 4	Up	Up
chr17	59876159	59877817	1658	Promoter (< =1kb)	17	59874539	59878762	4224	2	945	ENSG00000108423	TUBD1	tubulindelta 1	Up	Up
chr17	62625907	62632834	6927	Promoter (< =1 kb)	17	62627401	62693597	66197	1	0	ENSG000 00011028	MRC2	Mannose receptor C type 2	Up	Up
chr17	81512564	81514760	2196	Promoter (< =1 kb)	17	81509971	81512874	2904	2	0	ENSG00000184009	ACTG1	actingamma 1	Up	Up
chr17	81523347	81525549	2202	Promoter (< =1 kb)	17	81511346	81523847	12502	2	0	ENSG00000184009	ACTG1	Actin gamma 1	Up	Down
chr17	82092446	82094674	2228	Exon(uc002kd u.5/2194, exon 4 of 43)	17	82078340	82098187	19848	2	3513	ENSG000 00169710	FASN	fatty acid synthase	Up	Up
chr18	78978515	78981707	3192	Promoter (< =1 kb)	18	78980275	78997397	17123	1	0	ENSG000 00256463	SALL3	Spalt like transcripti on factor 3	Up	Up
chr20	21704172	21707170	2998	Promoter (< =1 kb)	20	21705659	21715982	10324	1	0	ENSG00000125813	PAX1	Paired box 1	Up	Up
chr20	57192385	57197334	4949	Intron(uc010gip.2/655, intron 3 of 6)	20	57177696	57187156	9461	2	−5230	ENSG000 00101144	BMP7	Bone morphogenetic protein 7	Up	Up
chr20	57215678	57220266	4588	Promoter (< =1 kb)	20	57214872	57215866	995	1	807	NA	BMP7- AS1	BMP7 antisense RNA 1	Up	Up
chr20	59298942	59301166	2224	Promoter (< =1 kb)	20	59300427	59325992	25566	1	0	ENSG00000124205	EDN3	Endothelin 3	Up	Up
chr20	61491753	61492584	831	Intron(uc032po k.1/1002, intron 2 of 15)	20	61570135	61937218	367084	1	−77551	ENSG000 00179242	CDH4	Cadherin 4	Down	No
chr20	63101496	63102157	661	Promoter (< =1 kb)	20	63102205	63104386	2182	1	−48	NA	HAR1A	Highly accelerated region 1A (non-protein coding)	Up	Up
chr21	10676553	10677170	617	Distal intergenic	21	10521575	10605716	84142	1	154979	ENSG000 00274391	TPTE	Transmembrane phosphatase with tensin homology	Up	Up
chr21	10698147	10699885	1738	Distal intergenic	21	10521575	10605716	84142	1	176573	ENSG000 00274391	TPTE	Transmembrane phosphatase with tensin homology	Up	Down
chr21	10699915	10727771	27856	Distal intergenic	21	10521575	10605716	84142	1	178341	ENSG000 00274391	TPTE	Transmembrane phosphatase with tensin homology	Up	Up
chr2	1.91E + 08	1.91E + 08	2186	Promoter (< =1 kb)	2	1.91E + 08	1.91E + 08	30381	2	0	ENSG000 00189362	NEMP2	Nuclear envelope integral membrane protein 2	Up	Up
chr4	49102367	49104228	1861	Distal intergenic	4	49028692	49031276	2585	1	73676	ENSG000 00109182	CWH43	Cell wall biogenesis 43 C-terminal homolog	Up	Up
chr4	49110432	49111161	729	Distal intergenic	4	49028692	49031276	2585	1	81741	ENSG000 00109182	CWH43	Cell wall biogenesis 43 C-terminal homolog	Up	Down
chr5	1.51E + 08	1.51E + 08	2147	Promoter (< =1 kb)	5	1.51E + 08	1.51E + 08	8555	1	0	ENSG000 00211445	GPX3	Glutathione peroxidase 3	Down	Up
chr6	2971675	2973873	2198	Promoter (< =1 kb)	6	2948159	2971856	23698	2	0	ENSG000 00124570	SERPIN B6	Serpin family B member 6	Up	Up
chr7	38025211	38027380	2169	Promoter (< =1 kb)	7	37909659	38025695	116037	2	0	ENSG000 00106483	SFRP4	Secreted frizzled related protein 4	No	Up
chr7	72446668	72448836	2168	Promoter (< =1 kb)	7	71779491	72447061	667571	2	0	ENSG00000183166	CALN1	Calneuron 1	Down	Up
chr8	1.24E + 08	1.24E + 08	801	Intron(uc003yq y.3/157376, intron 1 of 4)	8	1.24E + 08	1.24E + 08	125450	2	18854	ENSG000 00253868	FER1L6- AS2	FER1L6 antisense RNA 2	No	Up
chr8	1.41E + 08	1.41E + 08	649	Distal intergenic	8	1.41E + 08	1.41E + 08	2297	1	5267	ENSG000 00253595	LINC013 00	Long intergenic non-protein coding RNA 1300	Up	No
chr8	1.42E + 08	1.42E + 08	2264	Intron (uc003yw k.4/203062, intron 13 of 13)	8	1.42E + 08	1.42E + 08	10648	1	23195	ENSG000 00254008	LINC000 51	Long intergenic non-protein coding RNA 51	Up	Up
chr8	41308229	41311186	2957	Promoter (< =1 kb)	8	41261958	41309497	47540	2	0	ENSG000 00104332	SFRP1	Secreted frizzled related protein 1	No	Up
chr8	54456446	54459252	2806	Promoter (< =1 kb)	8	54457935	54460888	2954	1	0	ENSG00000164736	SOX17	SRY-box17	Up	Up
chr9	1.15E + 08	1.15E + 08	2194	Promoter (< =1 kb)	9	1.14E + 08	1.15E + 08	102890	2	0	ENSG00000095397	WHRN	Whirlin	Up	Up
chr9	1.22E + 08	1.22E + 08	593	Distal intergenic	9	1.22E + 08	1.22E + 08	8187	1	46250	ENSG000 00136848	DAB2IP	DAB2 interacting protein	Up	Up
chr9	96619350	96621535	2185	Promoter (< =1 kb)	9	96490241	96619830	129590	2	0	ENSG000 00081377	CDC14B	Cell division cycle 14B	Up	Up
chrY	10671476	10674208	2732	Distal intergenic	24	9910798	9911962	1165	1	760679	ENSG000 00239225	TTTY23	Testis-specific transcript, Y-linked 23 (non-protein coding)	Up	Up

Verification of the novel hydroxymethylated markers of cervical cancers

Finally, to find some novel potential biomarkers for cervical cancer diagnosis or prognosis, we drew the “vector picture” to clarify the sensitivity of variation trend of 5mC/5hmC (Figure 4). The 5mC (Y-axis) and 5hmC (X-axis) levels of the differential DMRs and DHMRs in three groups of cervicitis, I–IIa and IIb–IV of CSCC were documented, respectively (all the values of cervicitis were normalized as the original point). These three points of each gene were linked with the orientation from cervicitis to I–IIa to IIb–IV of CSCC. Here, two parameters could be described as biological meanings. One was the slopes of each vector between cervicitis to I–IIa or I–IIa to IIb–IV, which indicated the differential speed of variation between 5mC and 5hmC and the area A and C suggested the faster change of 5hmC while area B and D suggested the faster change of 5mC in the malignancy of CSCC. Also, the other was the length of each vector, which indicated the differential amplitude of variation, namely, the detectability of 5mC and 5hmC in cervical cancer. Here, we observed the DHMRs with a sensitively increasing 5hmC enrichment in ACTG1, SALL3, DNAJA3, SERPINB6 and CDC14B as well as a sensitively decreasing 5hmC enrichment in CALN1 compared to 5mC in cervical cancer. Likewise, VIM-AS1, LINC01069, CADM1, GPX3, SFRP4, AJAP1 and COL17A1 had the sensitive change of 5mC in CSCC. These results above suggested a list of potential hallmarks for CSCC diagnosis and prognosis.

Figure 4. Vector picture for clarification of the sensitivity of variation trend of 5mC/5hmC. The area of A and C mean the associated genes with 5hmC sensitive change while area B and D mean the associated genes with 5mC sensitive change in the malignancy of CSCC. The X/Y axis indicates the normalized log2-fold change of 5mC/5hmC.

Discussion

DNA methylation as a well-studied epigenetic hallmark has been determined to play a crucial role in cancers and other human diseases. In general, the presence of aberrant DNA hypermethylation is observed to cause the abnormal chromosome activity and gene silencing in most of the cancers [22]. As we know, DNA methyltransferases (DNMTs) contribute to the hypermethylation via the transfer of the methyl groups on the cytosines. Thus, DNMT inhibitors such as decitabine or azacytidine have been applied or combined with other drugs for the treatment of solid cancers [23]. Nevertheless, the persistently accumulated hypermethylation status of cancers also implies the other side that the process of DNA demethylation is blocking. For DNA demethylation, passive demethylation reduces the global 5mC level via suppression of DNA methylation information inheritance in daughter cells, while the active demethylation happens in situ to restore the 5mC to regular cytosine. It is widely acknowledged that 5mC is oxidized to 5hmC, 5foC and 5caC step by step and finally is decarboxylated to regular cytosine by thymine-DNA glycosylase (TDG). Here, we determine that the global 5mC is up-regulated while 5hmC was down-regulated (Figure 1), and 5hmC is negatively correlated with the malignancy of cervical cancer, which indicates that the active DNA demethylation is restrained in cervical cancer. We speculate that the active demethylation process is terminated and hydroxymethylation is blocked so that the overall hypermethylated DNA is abnormally accumulated in cervical cancer.

Since the overall 5mC/5hmC levels are investigated by dot plot and IHC, the genomic distribution of 5mC/5hmC are further addressed by (h)MeDIP. In general, 5mC are well studied in promoter regions for gene regulation, but our results show the pattern of the hypermethylated CpG loci or islands in gene body or UTR more than the promoter (Figure 2(D)). However, the hydroxymethylated CpG sites are majorly located at the promoter (Figure 2(C) and 3(B)), which suggests that most of the hydroxymethylation active zone for the occurrence and function of 5hmC are the gene promoter regions. Also, we find a significant difference of 5hmC levels at promoters rather than other regions in IIb–IV stage of CSCC compared to I–IIa (p=.037) and control (p=.012), but no difference of 5mC among three groups (Figure 2(C)). Moreover, we notice that most of the DMRs/DHMRs have a variation trend in the genetic regions (Figure 3(A)). Nevertheless, the up-regulated 5mC and down-regulated 5hmC are observed in CSCC using dot plot and IHC, which seems paradoxical with the presence of genetic distribution variation of 5mC and 5hmC. We speculate that 5hmC is supposed to be reduced in intergenic regions more than the elevation in genetic regions. Previous studies suggest that 5hmC enriches in intergenic regions besides introns and exons [24,25] in other disease models, which indicates that the regulatory way of 5hmC in non-genetic regions is different and may play a complicated role in gene regulation compared to DHMRs at the genetic regions. Consistently, the associated signaling pathways and biological functions related to DMRs and DHMRs are also different. DMRs enriches in the genes closely associated with Wnt signaling pathway (Figure 3(C)) while DHMRs may control ERK/MAPK, Rap1 and Hippo pathways (Figure 3(D)), which also implies the different roles of 5mC and 5hmC for gene regulation in CSCC.

Finally, we find that the 5hmC at the genes of ACTG1, SALL3, DNAJA3, SERPINB6, CDC14B and CALN1 are sensitive and specific to CSCC malignancy (Figure 4), which may provide a batch of novel 5hmC biomarkers for cervical cancer diagnosis. The epigenetic hallmarks such as DNA methylation/hydroxymethylation can usually reflect the earlier pathological events in cells compared to the transcriptome when diseases happen. In our case, plenty of regions with stable hypermethylation among the three groups are difficult to understand the actual reason, permanent hypermethylation with no change or the suppressive hydroxymethylation mediated hypermethylation. Thus, 5hmC of hallmarks can improve this issue. Secondly, CpG loci with 5hmC modification are extremely rare in genome compared to 5mC and can be accurately detected for evaluating cancer status. Altogether, our study describes the genomic landscape of 5hmC in cervical cancer and provides the significant 5hmC biomarkers of CSCC.

Authors’ contributions

J.W. and Y.S. performed the experiments and analyzed the data; Y.T. and Y.D. helped perform the experiments; X.W. designed the overall study and drafted the manuscript.

Disclosure statement

The authors declare no competing financial interests exist.