Investigation on the potential of circulating tumor DNA methylation patterns as prognostic biomarkers for lung squamous cell carcinoma.

BACKGROUND: Aberrant epigenetic modifications play a key role in lung tumorigenesis. In our study, we aimed to explore the clinical implications of baseline circulating tumor DNA (ctDNA) somatic and methylation profiles in patients with lung squamous cell carcinoma (LUSC).
METHODS: A total of 26 patients with LUSC of various stages were included in this study. Somatic mutations and methylation levels were profiled from the plasma-derived ctDNA obtained at the time of diagnosis using unique molecular identifier (UMI)-based targeted sequencing and bisulfite sequencing, respectively. The correlation between baseline ctDNA mutation and methylation profile, and overall survival (OS), were analyzed.
RESULTS: Somatic mutations were detected in 80.8% (20/26) of the patients. Patients harboring somatic mutations with maximum allelic fraction (maxAF) of >5% had significantly shorter OS compared to those with maxAF ≤5% (7.1 vs. 54.6 months; P=0.020). ctDNA methylation level was found to be strongly correlated with maxAF (Pearson correlation =0.934; P<0.001). Consistent with maxAF, higher methylation levels were also associated with poorer OS (hazard ratio =2.377; 95% CI: 1.283-4.405; P=0.006). Moreover, a total of 1,956 ctDNA methylation blocks were differentially methylated in patients with maxAF >0 (P<0.05). Least absolute shrinkage and selection operator (LASSO) regression analysis revealed a significant correlation between methylation signatures from 5 methylation blocks and OS (hazard ratio =183.20, 95% CI: 2.74-12,243.32; P=0.015). These 5 methylation blocks could serve as an alternative to maxAF and can be explored as prognostic biomarkers.
CONCLUSIONS: Our study identified several ctDNA methylation blocks that can potentially predict the prognosis of LUSC at the time of diagnosis. 2020 Translational Lung Cancer Research. All rights reserved.

Introduction

Lung cancer remains the leading cause of cancer-related death worldwide, including China (1-4). Among the non-small cell lung cancer (NSCLC) histological subtypes, lung squamous cell carcinoma (LUSC) is one of the most prevalent, accounting for approximately 15–30% of all lung cancers diagnosed in China (1,5-9). Histologically, LUSC is an epithelial malignancy characterized by keratinization and/or intercellular bridges and demonstrates positivity for squamous differentiation markers using immunohistochemistry (7). LUSC is strongly associated with exposure to tobacco and is predominantly diagnosed in males (1,8-10). The 5-year survival rate of LUSC is approximately 18% (1), with survival outcomes remaining unchanged over the last decade. The tumor-node-metastasis (TNM) staging system remains as the primary guide for the treatment and prognostication for NSCLC, including LUSC (11,12). As the standard of care, patients with early-stage LUSC are managed with surgery and adjuvant chemotherapy, with or without radiotherapy (13). Meanwhile, patients with advanced-stage LUSC are treated with cytotoxic chemotherapy as a first-line treatment, particularly with platinum-based doublet, and are recommended to undergo molecular testing for actionable mutations (13). However, due to the low frequency of actionable mutations typically found in lung adenocarcinoma, such as EGFR and ALK (14-16), the use of targeted therapy remains limited among patients with advanced LUSC. Late diagnosis and the lack of actionable mutations are major contributors to the poor prognosis of LUSC. Numerous efforts have been invested in identifying molecular factors with either prognostic or predictive value in LUSC; however, no biomarker has been validated for use in standard clinical practice (10,17,18). Thus, there is an urgent need to identify prognostic and predictive biomarkers to improve outcomes of patients with LUSC.

DNA methylation is a major epigenetic modification by which methyl groups are added to critical regions of DNA for the precise control of gene expression. The key biological processes under epigenetic control include the cell cycle, DNA repair, cell proliferation, and apoptosis, which are critical pathways involved in tumor development and progression (19,20). Aberrant DNA methylation, including genome-wide hypomethylation of repetitive elements and CpG-poor regions, and regional hypermethylation of CpG islands at gene promoters, has been implicated in the early stages of carcinogenesis in various solid tumors including LUSC (18,20-27). Numerous studies have explored the potential of epigenetic deregulation as a biomarker of early detection, cancer surveillance/disease monitoring, and prognosis of various human cancers (26,28-30).

Liquid biopsy samples are now being used as minimally-invasive alternative sources of tumor DNA for next-generation sequencing (NGS)-based genomic testing (31-34). The advances in molecular detection using liquid biopsy samples are particularly advantageous for patients who have tumors that are inaccessible or have a limited volume of tumor biopsy specimen. The presence of circulating tumor DNA (ctDNA) derived from blood samples, the most commonly used liquid biopsy specimen, is associated with minimal residual disease from resected early-stage lung cancer and serves as an early biomarker for cancer recurrence (29,34). With the possibility of longitudinal sampling, considerable attention has been invested towards exploring the utility of ctDNA methylation profiling for the early detection of cancer recurrence or relapse (30,33,35-43). Numerous genes have been identified to be differentially methylated in blood samples of patients with early-stage lung cancer, including SHOX2, RASSF1A, and p16, which could potentially serve as predictive biomarkers (30,33,35-42). As compared to lung adenocarcinoma, DNA hypermethylation of SHOX2 was highest in small cell and LUSC tumor tissue samples, indicating that methylation status is distinct among lung cancer histologies (40). However, investigations on ctDNA methylation patterns that enable the prediction of prognosis at the time of diagnosis were performed on lung cancer patients in general and independent investigation on LUSC is very limited. In our study, we aimed to identify differential methylation signatures in certain genomic CpG regions from blood samples at initial diagnosis that are associated with survival, and potentially predict survival outcomes in Chinese patients with LUSC. We present the following article in accordance with the REMARK reporting checklist (available at http://dx.doi.org/10.21037/tlcr-20-1070).

Methods

Patients

A total of 26 consecutive patients with stage IA–IV lung cancer with confirmed squamous cell histology diagnosed at our institution between March 2017 and September 2019 were prospectively included in this study. Clinical staging and histology were evaluated according to the eighth edition of the American Joint Committee on Cancer (11) and the 2015 World Health Organization classification of lung tumors (7), respectively. The main study inclusion criteria are as follows: (I) pathologically-confirmed squamous cell carcinoma; (II) provided consent to submit blood samples at diagnosis; (III) provided written informed consent to participate in the study. The main exclusion criteria are as follows: (I) non-squamous cell lung cancer histology; (II) history of other malignant tumors within the last 5 years; (III) failure to obtain blood samples prior to receiving treatment; (IV) refusal or withdrawal of consent. Baseline clinicopathological information was collected, including age at diagnosis, gender, smoking status, and treatment information, and is summarized in Table S1. The cohort was comprised of a majority of males (88.5%; 23/26), with a median age of 63 years (range, 38–78 years). All patients received therapeutic management based on their clinical stages following the guidelines of the National Comprehensive Cancer Network and the decision of the physician-in-charge. Baseline blood samples were obtained from each patient at the time of initial visit/diagnosis. The study was approved by our Institutional Review Board (approval number: NCC1576) and performed in compliance with the ethical standards of the institutional and/or national research committee. This study was also carried out in accordance with the Declaration of Helsinki (as revised in 2013). All patients provided written informed consent before inclusion in the study.

DNA isolation

DNA isolation and subsequent NGS library preparation were performed at Burning Rock Biotech, a College of American Pathologists (CAP)-accredited and Clinical Laboratory Improvement Amendments (CLIA)-certified clinical laboratory. In brief, circulating cell-free DNA (cfDNA) was isolated from plasma samples using QIAamp Circulating Nucleic Acid kit, according to the manufacturer’s standard protocol (Qiagen, Hilden, Germany). cfDNA quantity and quality were estimated using the Bioanalyzer 2100 instrument and high sensitivity DNA assay kit (Agilent Technologies, CA, USA). A minimum of 50 ng DNA was required for the NGS library construction.

Unique molecular identifier (UMI)-based targeted sequencing and data analysis

Targeted sequencing was performed as described previously (44). The cfDNA fragments were ligated with adapters containing UMIs. Target capture was performed using a panel consisting of 168 cancer-related genes (Burning Rock Biotech), spanning 273 kilobases of the human genome. The UMI-tagged indexed samples were subsequently sequenced on an Illumina NovaSeq 6000 (Illumina, San Diego, CA, USA) using 2×150 base pair (bp) cycles and a target sequencing depth of 30,000×. Sequencing data were analyzed using bioinformatics pipelines optimized to accurately distinguish between somatic mutations and sequencing artifacts.

Bisulfite targeted sequencing and data analysis

Whole-genome bisulfite sequencing library preparation was performed using the brELSA (Burning Rock Biotech) as described previously (35). Briefly, purified cfDNA were treated with sodium bisulfite (EZ-96 DNA Methylation-Lightning MagPrep, Zymo Research, Orange, CA, USA). The converted fragments underwent two series of ligation with appropriate adapters, amplification, and purification to obtain the whole-genome bisulfite sequencing library. Target enrichment was performed using custom-designed lung-cancer methylation profiling RNA baits covering 80,672 CpG sites, spanning 1.05 megabases of the human genome. The target libraries were finally quantified by real-time PCR (Kapa Biosciences, Wilmington, MA, USA) and sequenced on a NovaSeq 6000 (Illumina, San Diego, CA, USA) using 2×150 bp cycles. Bisulfite sequencing data analysis was performed using an optimized pipeline. Custom adapter sequences and low-quality bases were removed using Trimmomatic (v.0.32). Sequence alignment of paired-end reads to CtoT- and GtoA-transformed hg19 genome was performed using BWA-meth (v.0.2.2) (45). PCR duplicates from aligned sequences were marked with Samblaster (v.0.1.20) (46). Sambamba (v.0.4.7) (47) was then used to remove the sequencing reads with either low mapping quality (MAPQ <20) or improper pairing. To eliminate potential double-counting of methylation calls, overlapping reads from paired-read sequences were clipped and merged using in-house scripts.

Analysis of methylation patterns

Methylation levels in this study were analyzed and represented as scores to reflect methylation features of the sample. Using the methods described previously (35), the 80,672 CpG sites included in the panel were grouped into 8,312 methylation blocks to account for linkage disequilibrium. Methylation blocks were defined as the genomic region between the neighboring CpG sites with the r2 value calculated based on our modified correlation matrix. In brief, Pearson’s correlation analysis was performed to calculate the difference among the methylation frequencies of each pair of CpG sites, normalized by the difference in genomic distance and methylation level. Methylation frequencies for the entire set of methylation blocks were calculated as the ratio between the total number of Cs at all CpG sites within a methylation block, and the total number of C+Ts at the same methylation block.

Methylation pattern value (MPV1) was defined as the differential methylation pattern of CpG sites in a certain genomic region with an interval of 0–1 within a methylation block, and is represented by the following formula in Eq. [1]:

Wherein: n refers to the number of sequencing reads for the multiple CpG sites that cover the methylation block; L refers to the number of CpG sites within the methylation block that is covered by the i-th sequencing read; m refers to the number of methylation blocks with consecutively methylated CpG sites which appeared in the i-th sequencing read; l indicates the number of methylated CpG sites per methylation block.

Since each of the 8,312 methylation blocks had their individual MPV1 score, a methylation score was also calculated per patient sample using Eq. [2].

Wherein: n refers to the number of methylated blocks; T refers to the value corresponding to the i-th methylation block estimated for the patient circulating tumor DNA (ctDNA) sample; H refers to the value corresponding to the i-th methylation block estimated for a group of asymptomatic healthy samples, methylation data were derived from ctDNA mutation profiling of 198 asymptomatic healthy patients; N indicates the number of asymptomatic healthy patients included in the analysis; H refers to the value corresponding to the i-th methylation block for j-th sample from the asymptomatic healthy patients; refers to the mean value corresponding to the i-th methylation block for the remaining samples from the asymptomatic healthy patients.

Statistical analysis

Statistical analysis was performed using the R statistics package (R version 3.4.0; R: The R-Project for Statistical Computing, Vienna, Austria). Statistical significance was defined as P values <0.05. Differences between groups were calculated using a two-tailed Student’s t-test, as appropriate. Overall survival (OS) was defined from the date of diagnosis until death or the last day of follow-up. OS values corresponding to the last day of follow-up for patients who are still alive were censored and were indicated by tick marks in the survival curves. The data cut-off was October 1, 2019 with a median follow-up of 6 months (ranging from 1 month to 34 months). Survival differences between patients with certain maximum allelic fraction (maxAF) percentages were estimated using the Kaplan-Meier method and compared using the log-rank test. Univariate Cox proportional-hazards model was used to estimate the correlation between molecular features and survival outcome. Least absolute shrinkage and selection operator (LASSO) was implemented using the glmnet package to perform a regression analysis for identifying the smallest possible methylation block set that correlated with maxAF.

Results

Molecular profile of the cohort

To understand the baseline somatic mutation profile of the patients, targeted sequencing analysis was performed using blood samples obtained at the initial visit. Among the 26 patients with adequate baseline samples, the mutation detection rate was 80.8% [20/26], with the remaining 6 patients found to be wild-type for genes included in the panel. The most frequently mutated gene in the cohort was TP53, with a mutation rate of 73.0% [19/26]. Among the classic lung cancer driver mutations, PIK3CA (n=3), EGFR amplification (n=2), EGFR exon 19 deletion (n=1), KRAS Q61R (n=1), and MET amplification (n=1) were identified from our cohort. The mutation profile of the cohort is illustrated in Figure S1.

Correlation between maxAF and survival outcome

The maxAF represented the somatic mutation detected from the patient, and was defined as the highest allelic fraction observed among all the mutations detected from their blood sample, regardless of gene or mutation site. The median maxAF of the cohort was 2.06%, ranging between 0 and 52.86%. Among the 26 evaluable patients, 65.5% [17/26] of the patients had a maxAF of ≤5%, while the remaining 9 patients had a maxAF >5%. Patients with a maxAF >5% had significantly shorter OS compared to those with a maxAF ≤5% (P=0.020, ). The median OS for patients with a maxAF >5% was 7.1 months, while patients with a maxAF ≤5% had a median OS of 54.6 months.

Figure 1

Kaplan-Meier estimation of the overall survival (OS; expressed in days) of our cohort according to their maximum allelic fraction (maxAF). Red line represents patients with maxAF >5% (n=9), blue line represents patients with maxAF ≤5% (n=17). Tick marks at particular time points represent censored patients. The risk table below illustrates the number of patients included per time point.

Cox regression analysis further revealed an association between OS outcomes and maxAF [hazard ratio (HR) =54.12; 95% confidence interval (CI): 1.89–1,547.00; P=0.020], suggesting that ctDNA somatic mutations, represented by maxAF, hold potential as prognostic biomarkers for LUSC.

Correlation between methylation features and survival outcomes

We further analyzed the methylation profile of our cohort at the time of diagnosis to understand their epigenetic status at baseline. A strong correlation was found between maxAF and differential methylation patterns expressed as MPV1 (Pearson’s r=0.934, P<0.001; ). Consistently, Cox regression analysis also revealed an association between OS and differential methylation patterns (HR =2.377; 95% CI: 1.283–4.405; P=0.006), indicating that the higher the methylation levels, the worse the survival outcomes.

Figure 2

Methylation pattern value (MPV1) scores, reflecting the circulating tumor DNA (ctDNA) methylation pattern, were significantly correlated with maximum allelic fraction (maxAF). Dot plot illustrating the correlation. The gray area in the plot indicates the 95% upper and lower confidence intervals (CI).

Interestingly, in-house methylation data derived from 251 LUSC patient tissue samples at various stages (stage I–IV) revealed that a total of 402 methylation blocks correlated with OS. When considering only these 402 methylation blocks, the ctDNA methylation patterns of our cohort at baseline, expressed as MPV1, also demonstrated a strong correlation with OS (P=0.027, adjusted with maxAF, smoking, and stage).

These data suggest that ctDNA methylation status at the time of diagnosis could serve as a potential prognostic biomarker for LUSC.

Identifying maxAF-related methylation signatures

We further used LASSO to identify certain methylation blocks that were associated with maxAF and OS. A total of 1,956 methylation blocks were found to be differentially methylated from 21 samples with maxAF >0%, and were significantly correlated with maxAF (P<0.05). Using 7-fold cross-validation, the lambda.min and lambda.1se of the generalized linear model were used to determine the optimal methylation signature that was correlated with maxAF. illustrates the lambda and the corresponding mean-squared error (MSE) for various methylation block features generated from the LASSO model.

Figure 3

Cross-validated mean-squared error (MSE) generated by the least absolute shrinkage and selection operator (LASSO) model to predict maximum allelic fraction (maxAF). X-axis denotes the log of the lambda. Y-axis denotes the MSE. The values on the top of the graph indicate the number of methylation block features. Dotted line represents the lambda.min which reflects the minimum average cross-validated error, and lambda.1se which reflects the most normalized model whereby the error is within 1 standard error of the minimum.

When the lambda of the model equalled lambda.min yielding the minimum average cross-validated error of 3.5E-03, the model which included five methylation blocks had an MSE of 4.2E-03 (95% CI: 2.7E-03, 5.7E-3). When the lambda equalled lambda.1se, yielding the most normalized model whereby the error within 1 standard error of the minimum was 2.9E-02, the model which included 20 methylation blocks had an MSE of 5.6E-03 (95% CI: 2.7E-03, 8.5E-3). Using the model constructed with five methylation blocks on all of the 26 patient samples from the cohort including those with maxAF=0, the predicted methylation signature value was significantly associated with OS (P=0.015), and the corresponding HR was the largest compared to models constructed with 11, 15, and 20 methylation blocks (). lists the genomic position of the five differentially methylated blocks associated with prognosis. Table S2 lists the genomic position of all the methylation blocks identified by LASSO.

Table 1

Maximum allelic fraction (maxAF)-related methylation blocks associated with overall survival (OS)

Methylation block number	Hazard ratio	95% confidence intervals	P value
5	183.20	(2.74, 12,243.32)	0.015
11	102.22	(2.41, 4,338.23)	0.016
15	73.63	(2.21, 2,458.37)	0.016
20	62.23	(2.07, 1,931.42)	0.018

Table 2

The genes or regions covered by the 5 differentially methylated blocks identified by least absolute shrinkage and selection operator (LASSO) modeling to be associated with prognosis in patients with lung squamous cell carcinoma

Chromosome location	Covered regions/genes
1: 50886774-50887060	Exonic region of DMRTA2
1: 91182138-91182661	Promoter region of BARHL2
2: 223162695-223163049	Intronic region of PAX3 to promoter region of CCDC140
10: 101287805-101288002	Exonic region of LINC01475
10: 101293828-101293964	Promoter region of LINC01475 to intronic region of NKX2-3

Interestingly, MPV1 corresponding to the five methylation blocks was higher among the patients with maxAF >5% compared to those with maxAF ≤5% (). Table S3 summarizes the clinicopathological features, maxAF, and MPV1 scores for the 5 methylation blocks of the cohort. Patients with higher MPV1 and maxAF >5% had shorter OS, indicating that these differential methylation signatures from the five methylation blocks could serve as an alternative to maxAF, and potentially serve as prognostic biomarkers for LUSC.

Figure 4

Five maximum allelic fraction (maxAF)-related methylation blocks as potential prognostic markers of LUSC. Heat map illustrates the methylation patterns, expressed as methylation pattern value (MPV1), corresponding to the 5 methylation blocks identified by LASSO to be associated with maxAF. Blue denotes low MPV1 (MPV1 <0.1), red denotes high MPV1 (MPV1 >0.3) in the corresponding methylation block. X-axis denotes the genomic positions of the 5 methylation blocks. Y-axis denotes the maxAF of the 26 patients grouped according to the maxAF. Light green represents patients with maxAF ≤5%. Dark green represents patients with maxAF >5%. LUSC, lung squamous cell carcinoma.

Discussion

In addition to the clinical value of somatic mutation profiling in guiding therapeutic management for patients with lung cancer, mounting evidence supports the utility of DNA methylation signatures as diagnostic, prognostic, or predictive biomarkers (29-34). Compared to invasive tissue sampling, blood samples can be obtained in a minimally-invasive manner and sampled repeatedly during the treatment course. CtDNA, derived from DNA fragments shed into the plasma via necrosis, apoptosis, or active release of DNA from tumor cells, can reflect dynamic changes in the genetic and epigenetic status of the tumor without the intra-tumor heterogeneity associated with tissue biopsies (31-34). The concentration of ctDNA from the plasma-derived cfDNA is highly dependent on the tumor burden. Therefore, the detection of somatic mutations from ctDNA is severely limited in early-stage tumors (44). As compared with somatic mutations, DNA methylation in promoters of various tumor suppressor genes is more ubiquitous and occurs at the early stages of lung tumor development (18,20-27), which improves the sensitivity of detecting aberrant DNA methylation from cfDNA. Numerous reports have identified the promoter region of a single gene or a panel of genes that are differentially methylated in NSCLC; however, histology-based subcohort analysis on LUSC is limited (30,33,35-42). As compared to lung adenocarcinoma, LUSC has shown a histologically distinct methylation pattern in certain gene promoters (40); however, no study has yet reported on specific prognostic markers in LUSC.

In this study, we performed parallel somatic mutation and methylation profiling of ctDNA obtained from Chinese patients with various stages of LUSC to identify the molecular factors that can predict survival outcomes at the time of diagnosis. Our data demonstrated that maxAF >5%, derived from ctDNA somatic mutation profiling at baseline, was significantly correlated with shorter OS. Furthermore, we identified five genomic regions or genes that were differentially methylated at baseline and were significantly correlated with OS, suggesting their potential to predict prognosis even at the time of diagnosis. To the best of our knowledge, this study is the first to identify methylation signatures from the ctDNA of patients with LUSC that are significantly correlated with their prognosis. These methylation signatures are promising biomarkers for predicting survival outcome at the time of diagnosis for this subset of patients. Our findings highlight the potential clinical utility of ctDNA somatic mutation and methylation profiling at baseline for predicting the survival outcomes of patients diagnosed with various stages of LUSC. DNA methylation profiling of plasma samples could provide a minimally-invasive alternative to tissue biopsy procedures in the prognostication of patients with LUSC.

The ctDNA somatic mutation profile of our cohort was consistent with the mutation profile of LUSC patients in previous studies (14-16). TP53 was the most commonly mutated gene from our cohort, with an overall mutation rate of 73%. Only 8 patients from our cohort had potentially actionable mutations with available targeted agents. Stratifying the patients into two groups based on maxAF at baseline revealed a strong correlation between maxAF and OS.

ctDNA methylation profiling of our cohort revealed differential methylation patterns in various genomic regions, with 1,956 differentially methylated blocks at baseline correlated with maxAF. Among the 5 differentially methylated regions associated with maxAF and OS, DMRTA2 was also identified to be differentially methylated in lung tumors compared with normal lung samples and benign tumors (25). Moreover, PAX3 (22) and BARHL2 (23) were also identified as hypermethylated from primary LUSC tumors. However, no study has associated the differential methylation of these genes or other genomic regions with the prognosis of patients with LUSC.

Our study is limited by the inclusion of a small cohort of patients recruited from a single center, which might have introduced sampling bias. A prospective cohort study with a larger sample size in a multi-center setting is required to further establish the potential of differential methylation signatures from these genomic regions as prognostic biomarkers of patients with LUSC. Experimental studies on the individual genes that are correlated with prognosis would also advance the understanding of their regulation and functional role. It would also be interesting to analyze longitudinal ctDNA samples from patients with LUSC during their treatment course to establish methylation patterns that could predict treatment response, development of resistance, and disease relapse.

Conclusions

Our study identified molecular features that show promise as prognostic biomarkers for patients with various stages of LUSC at the time of their diagnosis. Our study provides an incremental contribution towards understanding the clinical value of ctDNA methylation and somatic mutation profiling at baseline in predicting the survival outcomes of patients with LUSC.

The article’s supplementary files as