Epigenomic analysis of 5-hydroxymethylcytosine (5hmC) reveals novel DNA methylation markers for lung cancers

Background DNA methylation at the fifth position of cytosine (5mC) is a common epigenetic alteration affecting a range of cellular processes. In recent years, 5-hydroxymethylcytosine (5hmC), an oxidized form of 5mC, has risen broad interests as a potential biomarker for lung cancer diagnosis and survival. Methods We analyzed the epigenome-wide 5hmC profiles of paired lung tumor and adjacent normal tissues, using the TET-Assisted Bisulfite (TAB) array – Infinium MethylationEPIC BeadChip (EPIC) approach. The differentially methylated CpG sites were identified, and the biological relevance of 5hmC was assessed by differential methylation regions (DMR) analysis and gene set analysis (GSA). Results We observed global hypomethylation of 5hmC comparing tumor to normal tissues, and hypermethylated 5hmC were enriched in CpG islands and gene upstream. Comparison of 5hmC and 5modC (total methylation: 5mC + 5hmC) profiling showed low correlation, and only a small proportion of the significant 5hmC loci overlapped with the significant total methylation loci. GSA analysis suggested that 5hmC was mainly involved in biological processes such as cellular process, biological regulation, and metabolic process. Conclusion This is the first study to analyze the epigenome-wide 5hmC profiles among paired lung tumor and normal tissues. We observed global hypomethylation of 5hmC in lung cancers, and hypermethylated 5hmC enriched in CpG islands and gene upstream. We found that the genome-wide 5hmC levels do not correlate with the total methylation, and the GSA suggested different biological functions of 5hmC compared to 5modC. Overall, our results demonstrate the potential of 5hmC as a novel biomarker for lung cancer.


Introduction
Lung cancer is the leading cancer-related cause of death in the world. There were an estimated 1.8 million lung cancer deaths globally [10,11] and 154 thousand deaths [45] in the United States in 2018 alone. The current five-year relative survival is about 19% [37], which is lower than most leading cancers and has improved only slightly over the past four decades. The five-year relative survival for localized tumors is 56%; however, only 16% of all lung cancer cases are diagnosed at an early stage and the five-year relative survival for distant tumors is only 5%. Therefore, early detection and intervention remain to be some of the most important strategies to improve long-term outcomes.
The rapid advancement of affordable`omics'-based technologies has allowed faster identification of putative biomarkers for improving lung cancer diagnosis and prognosis [48]. As one of the most important epigenetic mechanisms, DNA methylation at the fifth position of cytosine (5mC) plays a vital role in many cellular processes, such as embryonic development, gene expression, genomic imprinting, X-chromosome inactivation, and chromatin structure and stability [9,12,28,33,41]. DNA methylation abnormalities, such as global hypomethylation [13] and hypermethylation at promoters of tumor suppressor genes [18], are found in most tumors including lung cancer.
In recent years, 5-hydroxymethyl-cytosine (5hmC), an oxidized form of 5mC by Fe(II) and a-KG-dependent Ten-eleven translocation (TET) dioxygenases [24,47], has risen broad interests. It has been purported that 5hmC could act as an intermediate in DNA demethylation, either through active TET-assisted oxidation [20], or passive replication dependent loss [7,19]. Others also suggested that 5hmC could confer a distinct epigenetic function independent of DNA demethylation [7].
Early studies were made through global detection approaches, and global loss of 5hmC has been commonly observed across a spectrum of human tumors [15,22,27,30,51]. Mutations in IDH, SDH, and FH genes can control 5hmC genomic levels and trigger carcinogenesis by regulating TET activity [5,8,23,31].
Given the important epigenetic functions of 5hmC and its interplay with 5mC, it's essential to develop accurate and reliable methods to distinguish 5hmC from 5mC at single-base resolution, to better understand the biological mechanisms of 5hmC in cancer etiology and progression. However, the conventional bisulfite sequencing, which is the gold standard for methylation analysis, cannot distinguish 5hmC from 5mC. In the study, we analyzed the whole-genome 5hmC profiles of paired tumor and adjacent normal tissues from 8 patients with primary lung squamous cell carcinoma (LUSC), using a state-of-the-art technique of the TET-Assisted Bisulfite (TAB) array [52] with Infinium MethylationEPIC BeadChip (EPIC) [34]. We propose to use this as an exploratory study to assess the potential of using 5hmC as a novel epigenetic biomarker for lung cancer, and to lay the foundation for future clinical assay development and validation studies.

DNA samples
Genomic DNA was extracted for eight lung squamous cell carcinoma (LUSC) tissues and adjacent normal tissues. Tissues samples were obtained from lung cancer patients in the Harvard Lung Cancer Study (HLCS) who underwent surgical excision of tumors at Massachusetts General Hospital or Dana-Farber Cancer Institute. Manual microdissections of 5-u histopathologic sections have been conducted by one of our anatomic pathologists. Each specimen was evaluated for the amount and quality of tumor cells and histological classification was conducted by a pathologist using the WHO criteria. Tissue DNA was extracted from 50 mg of frozen tumor or normal tissue using the Maxwell 16 tissue DNA purification kit (Promega, USA). Tumor purity has been analyzed using LUMP analysis [4], which uses the mean normalized methylation level of 44 CpGs found to be unmethylated in blood cells.
5hmC and 5mC profiling Illumina MethylationEPIC with bisulfite and TET-assisted bisulfite treatment The 5mC and 5hmC profiling was completed at the Northwestern University Population Epigenetics Laboratory. A total of 1000 ng genomic DNA (gDNA) was split into two aliquots of 500 ng processed through either traditional bisulfite conversion (Bis) or TET-assisted bisulfite conversion (TAB) (WiseGene, Inc., Chicago, IL) [36]. As shown in Fig. 1, bisulfite treatment alone converted unmethylated cytosine into uracil that would be read as thymine (T) after PCR amplification, whereas both 5mC and 5hmC were resistant to conversion and would be read as cytosine (C). The TAB workflow involved protection of 5hmC with bglucosyltransferase, followed by TET-assisted oxidation of unprotected 5mC to 5caC, which would be converted to uracil by bisulfite treatment and be read as thymine just as unmethylated cytosine. Therefore, Bis workflow give readout of total methylation (5mC + 5hmC), while TAB workflow gives direct readout of 5hmC. Converted DNA of both workflows was input into the Illumina MethylationEPIC BeadChip (Illumina, Inc., San Diego, CA) [40].

Whole genome bisulfite sequencing with Illumina Methyl-Seq
To validate our results and to establish the feasibility of using next generation sequencing (NGS) to profile 5hmC in lung tissues, the total methylation profiling was repeated among 5 pairs of tumor and adjacent normal tissues using the targeted bisulfite sequencing (Methyl-Seq). Briefly, 1000 ng of gDNA was sheared to yield fragment DNA between 100-300 base pair. These fragment DNA was further end-repaired, 3 0end adenylated, ligated, and hybridized; and captured DNA was bisulfite-converted using the Illumina TruSeq Methyl Capture EPIC Library Prep Kit (Illumina, Inc., San Diego, CA) [32]. Last, DNA was PCR-amplified and sequenced on Illumina Hi-Seq2500 System.

Illumina MethylationEPIC Bis/TAB data
The raw Tab-array and EPIC array data in IDAT format were read directly into R (version 3.2.4) and were processed using the Chip Analysis Methylation Pipeline (ChAMP) [35].
Quality control was conducted separately for the Bis and TAB workflow using the same criteria. We excluded unqualified probes, including: 1) probes with detection p value less than 0.01 in one or more samples (N Bis = 1476, N TAB = 84,125); 2) probes with a bead count < 3 in at least 5% of samples (N Bis = 14,963, N TAB = 23,966); 3) Non CpG probes (N Bis = 2947, N TAB = 2526); 4) SNP-related probes [53] (N Bis = 96,512, N TAB = 87,925); and 5) multi-hit probes [38] (N Bis = 11, N TAB = 10). Note that the number of probes removed due to low detection p value is much higher for 5hmC as the genome-wide 5hmC level is much lower than that of 5mC. Because the p value filter was applied first, the number of probes deleted for annotation purposes were different for 5hmC and 5-modC accordingly. After filtering, 750,009 and 667,366 probes were included for Bis and TAB workflow respectively. Normalization was performed to minimize unwanted variation within and between samples and we chose peak-based correction (PBC) [11], which has been shown as the most effective method when comparing with the whole genome bisulfite sequencing results [49].
For each CpG site, the methylated intensity (M) and unmethylated intensity (U) was transformed to b ¼ M M þU . The beta values with Bis conversion represent the proportion of total methylation (5modC = 5mC + 5hmC), whereas the beta values with TAB conversion represent the proportion of 5hmC methylation for each CpG site. The M values were calculated as M ¼ log 2 ð M U Þ and were used for regression analyses. To detect the differential methylation patterns between tumor and normal tissues, differential methylation probes (DMP) analysis were conducted using the paired t-test with Limma [42]. The DMPs were grouped into clusters/regions and differential methylated regions (DMRs) were identified using the Bumphunter algorithm [21]. The functional relevance of 5hmC/5modC-containing genes was explored using Gene Set Analysis (GSA) [39] with 15,996 Gene Ontology Biological Process (GOBP) pathways [6], while accounting for the number of probes per gene [14]. False discovery rate (FDR) of <0.05 was used in all analyses to account for multiple comparisons.

Whole genome bisulfite sequencing data
Whole genome bisulfite sequencing data was analyzed in collaboration with Harvard Chan Bioinformatics Core. FASTQC (Version 0.11.8) [3] was used for quality control to examine sequence counts, sequence quality histograms, per sequence quality scores, per base sequence GC content, per base N content, sequence length distribution, sequence duplication levels, overrepresented sequences, and adapter content. We observed very low proportion of cytosine bases in R1s (forward reads or read-1 s of the read pairs), suggesting good bisulfite conversion. Reads were trimmed using Trim Galore (Version 0.6.0) [25] after assessment by M-bias plot [29] to remove low quality sequences, adapter sequences, and biased sequences. FASTQC was rerun after trimming to confirm efficacy. Bismark (version 0.21.0) [26] was used to align the reads to the Ensembl GRCh37 Homo sapiens assembly [2], and duplicates reads removed using Picard Tools (http://broadinstitute.github.io/picard/). Bismark also summarized the cytosine counts at cytosines followed by guanines (CpGs), cytosine followed by non-guanines followed by guanines (CHGs), and cytosines followed by at least two non-guanines (CHHs). We observed very low proportion (<0.2%) of CHG and CHH methylation, suggesting good bisulfite conversion.
CpG methylation percentages and read coverage of CpG generated from Bismark was assessed for adequate coverage (>10) and expected methylation percentage distributions in R using methylKit (version 1.8.1) [1]. The Bismark generated data was then imported into R using the bsseq Bioconductor package [16] and annotated with Ensembl GRCh37 Homo sapiens assembly using annotatr (version 1.8.0). Differential methylation analysis for the sequencing data was conducted using a general linear model accounting for tumor-normal pairing in DSS (version 2.30.1) [50], a new dispersion shrinkage method that estimates the dispersion parameter from beta-binomial distribution. Again, test statistics were calculated at each differentially methylated loci (DMLs) and significant DMLs were aggregated into DMRs.

Study subjects and sample characteristics
Tissue samples were obtained from 8 Caucasian males with early stage (stage 1 and 2) lung squamous cell carcinoma. The mean age at diagnosis is 69.2 years old (range = 63.1-76.9 years).
Using the LUMP method for tumor purity, the average tumor purity is around 40%, with the average normal tissue purity to be around 39%.

Genomic alterations of 5hmC among lung tissues
Using the TAB -EPIC and Bis -EPIC array, we profiled the 5hmC and 5modC levels in 8 pairs of tumor and adjacent normal samples ( Supplementary Fig. 1). While the distribution of total methylation demonstrated a typical two-mode distribution with peaks at both low and high b levels ( Supplementary Fig. 1B), the distribution of 5hmC showed a unimodal distribution with the b values of most probes to be less than 0.1 (Supplementary Fig. 1A).
We further examined the genomic 5hmC and 5modC alterations by CpG island (CGI) status as well as genomic features. With regards to CGI status, we observed similar patterns for both 5hmC and total methylation: there was higher likelihood of hypomethylation at opensea and shelf regions as well as hypermethylation at shore and CGI regions ( Supplementary Fig. 1C and D). For genomic features, 5hmC was hypermethylated at gene promoter regions (TSS1500 and TSS200), 5 0 UTR, and the 1st exon, but hypomethylated at gene bodies ( Supplementary Fig. 1E).

Comparison of 5hmC and 5modC profiling
We further compared the 5hmC and 5modC profiling measured by different methods. As 5hmC was considered as an intermediate of DNA demethylation, we examined the correlation between 5hmC and 5modC across CpG sites with significant differential hydroxymethylation (N = 25,611) (Fig. 3A) and across CpG sites with significant differential total methylation (N = 34,975) (Fig. 3B). The 5hmC levels at the significant CpG sites showed poorly correlation with the total methylation, suggesting potential distinct epigenetic function of 5hmC that is independent of 5mC. The average differences of 5hmC between tumor and normal tissues were also poorly correlated (r = 0.11) with the average difference of 5modC across genome (N = 656,768), though the average differences in 5modC measured by different methods (EPIC and MethylSeq) showed fairly high concordance (r = 0.74). The Venn diagram comparing the significant 5hmC loci (N = 25,611) and the significant 5modC loci (N = 34,975) found only 1441 (2.4%) overlapping CpG sites for both biomarkers, and among them 1288 loci (89.4%) were hypomethylated for both 5hmC and 5modC (Fig. 3E).
We conducted gene set analysis (GSA) to understand the biological processes that the differentially methylated probes/regions are involved in. In total, 430 and 189 GOBP pathways were significant for 5hmC and 5modC at FDR < 0.05 level, respectively. The top 10 significant pathways for 5hmC and 5modC are shown in Supplementary Table 2. We further explored the biological relevance of 5hmC by categorizing the GOBP pathways into 17 top-level categories (cellular process, biological regulation, metabolic process, multicellular organismal process, developmental process, response to stimulus, localization, cellular component organization or biogenesis, multi-organism process, signaling, immune system process, reproduction, locomotion, cell population proliferation, growth, biological adhesion, and behavior). For each category, we calculated the expected number of significant pathways and compared with the observed number of significant pathways in each category (Fig. 4). We found that the differentially methylated 5hmC genes were more likely to congregate at pathways involved in cellular process, biological regulation, and metabolic process, whereas the differentially methylated 5modC genes were more frequent among pathways related to localization, signaling, and multicellular organismal process.

Discussion
Though the incidence of lung cancer has decreased by 3% per year in men and 1.5% per year in women from 2011 to 2015 in the US [46], it remains to be the most common type of cancer and the leading cause of cancer-related death globally [10]. What's more, the current five-year relative survival of 19% is lower than most common cancers and has only improved slightly from 12% over the past three decades [37]. Therefore, sensitive and specific biomarkers for lung cancer are in urgent needs to improve diagnosis and long-term survival.
Due to the high tissue specificity and availability of high-throughput method for quantification, DNA methylation has been considered as a promising cancer biomarker and has been studied extensively in most cancers [17]. However, most studies on DNA methylation and cancer are based on bisulfite conversion methods, which cannot distinguish between 5hmC and 5mC. In this study, we applied a start-of-the-art technique of TET-Assisted Bisulfite (TAB) array with Illumina MethylationEPIC (EPIC) BeadChip, and quantified the novel 5hmC marker at 850,000 CpG sites across human genome.
In this study, we observed a global depletion of 5hmC with most CpG sites hypomethylated in tumors, which agrees with previous studies [15,22,27,30,51]. Notably, we found hypermethylation of 5hmC among tumors that were enriched in gene upstream (e.g. promoters), and previous in vitro studies reported that the presence of 5hmC at the promoter strongly suppressed transcription [43].
Though historically considered as an intermediate of DNA demethylation, the locus-specific 5hmC levels at differentially methylated CpG sites do not correlate with the 5modC levels, indicating that 5hmC may not simply be a demethylation intermediate. Moreover, comparison of the significant 5hmC loci with the significant 5modC loci found very few overlaps, and the differentially methylated 5hmC and 5modC genes were enriched in different biological systems, suggesting the potential of 5hmC as a new, independent biomarker for lung cancer diagnosis.

Conclusion
In this exploratory study, we analyzed the epigenome-wide 5hmC profiles among paired tumor and normal tissues using the TAB array -EPIC approach. We observed global hypomethylation of 5hmC in lung cancers, and hypermethylated 5hmC were enriched in CpG islands and gene upstream. We found that the genome-wide 5hmC levels do not correlate with total methylation, and gene set analysis revealed different biological functions of 5hmC. Overall, our results indicated the potential of 5hmC as a novel biomarker for lung cancer, and further work is needed to validate differentially methylated loci for cancer diagnosis and prognosis and to explore the functional role of 5hmC in lung cancer etiology.