# EPIGENETIC BIOMARKER AND PERSONALIZED PRECISION MEDICINE

EDITED BY : Jiucun Wang, Dongyi He, Momiao Xiong and Yun Liu PUBLISHED IN : Frontiers in Genetics and Frontiers in Cell and Developmental Biology

#### Frontiers eBook Copyright Statement

The copyright in the text of individual articles in this eBook is the property of their respective authors or their respective institutions or funders. The copyright in graphics and images within each article may be subject to copyright of other parties. In both cases this is subject to a license granted to Frontiers. The compilation of articles constituting this eBook is the property of Frontiers.

Each article within this eBook, and the eBook itself, are published under the most recent version of the Creative Commons CC-BY licence. The version current at the date of publication of this eBook is CC-BY 4.0. If the CC-BY licence is updated, the licence granted by Frontiers is automatically updated to the new version.

When exercising any right under the CC-BY licence, Frontiers must be attributed as the original publisher of the article or eBook, as applicable.

Authors have the responsibility of ensuring that any graphics or other materials which are the property of others may be included in the CC-BY licence, but this should be checked before relying on the CC-BY licence to reproduce those materials. Any copyright notices relating to those materials must be complied with.

Copyright and source acknowledgement notices may not be removed and must be displayed in any copy, derivative work or partial copy which includes the elements in question.

All copyright, and all rights therein, are protected by national and international copyright laws. The above represents a summary only. For further information please read Frontiers' Conditions for Website Use and Copyright Statement, and the applicable CC-BY licence.

ISSN 1664-8714 ISBN 978-2-88966-184-8 DOI 10.3389/978-2-88966-184-8

#### About Frontiers

Frontiers is more than just an open-access publisher of scholarly articles: it is a pioneering approach to the world of academia, radically improving the way scholarly research is managed. The grand vision of Frontiers is a world where all people have an equal opportunity to seek, share and generate knowledge. Frontiers provides immediate and permanent online open access to all its publications, but this alone is not enough to realize our grand goals.

#### Frontiers Journal Series

The Frontiers Journal Series is a multi-tier and interdisciplinary set of open-access, online journals, promising a paradigm shift from the current review, selection and dissemination processes in academic publishing. All Frontiers journals are driven by researchers for researchers; therefore, they constitute a service to the scholarly community. At the same time, the Frontiers Journal Series operates on a revolutionary invention, the tiered publishing system, initially addressing specific communities of scholars, and gradually climbing up to broader public understanding, thus serving the interests of the lay society, too.

#### Dedication to Quality

Each Frontiers article is a landmark of the highest quality, thanks to genuinely collaborative interactions between authors and review editors, who include some of the world's best academicians. Research must be certified by peers before entering a stream of knowledge that may eventually reach the public - and shape society; therefore, Frontiers only applies the most rigorous and unbiased reviews.

Frontiers revolutionizes research publishing by freely delivering the most outstanding research, evaluated with no bias from both the academic and social point of view. By applying the most advanced information technologies, Frontiers is catapulting scholarly publishing into a new generation.

#### What are Frontiers Research Topics?

Frontiers Research Topics are very popular trademarks of the Frontiers Journals Series: they are collections of at least ten articles, all centered on a particular subject. With their unique mix of varied contributions from Original Research to Review Articles, Frontiers Research Topics unify the most influential researchers, the latest key findings and historical advances in a hot research area! Find out more on how to host your own Frontiers Research Topic or contribute to one as an author by contacting the Frontiers Editorial Office<EMAIL_ADDRESS>
# EPIGENETIC BIOMARKER AND PERSONALIZED PRECISION MEDICINE

Topic Editors: Jiucun Wang, Fudan University, China Dongyi He, Shanghai Guanghua Rheumatology Hospital, China Momiao Xiong, University of Texas Health Science Center, United States Yun Liu, Fudan University, China

Citation: Wang, J., He, D., Xiong, M., Liu, Y., eds. (2020). Epigenetic Biomarker and Personalized Precision Medicine. Lausanne: Frontiers Media SA. doi: 10.3389/978-2-88966-184-8

# Table of Contents

*06 Prognostic and Predictive Value of Three DNA Methylation Signatures in Lung Adenocarcinoma*

Yanfang Wang, Haowen Deng, Shan Xin, Kai Zhang, Run Shi and Xuanwen Bao

Bruno Ramos-Molina, Lidia Sánchez-Alcoholado, Amanda Cabrera-Mulero, Raul Lopez-Dominguez, Pedro Carmona-Saez, Eduardo Garcia-Fuentes, Isabel Moreno-Indias and Francisco J. Tinahones

*91 DNA Methylation Biomarkers Predict Objective Responses to PD-1/PD-L1 Inhibition Blockade*

Gang Xue, Ze-Jia Cui, Xiong-Hui Zhou, Yue-Xing Zhu, Ying Chen, Feng-Ji Liang, Da-Nian Tang, Bing-Yang Huang, Hong-Yu Zhang, Zhi-Huang Hu, Xi-Yu Yuan and Jianghui Xiong

*103 A Hybrid Ensemble Approach for Identifying Robust Differentially Methylated Loci in Pan-Cancers*

Qi Tian, Jianxiao Zou, Yuan Fang, Zhongli Yu, Jianxiong Tang, Ying Song and Shicai Fan

*115 Maternal Smoking During Pregnancy Induces Persistent Epigenetic Changes Into Adolescence, Independent of Postnatal Smoke Exposure and is Associated With Cardiometabolic Risk*

Sebastian Rauschert, Phillip E. Melton, Graham Burdge, Jeffrey M. Craig, Keith M. Godfrey, Joanna D. Holbrook, Karen Lillycrop, Trevor A. Mori, Lawrence J. Beilin, Wendy H. Oddy, Craig Pennell and Rae-Chi Huang

*130 Predictive and Prognostic Value of Selected MicroRNAs in Luminal Breast Cancer*

Maria Amorim, João Lobo, Mário Fontes-Sousa, Helena Estevão-Pereira, Sofia Salta, Paula Lopes, Nuno Coimbra, Luís Antunes, Susana Palma de Sousa, Rui Henrique and Carmen Jerónimo

*146 Epigenetic Biomarkers in the Management of Ovarian Cancer: Current Prospectives*

Alka Singh, Sameer Gupta and Manisha Sachan

Dong-Mei Wu, Zheng-Kun Zhou, Shao-Hua Fan, Zi-Hui Zheng, Xin Wen, Xin-Rui Han, Shan Wang, Yong-Jian Wang, Zi-Feng Zhang, Qun Shan, Meng-Qiu Li, Bin Hu, Jun Lu, Gui-Quan Chen, Xiao-Wu Hong and Yuan-Lin Zheng

	- Jacob Peedicayil

Marcelo L. Ribeiro, Diana Reyes-Garau, Marc Armengol, Miranda Fernández-Serrano and Gaël Roué

*278 Association of Sperm Methylation at* LINE-1, *Four Candidate Genes, and Nicotine/Alcohol Exposure With the Risk of Infertility*

Wenjing Zhang, Min Li, Feng Sun, Xuting Xu, Zhaofeng Zhang, Junwei Liu, Xiaowei Sun, Aiping Zhang, Yupei Shen, Jianhua Xu, Maohua Miao, Bin Wu, Yao Yuan, Xianliang Huang, Huijuan Shi and Jing Du

*287 Perspectives on miRNAs as Epigenetic Markers in Osteoporosis and Bone Fracture Risk: A Step Forward in Personalized Diagnosis*

Michela Bottani, Giuseppe Banfi and Giovanni Lombardi

*312 From Genetics to Epigenetics, Roles of Epigenetics in Inflammatory Bowel Disease*

Zhen Zeng, Arjudeb Mukherjee and Hu Zhang

Victor G. Martinez, Ester Munera-Maravilla, Alejandra Bernardini, Carolina Rubio, Cristian Suarez-Cabrera, Cristina Segovia, Iris Lodewijk, Marta Dueñas, Mónica Martínez-Fernández and Jesus Maria Paramio

Clara Snijders, Julian Krauskopf, Ehsan Pishva, Lars Eijssen, Barbie Machiels, Jos Kleinjans, Gunter Kenis, Daniel van den Hove, Myeong Ok Kim, Marco P. M. Boks, Christiaan H. Vinkers, Eric Vermetten, Elbert Geuze, Bart P. F. Rutten and Laurence de Nijs

*427 E-Cadherin Downregulation is Mediated by Promoter Methylation in Canine Prostate Cancer*

Carlos Eduardo Fonseca-Alves, Priscila Emiko Kobayashi, Antonio Fernando Leis-Filho, Patricia de Faria Lainetti, Valeria Grieco, Hellen Kuasne, Silvia Regina Rogatto and Renee Laufer-Amorim

Xin Feng, Xubing Hao, Ruoyao Shi, Zhiqiang Xia, Lan Huang, Qiong Yu and Fengfeng Zhou

*471 New Analysis Framework Incorporating Mixed Mutual Information and Scalable Bayesian Networks for Multimodal High Dimensional Genomic and Epigenomic Cancer Data*

Xichun Wang, Sergio Branciamore, Grigoriy Gogoshin, Shuyu Ding and Andrei S. Rodin

Institute of Molecular Toxicology and Pharmacology,

Institute of Radiation

# Prognostic and Predictive Value of Three DNA Methylation Signatures in Lung Adenocarcinoma

Yanfang Wang1†, Haowen Deng2†, Shan Xin1,3, Kai Zhang4 , Run Shi 1 \* and Xuanwen Bao5,6 \*

1 Ludwig-Maximilians-Universität München, Munich, Germany, 2 Chair for Computer Aided Medical Procedures and

Cardiology, Sir Run Run Shaw Hospital, Zhejiang University School of Medicine, Hangzhou, China, 5

Helmholtz Center Munich, German Research Center for Environmental Health, Neuherberg, Germany, 4 Department of

Augmented Reality, Technical University Munich, Munich, Germany, 3

#### Edited by: Jiucun Wang,

Fudan University, China

#### Reviewed by:

Mariana Brait, Johns Hopkins University, United States David D. Eisenstat, University of Alberta, Canada Nejat Dalay, Istanbul University, Turkey Yi Wang, Fudan University, China

#### \*Correspondence:

Run Shi<EMAIL_ADDRESS>Xuanwen Bao<EMAIL_ADDRESS>
> †These authors have contributed equally to this work

#### Specialty section:

This article was submitted to Epigenomics and Epigenetics, a section of the journal Frontiers in Genetics

Received: 20 November 2018 Accepted: 01 April 2019 Published: 24 April 2019

#### Citation:

Wang Y, Deng H, Xin S, Zhang K, Shi R and Bao X (2019) Prognostic and Predictive Value of Three DNA Methylation Signatures in Lung Adenocarcinoma. Front. Genet. 10:349. doi: 10.3389/fgene.2019.00349 Biology, Helmholtz Center Munich, German Research Center for Environmental Health, Neuherberg, Germany, 6 Technical University Munich (TUM), Munich, Germany Background: Lung adenocarcinoma (LUAD) is the leading cause of cancer-related

mortality worldwide. Molecular characterization-based methods hold great promise for improving the diagnostic accuracy and for predicting treatment response. The DNA methylation patterns of LUAD display a great potential as a specific biomarker that will complement invasive biopsy, thus improving early detection.

Method: In this study, based on the whole-genome methylation datasets from The Cancer Genome Atlas (TCGA) and several machine learning methods, we evaluated the possibility of DNA methylation signatures for identifying lymph node metastasis of LUAD, differentiating between tumor tissue and normal tissue, and predicting the overall survival (OS) of LUAD patients. Using the regularized logistic regression, we built a classifier based on the 3616 CpG sites to identify the lymph node metastasis of LUAD. Furthermore, a classifier based on 14 CpG sites was established to differentiate between tumor and normal tissues. Using the Least Absolute Shrinkage and Selection Operator (LASSO) Cox regression, we built a 16-CpG-based model to predict the OS of LUAD patients.

Results: With the aid of 3616-CpG-based classifier, we were able to identify the lymph node metastatic status of patients directly by the methylation signature from the primary tumor tissues. The 14-CpG-based classifier could differentiate between tumor and normal tissues. The area under the receiver operating characteristic (ROC) curve (AUC) for both classifiers achieved values close to 1, demonstrating the robust classifier effect. The 16-CpG-based model showed independent prognostic value in LUAD patients.

Interpretation: These findings will not only facilitate future treatment decisions based on the DNA methylation signatures but also enable additional investigations into the utilization of LUAD DNA methylation pattern by different machine learning methods.

Keywords: LUAD, DNA methylation, regularized logistic regression, recursive feature elimination, LASSO Cox regression, metastasis

**6**

#### INTRODUCTION

Lung cancer is the leading cause of cancer-related mortality globally, causing over a million deaths a year (Genome Atlas Research Network., 2014; Jemal et al., 2018). There are two clinical types, one is the aggressive subtype small cell lung cancer and the other is non-small cell lung cancer (Hankey et al., 1999). Non-small cell lung cancer is histologically classified into four major subtypes by pathological and molecular characteristics: adenocarcinoma, large cell lung cancer, squamous cell lung cancer, and other types (Ettinger et al., 2010). Adenocarcinoma is the most common histological subtype of non-small cell lung cancer. Tobacco smoking is the major cause of lung adenocarcinoma (Toh et al., 2006). However, with the decrease in the number of smokers in many countries, the occurrence of LUAD in non-smokers has increased (Genome Atlas Research Network., 2014).

An accurate diagnosis of LUAD is one precondition to achieve a better treatment effect. Although the Mayo Clinic stage, size, grade, and necrosis (SSIGN) score, as well as the University of California Integrated Staging System can help improve the accuracy of the prognosis (Travis et al., 2011), the outcomes of patients with similar clinical characteristics or integrated systems scores still differ. Molecular characteristics may provide an indication for predicting the LUAD prognosis and response to therapy, thus offering great potential for improving individual treatment. Moreover, molecular characterization-based methods do not generally require bulk tissue samples, which can improve the patients' tolerance and reduce unnecessary operation steps. Among all the molecular characteristics, DNA methylation of CpG sites plays a crucial role in epigenetic regulation by reducing the activity of a DNA segment and repressing gene transcription (Jones, 2012; Du et al., 2015; Schübeler, 2015). DNA methylation is associated with carcinogenesis by repressing the expression of the tumor suppressor gene and promoting the expression of oncogenes (Herman et al., 1995; Schübeler, 2015; Vizoso et al., 2015; Klutstein et al., 2016). Hence, the cancer tissues have a distinct DNA methylation pattern compared to normal tissues. More importantly, unlike somatic genetic mutations in tumor tissues, DNA methylation patterns are inherently reversible changes and can therefore be promising targets for drug treatments (Ramchandani et al., 1999). Using DNA methylation signatures can help us make a better prognosis and predict the treatment response, thus prolonging the patients' survival.

Machine learning is a novel method to learn concept from data, which will help researchers discover the hidden insights. Based on DNA methylation patterns, machine learning techniques are developed and used to design models for precise classification and accurate prediction in medicine. In this study, we evaluated the possibility of DNA methylation signatures in identifying LUAD lymph node metastasis, differentiating between tumor tissue and normal tissue and predicting the OS of LUAD patients by applying TCGA whole-genome methylation datasets to several machine learning methods. Our results showed robust classifier effects with the AUC of both classifiers achieving values close to one for identifying lymph node metastasis and differentiating between tumor tissue and normal tissue. Cross-validation was applied to prevent overfitting. The LASSO Cox regression model was used to evaluate the patients' OS. Risk scores from the LASSO Cox model were combined with other clinicopathological risk factors to generate a nomogram to predict the prognosis and help the doctors to manage LUAD patients.

### METHODS

#### Data Source

The DNA methylation files and patients' information were obtained from Xena (https://xenabrowser.net/). Complete clinical, molecular, and histopathological data-sets are available at the TCGA website (https://portal.gdc.cancer.gov/).

### Feature Selection for DNA CpG Sites

We formulate critical methylation identification as a feature selection problem. Each CpG site is treated as a feature here and our goal is to find out which features are important for different tasks.

#### Variance Based Filtering

Variance is the squared deviation of the data from its mean, showing the spread of numbers. It is an important characteristic that reflects the distribution and discriminability of a feature. The variance σ of an observed sample sequence of a given feature {x1, x2, . . . , xi , . . . , xN} is computed by averaging across the squared difference of each value to the mean µ.

$$
\sigma^2 = \frac{\sum \left(\mathbf{x}\_i - \boldsymbol{\mu}\right)^2}{N}
$$

$$
\boldsymbol{\mu} = \frac{\sum \mathbf{x}\_i}{N}
$$

In general, a larger variance σ means a wider distributed and more separable feature space, which facilitates training a classifier to find class boundaries. On the other hand, variance σ is positively correlated with information entropy E, meaning that more information could be obtained with a larger variance σ. When σ is small, all the data is compressed and provides insufficient information for a classifier, so that we would avoid features with a small variance by setting a minimum threshold to filter out the indiscriminate features.

$$E = -\int P\left(\mathbf{x}\right) \log P\left(\mathbf{x}\right) d\mathbf{x}$$

#### Regularized Logistic Regression Model

Logistic regression is a widely applied and useful statistical, nonlinear model for predicting a binarized outcome based on a sequence of independent features. Assuming we have a general linear regression model y, which satisfies

$$\nu = \sum\_{i=0}^{N} \beta\_i \varkappa\_i$$

where xi stands for the i-th feature and βi is the correspondent coefficient. Since there is no constraint on the range of βi and xi , there is no maximum or minimum limit for y, i.e., y ∈ [−∞,∞]. Consider a standard logistic function

$$f\left(\mathbf{x}\right) = \frac{1}{1 + e^{-\mathbf{x}}}$$

This could map an input space from an infinite [−∞,∞] to a finite [0, 1]. By combining the linear regression model with the logistic function, we obtain the logistic regression model

$$\wp = \frac{1}{1 + e^{-\sum\_{i=1}^{N} \beta\_i x\_i}}$$

By thresholding y with threshold t, we obtain the binarized output.

$$o = \begin{cases} 1, y \ge t \\ 0, y < t \end{cases}$$

Note that, regularization term can be compounded with a logistic regression model, to force the learned coefficients to be sparser and more resistant to overfitting, which is highly beneficial for feature selection as well. We term the logistic regression model with regularization term as "Regularized logistic regression model."

#### Recursive Feature Elimination

Recursive feature elimination (RFE) adopts a brute-force and recursive way of undermining important features. Given a predefined model, which weighs all the features internally, RFE recursively uses the set of features to train the model and discard features that are the least important for the model (e.g., small weights) and repeats the training with the remaining features. This operation keeps recycling until certain expectations are reached, such as the maximum number of expected features Nexp. The process is described in **Algorithm 1**.

**Algorithm 1:** Recursive Feature Elimination.

**INPUT**: a set of features S = {f1, f2, . . . , fi , . . . , fn}, expected feature number Nexp **OUTPUT**: a set of kept features

Sbest = {fs1 , fs2 , . . . , fsi , . . . , fsn } **WHILE** size (S) > Nexp **DO**

#### **END WHILE**

Keep the final set of features as the set of most important features Sbest = S

#### Cox Regression

Cox regression, also called Proportional Hazards Regression, is a survival analysis model. It can be used to analyze relationships between different features and the survival time. The Cox model is based on the proportional hazards condition, which assumes that features have a proportional relationship to the exponential change of hazard. Thus, the model is formulated as a multiplication of a baseline hazard function with a sole time variable t, and an exponential function of the linear combination of all of the features as an input. Given a set of n samples {(**X**i , Yi ,si) | 0 ≤ i ≤ n, i ∈ **R**}, where Xi = (xi0, xi1, . . . , xik) and stands for the i-th sample of all the k features, Yi is the observation time and si is the survival status, the hazard function is

$$H\_i^{\;}\left(t\right) = H\_0\left(t\right)e^{X\_i^T \beta}$$

β= (β**0**,β1, . . . , βk ) is the coefficient vector weighing the contribution of the features. The partial likelihood of all the samples is

$$\begin{aligned} L\left(\mathcal{B}\right) &= \prod\_{i=1}^n L\_i\left(\mathcal{B}\right) \\ &= \prod\_{i=1}^n \frac{H\_i\left(Y\_i \mid X\_i\right)}{\sum\_{j:\,Y\_j \ge \, Y\_i} H\_i\left(Y\_i \mid X\_j\right)} \\ &= \prod\_{i=1}^n \frac{e^{X\_i^T \mathcal{B}}}{\sum\_{j:\,Y\_j \ge \, Y\_i} e^{X\_j^T \mathcal{B}}} \end{aligned}$$

By penalizing -log (L (β)), the optimal β could be uncovered.

#### LASSO Regularization

LASSO (Least Absolute Shrinkage and Selection Operator) is an important regularization in many regression analysis methods. The concept behind LASSO is that an L1-norm is used to penalize the weight of the model parameters. Assuming a model has a set of parameters {w0, w1, . . . , wn}, the LASSO regularization can be written as

$$\lambda \cdot \sum\_{i=0}^{n} \|\boldsymbol{w}\_i\|\_1$$

It can be also expressed as a constraint to the targeted objective function

$$\sum \left\| Y - Y^\* \right\|\_2 \text{, s.t. } \left\| \boldsymbol{w}\_i \right\|\_1 < t$$

An important property of the LASSO regularization term is that it can force the parameter values to be 0, thus generating a sparse parameter space, which is a desirable character for feature selection.

#### Workflow of the Coding Process

When it came to selecting the methylation features for the metastasis and tumor identification problems, we first used variance-based filtering to eliminate some of the least important CpG sites, and to decrease the computation for the following Regularized Logistic Regression Model and RFE. To avoid model overfitting and bias in the feature selection, cross validation was used in the following stages. The dataset was evenly divided into 5-folds, and further feature selections were conducted by applying Logistic Regression Model and RFE, following the standard pipeline of cross validation.

When predicting the OS of LUAD patients, we built the Cox proportional hazard regression model with LASSO regularization. 5-fold cross-validation was applied to avoid the overfitting. We plotted the plots in R software (R Foundation for Statistical Computing, Vienna, Austria. Version 3.4.3) and Python (Python Software Foundation. Python Language Reference, version 3.7).

### RESULTS

### Preparation of LUAD DNA Methylation Datasets

LUAD DNA methylation data and corresponding clinical data were downloaded from Xena (https://xenabrowser.net/) (Cline et al., 2013). After removing samples without a survival status and normalization, a total of 478 samples were analyzed in the present study (**Supplementary File 1**). The datasets included 409 samples for the recognition of metastasis, 428 samples for the recognition of tumor from normal tissue, and 446 samples for the prediction of OS [(**Supplementary Files 2**–(**4**].

### Identification of 3616-CpG-Based Signature for the Recognition of Metastasis

Variance-based selection was applied to filter features (methylation CpG sites). Features with small variances tend to be less discriminative, so we filtered out features with a standard variance smaller than 0.01 and 135,094 methylation signatures were selected. Regularized logistic regression and cross-validation were then applied to weigh the importance of each feature. The 428 LUAD samples were randomly assigned to a test set or a validation set by the cross-validation method. In short, five rounds of cross-validation were performed using different partitions and the validation results were combined over five rounds to overcome overfitting. By varying the value of the coefficient threshold, we obtained a different number of features that could be kept. When we used those kept features to regress the linear Logistic model by 5-fold cross-validation, the mean accuracy trend was as follows (**Figure 1A**). The number of kept features with regard to the different values of coefficient thresholds was shown in **Figure 1B**. The best performance was achieved at the threshold value 0.05 with 6,198 features kept with a 5-fold cross-validation. Recursive feature elimination with the same cross-validation configuration was tested and the result indicated that the kept features were the optimum minimal set of all the features (**Figure 1C**). The value of kept methylation CpG sites was shown in **Figure 1D**. We assessed the accuracy of the 3616-CpG-based classifier for detecting metastasis with a ROC analysis (**Figure 1E**) with the same cross-validation configuration, and averaged the weights of selected features across different set as the final coefficients. Furthermore, the metastatic probability of each sample were calculated by the coefficients of kept methylation CpG sites (**Figure 1F** and **Supplementary File 5**). The AUC for the classifier achieved values close to 1in all of the 5-fold cross-validation, indicating the robust classifier effect. The tumor tissues in total dataset were divided into high metastatic risk score and low metastatic risk score groups, respectively, using 0.5 as the cutoff. The patients in the low metastatic risk score group have a longer OS than those in the high metastatic risk score group in the total datasets (p < 0.0001, **Figure 1G**) as well as in the separated 5-fold training and validation sets (p < 0.05, **Supplementary Figure 1**). We assessed the prognostic accuracy of the 3616-CpG-based classifier metastatic classifier with a time-dependent ROC analysis at varying follow-up times (500, 1,000, 1,500, 2,000, 2,500, 3,000 days) (**Supplementary Figure 2**). The accuracy was all around 66%, indicating that the 3616-CpG-based classifier for identifying metastasis could also work well for predicting the OS of LUAD patients.

### Identification of 14-CpG-Based Signature to Recognize Tumor and Normal Tissues

134,015 features were kept by variance thresholding (0.01). Regularized Logistic regression and cross-validation were applied to weigh the importance of each feature as mentioned above. An accuracy of 100% can easily be achieved for the number of features range from 14 to 43,246 (**Figure 2A**). The number of kept features with regard to the different thresholding values was shown in **Figure 2B**. Recursive feature elimination with crossvalidation was tested and the result indicated that an accuracy of 100% can be achieved when the kept feature numbers reached 14 (**Figure 2C**). 14 CpG sites were kept: cg25774643, cg03502002, cg14789818, cg23479922, cg04864807, cg07915921, cg20146541, cg08862830, cg01016533, cg19191888, cg08094098, cg01912692, cg10707110, cg24103195. The value of kept methylation CpG sites was shown in **Figure 2D**. We then calculated the probability of being tumor for each sample by the coefficients of kept methylation CpG sites (**Supplementary File 6** and **Figure 2E**) in the same tradition as of the recognition of metastasis. The accuracy of the 14-CpG-based classifier was assessed by means of ROC analysis (**Figure 2F**). The results showed that the accuracy reached 100% in all 5-fold cross-validation, indicating the high sensitivity and specificity of the 14-CpG-based classifier in differentiating between LUAD tumor tissues and corresponding normal tissues. Furthermore, we applied the 14-CpG-based classifier on an external dataset to confirm the accuracy of the 14-CgG-based classifier (**Figure 2G**). The AUC value was 98.4% for differentiating the tumor and normal tissues (**Figure 2H**). The analysis before showed the regularized logistic model we applied worked well in different datasets.

### Identification of 16-CpG-Based Signature to Predict the OS of LUAD Patients

We used a LASSO Cox regression to build a prognostic model, which selected 16 methylation CpG sites from the CpG sites identified by the DNA methylation 450 k chip: cg00161124, cg01105229, cg03923535, cg10976778, cg12141052, cg12240358,

cg13297560, cg14139311, cg14184729, cg18140857, cg19410791, cg20268054, cg23146197, cg25229048, cg26709300, cg27018309 (**Figures 3A,B**). The values of the 16 methylation CpG sites for each patient were shown in **Figure 3B**. A formula was derived to calculate the risk score for every patient based on their individual 16 methylation β values (**Supplementary File 7**). The risk scores of tumor samples were calculated by the coefficients of the kept methylation CpG sites (**Figure 3C**). The patients were divided into high risk score and low risk score groups, respectively, with a cutoff of −0·54. Kaplan-Meier survival analysis (**Figure 3D**) showed that the survival probability of patients in lower risk score was significantly better than in high risk score group (logrank test, all p < 0.0001). We assessed the prediction accuracy of the 16-CpG-based model by means of time-dependent ROC analysis at varying follow-up times. The AUC values for 500, 1,000, 1,500, 2,000, 2,500, and 3,000 days were 0.688, 0.681, 0.697, 0.685, 0.738, and 0.758, respectively, which confirmed the effectiveness of the 16-CpG-based model to predict the OS of LUAD patients (**Figure 3E**).

According to their clinicopathological conditions, like epidermal growth factor receptor (EGFR) mutation, K-ras or Ki-ras (KRAS) mutation, lymph node metastatic (LNM) condition, and AJCC stage, LUAD patients were divided up into several subgroups to validate the independent diagnostic value of the methylation signature. EGFR mutation showed a striking correlation with LUAD patient characteristics, which were correlated with the clinical treatment response and then affected the OS of LUAD patients. The Kaplan-Meier curves regarding EGFR mutation and wildtype groups were shown in **Figures 4A,B**. Patients with low risk scores generally had significantly better survival than those with high risk scores in both groups (p < 0.0001). Similarly, patients with low risk scores had a significantly longer OS than those with high risk scores in both KRAS mutation and wildtype subgroup and both LNM positive and negative groups (**Figures 4C–F**, p < 0.0001). For the patients in AJCC stage I and AJCC stage II-IV, the survival probability of patients with low risk scores was higher than those with low risk scores (**Figures 4G,H**). The stratification analysis above revealed that the 16-CpG-based model could effectively predict the OS of patients regardless of the patients' clinicopathological properties, and provide prognostic power to complement the clinical stage and SSIGN scores.

Lastly, the risk scores were applied to the Cox regression model with the clinicopathological risk factors to perform multivariable survival analysis, thereby generating a nomogram to predict patients' survival probability for 3 and 5 years (**Figure 5A**). In the multivariable survival analysis, we included age, gender, EGFR status, AJCC stage, and risk scores from 16- CpG-based model. The nomogram was further verified with calibration plots (**Figure 5B**). The results showed that the nomogram fared well with the ideal mode for 3 and 5 years, indicating the nomogram worked well in predicting the OS of LUAD patients. According to the risk scores from the nomogram, patients were divided into high risk and low risk group. Kaplan-Meier survival analysis showed that the survival probability of patients with low risk score was significantly higher than those with high risk score (**Figure 5C**, log-rank test, p < 0.0001). The prognostic accuracy of the nomogram was further accessed by time-dependent ROC curves (**Figure 5D**). The results showed that the AUC values were all around 0.7 at varying follow-up times (500, 1,000, 1,500, 2,000, 2,500, 3,000 days), indicating the high effectiveness of the nomogram in predicting the prognostic OS of LUAD patients.

### DISCUSSION

The present study demonstrates the potential of using DNA methylation signatures to identify the lymph node metastasis of primary LUAD tissues, to differentiate between the LUAD tumor and normal tissues, and to predict the OS of LUAD patients. Invasive biopsy is the gold standard for the validation of tumor tissues and identification of histological subtypes. However, the collection of bulk tissue samples for immunohistochemical (IHC) staining may cause secondary damage to patients. An inadequate tissue yield or quality also creates barriers for the histological diagnosis. Besides, it may be difficult to identify lymph node metastasis during operation. Nowadays, molecular characterization methods provide new insights in pathological diagnosis (Tsou et al., 2007; Selamat et al., 2012; Zhang et al., 2013; Ogino et al., 2016). Since the global change of DNA methylation takes place at the beginning of carcinogenesis, DNA methylation has been considered a promising biomarker for the early detection and diagnosis of cancers (Franco et al., 2008; Hatano et al., 2015; Wu and Ni, 2015), which can complement the pathological IHC staining. Moreover, DNA methylation analysis does not require bulk tissue samples. Small amounts of tissue are enough for DNA extraction and methylation-chip or methylation-seq analysis, which will reduce the patients suffering. Hundreds of thousands of the DNA methylation CpG sites can be identified through genome-wide DNA methylation detection by DNA methylation chips or methylation-seq. Discovering a potential panel of DNA methylation-based biomarkers from the large DNA methylation files can be beneficial for the early diagnosis of cancer initiation and metastasis. Several research studies have shown the potential of utilizing DNA methylation profiles to help the diagnosis of different cancers (Diaz-Lagares et al., 2016; Zhang et al., 2017; Sandanger et al., 2018). One study applied an unsupervised clustering method on DNA methylation profiles to find potential subtypes of childhood B-cell acute lymphoblastic leukemia. The patients were allocated into two subgroups by the unsupervised hierarchical clustering of DNA methylation profiles, which showed a significant association between DNA methylation and disease-free survival (Sandoval et al., 2013a). Another study also utilized a similar strategy to find the association between DNA methylation signatures and the recurrence-free survival in non-small-cell lung cancer samples (Sandoval et al., 2013b). In our study, we applied supervised learning strategy (regularized logistic regression) to find the prognostic CpG cites in LUAD primary tissues. The RFE helped to eliminate the unnecessary features in regression, which constrained the numbers of key CpG sites for prognosis. Besides, LASSO Cox regression was useful to reduce the feature numbers in the COX survival analysis. One study built a

with regard to the different values of coefficient thresholds. (C) Recursive feature elimination with a cross-validation test. (D) Unsupervised hierarchical clustering and heat map associated with the methylation profile (according to the color scale shown) to differentiate between LUAD tumor tissues and corresponding normal tissues. (E) The probability of being tumor for each sample calculated by the coefficients of methylation signatures. (F) ROC curves showing the high sensitivity and specificity in differentiating between LUAD tumor tissues and normal tissues. (G) The workflow of model construction, internal validation and external validation. (H) ROC curve showing the high sensitivity and specificity in differentiating between LUAD tumor tissues and corresponding normal tissues on an external dataset.

calculated by the coefficients of methylation signatures from Lasso Cox analysis. (D) Kaplan-Meier curves of LUAD patients with a low or high risk of death, according to risk scores from the 16-CpG-based classifier. (E) Time-dependent ROC analysis at varying follow-up times (500, 1,000, 1,500, 2,000, 2,500, 3,000 days). We used AUC values at 500, 1,000, 1,500, 2,000, 2,500, 3000 days to assess the prognostic accuracy.

prognostic signature by LASSO Cox regression to predict the progression-free survival of LUAD patients and demonstrated the potential biological significance of DNA methylation in the etiology of LUAD (Bjaanæs et al., 2016). Another study built a mortality risk score by LASSO Cox regression (Zhang et al., 2017). The signature based on ten selected CpG sites exhibited strong association with all-cause mortality. Moreover, one recent study used blood-derived DNA methylation and gene expression profiles to identify CpG lung cancer markers prior to diagnosis. They emphasized the difference of prognostic CpG sites in smoking and non-smoking lung cancer patients (Sandanger et al., 2018). In this study, based on the methylation

FIGURE 4 | Kaplan-Meier survival analysis for LUAD patients according to the 16-CpG-based classifier. Patients were classified according to clinicopathological risk factors. (A,B) EGFR status; (C,D) KRAS status; (E,F) lymph node metastatic (LNM) status; (G,H) AJCC stage I and II-IV. The patients were divided into low-risk and high-risk groups. P-values were calculated using the log-rank test.

and 5 years, respectively. (C) Kaplan-Meier survival analysis for the OS of LUAD patients according to the risk scores from the nomogram. (D) Time-dependent ROC curves from the nomogram for overall survival in 3 and 5 years.

profiles of LUAD patients, we performed regularized logistic regression and LASSO Cox regression to identify the lymph node metastasis, to differentiate between tumor and normal tissues and to predict the OS of LUAD patients. From the primary LUAD tumor tissues, 3616 methylation CpG sites were kept to build a classifier to identify LUAD lymph node metastasis. ROC curves showed the high sensitivity and specificity of the 3616-CpGbased classifier in identifying lymph node metastasis from CpG sites of primary tumor tissues. All the samples came from the primary tumor tissues, which means that the metastatic behavior can be identified even without extracting tissues from lymph nodes. Therefore, it would work as a biomarker to predict the diagnosis of lymph node metastasis. Since the metastatic behavior of LUAD affects the OS of LUAD patients dramatically, we applied the metastatic classifier to check whether the model can be used to predict the OS of LUAD patients. The time-dependent ROC curves showed the effectiveness of the metastatic classifier in predicting the OS of LUAD patients at varying follow-up times. As expected, the patients in the high metastatic risk score group have a significantly worse OS than those in the low metastatic risk score group.

 More than 70% of the tumor tissues are cancer cells. The heterogeneity of tumor tissues may influence the accuracy of the diagnosis. We compared the heterogeneous tumor tissues with the normal tissues, which also include the vascular epithelial cells and other cell types. Considering the heterogeneous tumor tissues and heterogeneous normal tissues as a whole for each, we tried to eliminate the influence brought about by heterogeneity (Li et al., 2014). In this study, to differentiate between tumor and normal tissues, we concluded that 14 CpG methylation sites were enough for the diagnosis. To check the overfitting potential, we applied 5-fold cross-validation. The efficiency of the model above was tested by a ROC curves in five different training and validation datasets, which showed the high efficiency and specificity of the 14-CpG-based classifier in differentiating between LUAD tumor tissues and the normal tissues. Furthermore, we also validated our regression model on the external dataset from one study (Bjaanæs et al., 2016). Results showed an AUC value of 98.4% to differentiate the tumor and normal tissues by the ROC analysis. The external dataset further confirmed the accuracy of the regularized logistic model which we applied to build the both classifiers above. From the two classifiers, we obtained an overlap cluster of CpG sites: cg03502002 and cg07915921. The information for cg07915921 is not clear. cg03502002 is on the CpG island of the promotor region of the GALR1 gene. The methylation status of the GALR1 promoter and the level of GALR1 gene expression have been correlated in a large number of head and neck squamous tumor specimens (Misawa et al., 2008). Ectopic expression of GALR1 suppresses tumor cell proliferation through Erk1/2-mediated regulation of cyclin-dependent kinase inhibitors and cyclin D1 (Kanazawa et al., 2009). One study revealed that hypermethylated GALR1 plays important roles in smoking-associated LUAD (Tan et al., 2013).

We also built a model to predict the OS of LUAD patients by means of methylation CpG sites. The LASSO Cox regression model generated risk score for each patient. When we assessed the survival status and distribution of risk scores, patients with low risk scores generally had a better OS than those with high risk scores. The model will help guide individualized follow-up schedules for LUAD patients. The high-risk patients have poor OS prediction. This could be the basis of a future clinical trial. The LASSO Cox regression results were further confirmed by the time-dependent ROC analysis. When we compared the timedependent ROC from the OS-prediction model and metastasisprediction classifier, the OS-prediction model turned out to be more precise in the long-term survival prediction while the metastasis-prediction classifier worked better in the short-term survival prediction. One explanation could be that when the LUAD patients were accompanied by lymph node metastasis, the tumor progressed and the patients had a poorer prognosis. The OS expectation of patients with lymph node metastasis was shorter than those without lymph node metastasis. Hence, the metastasis-prediction classifier would work better for the short-term prediction.

To further utilize the risk scores from the Cox regression model, we classified patients into several subgroups according to the clinicopathological risk factors (EGFR mutation, KRAS mutation, LNM status and AJCC stages). The 16-CpG-based classifier still showed clinical and statistical significance regardless of the clinicopathological status of LUAD patients.

The independent prognostic values of the 16-CpG-based model were validated by multivariable survival analysis, which integrated other clinicopathological risk factors for the OS of LUAD patients. The Cox regression risk scores were applied together with age, gender, EGFR status, AJCC stages as indicators to generate a nomogram to predict the 3- and 5-year survival probability. We verified the performance of the nomogram by calibration plots. The predicted OS of LUAD patients by the nomogram was highly consistent with the observed 3- and 5-year OS of LUAD patients. Log-rank test and time-dependent ROC curves at vary follow-up times further confirmed the nomogram. Thus, the nomogram could provide an accurate and simple prognostic prediction for LUAD patients.

In previous studies, mRNA expression files (Beer et al., 2002), the mutation of key genes (Takano et al., 2008; Kosaka et al., 2009), long no-coding RNA expressions (Kosaka et al., 2009; Huarte, 2015; Zhou et al., 2016), and histone modifications (Seligson et al., 2009; Zhou et al., 2016) showed the prognostic potential for different types of cancer. Here, we emphasized that the methylation patterns could also be a meaningful tool for the prognosis of LUAD patients. Some studies have identified that multiple CpG sites are differentially methylated in lung cancer compared to normal tissues (Genome Atlas Research Network., 2014; Poirier et al., 2015; Hao et al., 2017). The key for methylation pattern-based early diagnosis is the identification of crucial CpG sites in LUAD. The use of supervised machine learning methods allowed us to integrate all methylation CpG sites identified by the methylation chip into one model, which improved the prognostic accuracy over that of a single CpG site alone. Our findings show that three DNA CpG signature-based models can effectively identify lymph node metastasis by the CpG sites from primary tumor tissues, differentiate between tumor and normal tissues, and predict the OS of LUAD patients. The tissues would be collected by preoperative biopsy or at surgery. The classifiers for identifying lymph node metastasis and differentiation between tumor and normal tissues would help the preoperative diagnosis. The Lasso Cox model would be helpful for adjuvant treatment and prognostic planning. Therefore, the 3 methylation signatures could be of great value in assessing the status, predicting prognosis and achieving individualized treatments of LUAD patients.

The limitations of our study should be mentioned. The methylation 450 k chip did not identify as many CpG sites as the methylation 850 k chip or methylation sequencing. The methylation CpG site candidates identified here did not represent the complete CpG sites in the genome of LUAD patients.

In conclusion, we built three DNA CpG signature-based models to identify LUAD lymph node metastasis by the CpG sites from primary tumor tissues, differentiate between tumor tissue and normal tissue, and predict the OS of LUAD patients, which highlight the relationship between clinical results (metastasis, survival) and methylation biomarkers in LUAD patients. The nomogram comprising LASSO Cox risk scores and clinicopathological factors may help predict the OS of LUAD patients and help individualized treatment of LUAD patients.

#### AUTHOR CONTRIBUTIONS

XB and RS conceived and designed the experiments. XB and HD wrote the code (The Python code is available due to request). XB, YW, and SX wrote the paper. RS and KZ reviewed the manuscript. All authors read and approved the final manuscript.

#### ACKNOWLEDGMENTS

We greatly thank Dr. Michael Rosemann and Prof. Dr. Michael J. Atkinson for helpful discussions and suggestions. We sincerely thank Sino-German (CSC-DAAD) Postdoc Scholarship Program for supporting the research of YW, and greatly thank the China

#### REFERENCES

Scholarship Council (CSC) for supporting the research and work of SX, RS, and XB.

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene. 2019.00349/full#supplementary-material

Supplementary Figure 1 | Kaplan-Meier survival according to risk scores from the 3616-CpG-based classifier in the training, validation sets for 5-fold cross-validation. (A, C, E, G, I) The training sets 1-5. (B, D, F, H, J) The validation sets 1-5.

Supplementary Figure 2 | Time-dependent ROC analysis at varying follow-up times (500, 1,000, 1,500, 2,000, 2,500, 3,000 days) according to risk scores from the 3616-CpG-based classifier.

adenocarcinoma and integration with mRNA expression. Genome Res. 22, 1197–1211. doi: 10.1101/gr.132662.111

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Wang, Deng, Xin, Zhang, Shi and Bao. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Genome-Wide Analysis of Lung Adenocarcinoma Identifies Novel Prognostic Factors and a Prognostic Score

Donglai Chen1† , Yueqiang Song2† , Fuquan Zhang3† , Xiaofan Wang3 , Erjia Zhu1 , Xi Zhang2 , Gening Jiang1 , Siguang Li2 \*, Chang Chen1 \* and Yongbing Chen3 \*

 Department of Thoracic Surgery, Shanghai Pulmonary Hospital, Tongji University School of Medicine, Shanghai, China, Department of Regenerative Medicine, Stem Cell Center, Tongji University School of Medicine, Shanghai, China, Department of Thoracic Surgery, The Second Affiliated Hospital of Soochow University, Medical College of Soochow University, Suzhou, China

#### Edited by:

Dongyi He, Shanghai Guanghua Rheumatology Hospital, China

#### Reviewed by:

Steven G. Gray, St. James's Hospital, Ireland David D. Eisenstat, University of Alberta, Canada

#### \*Correspondence:

Siguang Li<EMAIL_ADDRESS>Chang Chen<EMAIL_ADDRESS>Yongbing Chen<EMAIL_ADDRESS>†These authors have contributed

equally to this work

#### Specialty section:

This article was submitted to Epigenomics and Epigenetics, a section of the journal Frontiers in Genetics

Received: 22 January 2019 Accepted: 06 May 2019 Published: 22 May 2019

#### Citation:

Chen D, Song Y, Zhang F, Wang X, Zhu E, Zhang X, Jiang G, Li S, Chen C and Chen Y (2019) Genome-Wide Analysis of Lung Adenocarcinoma Identifies Novel Prognostic Factors and a Prognostic Score. Front. Genet. 10:493. doi: 10.3389/fgene.2019.00493 Background and Objective: Lung adenocarcinoma (LUAD) is the most common histological type of all lung cancers and is associated with genetic and epigenetic aberrations. The tumor, node, and metastasis (TNM) stage is the most authoritative indicator of the clinical outcome in LUAD patients in current clinical practice. In this study, we attempted to identify novel genetic and epigenetic modifications and integrate them as a predictor of the prognosis for LUAD, to supplement the TNM stage with additional information.

Methods: A dataset of 445 patients with LUAD was obtained from The Cancer Genome Atlas database. Both genetic and epigenetic aberrations were screened for their prognostic impact on overall survival (OS). A prognostic score (PS) integrating all the candidate prognostic factors was then developed and its prognostic value validated.

Results: A total of two micro-RNAs, two mRNAs and two DNA methylation sites were identified as prognostic factors associated with OS. The low- and high-risk patient groups, divided by their PS level, showed significantly different OS (p < 0.001) and recurrence-free survival (RFS; p = 0.005). Patients in the early stages (stages I/II) and advanced stages (stages III/IV) of LUAD could be further subdivided by PS into four subgroups. PS remained efficient in stratifying patients into different OS (p < 0.001) and RFS (p = 0.005) when the low- and high-risk subgroups were in the early stages of the disease. However, there was only a significant difference in OS (p = 0.04) but not RFS (p = 0.2), between the low-risk and high-risk subgroups when both were in advanced stages.

Conclusion: PS, in combination with the TNM stage, provides additional precision in stratifying patients with significantly different OS and RFS prognoses. Further studies are warranted to assess the efficiency of PS and to explain the effects of the genetic and epigenetic aberrations observed in LUAD.

Keywords: lung adenocarcinoma, genome-wide, prognostic factor, survival prediction, TCGA

### INTRODUCTION

fgene-10-00493 May 21, 2019 Time: 18:26 # 2

Lung cancer is the leading cause of global cancer-related mortality, and ranks second in the estimated new cases of cancer in both sexes in the United States (Siegel et al., 2017). Lung adenocarcinoma (LUAD) is the most common histological type of lung cancer, accounting for approximately 50% (Shedden et al., 2008; Warth et al., 2012; Cancer Genome Atlas Research Network, 2014). Currently, the tumor, node, and metastasis (TNM) stage is the most accepted system for estimating the prognosis of patients with LUAD in clinical practice (Warth et al., 2012). However, prognoses of LUAD patients who share the same pathological stage vary considerably (Tsao et al., 2015; Zhang et al., 2015; Liang et al., 2017; Dalwadi et al., 2018). Therefore, a more accurate system is in demand to predict the outcomes of patients with LUAD that can add further valuable information to the TNM stage.

Aberrant genetic and epigenetic modifications of oncogenes and tumor suppressors contribute to the tumorigenesis and progression of LUAD (Khalil et al., 2018; Rowbotham et al., 2018; Tessema et al., 2018; Toyokawa et al., 2018; Wang et al., 2018). Genetic and epigenetic abnormalities have been associated with LUAD patient survival, especially the aberrant expression of cancer-related genes and DNA methylation at specific sites (Uruga et al., 2017; Zhang et al., 2017; Gonzalez-Vallinas et al., 2018; Wang et al., 2018). For instance, using genome-scale DNA methylation profiling, a study identified 164 hypermethylated genes and 57 hypomethylated genes involved in cell differentiation and the epithelial-tomesenchymal transition in LUAD (Selamat et al., 2012). Notably, DNA methylation also accounts for the alteration of gene expression in LUAD (Zhang et al., 2017; Gao et al., 2018; He et al., 2018), and may thus indirectly affect the biological behaviors and processes of LUAD. Specifically, He et al. (2018) identified an association between aberrant CpGmethylation and the prognostic value of the corresponding gene expression based on 1095 LUAD samples, and identified 10 aberrantly methylated and dysregulated genes with independent prognostic value.

In recent years, a class of small non-coding RNA molecules, called microRNA (miRNA), has been increasingly investigated (Fu et al., 2017; Greenawalt et al., 2018; Othman and Nagoor, 2019). miRNAs can regulate the expression of protein-coding genes by base pairing with the target mRNAs, inducing the degradation or translational repression of the bound mRNAs (Ha and Kim, 2014; Hou et al., 2018). The prognostic significance of miRNAs has also been investigated and confirmed in many studies (Zhang et al., 2017; Gonzalez-Vallinas et al., 2018; Xu et al., 2018). For instance, mir-486 was shown to be a miRNA that is differentially expressed in LUAD and potentially interacts with ITGA11, a cancer-promoting gene (Zhang et al., 2017). Gonzalez-Vallinas et al also reported a significant association between mir-539, mir-323b, and mir-487a upregulation and worse disease-free survival in non-smoker patients with LUAD (Gonzalez-Vallinas et al., 2018).

So far, many studies have established panels of prognostic factors that predict the outcomes of patients with LUAD, based on multiple lines of evidence. However, most studies were conducted without integration of the network constituted by dysregulations at different levels. Because LUAD represents a set of heterogeneous diseases in which aberrations can exist at genome and epigenome levels, we performed a genome-wide analysis, which should provide more comprehensive insight into survival prediction. Using the data of 445 patients from The Cancer Genome Atlas (TCGA) database, we identified prognostic value of two miRNAs, two mRNAs and two methylation sites. A prognostic score (PS) was developed by integrating these factors to stratify LUAD patients with different lengths of survival into subgroups. From our data, combining PS and the TNM stage achieved greater accuracy in predicting the prognoses of patients with LUAD, indicating that PS is a promising system for personalized and precise medicine.

### MATERIALS AND METHODS

### Data Extraction and Prepossessing

The genome-wide data for 706 LUAD patients was downloaded from TCGA database1 , including the expression levels of 20530 mRNAs, 2228 miRNAs and 485577 DNA methylation sites, together with the outcomes of 630 patients. The exclusion criteria were listed as follows: (1) Patients whose genomic or epigenomic information was absent; (2) Genes lacking information on either their transcript (mRNA or miRNA) or DNA methylation levels in more than half the LUAD samples; (3) Patients whose survival information was unavailable. Ultimately, a total of 445 LUAD patients were included in the study, together with 16928 mRNAs, 453 miRNAs and 395963 DNA methylation sites.

### Identification of Survival-Associated Transcripts and DNA Methylation Sites

A Cox regression model was used to evaluate the association of gene transcripts (mRNA or miRNA) and DNA methylation sites with lengths of overall survival (OS). A univariate Cox regression analysis was initially used, followed by the screening of included potential factors with a p ≤ 0.1 for further analysis.

Afterward, considering the remaining large numbers of gene transcripts and methylated sites, we performed a Lasso-Cox analysis to screen and shrink the data. We then used multivariate Cox regression to further analyze the association between the gene transcripts or DNA methylation sites with OS, while adjusting for other clinicopathological factors.

### Identifying and Screening Potential miRNA Targets

We retrieved the potential target genes of miRNAs that had already been shown to be significantly associated with OS from miRTarBase (the experimentally validated microRNAtarget interactions database, release 7.0) (Hou et al., 2018). Lasso-Cox regression was then used to screen the mRNAs of

1https://xenabrowser.net

genes which were identified as targets of miRNAs with highconfidence (p ≤ 0.05).

### Calculation of Spearman's Correlation Coefficients

The direction of association among the transcripts and methylation sites were calculated with Spearman's correlation in the 445 LUAD tissues. If an mRNA tended to increase when miRNA or methylation increased, the Spearman's correlation coefficient was positive. If an mRNA tended to decrease when miRNA or methylation increased, the Spearman's correlation coefficient was negative. We set a threshold of 0 with which to assess the candidate miRNAs, mRNAs and methylation sites. Any pair with a correlation coefficient value < 0 was considered to be negatively correlated, whereas any pair with a correlation coefficient value > 0 as positively correlated.

### External Validation of Identified Transcripts and Methylation Sites

We validated the prognostic value of our candidate transcripts in KM-Plot2 . The impact of the candidate methylation sites on survival was confirmed in MethSurv3 (Modhukur et al., 2018). We used Jetset to select the corresponding probe sets for the candidate mRNAs and miRNAs because a given gene may be detected by multiple probe sets, which may lead to inconsistent or even contradictory measurements (Li et al., 2011).

### Construction and Validation of PS

To further assess the predictive ability of the significant factors identified, we constructed a PS as an integrated predictor. PS was calculated as a weighted sum of the expression levels of the transcripts and DNA methylation sites present in a given sample (Hou et al., 2018). For specimen i the calculation formula for PS was shown as follows:

$$\text{PS} = \sum\_{i=1}^{n} \beta i \mathbf{x} i$$

The weight of each variable is represented by the Cox regression coefficient β, and the expression level is denoted by x. A greater value of PS indicates a worse prognosis.

We divided the patients into either the high-risk or low-risk group according to the median value of PS. Each group was subdivided into the early-stage (stages I–II) and advanced-stage (stages III–IV) subgroups based on the pathological stage. The Kaplan–Meier method and log-rank tests were used to assess the differences in OS and RFS between in the subgroups.

#### Statistical Analysis

All statistical analyses were performed with R version 3.4.4 (packages glmnet\_2.0-16, survival\_2.4-3; Institute for Statistics and Mathematics, Vienna, Austria). A two-tailed p < 0.05 was considered statistically significant.

#### RESULTS

### General Information on Patients With LUAD

The clinicopathological characteristics of the LUAD patients in our study are shown in **Table 1**. Of the 445 patients, 210 (47.2%)

TABLE 1 | Distributions of the demographic and clinical variables of 445 patients with lung adenocarcinoma patients.

Adc, adenocarcinoma; RUL, right upper lobe; RML, right middle lobe; RLL, right lower lobe; LUL, left upper lobe; LLL, left lower lobe.

2http://kmplot.com/analysis/index.php?p=background

3https://biit.cs.ut.ee/methsurv/

Adc, adenocarcinoma; HR, hazard ratio; CI, confidence interval; NA, not available.

were male and 235 (52.8%) were female. The median age was 66 years (ranging from 39 to 88). Patients with early-stage LUAD constituted the majority of our cohort. The primary tumor mainly was mainly located in the upper lobe on either side.

As shown in **Table 2**, we examined the association between each clinicopathological characteristic and OS. A univariate Cox regression analysis indicated that a higher TNM stage was significantly associated with poorer OS (**Table 2**). Meanwhile, a Kaplan-Meier survival analysis showed significantly different OS among the patients with different TNM stages (**Figure 1C**), but not among those who differed in age or sex (**Figures 1A,B**). Interestingly, a trend toward different OS among patients with different histologic subtypes was observed, but the p value was only marginally significant (**Figure 1D**). In the multivariate regression analysis, only a higher TNM stage remained a significant risk factor for OS (**Table 2**).

#### Identification of Transcripts and DNA Methylation Sites as Prognostic Factors

From 16928 mRNAs, 453 miRNA and 395963 DNA methylation sites, a total of 26 miRNAs, 15 mRNAs and 11 DNA methylation sites were initially identified as factors associated with OS using univariate Cox and Lasso-Cox analyses (**Supplementary Table S1**). Next, a list of 2882 genes was then retrieved from the miRTarBase database, of which 21 were identified as highconfidence (p ≤ 0.05) targets of the 26 survival-associated miRNAs. A Lasso-Cox analysis was used to select the mRNAs of the 21 genes that interacted with the corresponding survivalrelated miRNAs. After a multivariate Cox regression analysis of the 73 potential prognostic factors, the overexpression of two miRNAs (MIMAT0002890 and MIMAT0000426) and the hypermethylation of two sites (cg12141052 and cg16404170) were confirmed as significant predictors of worse prognosis (**Table 3**), and the higher expression level of two mRNAs (CDADC1, FAHD2B) was significantly associated with a better prognosis (**Table 3**). Therefore, the final list of candidate prognostic factors for LUAD contained 6 biomarkers, including two miRNAs, two mRNAs and two methylation sites.

The Spearman's rank correlation coefficients for these candidate transcripts and methylation levels were then calculated for the LUAD cohort of 445 patients (**Table 4** and **Supplementary Figure S1**).

TABLE 3 | Genome-wide prognostic factors identified in our study.

HR, hazard ratio; CI, confidence interval; SE, standard error; z value, Wald z-statistic value.

### External Validation of Candidate Transcripts and Methylation Sites

As shown in **Supplementary Figures S2**–**S4**, a univariate Cox proportional hazards regression analysis showed that the six candidate factors identified from either genome or epigenome of LUAD were significantly associated with the survival of patient cohorts in other databases. Moreover, the relationships between their expression levels and the survival rate of LUAD patients were consistent with our findings.

TABLE 4 | Spearman's correlation coefficients among the prognostic factors identified in the study.

### Validation of the Integrated Prognostic Factors

To further assess the predictive capacity of all the candidate prognostic factors, PS was established as an integrated prognostic predictor. To verify the efficiency of PS, the 445 LUAD patients were divided into two groups stratified by the median PS. The high-risk group (PS > 1.88) included 223 patients and the low-risk group (PS < 1.88) included 222 patients. As shown in **Figure 2**, the patients in the low- and high-risk groups displayed significantly different median OS (1070.8 vs. 753.9 days, p < 0.001) and RFS (900.8 vs. 668.2, p = 0.005). As shown in **Figures 3A**, **4A**, the Kaplan-Meier curves and logrank tests indicated significant differences in the OS [hazard ratio (HR): 2.861, 95% confidence interval (CI): 2.052–3.988, p < 0.001] and RFS (HR: 1.77, 95% CI: 1.255–2.497, p = 0.001) between two groups.

Significantly different OS and RFS were also observed among the subgroups in further analyses (**Figures 3B**, **4B**). On the one hand, PS remained efficient in stratifying the patients into different OS (HR: 3.177, 95%CI: 2.110–4.783, p < 0.001) and RFS (HR: 1.752, 95% CI: 1.184–2.595, p = 0.005) when the lowrisk and high-risk subgroups were in the early stages of the disease (**Figures 3C**, **4C**). However, there was only a significant difference in OS (HR: 1.806, 95% CI: 1.019–3.200, p = 0.04) but not RFS (HR: 1.594, 95% CI: 0.763–3.333, p = 0.2) between the two subgroups when both were in the advanced stages of the disease (**Figures 3D**, **4D**). On the other hand, the pathological stage could distinguish significantly different OS (low-risk group: HR: 3.341, 95% CI: 1.888–5.912, p < 0.001; high-risk group: HR: 1.955, 95% CI: 1.305–2.929, p < 0.001) but not RFS (lowrisk group: HR: 1.472, 95% CI: 0.878–2.467, p = 0.1; high-risk group: HR: 1.604, 95% CI: 0.829–3.104, p = 0.2) in the lowrisk and high-risk groups. Thus, PS was proved to be a useful prognostic indicator that can supplement additional information to the TNM stage, especially for LUAD patients in the early stages of the disease. Our study suggests that the combination of the TNM stage and PS increases the accuracy in predicting the outcomes of patients with LUAD.

FIGURE 3 | Kaplan Meier curve showing the overall survival (OS) of the patient cohort grouped by (A) prognostic score (PS), and (B) PS plus pathological stage. OS of patients stratified by PS in subgroups with (C) early-stage tumors and (D) advanced- stage tumors.

#### DISCUSSION

The past decade has witnessed rapid progress in next-generation sequencing and its increasing application in preclinical practice. In recent years, several studies have attempted to associate the transcriptome or epigenome with the clinical outcomes of patients with LUAD (Selamat et al., 2012; Zhang et al., 2017; Gao et al., 2018; He et al., 2018). Zhang et al. analyzed and validated the expression profiles and prognostic values of the mRNAs of five differentially expressed genes associated with DNA methylation in LUAD (Zhang et al., 2017), increasing the likelihood that altered signature genes will become useful biomarkers. Using a TCGA dataset, He et al. (2018) disentangled the relationships between aberrant CpG-methylation and gene expression to identify 10 aberrantly methylated and dysregulated genes. However, their study only focused on the ability of individual genes to predict OS. Another TCGA-based study examined the feasibility of integrating prognosis-related methylation-driven genes into a risk model to predict the OS of patients with LUAD, which also involved a joint survival analysis based on methylation sites and gene expression (Gao et al., 2018). Nevertheless, it remained unclear whether a risk model could improve the accuracy of the TNM stage for survival estimation. Furthermore, no information was given on the predictive value of their method in distinguishing RFS in LUAD patients. None of these studies included the histologic

FIGURE 4 | Kaplan Meier curve showing recurrence-free survival (RFS) of the patient cohort grouped by (A) prognostic score (PS), and (B) PS plus pathological stage. RFS of the patients stratified by PS in subgroups with (C) early-stage tumors and (D) advanced-stage tumors.

subtypes proposed by the International Association for the Study of Lung Cancer/American Thoracic Society/European Respiratory Society (IASLC/ATS/ERS) (Travis et al., 2011; Warth et al., 2012; Hung et al., 2014) as an independent prognostic factor.

To the best of our knowledge, this is the first study to integrate genetic and epigenetic modifications for survival prediction in LUAD patients using TCGA samples. With a comprehensive analysis and screening of mRNA expression, miRNAs and DNA methylation sites based on samples from 445 patients, we identified a set of prognostic factors from both the transcriptome and epigenome. Notably, we included the histologic subtypes and the TNM stages in our initial survival analysis. In this way, we developed a novel subgrouping system that integrates PS and the TNM stage to predict the survival of patients with LUAD.

We started by identifying the clinicopathological characteristics associated with the OS of patients with LUAD. Both Cox regression and Kaplan-Meier survival analyses confirmed the significant prognostic impact of the TNM stage (**Table 2** and **Figure 1**). Further screening of genetic and epigenetic aberrations identified a collection of 26 miRNAs, 15 mRNAs and 11 DNA methylation sites whose expression or methylation levels were significantly associated with OS (**Supplementary Table S1**). Since miRNAs exert their function by regulating the expression of their target mRNAs, we retrieved the potential targets of these 26 miRNAs and performed a

LASSO-Cox analysis to select 21 mRNAs as high-confidence miRNA targets. This provided clues to the potential molecular interactions by which these miRNAs affect the clinical outcomes. Considering the interactions between these candidate prognostic factors into consideration, we performed a multivariate Cox regression to finally identify a list of six survival-related biomarkers. From our data, the expression levels of two mRNAs (CDADC1 and FAHD2B), two miRNAs (MIMAT0002890 and MIMAT0000426) and methylation of two sites (cg12141052 and cg16404170) were strongly associated with the clinical outcomes. PS was then computed as a predictor that integrated these candidate biomarkers and stratified the patients into low-risk (PS < 1.88) and high-risk groups (PS > 1.88). The efficiency of PS was confirmed by our success in distinguishing the OS and RFS of LUAD patients (**Figures 3**, **4**). A subgroup analysis further demonstrated that a more precise prediction of survival could be achieved for patients with LUAD by combining PS with the TNM stage, which should allow more timely therapeutic interventions.

To be noted, Targetscan4 was preferentially considered for the validation of our candidate miRNAs, however, the small number of miRNA targets shared between miRTarBase and Targetscan limited its use (**Supplementary Figure S5**).

Ten survival-associated genes, whose aberrant expression was affected by methylation, have been identified previously by He et al. (2018) from the TCGA data portal5 . Therefore we attempted to include the mRNAs of these 10 genes in our transcripts for further screening. However, as shown in **Supplementary Table S2** and **Supplementary Figures S6**–**S8**, integrating the mRNA of BLK which was identified in the Cox regression analysis into PS did not improve its predictive ability.

In terms of the disproportionate number of non-smokers in the selected patient cohort, secondary analyses were therefore performed to assess the potential value of PS in predicting the survival of the non-smokers and smokers in our cohort. As shown in **Supplementary Figure S9A**, Kaplan-Meier curves and the log-rank test indicated a significant difference in OS between two groups (HR: 2.785, 95% CI: 1.071–7.24, p = 0.03). Significantly different OS was also observed among subgroups stratified by PS plus the TNM stage (**Supplementary Figure S9B**). However, the performance of PS was not satisfactory for the non-smokers (**Supplementary Figures S9C,D**), especially in stratifying patients with advanced-stage LUAD into subgroups with different OS, which might be attributed to the limited number of non-smokers (n = 65). On the contrary, PS remained consistently efficient in stratifying OS in the smokers (n = 359) (**Supplementary Figure S10**).

There were several limitations to our study. For instance, risk factors such as packages of cigarettes and adjuvant therapy were not included in our analysis because of their interpatient heterogeneity. Moreover, the histologic subtypes was unsatisfactory in distinguishing prognoses in the multivariate analysis which could be explained by the missing histologic information for almost half the patients. It is noteworthy that the failure of PS to distinguish RFS in the advanced-stage subgroups

4http://www.targetscan.org/vert\_72/

5https://cancergenome.nih.gov

(p = 0.2) should be possibly attributed to the limited number of LUAD patients with advanced-stage disease. Last but not least, the clinical utility of PS identified here may be limited in patients with small-sized lesions because of the difficulty in extracting sufficient RNA and protein. More studies are warranted to assess the roles of these candidate prognostic factors in LUAD.

### CONCLUSION

In conclusion, using a TCGA dataset of 445 LUAD patients, we identified six prognostic factors (two mRNAs, two miRNAs and two DNA methylation sites) for LUAD from the genome and epigenome, and developed PS from them. Combining the TNM stage and PS provided additional precision in stratifying patients into significantly different OS and RFS subgroups. Further studies are warranted to assess the efficiency of PS and to explain the effects of these observed genetic and epigenetic aberrations in LUAD.

### AUTHOR CONTRIBUTIONS

DC, YS, SL, CC, and YC designed the study. DC, YS, FZ, and EZ performed the data analysis. DC, YS, FZ, XW, and XZ wrote the manuscript. GJ, SL, CC, and YC reviewed and edited the manuscript. YC and FZ afforded to conduct our study.

### FUNDING

The work was supported from the National Natural Science Foundation pre Research Fund Project (SDFEYGJ1704), Jiangsu Provincial Commission of Health and Family Planning (Grant H201521), the Natural Science Foundation of Jiangsu Province (Grant BK20161224), Suzhou Key Discipline for Medicine (SZXK201803), and the Science and Technology Research Foundation of Suzhou Municipality (SYS2018063).

#### ACKNOWLEDGMENTS

We thank International Science Editing (http://www. internationalscienceediting.com) for editing this manuscript.

### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2019. 00493/full#supplementary-material

FIGURE S1 | Plot of Spearman's rank correlation coefficients among candidate transcripts and methylation levels in LUAD (n = 445).

FIGURE S2 | Differential expression and prognostic impact of (A) FADH2B and (B) CDADC1 (two candidate mRNAs) in LUAD patients. Kaplan-Meier curves of 720 LUAD patients, who were separated into high-expression and low-expression groups, using as cutoffs the best-performing thresholds of the different genes. All values were significant (p < 0.05).

FIGURE S3 | Differential expression and prognostic impact of (A) MIMAT0002890 and (B) MIMAT0000426 (two candidate miRNAs) in LUAD patients. Kaplan-Meier curves of 720 LUAD patients, who were separated into high-expression and low-expression groups, using as cutoffs the best-performing thresholds of the different miRNAs. All values were significant (p < 0.05).

FIGURE S4 | Differential expression and prognostic impact of (A) cg12141052 and (B) 16404170 (two candidate methylation sites) in LUAD patients. Kaplan-Meier curves of 720 LUAD patients, who were separated into high-expression and low-expression groups, using as cutoffs the best-performing thresholds. All values were significant (p < 0.05).

FIGURE S5 | Venn diagrams representing 1161 miRNA targets that overlapped between miRTarBase and Targetscan.

FIGURE S6 | Median overall survival (A) and recurrence-free survival (B) of patients in the high-risk (PS > 2.78) and low-risk groups (PS < 2.78). Patients in the high-risk group showed a significantly shorter survival than those in the low-risk group (∗∗∗p < 0.001).

FIGURE S7 | Kaplan Meier curve showing the overall survival (OS) of the patient cohort grouped by (A) recombinant prognostic score (PS), and (B) recombinant

### REFERENCES

PS plus pathological stage. OS of the patients stratified by PS in subgroups with (C) early-stage disease and (D) advanced-stage disease.

FIGURE S8 | Kaplan Meier curve showing recurrence-free survival (RFS) of the patient cohort grouped by (A) recombinant prognostic score (PS), and (B) recombinant PS plus pathological stage. RFS of the patients stratified by PS in subgroups with (C) early-stage disease and (D) advanced-stage disease.

FIGURE S9 | Kaplan Meier curve showing the overall survival (OS) of the non-smokers in our cohort grouped by (A) prognostic score (PS), and (B) PS plus pathological stage. OS of the patients stratified by PS in subgroups with (C) early-stage disease and (D) advanced-stage disease.

FIGURE S10 | Kaplan Meier curve showing the overall survival (OS) of the smokers in our cohort grouped by (A) prognostic score (PS), and (B) PS plus pathological stage. OS of the patients stratified by PS in subgroups with (C) early-stage disease and (D) advanced-stage disease.

TABLE S1 | Transcripts and DNA methylation sites whose expression levels showed significant association with overall survival.

TABLE S2 | Integrated genome-wide prognostic factors in our study.

in stage II and III lung adenocarcinomas and nodal metastases. J. Thorac. Oncol. 12, 458–466. doi: 10.1016/j.jtho.2016.10.015

cancer (NSCLC): an analysis of the surveillance, epidemiology, and end results (SEER) registry. J. Thorac. Oncol. 10, 682–690. doi: 10.1097/jto.000000000000 0456

Zhang, Y., Zhao, W., and Zhang, J. (2017). Comprehensive epigenetic analysis of the signature genes in lung adenocarcinoma. Epigenomics 9, 1161–1173. doi: 10.2217/epi-2017-0023

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Chen, Song, Zhang, Wang, Zhu, Zhang, Jiang, Li, Chen and Chen. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# F8 Inversions at Xq28 Causing Hemophilia A Are Associated With Specific Methylation Changes: Implication for Molecular Epigenetic Diagnosis

Muhammad Ahmer Jamil1† , Amit Sharma1† , Nicole Nuesgen1 , Behnaz Pezeshkpoor1 , André Heimbach2 , Anne Pavlova1 , Johannes Oldenburg1 and Osman El-Maarri1 \*

#### Edited by:

Dongyi He, Shanghai Guanghua Rheumatology Hospital, China

#### Reviewed by:

Apiwat Mutirangura, Chulalongkorn University, Thailand Hehuang Xie, Virginia Tech, United States

\*Correspondence:

Osman El-Maarri<EMAIL_ADDRESS>†These authors have contributed equally to this work as first authors

#### Specialty section:

This article was submitted to Epigenomics and Epigenetics, a section of the journal Frontiers in Genetics

Received: 09 February 2019 Accepted: 09 May 2019 Published: 29 May 2019

#### Citation:

Jamil MA, Sharma A, Nuesgen N, Pezeshkpoor B, Heimbach A, Pavlova A, Oldenburg J and El-Maarri O (2019) F8 Inversions at Xq28 Causing Hemophilia A Are Associated With Specific Methylation Changes: Implication for Molecular Epigenetic Diagnosis. Front. Genet. 10:508. doi: 10.3389/fgene.2019.00508 1 Institute of Experimental Hematology and Transfusion Medicine, University of Bonn, Bonn, Germany, 2 Institute of Human Genetics, School of Medicine, University of Bonn – University Hospital Bonn, Bonn, Germany

Diverse DNA structural variations (SVs) in human cancers and several other diseases are well documented. For genomic inversions in particular, the disease causing mechanism may not be clear, especially if the inversion border does not cross a coding sequence. Understanding about the molecular processes of these inverted genomic sequences, in a mainly epigenetic context, may provide additional information regarding sequencespecific regulation of gene expression in human diseases. Herein, we study one such inversion hotspot at Xq28, which leads to the disruption of F8 gene and results in hemophilia A phenotype. To determine the epigenetic consequence of this rearrangement, we evaluated DNA methylation levels of 12 CpG rich regions with the coverage of 550 kb by using bisulfite-pyrosequencing and next-generation sequencing (NGS)-based bisulfite re-sequencing enrichment assay. Our results show that this inversion prone area harbors widespread methylation changes at the studied regions. However, only 5/12 regions showed significant methylation changes, specifically in case of intron 1 inversion (two regions), intron 22 inversion (two regions) and one common region in both inversions. Interestingly, these aberrant methylated regions were found to be overlapping with the inversion proximities. In addition, two CpG sites reached 100% sensitivity and specificity to discriminate wild type from intron 22 and intron 1 inversion samples. While we found age to be an influencing factor on methylation levels at some regions, covariate analysis still confirms the differential methylation induced by inversion, regardless of age. The hemophilia A methylation inversion "HAMI" assay provides an advantage over conventional PCR-based methods, which may not detect novel rare genomic rearrangements. Taken together, we showed that genomic inversions in the F8 (Xq28) region are associated with detectable changes in methylation levels and can be used as an epigenetic diagnostic marker.

Keywords: structural rearrangement, epigenetic, DNA methylation, inversion, hemophilia, molecular diagnosis

## INTRODUCTION

fgene-10-00508 May 27, 2019 Time: 14:39 # 2

Implications of human DNA sequence variations have received considerable attention in recent years and structural variants (SVs) are considered an important contributor among them. SVs, particularly inversions, can vary in size from few nucleotides to large-scale chromosomal rearrangements. The inversions can have functional consequences by truncating a given gene (or genes) or by rearranging the regulatory element in the local proximity, both having a disproportionate impact on gene expression and transcriptional variability (Puig et al., 2015).

Hemophilia A (OMIM #306700), an inherited bleeding disorder, harbors two such rearrangements at chromosome X (Xq28) involving the coagulation factor VIII (F8) gene. F8 (∼186 kb; 26 exons) is located at the telomeric end of the X- chromosome and contains regions with high GC content, which makes it more susceptible toward the methylated cytosine deamination mutations [2,537 mutations, reported by CDC Hemophilia A Mutation Project (CHAMP) (Payne et al., 2013)]. In addition, two hotspot inversions (known as intron 1 and intron 22 inversions) are reported accounting for 40– 50% of patients with severe hemophilia A (Andrikovics et al., 2003; Oldenburg and El-Maarri, 2006; Zimmermann et al., 2011). These hotspot recurring inversions are caused by intrachromosomal homologous recombination between identical inverted repeats: two long repeats located within the F8 locus: the Int22h-1 in intron 22 and the Int1h-1 in intron 1 (Lakich et al., 1993; Naylor et al., 1995). The former is 9.1 kb in length and has two additional homologs (Int22h-2 and 3) at about 500– 580 kb distance toward the telomere, while the latter is about 1 kb and has one homolog located 141 kb toward the telomere (**Figure 1**) (UCSC genome browser). Both repeats are prone to intra-chromosomal homologous recombination leading to an inversion of the intervening sequence, thus leaving the F8 split into two parts of opposite transcriptional direction (Bagnall et al., 2002, 2006). The clinical result of such inversions is a severe hemophilia A (HA) phenotype with no functional FVIII protein. Inversion events leading to human diseases are not limited to F8 gene, other genes, such as IDS gene (Hunter syndrome), MSH2 gene (Lynch syndrome), EML4-ALK rearrangement in non-small cell lung cancer (NSCLC), AP3B1 (Hermansky-Pudlak syndrome type 2), have been previously implicated (Bondeson et al., 1995; Soda et al., 2007; Jones et al., 2013; Rhees et al., 2014).

However, the effect of a given DNA inversion may not be as clear as the above examples. It is not known whether SVs without a clear gene-destruction effect are still benign in nature. For instance, it is likely that an inversion can disturb normal chromatin architecture and this could be translated into interchanging of hetero- and euchromatic states. This would lead to abnormal methylation patterns, and ultimately to alterations in gene expression that may have a phenotypic impact. Thus, a given gene could cross the borders between an actively transcribed and a non-active region as a result of the inversion. A clear example is what has been observed in Drosophila position-effect variegation, where an inversion of DNA shifts the w+ and rst+ genes from an euchromatin to a heterochromatin domain, thus resulting in white color eyes (Schotta et al., 2003).

In humans, indications for the none gene breaking effects of inversions came from Gonzalez et al. (2014) who reported on the effect of a common 0.45 Mb inversion at 16p11.2 on local gene expression and found that inverted alleles strongly correlated to neighboring gene expression. Expression effects were seen on single copy genes within the inverted regions as well as on genes flanking the duplicated regions (where the inversion breakpoints occur). Additionally, the multiple copy genes located in the duplications were also affected. Some genes are over-expressed, while others are under-expressed in the inverted allele. However, a large proportion remained unaffected. Although the molecular mechanism behind this set of observations goes beyond the scope of this particular study, it provides an indication for a cause-effect relationship between common human inversions and gene expression and its link to a disease phenotype: the joint susceptibility to asthma and obesity.

To date, there is no molecular mechanism that explains the biogenesis for the inversion effect on expression. Furthermore, it is still not possible to predict the effects of a given inversion. It has been hypothesized that changes in the chromatin structure comprise a possible underlying reason, but a clear model detailing the interplay between a given methylation and the different parameters that affect the gene expression, such as histones modifications, DNA methylation, nucleosome occupancy and three-dimensional chromatin structure, remains elusive. The above-described F8 inversions are well characterized and their breakpoints are within defined unique repeats regions. Therefore, these two inversions are a suitable model for investigating the effect of inversions on gene expression as well as chromatin structure and epigenetic modifications.

In this study, we took advantage of the F8 gene inversions model to analyze DNA methylation levels of CpG rich regions within and flanking the inverted DNA regions in wild type and inverted DNA (with intron 1 or intron 22 inversions). In summary, our results show clear detectable DNA methylation changes associated with inversions that are flanking the inverted regions. Therefore, methylation aberrations are a useful diagnostic tool to identify inversion structural variations.

### MATERIALS AND METHODS

#### DNA Samples

DNA samples corresponding to healthy controls (21 nonhemophilic males) and to male hemophilia patients with known intron 1 (16 samples) or intron 22 inversions (19 samples) were obtained from the hemophilia center at the Institute of Experimental Hematology and Transfusion Medicine (University Clinic Bonn, Germany) and from the Institute for Human Genetics (University of Wuerzburg, Germany). The samples used are derived from DNA collected for molecular diagnostic purposes. All blood samples from patients and healthy controls were obtained upon written informed consent. The Ethics Committee of the University Clinic Bonn authorized the use of pre-collected DNA samples for research purposes (approval number 091/09 date 05/06/2009).

FIGURE 1 | Pyrosequencing methylation data on 12 selected regions from intron 22 and intron 1-inversion samples as well as healthy male controls. (A) Detailed map on X chromosome (Chr X: 154,027,275-154,751,861:hg19) showing F8, the three Int22h and the two Int1h repeats involved in the inversion mutations. The positions of the studied regions are indicated in the middle region by capital letters. The inversion prone regions are labeled with red and blue horizontal lines for intron-22 and intron-1 inversions, respectively. (B) Methylation data represented by heatmaps, sample PCA and variable PCA plots. (C) Detailed data of the methylation values for individual samples at the two best regions (H and F) that clearly distinguish between inverted and non-inverted control samples.

## Methylation Analysis and Pyrosequencing Assays

fgene-10-00508 May 27, 2019 Time: 14:39 # 4

CpG rich regions within and around the inverted sequences were identified using the UCSC website. Feasible regions for primer designs were selected for methylation analysis. Primers for PCR amplifications as well as pyrosequencing primers were designed using the PyroMark assay-design Q24 software (Qiagen, Germany). Bisulfite treatment was done using the EZ 96-DNA methylation kit (Zymo Research, Irvine, CA, United States) following manufacturer's protocol. Bisulfite PCRs were done using HOT FIREPol (Solis Biodyne, Tartu, Estonia). Pyrosequencing was done on a PyroMark Q24 or Q96 machine (Qiagen, Germany). Primers used for amplification are listed in **Supplementary Table S1**. In total, 12 different regions were studied and designated as regions A to O.

### Verification of Results by Bis-Seq NGS Based Assay

Since the pyrosequencing assay is restricted to check only few individual CpGs and provides an estimated average of methylation the results had to be verified by covering relatively larger regions and the spatial relationships (phase) between different CpGs in the same region had to be revealed. Such data could be provided by NGS-based resequencing assays. For this purpose, we chose the SeqCap Epi Enrichment system from NimbleGen (Roche, Switzerland). Using this system, we targeted the F8 region: chrX: 154,027,275-154,751,861. Samples included four intron 1, six intron 22 inversions and four wild type controls. After obtaining the data we filtered for the overlapping reads with our pyrosequencing assays. All data are submitted to EBI as a mapped "BAM" file under study accession number "ERP113762."

### Next Generation Sequencing Analysis

Sequencing data was generated using Illumina HiSeq 2500 v4 with read length of 2 × 125 bp. Reads were generated in fastq file format. Reads were pre-filtered for any adapters' sequence. Reads quality was tested using fastqc1 and all reads were passed for the quality cut-off of 10. Reads were mapped using BSMAP (Xi and Li, 2009) program to HG38 genome downloaded from UCSC with parameter settings to WGBS mode (−s = 16), all four strands mapping (−n = 1) and with parallel computing of four processor cores (−p = 4). Mapped reads were split into top and bottom strand using bamtools (Barnett et al., 2011) to separately remove duplicates for both strands. Duplicates were removed using the "MarkDuplicates" function in picard tool2 . Removed duplicate removed strands were merged together into single mapped file using bamtools. Filtered reads were filtered again for the properly paired reads using bamtools filter with parameters of "-isMapped true," "isPaired true" and "isProperPair true." Properly paired reads were further processed using "clipOverlap" function in "bamUtil" (Jun et al., 2015) to clip overlapping paired-end reads to correct bias for methylation calculation. Methylation percentages were determined using the "methRatio.py" function in BSMAP with the parameters of minimum number of reads per CpGs set to 1 (m = 1) and report to zero methylation (−z). A final methylation table with number of Cs, Ts and coverage for every CpG was created by removing the uncovered region via NimbleGen. Methylation analyses were further carried out in R using the "methylkit" (Akalin et al., 2012) package. Fisher's exact test was performed to calculate the p-value between samples for every CpG site.

### Statistical Analysis and Data Visualization

Statistical analysis was done using R or Prism (GraphPad software). Additional data analysis and visualization were done using Qlucore Omics Explorer (Sweden) and ProFit software from Quatum Soft (Switzerland). Regression analyses using R were performed to understand the effect of covariates (age). Formula for regression analysis used were "aov(lm(MethDiff∼CaseControl+Age+CaseControl∗Age)."

## RESULTS

### CpGs Regions at the Border of the Inverted DNA Are Prone to Significant Differential Methylation

The main aim of this study was to detect differential methylation region(s) that could serve as markers for identification of F8 inversions rearrangements. We initially designed and selected the regions based on (1) feasibility of reading methylation of at least three CpGs, (2) their presence in a region between the three prime regions of F8 and the Int22h3 repeat regions, and (3) their presence in non-repetitive regions (like L1 and Alu). Next, we could retain 12 regions whose methylation was neither constantly 0 nor 100% for all samples: i.e., variable methylation. We then studied three groups of samples: int22 and int1 inversions and healthy controls. Three of the regions failed to show significant statistical differences when applying statistical test to compare between the groups, namely regions G, E, and O (**Figure 1**).

The rest of the nine regions showed statistical significance for at least one CpG at one of three comparisons (**Figure 1** and **Supplementary Table S2**). For intron 22 inversion samples, eight individual CpGs in six regions (regions H, L, A, N, J, and I) were

1http://www.bioinformatics.babraham.ac.uk/projects/fastqc

2http://broadinstitute.github.io/picard/

statistically different compared to healthy controls. The most significant region was region H (average meth. diff. = 6% at CpG2; t-test p < 0.0001) embedded within the Int22h repeats, followed by region L in exon 14 of factor 8 (average meth. diff. = 4% at CpG2; t-test p = 0.0005).

For intron 1 inversion samples, eight individual CpGs were statistically also significantly differentially methylated in comparison to healthy controls (**Figure 1** and **Supplementary Table S2**), covering five different regions: A, C, F, J, and I. The most significant region was region F embedded in the Int1h repeat (average meth. diff. = 23.7% at CpG3; t-test p < 0.0001), followed by region C (average meth. diff. = 2.9% at CpG1; t-test p = 0.0004). However, region I showed higher average differential methylation reaching 12.6%, but a p-value of 0.0011.

Of significance, the regions that showed the highest differential methylation were situated either within the repeats involved in the homologous recombination leading to the inversion (region H: intron 22 inversion and region F: intron 1 inversion) or close to that border (region C, L, and I).

### Two Regions Show Promising Biomarkers Properties: High Sensitivity and Specificity, Making Them Eligible as Diagnostic Markers

In order to use methylation at a given CpG as a diagnostic marker to detect inverted DNA we calculated sensitivity and specificity for each CpG that showed a statistical significance difference between inversions and healthy controls. For this purpose, we defined sensitivity as the fraction of the inverted DNA sample that is identified as differentially methylated in comparison to the wild type controls. Whereas specificity is defined as the fraction of healthy samples within the normal range of methylation and not overlapping with inverted DNA. For intron 1 and intron 22 inversions, a sensitivity and specificity of 1 were reached for region F CpG3 and region H CpG2 (**Supplementary Table S2**).

### Investigation of Factors Influencing DNA Methylation: Age and DNA Polymorphism

#### Age Effect: Healthy Group Shows Statistically Significant Linear Correlation Between Age and Methylation at Regions F and I and a Clear Tendency at Region L

A significant correlation between age and methylation difference was observed for some CpG sites. In order to understand whether the difference is due to age or rearrangements of the inversion region, we performed rigorous regression analysis between inversion samples and controls with age as a covariate (**Supplementary Figure S2**). Regression analysis revealed that some CpGs sites, i.e., F-CpG1, F-CpG2, F-CpG3, I-CpG1, and I-CpG2, showed statistical significance between intron-1 inversion and control in the difference between age and methylation and the difference between phenotype and methylation, while the difference between age and phenotype was not found to be significant. Thus, the difference in methylation due to intron-1 inversion will be statistically significant at any age range (**Supplementary Figure S2A**). Regarding intron-22 inversion, we found no statistically significant difference between age to methylation or phenotype to methylation (**Supplementary Figure S2B**).

In order to re-emphasize the age effect on region F and to exclude an effect on the ability to discriminate wild type from inversions at any age group, we calculated observed – expected – methylation levels for all samples in the healthy and the intron 1 inversion groups. For this purpose, we calculated the expected methylation values according to an equation of best fit linear regression model of healthy samples for each of the three CpGs and the average of the three CpGs in region F (**Supplementary Figure S3A**). Comparison of observed and expected levels showed a highly significant difference only at intron1 samples (at all three CpGs), which indicates that the observed differences between intron 1 inversions and healthy controls are not solely due to an age effect (**Supplementary Figure S3B**). Moreover, observed methylation values minus calculated expected methylation values (according to age using the linear regression fitting equation of the healthy samples) revealed high significance between healthy controls and intron 1 inversion compared to intron 22 inversions (**Supplementary Figure S3C**). This once more indicates that the differences in comparison to healthy controls are largely due to the intron 1 inversion of DNA.

### DNA Polymorphism Effect

As the DNA polymorphism may affect the level of methylation at neighboring CpGs, we searched the UCSC databases for occurrences of polymorphisms in a window of 1 Kb surrounding each investigated CpG. The results are displayed in **Supplementary Table S3**. While we found some SNPs with minimum allele frequency up to 0.21 in European populations, especially in the two regions with high discrimination power to detect intron 1 (region F) and intron 22 (region H), no reported SNPs with MAF > 0.05 have been reported. Therefore, we could largely exclude a broader effect of polymorphism on the level of methylation at the two relevant regions H and F.

### Methylation Correlation Between Different Regions Suggests Stochastic Random Effect, While Top Differentially Methylated Regions Are Indeed Correlated

In intron 1 and intron 22 inversions, we observed abnormal methylation at several CpG sites. Therefore, we queried whether these changes are coordinated and if they are correlated. In other words, are these changes in methylation in parallel at two or more altered regions for a given inversion type. If this is the case, a statistically significant correlation should be observed. Indeed, we calculated all correlations pairwise for all 22 CpG sites for every group (intron 22, intron 1 and healthy controls as separate groups). While we observed 15, 15, and 16 inter CpGs correlations in intron 22, intron 1 and healthy controls, we found little overlaps between all three groups. This was mainly observed at the intra-CpG correlations within the regions N and F

triangle). The CpG rich region names are labeled with capital letters, while the individual CpGs are labeled with numbers, whereby red ones represent statistically significant ones. When two significant CpGs from two regions are correlated they are highlighted with a blue circle. (B) Correlation graphs of the circled ones of part A. The best fit linear curves as well as the 95% confidence intervals are shown in red.

(**Figure 2A**). The absence of overlap suggests a change in the nature of epigenetic marks from the normal non-inverted to an inverted DNA.

In this context, three inter region correlations were observed in inversion groups that involve regions that are differentially methylated between inversions and controls. Two of these are observed in intron 22 only and are not present in controls, namely H-CpG1 vs. I-CpG1 and L-CpG2 vs. N-CpG2 (**Figure 2**). Possibly, this is specific for the inversion samples and is induced by the rearrangement. This hypothesis is supported by two arguments: (1) all four involved CpGs are at the top differentially methylated between intron 22 inversions and controls and (2) such correlation is absent in normal samples.

Yet, the third correlation was observed in intron 1 samples between F-CpG1 and I-CpG1. These CpGs also showed a significant methylation difference between intron 1 inversion samples and controls. In fact, these two CpGs are the highest two differentially methylated CpGs with an average difference between intron 1 samples and controls of 23.7 and 12.6% for F-CpG1 and I-CpG1, respectively. However, this correlation is induced by the age effect on methylation as this involves two CpGs that show high correlation between age and methylation (**Supplementary Figure S2**). In addition, this correlation is observed in healthy controls, re-emphasizing that the correlation between age and methylation is the driving force behind this correlation between the CpG1 at region F and CpG1 at region I.

FIGURE 3 | NGS results of the studied regions shown in Figure 1A. (A) Each graph represents one region; regions A and N have no enough coverage and are absent. The number of reads for each CpG is shown below the corresponding CpG, the p-value of Fisher's exact test is shown when significant (marked by X) between healthy samples (green) and intron 1-inversion samples (blue) or intron 22-inversion samples (red). The corresponding pyrosequencing CpGs are in red and underlined. (B) Correlation graphs between the pyrosequencing and the NGS methylation levels results.

### Targeted Bisulfite Re-sequencing Largely Confirms the Pyrosequencing Results

#### Confirmation of Bisulfite Pyrosequencing

In order to further confirm the above results via an alternative method, we performed targeted bisulfite re-sequencing with the SeqCap Epi Enrichment system from NimbleGen (Roche, Switzerland) to capture a region containing the F8 and extending up to the extragenic Int22h repeats (hg19, ChrX: 154,027,275– 154,751,861; **Figures 3**, **4**). Ten of the studied pyrosequencing regions could be covered, while two were insufficiently covered with low read counts (regions A and N). In order to increase coverage and to decrease the effect of inter-individual differences, we merged the reads that belonged to the same group of samples. This resulted in a pool of reads of three groups including six, two and four individual DNA samples for intron 22 inversion, intron 1 inversion and healthy controls, respectively. Moreover, this merging approach increased the read numbers at each CpG site with ranges 23-717, 4-227, and 9-466 for intron 22 inversion, intron-1 inversion and healthy controls, respectively.

Nineteen individual CpGs were overlapping between the pyrosequencing and the NGS enrichment approach, of which 12 were showing complete concordance in results of significance (**Figure 3** and **Supplementary Figure S1**). This is also reflected by the correlation of average methylation between both approaches across the three samples cohorts (**Figure 3B**). Of particular interest are the two highly differentially methylated regions that showed high sensitivity and high specificity for distinguishing inversions from non-inversions, namely regions F and H for detecting intron 1 and intron 22, respectively. Both showed high significance in correlation values and overlapping, confirming results in both methods. The bisulfite targeted enrichment analysis provided additional confidence in the inversion-induced methylation aberrations and in the ability of such methylation assays to detect the inversions.

### Trend Line of Methylation Changes Over the F8 Till Intr22h3 Covered Region

Next, a global trend line was drawn of the methylation differences including all CpGs captured by the enrichment protocol (i.e., not only overlaps with pyrosequencing results as presented in the previous section). For this approach, we filtered the data to exclude any CpG overlapping with a repeat or a known SNP. Additionally, we excluded data for any CpG that had less than 30 reads in one of the two compared categories. In a next step, a trend line of difference in methylation to the healthy male controls was plotted in a map showing relative position to the studied pyrosequencing regions (**Figure 4A**). As expected, this approach indicated a major hypermethylated domain overlapping with the regions F and H for intron 1 and intron 22, respectively. However, we also noticed that the inversion breakpoints (shown as blue and red stars in **Figure 4A**) are lying in "methylation-disturbed" domains. The largest methylation disturbance in both magnitude and length of the domain appear to be overlapping with the inversion junctions. All of the above suggest that the observed methylation alterations are indeed reflection of new genomic architecture caused by the DNA inversion.

### Characteristics of CpGs Showing Differential Methylation

We investigated the relationship between the degree of CpG methylation difference and the density of CpGs in a window of 50 bp where the CpG in question is in the center (**Figure 4B**). Using all data for all CpGs (regardless of statistical significance) we found that a clear and highly significant relationship between the methylation differences and CpG density exists where relatively CpG dense regions are more stable and show smaller methylation differences. This applies for the comparisons healthy vs. intron 1 (r = −0.909, Fisher's exact test p < 0.0001) and healthy vs. intron 22 (r = −0.904, Fisher's exact test p < 0.0001) (**Figure 4B**). It is our opinion that this is a general phenomenon of variability of methylation at "stand-alone" CpGs where they are more prone to uncontrolled "natural" fluctuations. However, the CpG methylation at significantly differentially methylated CpGs failed to show this correlation indicating that the latter are the result of induced aberrant methylation due to DNA rearrangement. From this analysis, we conclude that statistically significant methylation changes are more likely to occur at CpG rich regions or at clusters rather than at sole dispersed (nonclustered) CpGs.

### DISCUSSION

The human genome shows significant variability between individuals (Auton et al., 2015). This variability is caused by single nucleotide polymorphisms, deletions, duplications, translocations and inversions. The effect of which may either be detectable as a change in the phenotype (which include disease manifestation) or be benign without observable phenotype. The molecular mechanism for the former can be explained for SNPs, deletions or duplications by virtue of possible changes in the DNA sequences leading to altered gene expression or protein structure. However, in the cases of translocations or inversions, there is no net gain or loss in DNA. Therefore, the association to a phenotype is difficult to explain by DNA changes unless the breakpoints disrupt a coding sequence or an expression-regulatory element (like a promoter or enhancer) (Harewood and Fraser, 2014). However, an additional scenario could be responsible to cause a

FIGURE 4 | Global visualization of NGS data in the F8 region (hg19: Chr X 154,027,275–154,751,861). (A) Upper panel shows the relative positions of the studied pyrosequencing regions, the middle panel shows the NGS data for intron 1 inversion samples and the lower panel the intron 22 inversion samples. The covered individual CpG methylation data are represented by a gray dot, while additionally, the data is represented by a smooth curve representing the trend of changes between the inverted and the control samples. CpG sites with less than 30 coverage or overlapping with known SNPs or repeats were excluded. Red and blue stars indicate the DNA inversion junctions. (B) Correlation between the methylation differences at a given CpG and the density of CpG within 50 bp flanking region. Left and right side include all CpG data and only significant data (Fisher's exact test), respectively.

phenotype: a shift in chromatin structure or – as it is also known – a position effect variegation (PEV).

Position effect variegation is one of similar phenomena which occur due to relocation of a genomic segment from one region to another and it has been extensively studied in Drosophila, yeast, mice and cultured human cells (Tham and Zakian, 2002; Pedram et al., 2006; Elgin and Reuter, 2013; Tchasovnikarova et al., 2015). Inversion prone position effects are not only limited to other species, it has also been reported in some human disease conditions, such as aniridia (PAX6), campomelic dysplasia (SOX9), familial adenomatous polyposis (APC) and Saethre-Chotzen syndrome (TWIST1) (Fantes et al., 1995; de Chadarevian et al., 2002; Cai et al., 2003; Velagaleti et al., 2005). Of note, some inversion variants can also act as risk factor for the offspring in microdeletion syndromes, such as Williams–Beuren syndrome, Angelman syndrome and Sotos syndrome (Osborne et al., 2001; Gimelli et al., 2003; Visser et al., 2005).

The above would lead to recreation of chromatin domains that result in local and regional epigenetic changes like DNA methylation aberrations. In this study, we used the two inversion hotspots in the F8 gene at Xq28 as a model to investigate the global methylation aberration. Indeed, we found specific changes associated with each of the two inversions. With one specific region for each of the inversions showing high sensitivity and high specificity, our results pave the way for the use of methylation-based assay to detect the inversion. The hemophilia A methylation inversion "HAMI" assay will have several advantages over traditional assays. It is noteworthy to mention that repetitive elements also play an important role in generating structural variants (SVs) in humans (Xing et al., 2009). Among all mobile element types, long interspersed element-1 (LINE-1, or L1) has been previously investigated for DNA methylation-related changes in diseased conditions (Nusgen et al., 2015; Sharma et al., 2019). In this particular study, we took advantage of one such full length L1 repeat (region O) located in the vicinity of the F8 gene and evaluated the methylation status of this repeat in both inversion type patients. However, no differences were found between inversion and wild type.

Currently, the gold standard molecular diagnostic assay to detect the inversion is the inverse based PCR assay (Rossetti et al., 2005), a procedure that needs up to 2–3 working days to complete and requires a skilled technician to perform a critical ligation step. In comparison, the HAMI assay includes three fail-free steps: (1) bisulfite conversion, (2) PCR and (3) quantitative pyrosequencing, all of which could be performed in 1 day. An additional advantage for HAMI is that it does not detect a specific DNA junction. Therefore, no specific amplification primers to detect only known inversions are required, while any rearrangement that could still be missed by specific amplification across known rearrangement junctions will be detected. However, disadvantages and limitations of such an assay include establishment of controls to define the relative borders (cut-off) of normal levels, as this could be population- or ethnicity-specific.

Overall, we could determine the methylation levels at multiple regions surrounding/overlapping F8 associated genomic inversions at Xq28 region. Further evaluations are required to establish whether these epigenetic changes are cause or consequence of these inversion events.

### DATA AVAILABILITY

The datasets generated for this study can be found in the ENA, https://www.ebi.ac.uk/ena/data/view/PRJEB31235.

### ETHICS STATEMENT

The Ethic committee of the University Clinic in Bonn authorized the use of pre-collected DNA samples for research purposes, approval number 091/09 date 05/06/2009.

### AUTHOR CONTRIBUTIONS

MJ, AS, and OE-M wrote the manuscript. MJ, AS, BP, and OE-M analyzed the results. MJ and OE-M performed the bioinformatics analysis. AS and NN performed the experiments. BP, AH, AP, and JO provided the samples and infrastructure and commented on the manuscript. OE-M designed the study.

## FUNDING

This study was supported by an ASPIRE hemophilia research award (2014) from Pfizer (grant reference WI193463) and institutional fund of Institute of Experimental Hematology and Transfusion Medicine, University of Bonn, Germany.

## ACKNOWLEDGMENTS

We thank the patients for donating DNA material for research purposes.

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2019. 00508/full#supplementary-material

FIGURE S1 | Detailed pyrosequencing methylation data on 12 selected regions from intron 22 and intron 1-inversion samples as well as healthy male controls. (A) Detailed map of X chromosome (154,027,275–154,751,861:hg19) showing the F8, the three Int22h and the two Int1h repeats involved in the inversion mutations. The positions of the studied regions are represented in the middle by capital letters, below which is the max-observed difference in methylation average. The inversion prone regions are labeled with red and blue horizontal lines for the intron-22 and intron-1 inversions, respectively. (B) Vertical scatter plots represent the detail methylation levels of individual studied samples for all individual CpG sites. The significant t-test is shown in the figure.

FIGURE S2 | Age covariate regression analysis showing correlation between age and methylation levels for healthy controls in comparison to intron 1-inversion samples (A) and to intron 22-inversion samples. (B) Every plot shows methylation data vs. age. Above the individual plots are Pearson correlation p and rho values, while below the p-values of age-covariate analysis are shown. All significant p-values are written in bold. In case of significant Pearson correlation, the values are labeled with solid transparent red, green or blue rectangles. The plots corresponding to significant differences between cases and controls, even after

### REFERENCES

considering age as covariates (P-value. MethDiff∼CaseControls), are indicated by a red cadre.

FIGURE S3 | Calculation of observed vs. expected methylations values according to predicted linear regression formula of methylation vs. age of healthy group. (A) Linear regression curves of methylation vs. age for the three groups of samples (intron 1-inversion samples, intron 22-inversion samples and healthy controls). Also, equations are shown for every CpG (red, green, and blue for CpGs 1, 2, and 3, respectively) and for the average methylation of three CpGs (in black). (B) Comparison between observed and calculated expected values according to the linear regression equation of healthy controls. T-test p-values showed significance of all CpGs and their average only in the intron 1 inversion group. (C) Comparisons of observed-expected methylation values between inversion groups and healthy controls.

TABLE S1 | Primers list used in this study.

TABLE S2 | Summary values of data presented in Figure 1 together with calculated sensitivity and specificity.

TABLE S3 | Common SNPs around the studied CpG site.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Jamil, Sharma, Nuesgen, Pezeshkpoor, Heimbach, Pavlova, Oldenburg and El-Maarri. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Epigenetic Changes in the Pathogenesis of Rheumatoid Arthritis

Marina V. Nemtsova1,2 \*, Dmitry V. Zaletaev1,2, Irina V. Bure1 , Dmitry S. Mikhaylenko1,2 , Ekaterina B. Kuznetsova1,2, Ekaterina A. Alekseeva1,2, Marina I. Beloukhova1 , Andrei A. Deviatkin1 , Alexander N. Lukashev1,3 and Andrey A. Zamyatnin Jr.1,4 \*

1 Institute of Molecular Medicine, I.M. Sechenov First Moscow State Medical University (Sechenov University), Moscow, Russia, 2 Laboratory of Epigenetics, Research Centre for Medical Genetics, Moscow, Russia, 3 Martsinovsky Institute of Medical Parasitology, Tropical and Vector Borne Diseases, I.M. Sechenov First Moscow State Medical University (Sechenov University), Moscow, Russia, 4 A.N. Belozersky Institute of Physico-Chemical Biology, Lomonosov Moscow State University, Moscow, Russia

#### Edited by:

Yun Liu, Fudan University, China

#### Reviewed by:

Marco Magistri, University of Miami, United States Rowan Hardy, University of Birmingham, United Kingdom

\*Correspondence:

Marina V. Nemtsova<EMAIL_ADDRESS>Andrey A. Zamyatnin Jr<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Epigenomics and Epigenetics, a section of the journal Frontiers in Genetics

Received: 26 March 2019 Accepted: 31 May 2019 Published: 14 June 2019

#### Citation:

Nemtsova MV, Zaletaev DV, Bure IV, Mikhaylenko DS, Kuznetsova EB, Alekseeva EA, Beloukhova MI, Deviatkin AA, Lukashev AN and Zamyatnin AA Jr (2019) Epigenetic Changes in the Pathogenesis of Rheumatoid Arthritis. Front. Genet. 10:570. doi: 10.3389/fgene.2019.00570 Rheumatoid arthritis (RA) is a systemic autoimmune disease that affects about 1% of the world's population. The etiology of RA remains unknown. It is considered to occur in the presence of genetic and environmental factors. An increasing body of evidence pinpoints that epigenetic modifications play an important role in the regulation of RA pathogenesis. Epigenetics causes heritable phenotype changes that are not determined by changes in the DNA sequence. The major epigenetic mechanisms include DNA methylation, histone proteins modifications and changes in gene expression caused by microRNAs and other non-coding RNAs. These modifications are reversible and could be modulated by diet, drugs, and other environmental factors. Specific changes in DNA methylation, histone modifications and abnormal expression of non-coding RNAs associated with RA have already been identified. This review focuses on the role of these multiple epigenetic factors in the pathogenesis and progression of the disease, not only in synovial fibroblasts, immune cells, but also in the peripheral blood of patients with RA, which clearly shows their high diagnostic potential and promising targets for therapy in the future.

Keywords: rheumatoid arthritis, epigenetics, DNA methylation, miRNA, histone modifications, circRNA

### INTRODUCTION

Rheumatoid arthritis (RA) is a chronic auto-inflammatory disease of connective tissue with progressive joint damage and systemic disorders that affects around 1% of the world's population (Cribbs A. et al., 2015). RA can cause various symptoms, clinical forms and prognoses. The incidence of RA begins to increase at the age of 25 years; at the age of 55, the incidence of RA is reaching a plateau (Gabriel, 2001). For example, RA is more than six times greater in 60- to 64 year-old women compared to 18- to 29-year-old women (Melorose et al., 2015). The prevalence of RA varies in different ethnic groups. For example, the incidence of RA among American Indians is 7%, while for some other nations it is 0.2–0.4% (Ferucci et al., 2005). As in most other autoimmune diseases, RA is more common among women than among men in a ratio of 2–3 to 1

(van Vollenhoven, 2009). Based on this fact, there are assumptions that estrogens are actively involved in the pathogenesis of the disease (Wluka et al., 2000).

The etiology of RA remains unknown. It is considered to occur in the presence of genetic predispositions and provoking environmental factors. The heritability of RA has been shown from twin studies to be 60% (Yarwood et al., 2016). Earlier genealogical studies and modern molecular-genetic investigations confirm the involvement of genetic factors in RA development. Accumulation of the disease cases was revealed within families along with an increased risk of RA among firstdegree relatives of the patients (Sparks and Costenbader, 2014).

The pathological process in RA represents an autoimmune inflammation of the synovial membrane of joints with synovial cells proliferation and pannus formation. This tumor-like aggressive granulation tissue promotes articular cartilage erosion and bones destruction. Synovial tissue dysfunction allows macrophages, fibroblasts and activated lymphocytes to penetrate into it. T-lymphocytes produce a variety of proinflammatory cytokines, predominantly belonging to tumor necrosis factor (TNF) and interleukin (IL) superfamilies as well as growth factors (Firenstein et al., 2013). B-lymphocytes are involved in the production of autoantibodies such as rheumatoid factor (RF) and antibodies against cyclic citrullinated peptide (anti−CCP). Differencies in expression of anti−CCP and RF, rate of disease manifestation and variability of response to therapy cause heterogeneity of RA patients indicating different pathophysiological mechanisms implication in the disease development and progression.

Genetic heterogeneity does not explain all the features of RA (Viatte et al., 2013). Thus, investigation of epigenetic factors and mechanisms associated with the progression of the disease and response to treatment is increasingly important. Investigation of the epigenetic landscape can provide novel therapeutic targets (Glant et al., 2014).

Different levels of DNA organization and chromatin packing in the eukaryotic cell's nucleus are points of application for the epigenetic regulation. Epigenetic mechanisms regulate chromatin structure and create clear patterns of gene expression during cell differentiation.

The chromatin structure regulates gene transcription by altering DNA regulatory regions (such as promoters and enhancers) availability for transcription factors (TF). An open chromatin structure, euchromatin, enables DNA-binding proteins and TF to interact with regulatory DNA sequences, leading to active gene transcription. Conversely, heterochromatin is a closed condensed chromatin state, where DNA is tightly bound to protein complexes forming a superspiralized structure. It prevents TF interaction with regulatory sequences, thus inactivating gene expression and leading to its silencing. Transcription factors, non-coding RNAs (ncRNAs), DNA methylation, histone modification and microRNAs (miRNAs) affect gene transcription without changing the DNA sequence itself (reviewed by Golbabapour et al., 2011). Specific epigenetic landscape of chromatin determines a differential gene expression and regulates various cellular processes in physiological and pathological conditions.

Epigenetic changes in RA have been studied both in mononuclear cells of peripheral blood and in different types of immune cells such as monocytes, T-cells and B-cells (Ospelt, 2016). At the same time, the epigenetic modifications in the rheumatoid arthritis synovial fibroblasts (RASFs) are of particular interest because of their aggressive phenotype, which is stable for several passages in cell culture (Hardy et al., 2013). RASFs are clue cells of joint damage and inflammation development due to pro-inflammatory and catabolic molecules synthesis, promoting abnormal proliferation and invasiveness. Implantation of RASFs together with normal human cartilage to immunodeficient mice revealed cell attachment and cartilage destruction without any proinflammatory stimuli. Such behavior was not observed in osteoarthritis (OA) synovial fibroblasts and is presumably related to epigenetic changes in these cells due to the specific pathology only (Lefèvre et al., 2009).

### ABERRANT DNA METHYLATION IN IMMUNE CELLS AND PERIPHERAL BLOOD CELLS IN RA