To understand the map pathways and the genes that were differentially expressed in the OC patient cell lines, we used MetacoreTM software (Calrivate Analytics) to perform enrichment analysis (EA) using a widely known database for protein–protein signaling. EA identified the gene IDs of the potential targets from the DEG sets with gene IDs via the functional ontology function in MetaCore. The possibility of a random intersection of a gene set and the corresponding ontological entities was determined according to the hypergeometric intersection p-value. A reduced p-value suggested that the object was more relevant to the dataset, which indicated that it had a higher rating. Comparative enrichment analysis of the DEG dataset identified the top 10 enriched pathways, and the maps, GO cellular processes, networks, and biomarkers (by disease) are shown in **Figures 5A–D**. These are the most statistically significant data for each category based on a very low p-value. The pictorial representation of the top-scored pathway map (lowest p-value) is based on the distribution of gene enrichment, as shown in **Figure 6A**; similarly, the second scored map (second-lowest p-value) is shown in **Figure 6B**. In **Figures 6A,B**, the well-characterized proteins or protein complexes are displayed as individual symbols; the data from all experiments are shown and linked on the maps as thermometer-like symbols. A red upwardfacing thermometer indicates the upregulated genes, and a blue thermometer indicates the expression level of a downregulated gene. The AN algorithm in GeneGo was used to identify the most relevant biological networks by prioritizing the number of fragments of the canonical pathways in the network, as shown in **Table 3**. The top regulated network processes are presented in **Supplementary Figures 1A,B**, illustrating the two major pathways involving DEGs that were commonly affected in both OC groups. We identified several crucial hub genes, including TCF4, frizzled family proteins (FZD2, FZD8, and FZD6), RUNX2, CDC25 family protein (CDC25A), protein kinase family proteins (CDK2), BRCA1, ATM, and RBBP8. The selected hub genes were mainly involved in the regulation of the canonical Wnt signaling pathway, cell–cell signaling mediated by Wnt, cell cycle phase transition, and the positive regulation of the cell cycle (**Figure 6**, **Supplementary Figures 1A,B**).

## Survival Analysis and Expression Levels of Hub Genes

GEPIA survival assessment was used to investigate the overall association with survival of 10 hub genes from both the lowand high-expression OC groups. As a result, we noticed that the high expression of FZD2 (HR = 0.93) (**Figure 7C**), FZD8 (HR = 0.88) (**Figure 7D**), CDC25A (HR = 0.83) (**Figure 7F**), CDK2 (HR = 0.86) (**Figure 7G**), and RBBP8 (HR = 0.95) (**Figure 7J**) were associated with improved overall survival in the OC cell line. However, the high expression of TCF4 (HR = 1) (**Figure 7A**), FZD6 (HR = 1) (**Figure 7B**), RUNX2 (HR = 1) (**Figure 7E**), BRCA1 (HR = 1.1) (**Figure 7H**), and ATM (HR = 1.2) (**Figure 7I**) were linked with worse overall survival in the OC cell line. Taken together, the results show that FZD6, FZD8, CDK2, and RBBP8 function as core genes that have a close relationship with OC. Furthermore, the GEPIA box plot analysis investigated the level of expression of the core genes in 426 OC tissue samples and 88 normal tissue samples. The boxplot in **Figures 8A–D** shows a considerable increase in the level of core

gene expression (FZD6, FZD8, CDK2, and RBBP8) in the OC cell line.

## DISCUSSION

Microarray technology is one of the most important approaches used by many researchers worldwide to explore the expression levels of genes involved in complex disorders (Russo et al., 2003; Babu, 2004; Perez-Diez et al., 2013). Therefore, studying the expression profiles of DEGs and predicting the target genes of OC is of the utmost importance. In this study, data from a total of three OC cell lines and three OC cell lines with SINE resistance were obtained from the GEO database (GSE126519). A total of 2708 DEGs were screened, including 809 upregulated and 700 downregulated genes. In silico methods have typically shown good efficiency, and networks have been demonstrated to be a reliable way to depict genomic data. The topological interpretation of upregulated and downregulated genes is required for large PPI networks and is thus substantially based on integrated local components, such as the distribution node of the degree, the topological coefficient, the shortest path length distribution, and the betweenness and closeness

centralities (Assenov et al., 2008). These parameters were used to analyze the nodes in individual PPI networks of the DEG dataset to deduce their significance in networks with different characteristics. Furthermore, we implemented GO and KEGG pathway analyses to determine MF, BP, CC, and pathways involving the DEGs using the DAVID online tool. The GO BP terms and KEGG assessment indicated that the upregulated DEGs were enriched primarily in the cell cycle, mitotic cell cycle process, mitotic nuclear division, DNA replication, cell cycle, nucleotide excision repair, DNA mismatch repair, and Fanconi anemia pathways. Interestingly, mutations in mismatch repair (MMR) and Fanconi anemia pathway-related genes in women have been shown to be one of the primary causes of hereditary OC (Norquist et al., 2016). Therefore, our observed results are consistent with the role of upregulated genes in pathways that cause OC. Similarly, the downregulated DEGs were mainly enriched in organ morphogenesis, cell development, cell

organized into a hierarchical tree.

morphogenesis, mesenchymal development and the interaction of the ECM receptor, PI3K-Akt signaling, and arginine and proline metabolism pathways. In line with this, a significant cause of cancer would appear to be the abnormal functioning of the cell cycle and mitosis (Kastan and Bartek, 2004; Malumbres and Barbacid, 2009).

The findings from the STRING, Cytoscape, GO, and KEGG analyses indicated that many pathways were primarily affected in OC. Several studies have used Cytoscape plugins such as MCODE, cytoHubba, CytoCluster, CytoKegg, and CytoNCA to elucidate the core interactions in PPI networks (Lan et al., 2015; Villaveces et al., 2015; Sriroopreddy and Sudandiradoss, 2018; Zhang et al., 2019). To delineate the molecular interactions and pathways identified from the STRING, GO, and KEGG analyses, we utilized GeneGo Metacore, which has a massive amount of information about regulatory and metabolic pathways and contains precisely curated biological networks. To obtain a

FIGURE 6 | (A) The cell cycle metaphase checkpoint pathway. (B) The APC cell cycle regulation pathway. The two top-scored regulated pathways were activated in the OC cell lines. The pathway images were generated by GeneGo MetacoreTM enrichment analysis. Well-characterized proteins or protein complexes are shown as individual symbols within the image; experimental data from all the records are connected and depicted as thermometer-like figures on the maps. Upward-facing thermometers are shown in red and indicate upregulated gene transcripts. The linkage of proteins by arrows depicts the stimulatory and inhibitory effects or interaction of the encoded protein on the desired protein. The hub genes (protein families) that were involved in the canonical signaling pathways are marked in a circle (scarlet). Further explanations are provided at https://portal.genego.com/help/MC\_legend.pdf.

TABLE 3 | Most relevant biological networks were generated using GeneGo Analyze Networks (AN) algorithm.

*In this workflow, the networks were prioritized based on the number of fragments of canonical pathways on the network.*

comprehensive picture of the DEGs involved in OC, we used GeneGo MetacoreTM software to identify the most significant genes and signaling pathways based on the calculated pvalues. Among the top 10 enriched pathways, the cell cycle metaphase checkpoint, APC in cell cycle regulation, chromosome separation, spindle assembly, and DNA damage/ATR regulation of G2/M phase checkpoint/ATM were highly significant in the DEG datasets from both groups (**Figure 5A**). The GO cellular processes showed that the DEGs were enriched in a variety of cellular processes (**Figure 5B**), and these processes are mainly utilized in the enrichment analysis and to prioritize the genes in the constructed networks. The GO process networks were enriched in various groups. Among the top 10 process networks, we selected the four that were the most significant based on the calculated p-values, which included the mitotic cell cycle, the S phase in the cell cycle, and cytoskeletonspindle microtubules (**Figure 5C**). The biomarkers of the diseases distinctly showed that the DEGs with the highest representation in the dataset were also known to contribute to other cancer types (**Figure 5D**). Additionally, there were two top-scored regulated pathways that were activated in the OC cell line that were involved in the cell cycle: the metaphase checkpoint and APC cell cycle regulation pathways (**Figure 6**). Furthermore, we analyzed the biological network of the upregulated DEGs in the signaling pathways by utilizing the AN algorithm in GeneGo. As a result, we determined the two most significant networks that were commonly affected in both of the OC groups. The components of these networks included several crucial hub genes, including TCF4, the frizzled family proteins (FZD2, FZD8, and FZD6), RUNX2, a CDC25 family protein (CDC25A), the protein kinase family proteins (CDK2), BRCA1, ATM, and RBBP8. Among them, the TCF4, frizzled, and RUNX2 genes are primarily involved in the regulation of the canonical WNT signaling pathway (58%) and cell–cell signaling mediated by WNT (56%). Genes such as CDC25A, CDK2, BRCA1, ATM, and RBBP8 are mostly involved in the mitotic cell cycle process (53.1%), mitotic cell cycle phase transition (42.9%), positive regulation of the cell cycle (46.9%), and cell cycle phase transition (42.9%). Finally, the GEPIA web server was used to assess the association between hub gene expression and OC prognosis. The overall survival analysis indicated that high expression of FZD2, FZD8, CDC25A, CDK2, and RBBP8 were associated with better survival, and high expression of TCF4, FZD6, RUNX2, BRCA1, and ATM were associated with decreased survival in the OC cell line. Collectively, FZD6, FZD8, CDK2, and RBBP8 were identified as core genes that were strongly associated with overall survival in OC. Therefore, these four genes could contribute to OC metastasis.

**Supplementary Figure 1A** shows that the frizzled family of proteins is involved in canonical Wnt signaling pathway regulation. FZD6, also known as frizzled class receptor 6, is a member of the "frizzled" gene family, which consists of 7 transmembrane domain proteins that are Wnt signaling protein receptors. Many studies have observed through mutagenesis experiments that several residues in the intracellular loops and the C-terminus of FZD play a prominent role in signaling (Cong et al., 2004; Wallingford and Habas, 2005). Kim et al. (2015) found that the expression of FZD6 was increased in colorectal cancer (CRC) patients when compared to that in nontumor tissues. Furthermore, they discovered that FZD6 expression in CRC cells was negatively regulated by miR199a-5p (Kim et al., 2015). In recent research, Corda et al. (2017) observed that the Wnt receptor-encoding gene FZD6 is often duplicated in breast cancer and confers a higher risk of triple-negative breast cancer. For the assembly of the fibronectin matrix, FZD6 signaling is intrinsically required and interferes with actin cytoskeletal organization. The researchers concluded that in highly metastatic forms of breast cancer, such as TNBC, the FZD6-fibronectin actin axis could be targeted for drug treatment (Corda et al., 2017). In our study, we observed the overexpression of frizzled class receptor 6 in OC cell lines, and the overexpression of FZD6, which acts as an adverse prognostic factor, was associated with decreased survival in OC patients.

Frizzled class receptor 8 is also a "frizzled" gene family member that serves as a Wnt signaling protein receptor (Bhanot et al., 1996). Most frizzled receptors are also associated with the canonical signaling pathway of beta-catenin (Dann et al., 2001). Li et al. (2017) reported a higher level of FZD8 expression in

bone metastases in prostate cancer (PCa), which is frequently diagnosed among men. This research group also found that the silencing of FZD8 suppressed the migration and invasion of cells and the occurrence of PCa bone metastasis in vitro and in vivo by activating the canonical β-catenin/Wnt signaling pathway, and the data suggest that FZD8 could be a potential therapeutic target for the treatment of bone metastasis in PCa (Li et al., 2017). Chakravarthi et al. (2018) reported that ETSrelated gene (ERS) specifically targets and activates FZD8 directly by binding to its promoter region rather than ETV1 and suggested that the overexpression of ERG in PCa leads to FZD8 induction and the activation of the Wnt pathway (Chakravarthi et al., 2018). The research group led by He et al. recently found that miR-520b overexpression results in the inhibition of cell proliferation, migration, and invasion in human spinal osteosarcoma (OS) tissues and cell lines by inactivating the Wnt/β-catenin signaling pathway through the downregulation of FZD8 and thus provides a new spinal OS therapeutic target (Wang et al., 2017). Similarly, Liu et al. (2019) reported a reduced level of miR-99b-5p in non-small cell lung cancer (NSCLC) cell lines. They validated FZD8 as a specific target of miR-99b-5p and found that increased expression of miR-99b-5p inhibited NSCLC proliferation, migration, and invasion in vitro (Liu et al., 2019). The findings of our study suggest that the overexpression of FZDs in OC results in the anomalous activation of the canonical Wnt signaling pathway and may increase their function during the development of OC.

As seen in **Figure 6B**, CDK2 is mainly involved in the APC cell cycle regulation pathway, and the overexpression of CDK2 results in the upregulation of the G1/S phase transition, resulting in cancer cell proliferation. CDK2 and other relevant genes that are upregulated in OC (red circles in **Supplementary Figure 1B**) are shown in the closed network of the APC cell cycle regulation pathway; because of the increase in CDKs, APC failed to inactivate the CDK complexes by inducing their degradation. Liu et al. (2011) revealed that CDK2 expression was significantly higher in laryngeal squamous cell cancer tissues when compared to that in paired adjacent normal laryngeal tissues (Liu et al., 2011). Duong et al. (2012) reported that low-molecular-weight cyclin E (LMW-E) required kinase activity associated with CDK2 to induce the formation of mammary tumors by disrupting the growth of acinar cells. They used a combination of therapy with a CDK inhibitor (roscovitine) plus a b-Raf-targeting pan-kinase inhibitor (sorafenib) or an mTOR inhibitor (rapamycin) to arrest

the G1/S cell cycle in breast cancer cells; thus, the b-Raf-ERK1/2 mTOR signaling pathway could be suppressed (Duong et al., 2012). Kanwal et al. (2016) found that the expression of CDK2 is significantly increased in OC tissues when compared to that in normal ovarian tissues (Kanwal et al., 2016). The pathways involving cyclin-dependent kinase (CDK) are significant and well-established cancer treatment targets. The role of CDK2 remains controversial in several cancer types (McCurdy et al., 2017). Many studies have suggested that CDK2 could be a crucial factor in the progression of cancer by regulating several pathways and might be a prospective biomarker and indicator of prognosis. Therefore, CDK2 and its cyclin binding partners are possible therapeutic targets for future cancer treatments (Yin et al., 2018; Zhang et al., 2018; Wood et al., 2019).

RBBP8, also known as retinoblastoma-binding protein 8, aids in regulating cell proliferation and DNA repair by homologous recombination (Fusco et al., 1998).  The research group led by Rose et al. suggested that RBBP8 was significantly hypermethylated in bladder cancer (BLCA) and was associated with more prolonged overall survival, and they indicated that it may be used as a complementary marker for the detection of BLCA in urine (Mijnes et al., 2018). Miao et al. (2019) reported that the downregulation of the long noncoding RNA (lncRNA) cancer susceptibility candidate 2 (CASC2, enhanced tumor development, increased miR-18a-5p levels, and reduced the expression of RBBP8 in nasopharyngeal carcinoma (NPC). The upregulation of CASC2 resulted in decreased proliferation and increased apoptotic cell death in vivo (Miao et al., 2019). Our data clearly showed that RBBP8 was differentially expressed and involved in cell cycle regulation (**Supplementary Figure 1B**). Additionally, it contributes to the development of OC in both groups. However, the role of RBBP8 in OC is unclear and requires further research.

Furthermore, we conducted an interrelation analysis of the identified hub genes to elucidate the interactions between them, primarily among the genes that interacted with one another directly or indirectly. As shown in **Figure 9**, the mitotic cell cycle phase transition pathway interacts with the regulation of G1/S phase transition and APC in the cell cycle regulating pathway via the essential genes RBBP8, BRCA1, CDC25A, ATM, and CDK2 (D'Andrilli et al., 2004; Soria-Bretones et al., 2013; Xiao et al., 2019). In contrast, frizzled family proteins (FZD2, FZD6, and FZD8) are directly involved in the Wnt signaling pathways because they are receptors of Wnt proteins (Janda et al., 2012). The TCF4 and RUNX2 genes are involved directly or indirectly in the Wnt signaling network, resulting in tumorigenesis (Gaur et al., 2005; Hrckulak et al., 2018; Komori, 2019). Taken together, these findings showed that node genes involved in the development of OC could be significant factors in cell cycle regulation and the Wnt signaling pathway.

Overall, our systematic bioinformatics assessment demonstrated that DEGs might play a pivotal role in the incidence, prognosis, growth, and development of OC. In this study, a total of 2708 DEGs and 10 hub genes were identified, and FZD6, FZD8, CDK2, and RBBP8 could be the core genes involved in OC and SINE-resistant OC. Expression analysis and the correlation of the multiple genes will undoubtedly aid in the understanding of the roles of such genes in the growth and development of OC. Several research groups have demonstrated that preclinical models have showed some success in reducing tumor growth and decreasing the side effects of existing chemotherapy drugs (Cicenas et al., 2015; Whittaker

were involved in the mitotic cell cycle phase transition pathway; and *CDC25A* and *CDK2* were involved in the APC cell cycle regulation pathway.

et al., 2017; Xia et al., 2018). We need to conduct a series of experimental studies to prove this hypothesis to obtain more precise correlation reports. However, the findings from this study could enhance the understanding of the molecular pathogenesis of OC. Furthermore, the core genes and pathways might be potential biomarkers that could be used for the detection and targeting of OC and SINE-resistant OC cells for therapy.

## DATA AVAILABILITY STATEMENT

All datasets generated for this study are included in the article/**Supplementary Material**.

## AUTHOR CONTRIBUTIONS

SK, DK, and CD were involved in the design of the study and the acquisition, analysis, and interpretation of the data. SK, DK, CD, and HZ were involved in the interpretation of the data and drafting the manuscript. CD, RS, and HZ supervised the entire study and were involved in study design,

## REFERENCES

the acquisition, analysis, interpretation of the data, and drafting the manuscript. The manuscript was reviewed and approved by all the authors.

## FUNDING

This work was supported by Qatar University Grant# QUST-2- CHS-2019-3.

## ACKNOWLEDGMENTS

The authors would like to take this opportunity to thank the management of Vellore Institute of Technology for providing the necessary facilities and encouragement to carry out this work.

## SUPPLEMENTARY MATERIAL

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

**Conflict of Interest:** 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 Kumar, Kumar, Siva, Doss and Zayed. 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.

# Identification of Biological Pathways Contributing to Marbling in Skeletal Muscle to Improve Beef Cattle Breeding

Zahra Roudbari 1,2\*, Susan L. Coort <sup>2</sup> , Martina Kutmon2,3, Lars Eijssen2 , Jonathan Melius <sup>2</sup> , Tomasz Sadkowski <sup>4</sup> and Chris T. Evelo2,3

<sup>1</sup> Department of Animal Science, Faculty of Agriculture, University of Jiroft, Jiroft, Iran, <sup>2</sup> Department of Bioinformatics-BiGCaT, NUTRIM School of Nutrition and Translational Research in Metabolism, Maastricht University, Maastricht, Netherlands, <sup>3</sup> Maastricht Centre for Systems Biology (MaCSBio), Maastricht University, Maastricht, Netherlands, <sup>4</sup> Department of Physiological Sciences, Faculty of Veterinary Medicine, Warsaw University of Life Sciences - SGGW, Warsaw, Poland

### Edited by:

Jyoti Sharma, Institute of Bioinformatics (IOB), India

### Reviewed by:

Duo Xu, Cornell University, United States Abhishek Kumar, University of Kiel, Germany Bin Tong, Inner Mongolia University, China

\*Correspondence: Zahra Roudbari<EMAIL_ADDRESS>
### Specialty section:

This article was submitted to Bioinformatics and Computational Biology, a section of the journal Frontiers in Genetics

Received: 04 September 2019 Accepted: 16 December 2019 Published: 07 February 2020

### Citation:

Roudbari Z, Coort SL, Kutmon M, Eijssen L, Melius J, Sadkowski T and Evelo CT (2020) Identification of Biological Pathways Contributing to Marbling in Skeletal Muscle to Improve Beef Cattle Breeding. Front. Genet. 10:1370. doi: 10.3389/fgene.2019.01370 Red meat is an important dietary source that provides part of the nutritional requirements. Intramuscular fat, known as marbling, is located throughout skeletal muscle. Marbling is a trait of major economic relevance that positively influences sensory quality aspects. The aim of the present study was to identify and better understand biological pathways defining marbling in beef cattle. Pathway analysis was performed in PathVisio with publicly available transcriptomic data from semitendinosus muscle of well-marbled and lean-marbled beef. Moreover, for Bos taurus we created a gene identifier mapping database with bridgeDb and a pathway collection in WikiPathways. The regulation of marbling is possibly the result of the interplay between signaling pathways in muscle, fat, and intramuscular connective tissue. Pathway analysis revealed 17 pathways that were significantly different between well-marbled and lean-marbled beef. The MAPK signaling pathway was enriched, and the signaling pathways that play a role in tissue development were also affected. Interestingly, pathways related to immune response and insulin signaling were enriched.

Keywords: marbling, curation pathway for cow, signaling pathway, improve breeding selection, transcriptomics profiling

## INTRODUCTION

Red meat is as an important dietary source that provides part of the nutritional requirements such as proteins, minerals, B-complex vitamins, and essential fatty acids (McAfee et al., 2010). Control of meat quality is very important for meat producers and meat sellers to satisfy customer's preferences (Bernard et al., 2007). Marbling, a trait that describes the presence intramuscular fat, is of major economic relevance for beef producing cattle that has a positive impact on sensory quality traits, such as flavor, juiciness, and tenderness of meat. Studies have shown that marbling depends on factors such as breed, genotype, age, diet, husbandry, and growth stages. Although in marbling the environmental factors play an important role, the genetic background of the animals is the major factor defining the marbling status (Yamada et al., 2006). O'Connor et al. studied the effect of breedtype on marbling, their results demonstrated that increase in meat marbling from Bos taurus cattle (Hereford, Red Angus, Angus, and Tarentaise breeds) can increase the tenderness, more than the Bos

**48**

indicus cattle (Braford, Red Brangus, and Simbrah breeds) (O'Connor et al., 1997). Although, Shackelford et al., found that, the lower tenderness of meat from Bos indicus cattle is mainly because of decreased postmortem proteolysis which result from elevated calpastatin activity that the possibility existed for an interaction between breed and the influences of marbling score on tenderness (Shackelford et al., 1991). Also, Wulf et al., reported that, in Charolais and Limousin breeds, marbling correlated with calpastatin activity and shear force. He suggested that selection for increased marbling based on these genetic effects, in those two breeds, might be effective for enhancing beef tenderness (Wulf et al., 1996). There is substantial evidence from transcriptomics studies that gene expression profiles affect phenotypic variation for marbling (Cesar et al., 2015). Understanding the signaling pathways that make up the regulatory network in the marbling process can help steer the breeding process (Thaller et al., 2003). Therefore, animal breeding specialists have attempted system-oriented approaches to investigate major economic traits (Lee et al., 2010). It has been shown that marbling differences may be a function of a number of complex interactions among biological pathways. Therefore, a pathway analysis with differential gene expression patterns can result in a better understanding of muscle physiological states and their influence on beef quality and animal welfare (Hocquette et al., 2012). The public availability of transcriptomics data from beef producing cattle, provides new opportunities to explore the global gene expression in muscle to investigate physiological processes and their influence on meat sensory quality traits (Lee et al., 2010; Hocquette et al., 2012).

Within the genomic region of marbling there are several genes considered as parts of QTLs such as EDGPR1, Titin, Akirin 2, and RPL27 (Takasuga et al., 2007) which were mapped in a half-sib family of Japanese Black cattle (Yamada et al., 2006). Thus, these genes were considered as positional functional candidates for the genes responsible for marbling. This study aims at identifying genes and biological pathways regulating marbling of muscle tissue in beef cattle based on publicly available transcriptomics data obtained from a study by Sadkowski and coworkers (2014). We updated and extended the pathway collection for B. taurus at WikiPathways (Slenter et al., 2017) an online pathway repository, and a B. taurus gene product identifier mapping BridgeDb database was created to allow mapping of expression data to the gene databases identifiers used in the pathways (van Iersel et al., 2010). Sadkowski et al., 2014 measured global gene expression in skeletal muscle of three cattle breeds, i.e., Limousin, Holstein-Friesian, and Hereford, using Agilent microarray chips. Pathway and network analysis were performed to select the important biological pathways involved in marbling and their interactions.

## MATERIALS AND METHODS

## Transcriptomics Data Set

The study by Sadkowski et al., 2014 compared gene expression in semitendinosus skeletal muscle of well-marbled beef (HolsteinFriesian and Hereford) versus lean-marbled beef (Limousin). Their publicly available microarray data set was used in the present study (NCBI GEO GSE46411). The Holstein-Friesian, Hereford, and Limousin groups consisted of four animals each. Samples for total RNA isolation were taken instantaneously after slaughter from semitendinosus muscle and were kept in liquid nitrogen for transportation and then at −80°C until analyzed. Quality of RNA samples was evaluated using Bioanalyzer 2100 (Agilent Technologies, USA). Only samples with RIN ≥ 8 were further analyzed (Sadkowski et al., 2014).

## Agilent Microarray Data Analysis

Global gene expression was measured with Agilent Two-Color Mi Bovine (V2) 4 x 44K Gene Expression Microarray oligonucleotide slides (Agilent, USA). Sadkowski and coworkers checked the quality of the data and performed LOWESS normalization. The normalized transcriptomic data compared well-marbled beef Holstein-Friesian (n = 4) or Hereford (n = 4) to lean-marbled beef (n = 4). The four log10 fold change (log10FC) values for each group comparison were averaged to obtain an estimate of the 10logFC between the entire groups. Furthermore, a one-sample t-test was performed on both sets of four values, comparing those to 0 (giving a p-value indicating the significance of these values being different from 0 = no change).

## B. taurus Pathway Collection

The online biological pathway repository, WikiPathways (Slenter et al., 2017), contains pathways of different species, however a B. taurus collection was missing. We updated and extended the pathway collection for B. taurus. We also created a B. taurus gene identifier (ID) mapping database based on mappings present in the Ensembl-based BridgeDb framework (van Iersel et al., 2010). The newly created B. taurus ID mapping database was used to annotate genes and proteins in pathways from WikiPathways and to perform pathways analysis. A online and freely available version of the database for the Ensembl build 85 is accessible at (http://bridgedb.org/data/gene\_database/archive/r85/Bt\_Derby\_ Ensembl\_85.bridge.zip). Second, the WikiPathways homology based the homology mapper which is available at GitHub (https://github.com/PathVisio/homology.mapper) was updated to improve homology coverage for gene products that were annotated with different data sources. The pathways were converted from human pathways, with a required minimum successful conversion of at least 50% of the original human genes. Third, we manually curated all converted pathways to check whether the genes were correctly annotated and pathways are relevant in B. taurus. Finally, new pathways directly derived from cow breeding literature and not present in the WikiPathways collection were designed in PathVisio (v3.2.0) (Kutmon et al., 2015), the pathway creation, visualization, and analysis tool. All pathways were uploaded in gpml format to WikiPathways using the WikiPathways plugin (https://www. pathvisio.org/plugin/wikipathways-plugin/) for PathVisio.

## Pathway-Based Over-representation Analysis

To analyze and visualize the molecular changes in marbling at biological process level a pathway-based over-representation analysis was performed in PathVisio (v3.2.0) (Kutmon et al., 2015). The B. taurus WikiPathways pathway collection, containing 286 pathways (6/30/2015), and the B. taurus ID mapping database, was used in the analysis. The pathways are ranked based on a standardized difference score (Z-score) based on the expected value and standard deviation of the number of significantly (p ≤ 0.05) and differentially (absolute FC ≥ 1.3) expressed genes in a pathway. Biological pathways significantly changed when (i) Z-score > 1.96, (ii) permuted p-value < 0.05 and (iii) minimum number of changed genes is 3. Additionally, alterations in gene expression (log10FC and p value) when comparing Hereford to Limousin were visualized on the B. taurus pathways with PathVisio.

## Gene Ontology Overrepresentation Analysis

To find the biological processes in which differentially expressed genes were over represented while no pathways for these processes were present in the B. taurus collection at WikiPathways we performed Gene Ontology (GO) analysis via the GO-Elite webinterface (Zambon et al., 2012). GO-Elite is a flexible tool for GObased over-representation analysis. To identify GO processes the following settings in GO-Elite were used: (i) 2000 permutations, (ii) Z-score GO pruning algorithm, (iii) Z-score threshold >1.96, (iv) p-value threshold <0.05 and (v) minimum number of changed genes is 3 (apart from the method specific permutations those are the same criteria as used for the pathway analysis). This approach not only helps to unify the characteristics and functions of the genes but also to attain a broader perspective of the muscle physiological processes and their influence on meat quality.

## Integrated Network Analysis

To visualize the pathway and GO analysis results and their interactions the network analysis and visualization tool, Cytoscape (version 3.2.0), was used (Shannon et al., 2003). First, all enriched pathways and the differentially expressed genes present in these pathways were selected.  Third, both results were combined into one network showing the interaction between pathways and GO classes based on corresponding differentially expressed genes. Finally, differences in gene expression between well-marbled and lean-marbled skeletal muscle were visualized in the network.

## RESULTS

## Identification of Differentially Expressed Genes Between Well-Marbled and Lean-Marbled Skeletal Muscle

In the selected transcriptomic data set of beef marbling 42,990 microarray reporters were measured in both lean marbling beef (Limousin) and well marbling beef (Hereford and Holstein-Friesian) animals. Statistical analysis was performed on 29,677 reported genes that remained from the 42,990 reporters after quality control and annotation with Ensembl gene IDs. In the Hereford breed compared to the Limousin breed, 1,513 were higher expressed and 1,556 lower expressed (absolute log10FC >0.11 and p-value <0.05). When comparing the Holstein-Friesian breed to the Limousin breed, 1,772 genes were higher expressed and 2,458 lower expressed in the Holstein Friesian breed. The genes that met these criteria were used for further analysis.

## Creating of B. taurus Pathway Collection and Pathway Design

In total, 282 human pathways were converted from human pathways to cattle pathways. All these pathways were manually checked and are available at (https://www.pathvisio.org/ downloads/download-pathways/). Moreover, 4 pathways were newly created based on the bovine breeding literature: Growth hormone signaling (WP2890) (Roudbari Z and Kutmon M: Growth Hormone (GH) Signaling (B. taurus); (https://www. wikipathways.org/instance/WP2890), Growth hormone receptor signaling (WP2891) (Roudbari Z, Hanspers K, Evelo C, Kutmon M: Growth Hormone Receptor (GHR) Signaling (B. taurus) (https://www.wikipathways.org/instance/WP2891), IGF1-signaling (WP2892) (Roudbari Z, Evelo C, Willighagen E, Mélius J, Hanspers K, Kutmon M: IGF1-signaling (B. taurus) (https://www.wikipathways.org/instance/WP2892), and Gonadotropin- releasing hormone signaling (WP2901) (Roudbari Z, Kutmon M, Pico A, Willighagen E, Mélius J: Gonadotropin-releasing hormone (GNRH) signaling (B. taurus) (https://www.wikipathways.org/instance/WP2901). As an example, the newly designed GNRH signaling pathway is shown in Figure 1. The elements of this process are the key factors stimulating gonadotropin release from the pituitary, which controls the release of luteinizing hormone and folliclestimulating hormone, and reproductive development in mammals.

## Pathway Analysis

When comparing Hereford with Limousin breed ten biological pathways that were formerly known to be involved in marbling (Cui et al., 2012; Lim et al., 2013; Silva-Vignato et al., 2017) were found to be significantly enriched in differentially expressed genes (z-score > 1.96) (Table 1). Four biological processes such as: The Hypertrophy Model, P38 MAPK signaling, IL-1 signaling, and insulin signaling pathways, which are known to be important in marbling development are described in more detail, and the pathways are shown in Results section. Interestingly, some pathways not yet known to play a role in marbling were also found to be enriched in differentially expressed genes when comparing Hereford with Limousin breed. These included histone modifications and vitamin D metabolism pathways, in addition to Hereford breed that is used for meat production, the Holstein-Friesian breed that is a dairy cattle, was also

compared with the Limousin breed. Pathway analysis revealed that some but not all of the marbling related pathways found for the comparison between Hereford breed and Limousin breed were also found when comparing Holstein-Friesian breed with Limousin. Examples of consistently affected pathways are the P38 MAPK signaling and the Hypertrophy Model pathways (Table 2).

## Hypertrophy Model Pathway

Muscle hypertrophy is known to increase the muscle mass, and is determined by increased protein mass per fiber which results from an increase of protein synthesis (Glass, 2005). In the B. taurus Hypertrophy model (http://www.wikipathways. org/instance/WP982) the overall gene expression was higher in skeletal muscle of Hereford and Holstein-Friesian compared to Limousin (Figure 2). The expression of the Il1a, Ifrd1, Cyr61, ATF3, and Ankrd1 genes were significantly higher in Hereford in this pathway and the Il18, Eif4ebp1, and Il1r1 genes were significantly lower in the model (p-value < 0.05). Among them was IL-1 which plays a significant role in lipid metabolism by regulating insulin levels under physiological conditions (Matsuki et al., 2003); Atf3 which works together with p38c in a common pathway in the intestine to regulate lipid metabolism and immune homeostasis (Chakrabarti et al., 2014); and Frd1, Cyr61, and Ankrd1 genes. Two of the seven genes were significantly lower expressed including Ef4ebp1 contributing to the development of obesity through increased adipogenesis and fat metabolism alterations (Le Bacquer et al., 2007) and Il18 gene. TNF and interleukin (IL)-1 may cause negative inotropic effects indirectly through activation or release of IL-18 (Mann, 2015).

## p38 MAPK Signaling Pathway

The p38 Mitogen-activated protein kinase (p38 MAPK) signaling pathway has found to be responsible for transduction of extracellular signals to their intracellular targets in different types of cells, including skeletal muscle cells and which leads to several biological effects for example proliferation, differentiation, migration, growth, apoptosis, and more specifically to muscle cells, hypertrophy (Yu et al., 2010; Silva-Vignato et al., 2017). The p38 MAPK is one intracellular signaling pathway activated during the differentiation of myogenic cell lines and this pathway is a chief regulator of skeletal muscle development (Keren et al., 2006). The p38 MAPK signaling pathway is a well-known pathway that affects lipid metabolism (Zhang and Liu, 2002). In the B. taurus p38 MAPK signal pathway five of the seven genes present were significantly higher expressed in Hereford and Holstein-Friesian compared to Limousin (Figure 3).

## IL-1 Signaling Pathway

The IL-1 signal pathway is a major mediator of innate immune reactions.This pathway regulates extracellular and intracellular signaling of IL-1a or IL-1b including positive and negativefeedback mechanisms which strengthen or terminate the IL-1 response. In reply to ligand binding of the receptor, a complicated sequence of combinatorial phosphorylation and

\*Pathways previously known to be related to marbling.

The pathways are ranked based on Z score. Per pathway the following is listed; Positive = amount of genes differentially expressed, Measured = the amount of genes measured in the study, Total = the amount of genes in the pathway and P-value = the significance level.

ubiquitination events lead to activation of nuclear factor kB signaling and the JNK and p38 mitogen-activated protein kinase pathways (Weber et al., 2010). The members of the B. taurus IL-1 signaling pathway (http://www.wikipathways.org/instance/ WP3271), such as IL-1a, IL-1b, MAP3K1, UBE2N, MAPK14, REL, ATF2, and JUN were significantly up-regulated in Hereford breed versus Limousin breed (Figure 4). Among them IL-1a which was found to play a role as an inhibitor of the expression of peroxisome proliferator-activated receptor gamma (PPARG), a key transcriptional factor for adipocytes differentiation, (Um et al., 2011). IL-1b has been reported to inhibit adipocyte differentiation from preadipocytes and to reduce the lipid content in mature adipocytes (Simons et al., 2005). Some of aforementioned genes: MAP3K1, UBE2N, MAPK14, REL, ATF2 can directly bind to the peroxisome proliferator-activated receptor promoter and activate transcription to regulate adipocyte differentiation (Maekawa et al., 2010). The significantly down-regulated genes are: RELA, MAPK1, IKBKG, MAP2K4, SQSTM1, and IL1R1 (Figure 4). Some of them are known to participate in lipid metabolism processes; activation of p62/SQSTM1 and peroxisome proliferatoractivated receptor gamma is induced by palmitate

\*Pathways previously known to be related to marbling.

The pathways are ranked based on Z score. Per pathway the following is listed; Positive = amount of genes differentially expressed, measured = the amount of genes measured in the study, Total = the amount of genes in the pathway and P-value = the significance level.

internalization, which triggers lipid metabolism and limits inflammation (Krausgruber et al., 2011).

## Insulin Signaling Pathway

Genes engaged in the insulin signaling pathway regulate several aspects of cellular function, including most notably the regulation of cellular growth and maintaining glucose homeostasis (DeBosch and Muslin, 2008). For Hereford vs. Limousine comparison, twenty six genes present in the B. taurus Insulin signaling pathway (http://www.wikipathways.org/instance/WP966) showed significant expression differences. Fourteen genes were identified as up-regulated in Hereford breed (SOS2, PIK3R3, PIK3CA, PIK3C2A, CBLB, CBLC, SNAP25, JUN, EGR1, MAP3K1, MAPK14. ENPP1, and XBP1), and twelve were down-regulated in Hereford breed (PFKM, PFKL, ARF1, STXBP2, EIF4EBP1, PIK3CD, GAB1, IGF1R, MAPK1, MAPK13, MAP2K4, and ELK1) when compared to Limousin breed (Figure 5).

FIGURE 2 | Skeletal muscle gene expression in Hereford and Holstein-Friesian vs Limousin visualized on the Hypertrophy model. In the hypertrophy model from WikiPathways (WP982) the changes in gene expression between Hereford and Holstein-Friesian with Limousin in skeletal muscle are visualized. The logFC (Hereford and Holstein-Friesian vs Limousin) is indicated with a color gradient (blue to red over white), i.e., blue represents a negative value (= lower expressed in Hereford and Holstein-Friesian) and red a positive value (= higher expressed in Hereford and Holstein-Friesian). The p-value is colored based on a rule, i.e. p-value < = 0.05 (= significant) is shown in green and p-value > 0.05 in white.

FIGURE 3 | Skeletal muscle gene expression in Hereford and Holstein-Friesian vs Limousin visualized on the p38 MAPK signaling. In the p38 MAPK signaling pathway from WikiPathways (WP1037) the changes in gene expression between Hereford and Holstein-Friesian with Limousin in skeletal muscle are visualized. The logFC (Hereford and Holstein-Friesian vs Limousin) is indicated with a color gradient (blue to red over white), i.e., blue represents a negative value (= lower expressed in Hereford and Holstein-Friesian) and red a positive value (= higher expressed in Hereford and Holstein-Friesian). The p-value is colored based on a rule, i.e., p-value < = 0.05 (= significant) is shown in green and p-value > 0.05 in white.

FIGURE 4 | Skeletal muscle gene expression in Hereford and Holstein-Friesian vs Limousin visualized on IL-1 signaling. In the IL-1 signaling pathway from WikiPathways (WP3271) the changes in gene expression between Hereford and Holstein-Friesian with Limousin in skeletal muscle are visualized. The logFC (Hereford and Holstein-Friesian vs Limousin) is indicated with a color gradient (blue to red over white), i.e., blue represents a negative value (= lower expressed in Hereford and Holstein-Friesian) and red a positive value (= higher expressed in Hereford and Holstein-Friesian). The p-value is colored based on a rule, i.e., p-value < = 0.05 (= significant) is shown in green and p-value > 0.05 in white.

## Gene Ontology Analysis

Pathway analysis gave an insight in the biological processes involved in marbling. However, only 69% of all measured genes are present in the investigated pathways from the B. taurus WikiPathways collection. In order to obtain a better insight in the biological role of the differentially expressed genes not present in WikiPathways a GO analysis was performed. (Table 3). This approach not only helps to unify the characteristics and functions of the genes but also to attain a broader perspective of the muscle physiological processes and their influence on meat quality related to marbling.

## Integrated Network of Altered Pathways With GO-Terms

The significant pathways and GO terms were merged together and shown in Figure 6. Some of the highly connected nodes are IL1a,

TGFB2, PAK1, (as Pak1 deficiency led to upregulation of reverse cholesterol transporters in ApoE−/− mice in response to Western diet feeding, it might be suggested that Pak1 exerts a negative modulatory influence on these transporters and thereby might promote lipid retention in inflamed arteries which cause atherogenesis (Singh et al., 2015), TGFB2 (that TGF-b2 might control adipocyte differentiation in bone marrow stromal cells in vivo by inducing PPARg phosphorylation. Whether Smad activation induced by TGFb2 might play a role with MAPK in the inhibition of adipocyte differentiation induced by TGFb<sup>2</sup> in vivo requires more investigation (Ahdjoudj et al., 2005), IL1b, UBE2N, IRAK2, RCAN1, LRP6 (Treatment of LRP6 knockdown-Human mesenchymal stem cells with adipogenic supplements led to the accumulation of fat vacuoles, which was demonstrated by Oil Red O staining (Peröbner et al., 2012), IKBKG, TGFB3, RELA, BCL10, FGF2 (FGF-2 treatment of human preadipocytes also

TABLE 3 | The enriched processes found by GO-Analysis. A description of the process together with the positive gene number, Z Score, and P-value are given.

resulted in increased adipocyte differentiation, suggesting that this feature might be common to members of the fibroblast growth factor family, although FGF-1 was consistently the more potent adipogenic agent, particularly in cells from subcutaneous depots (Hutley et al., 2004), IGF1R, and TLR4 (TLR4 knockdown in H9C2 cardio myocytes decreases fatty acid-induced lipid accumulation (Dong et al., 2012) which are present in at least two different pathways and GO terms.

## DISCUSSION

The aim of this study was to determine transcriptional profiles of high and low marbled beef with a focus on pathways of muscle cell origin that might play a role in the regulation of marbling development. The regulation of marbling is suggested to be the result of interaction of signaling pathways in muscle, fat, and intramuscular connective tissue (Hocquette, 2010). Identifying these processes with pathway analysis can help to decipher the key processes involved in marbling development. Pathway analysis revealed 17 pathways that were significantly different (z-score > 1.96) between well-marbled and lean marbled breeds. P38 MAPK signaling pathway well known to affects lipid metabolism and muscle development, was enriched when we compared gene expression in well and low marbling breeds. In addition, the signaling pathways "Hypertrophy Model", "MicroRNAs in cardiomyocyte hypertrophy" and "Physiological and pathological hypertrophy of the heart" that play a role in tissue development were affected. Interestingly, the analyses also demonstrated that pathways related to immune response (IL signaling, TCR signaling, and Toll-like receptor signaling pathways) and insulin signaling, mitochondrial gene expression and vitamin D metabolism were enriched and might act together with pathways related to lipid metabolism. We explored regulatory pathways that control gene expression in bovine muscle and the relationships between gene expression and the marbling trait to identify markers that effect on marbling. A similar study done by (Hong et al., 2014) investigated the biological characteristics of differentially expressed genes in high marbled muscle in pig compared to a low marbled muscle. They indicated that the differentially expressed genes were clustered to three group related to energy metabolism, protein synthesis, and immune response in high marbling pigs. These finding suggested that the genes related to energy metabolism, protein synthesis, and immune response contribute to growth performance and meat quality. Our results also showed differentially expressed genes take part in these processes. The hypertrophy model pathway was found to be enriched with the highest Z score in the present study and during muscle hypertrophy there is an equilibrium between protein synthesis and degradation that might bring about protein deposition, and hence muscle growth. Together, these processes will lead to differences in muscle and fat deposition, and for this reason animals have different proportions of ribeye area and back fat thickness (Silva-Vignato et al.,2017). Mitogenactivated protein kinase (MAPK) signals have been shown to play a significant role in intracellular signaling associated with a variety of cellular activities including cell proliferation, differentiation, survival, and death (Yu et al., 2010). In mammalian cells, three MAPK families have been characterized: classical MAPK (also known as ERK), C-Jun Nterminal kinase/stress activated protein kinase and p38 kinase pathways (Zhang and Liu, 2002). Each mammalian MAPK signaling route comprises at least three components: a MAPK kinase kinase (MAP3K), a MAPK kinase (MAP2K), and a MAPK. Activated MAPKs phosphorylate various substrate proteins including transcription factors such as ATF2 and Jun (Kim & Choi, 2010). Philip and coworkers (Philip et al., 2005) discovered that the p38 MAPK played a key role in GDF-8 induced inhibition of proliferation and upregulation of the cyclin kinase inhibitor p21. In addition, their results showed a functional link between the p38 MAPK and GDF-8-activated Smad pathways, and identify an important role for the p38 MAPK in GDF-8's function as a negative regulator of muscle growth (Philip et al., 2005). In comparative muscle transcriptome associated with carcass traits of Nellore cattle, Silva-Vignato and colleagues indicated that MAPK signaling pathway involved in muscle and fat deposition, which are economically important carcass traits for beef production (Silva-Vignato et al., 2017). The third pathway found in the present study was IL-1, the IL-1 family of cytokines includes 11 proteins encoded by 11 different genes and gene regulation of IL-1 signal is activation of MKK4, MKK3, and MKK6 gene that activate NF-kB and p38 MAPK pathways (Weber et al., 2010). These two signaling pathway are needed to upregulate the expression of the key E3 ligases, MuRF1, which mediate the inhibition of protein synthesis (Clarke et al., 2007). Moreover, the insulin-signal transduction pathway, which was another pathway identified in the present study, is a highly conserved pathway that regulates cellular growth and when insulin binding to its cell-surface receptor, insulin receptor, activates a complex intracellular signaling network through insulin substrate proteins and the canonical PI3K and ERK cascades (Hocquette et al., 2010). Interestingly, insulin signaling is one of important factors involved in muscle development since stimulation of glucose utilization in fat and muscle cells in calves is occurring by enhancing insulin intracellular signaling (Jovanović et al., 2017).

Currently, systems biology approaches have become one of the most effective manners to accelerate the genetic improvement of beef and dairy cattle herds (Kadarmideen, 2014). It allows the selection of desired characteristics through the use of transcriptome profiles. The development of high throughput data and bioinformatics tools allow the selection of superior breeds without wasting time and money, contributing to the widespread use of transcriptome analysis in beef cattle operations. Our study shows in cattle that integration of pathway expression profiles in a systems biology approach will contribute to a better understanding of the genes and regulatory processes involved in marbling. These novel insights can be used in the future to take into account when improving the meat quality in beef cattle. The molecular mechanisms which underlie fat content in muscle can provide vital information for the production of healthier beef for human consumption.

## CONCLUSIONS

The outcome of our research is the identification of biological pathways where we highlighted changed genes which are related with marbling in beef cattle. These results give a better understanding of mechanisms involving marbling in beef cattle, which is economically important carcasses trait for meat quality. Moreover, the genes involved in the highlighted pathways can potentially be utilized as an early biological marker for marbling fat content in breed-specific differences in growth performance and meat quality.

## DATA AVAILABILITY STATEMENT

The datasets analyzed for this study can be found in the NCBI Gene Expression Omnibus (GEO) under the accession number GSE46411 and in the wikipathways.org knowledge base under the accession numbers: WikiPathways: WP2890, WikiPathways: WP2891, WikiPathways: WP2891, and WikiPathways: WP2901.

## ETHICS STATEMENT

Ethical review and approval was not required for the animal study because this study is just analysis.

## AUTHOR CONTRIBUTIONS

ZR, SC, and CE conceptualized and designed the study. ZR, SC, MK, LE, and JM analyzed the data. ZR collected and assembled the data and wrote the manuscript and SC, CE, and TS contributed to reviewing and editing the manuscript and also read and approved the final manuscript.

## ACKNOWLEDGMENTS

The publication of the manuscript was financed by the KNOW (Leading National Research Centre) Scientific Consortium "Healthy Animal–Safe Food," under Ministry of Science and Higher Education Decision No. 05-1/KNOW2/2015 (Poland).

## REFERENCES

and interleukin-1b-treated human preadipocytes are potent leptin producers. Cytokine 32 (2), 94–103. doi: 10.1016/j.cyto.2005.08.003

and pathways interacting to modulate divergent growth in cattle. BMC Genomics, 14 (1), 798. doi: 10.1186/1471-2164-14-798

Conflict of Interest: 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 © 2020 Roudbari, Coort, Kutmon, Eijssen, Melius, Sadkowski and Evelo. 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 Identification and Characterization of DNA Methylation and Long Non-Coding RNA Expression in Gastric Cancer

Peng Song1† , Lei Wu2† and Wenxian Guan1\*

<sup>1</sup> Department of General Surgery, Nanjing Drum Tower Hospital, The Affiliated Hospital of Nanjing University Medical School, Nanjing, China, <sup>2</sup> Department of Laboratory Medicine, The First Affiliated Hospital of Nanjing Medical University, Nanjing, China

### Edited by:

Manoj Kumar Kashyap, Amity University Gurgaon, India

### Reviewed by:

Arivusudar Everad John, Mazumdar Shaw Medical Centre, India Yashwanth Subbannayya, Norwegian University of Science and Technology, Norway

> \*Correspondence: Wenxian Guan

<EMAIL_ADDRESS>
† These authors have contributed equally to this work

### Specialty section:

This article was submitted to Bioinformatics and Computational Biology, a section of the journal Frontiers in Genetics

Received: 26 November 2019 Accepted: 27 January 2020 Published: 27 February 2020

### Citation:

Song P, Wu L and Guan W (2020) Genome-Wide Identification and Characterization of DNA Methylation and Long Non-Coding RNA Expression in Gastric Cancer. Front. Genet. 11:91. doi: 10.3389/fgene.2020.00091 Abnormal DNA methylation, an epigenetic modification, has increasingly been linked to the pathogenesis of many human cancers. However, there has been little focus on the DNA methylation patterns of genes encoding long noncoding RNAs (lncRNAs) in gastric cancer (GC). This study comprehensively determined DNA methylation and lncRNA expression profiles in GC through genome-wide analysis. Differentially methylated loci and lncRNAs were identified by integrating multi-omics data. In total, 548 differentially methylated CpG sites in lncRNA promoters and 2,399 differentially expressed lncRNAs were screened that were capable of distinguishing GC from normal tissues. Among them, 22 differentially methylation sites in 17 lncRNAs were inversely related to expression levels. Further analysis of DNA methylation status and gene expression level in GC revealed that three CpG sites (cg01550148, cg22497867, and cg20001829) and two lncRNAs (RP11- 366F6.2 and RP5-881L22.5) were significantly associated with GC patient overall survival. Molecular function analysis showed that these abnormally methylated lncRNAs were mainly involved in transcriptional activator activity. Our study identified several lncRNAs regulated by aberrant DNA methylation that have clinical utility as novel prognostic biomarkers in GC. These findings help improve the understanding of methylated patterns of lncRNAs and further our knowledge of the role of epigenetics in cancer development.

Keywords: DNA methylation, long non-coding RNA, epigenetics, prognosis, gastric cancer

## INTRODUCTION

Gastric carcinoma (GC) is the fourth most prevalent malignancy and third leading cause of cancer death worldwide (Torre et al., 2015). Histologically, GC demonstrates marked heterogeneity at the cytologic level, resulting in the classification of tumor subtypes. Distinct molecular genetic profiles, morphology, and expression of specific markers have been used to investigate the diversity and characteristics of GC (Zouridis et al., 2012; Cancer Genome Atlas Research Network, 2014). Therefore, identifying potential biomarkers to further understand the pathogenesis of GC is critical.

**60**

Long noncoding RNAs (lncRNAs) are loosely defined as RNAs more than 200 bases in length with no apparent coding capacity (Mattick and Rinn, 2015). LncRNAs regulate gene expression at transcriptional and post-transcriptional levels and thus are involved in diverse biological functions. Furthermore, recent studies have demonstrated a role for lncRNAs in carcinogenesis (Spizzo et al., 2012; Zhuo et al., 2019). DNA methylation, a key epigenetic mechanism, plays a crucial role in the regulation of gene expression, genomic imprinting, genome stabilization, and chromatin modification. Aberrant DNA methylation has been reported to be involved in the formation and progressions of diseases, especially cancers (Guo et al., 2018; Xu et al., 2019). Recent studies have showed that expression alterations of lncRNA-encoding genes mediated by changes in methylation can subsequently affect their downstream targets. For instance, the lncRNA C5orf66-AS1 functions as a tumor suppressor gene in GC, and aberrant hypermethylation of the regions around its transcription start site (TSS) is associated with its expression and is cancer-specific (Guo et al., 2018). This study indicated that hypermethylation of the C5orf66-AS1 promoter may serve as a potential prognostic marker in predicting GC patient survival. Shahabi et al. identified an epigenetically deregulated lncRNA linc00261, whose expression was lost in lung adenocarcinoma through DNA methylation silencing. The authors found that linc00261 acted upstream of ATM activation to facilitate DNA damage response activation and its loss resulted in malignant phenotypes and predisposed lung cells to cancer development (Shahabi et al., 2019). Additionally, lncRNAs showing aberrant DNA methylation may serve as potential epigenetically-based diagnostic factors. Silencing from CpG-island methylation of promoter-induced transcribed ultraconserved regions (T-UCRs) is common in many tumors and linked to colorectal cancer diagnosis (Kottorou et al., 2016; Honma et al., 2017). Therefore, elucidating the relationship between DNA methylation and lncRNA expression is essential for understanding GC development and potentially identifying new prognostic or diagnostic markers.

Here, we employed multigenomic data from The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) datasets to systematically characterize global DNA methylation levels, lncRNA expression profiles, and clinical features in GC. Our results decode the landscape of DNA methylation-mediated regulation for lncRNAs and provide promising biomarkers in the diagnosis and treatment of GC.

## METHODS

## DNA Methylation and Gene Expression Data

The DNA methylation array data (Illumina Infinium Human Methylation27, 450 BeadChip) were downloaded from the UCSC Xena browser (https://xenabrowser.net/). A Human Methylation27 BeadChip array of GC (GSE30601) was obtained from the GEO database (https://www.ncbi.nlm.nih. gov/geo/). Level-3 RNA-sequencing data (HTSeq-Counts and HTSeq-FPKM-UQ) and the clinicopathological and survival data of GC patients were also downloaded from the Xena website.

## Analysis of DNA Methylation Data

Differentially methylated CpG sites (DMCs) and differentially methylated regions (DMRs) between GC samples and adjacent tissues were identified using the minfi package (Version: 1.32.0; http://www.bioconductor.org/packages/release/bioc/html/minfi. html) (Fortin et al., 2017). Bump hunting method was applied to identify DMRs. False discovery rate (FDR) was calculated from multiple testing corrections of raw P-value by the Benjamini and Hochberg method (Benjamini et al., 2001). The genomic annotation of each CpG site was conducted using the hm27.hg38.manifest file (http://zwdzwd.io/InfiniumAnnotation). The coordinates of the individual lncRNA were extracted from GENCODE v22 (https://www.gencodegenes.org/human/release\_ 22.html). After the preprocessing the coordinates of CpG sites and lncRNAs, we further integrated both information based on the genomic location, considering differentially methylated loci within promoter regions (DNA sequences between –2,500 and 1,000 bp relative to the putative TSS). Manhattan plot was constructed to depict the distribution of CpG sites according to FDR via qqman package (Version: 0.1.4; https://cran.r-project.org/web/packages/ qqman/index.html) (Turner, 2018).

## Differential Long Non-Coding RNA Expression Analysis

Read count tables were imported into the edgeR package for identifying differentially expressed transcripts (Version: 3.7, https://bioconductor.org/packages/release/bioc/html/edgeR. html) (McCarthy et al., 2012). LncRNA catalogue was retrieved from GENCODE v22. Genes with FDR < 0.05 and absolute fold change (FC) > 2 were considered differentially expressed lncRNAs (DElncs).

## Integrated Analysis of DNA Methylation and Long Non-Coding RNA Expression

The correlation analysis between DMCs and DElncs was calculated and those with |coefficient of correlation| > 0.3 and P-value < 0.05 were considered significant. The visualization of the potential regulation of CpG sites to genes was constructed in Cytoscape 3.7.1 (Shannon et al., 2003).

## Functional Annotation and Enrichment Analysis for Long Non-Coding RNAs

ClusterProfiler tool (Version: 3.8.1, https://bioconductor.org/ packages/release/bioc/html/clusterProfiler.html) was used to perform Gene Ontology (GO) function and Gene Set Enrichment Analysis (GSEA) for DElncs with the DMCs (Yu et al., 2012). Spearman's correlation coefficients of expression levels between DElncs and protein-coding RNAs were calculated. The deregulated protein-coding genes were considered for GO analysis, setting parameters as "pAdjustMethod" = "BH," "pvalueCutoff" = 0.05, and "qvalueCutoff" = 0.05 for multiple comparisons. Terms from GO (molecular function) database slice were tested for enrichment. Mappings between GO terms and Entrez Gene IDs relied on the regularly updated R package org.Hs.eg.db (Version 3.10.0; https://bioconductor.org/packages/ release/data/annotation/html/org. Hs.eg.db.html). The value of log2 (FC) calculated by edgeR package was used as ranking metric for GSEA. We used the canonical pathways sub-collection of the C2 collection in the Molecular Signatures Database as the gene sets in the analysis. The leading-edge subset of genes in an enriched gene set are defined as those that appear in the ranked list before the point at which the running sum reaches its maximum deviation from zero.

## Statistical Analysis

Cox proportional hazard regression analyses were carried out to compare clinical features, DNA methylation, DElncs expression, and GC patients' prognosis. Survival curves were compared using Kaplan-Meier Plotter with log-rank test. All statistical analysis was two-sided and P < 0.05 was defined as statistically significant. Statistical analysis was performed using R programming language v.3.5.3.

## RESULTS

## Characteristics of the DNA Methylation Pattern in Gastric Cancer

Because the TCGA-450k set contained only two normal samples and the number was too small to reach statistical significance with respect to determining the DNA methylation profile of GC, we used the TCGA-27k set to identify DMCs and DMRs. We obtained 6,404 CpG sites with FDR < 0.05 between 48 GC and 25 non-tumor samples and identified 1,078 DMCs with a delta-beta value > 0.2. A total of 103 DMRs were identified based on the following parameters: resamples = 100, cut off = 0.2, and probe number ≥ 2. The 103 DMRs included 65 hypermethylated regions and 38 hypomethylated regions.

To identify DNA methylation alterations in lncRNA promoter regions, 3,010 CpG sites were examined. The methylation distribution in lncRNAs showed a V-shaped curve around the TSSs, indicating a relative reduction of the methylation density at the TSS (Figure 1A). As shown by Manhattan plot (Figure 1B), the CpG sites were distributed in all chromosomes, and 698 probes were found using the threshold of FDR < 0.05. We subsequently validated the DNA methylation patterns of the CpG sites of interest in an independent cohort (GSE30601) and found that 548 probes overlapped with the 698 sites reported in TCGA-27k set (Supplementary Table S1).

## Characteristics of Long Non-Coding RNA Expression in Gastric Cancer

To determine the lncRNA expression profile in GC, RNA-seq data of 375 GC tumors and 32 normal tissues were retrieved from TCGA. Among the 6,820 lncRNAs, we identified 2,399 DElncs, including upregulated 1,830 lncRNAs and 569 downregulated lncRNAs, using the criteria of FDR < 0.05 and absolute FC > 2 (Figure 2A, Supplementary Table S2). We then analyzed the categories of the 2,399 DElncs, as shown in Figure 2B. Long intergenic non-coding RNAs (lincRNAs) accounted for 54.3% of all DElncs, followed by antisense transcripts (30.3%). The remaining non-coding transcript types were sense\_intronic transcripts (4.5%), processed\_transcripts (3.7%), and sense\_overlapping transcripts (1.3%).

## Integrated Analysis of Differential Methylation and Long Non-Coding RNA Expression Data

After the profiles of DNA methylation and lncRNA expression were preprocessed, we combined the two omics data for further analysis. By associating the 548 DMCs to 2,399 DElncs, 31

lncRNAs (absolute FC >2 and FDR < 0.05). (B) Pie chart shows the number of differentially expressed lncRNAs in each category.

negative correlated pairs and 1 positive correlated pair were obtained in TCGA-27k set. DNA methylation in promoters is well known to negatively correlate with corresponding gene expression (Mosquera Orgueira, 2015). The 31 negative correlated pairs were validated in TCGA-450k set, and 22 probes showed a significantly inverse correlation with the promoter methylation of 17 aberrantly expressed lncRNAs (Supplementary Table S3). We used the negatively correlated pairs to construct a DNA methylation-regulated network that was composed of 39 nodes, including 6 hypermethylated DMCs, 16 hypomethylated DMCs, 13 upregulated lncRNAs, and 4 downregulated lncRNAs (Figure 3A). As shown in Figure 3B,

FIGURE 3 | Relation between DNA methylation and long noncoding RNA (lncRNA) expression. (A) Correlation between differentially methylated CpGs (DMCs) and lncRNAs. Circles and rectangles represent lncRNAs and DMCs, respectively. Red color indicates upregulated or hypermethylated, and blue indicates downregulated or hypomethylated. (B) Correlation (P values derive from Spearman's correlation) between DNA methylation and the expression of HOX family genes associated with five lncRNAs in matched samples.

the expression of 5 lncRNAs (HOTAIR, HOTTIP, HOXA11-AS, HOXB-AS4, and HOXC-AS3) generated in HOX family genes were negatively correlated with their methylation levels. The functions of these 17 lncRNAs are listed in Supplementary Table S4; only 6 of the lncRNAs have been reported to function GC (Zhang et al., 2014; Liu et al., 2015; Sun et al., 2016; Wu et al., 2017; Liu et al., 2018; Zhang et al., 2018; Song et al., 2019).

## Impact of DNA Methylation and Long Non-Coding RNA Expression on Gastric Cancer Survival

Univariate Cox regression was used to evaluate the association of the 22 probes and 17 lncRNAs with overall survival in GC and the results identified three CpG sites (cg01550148, cg22497867, and cg20001829) and two lncRNAs (RP11-366F6.2 and RP5- 881L22.5) with P < 0.05. Forest plot demonstrated that the methylation of the three probes and two lncRNAs were associated with the overall survival time of GC patients (Figure 4A). Figure 4B shows the Kaplan-Meier curves for survival in GC patients according to RP11-366F6.2 and RP5- 881L22.5 expression. High expression of RP11-366F6.2 was significantly correlated with poor survival compared with low expression. Additionally, poor survival was observed for patients with low expression of RP5-881L22.5 compared with patients with high levels (P = 0.039). In the multivariate Cox analysis, even after adjustment by tumor stage and other covariates, the expression of the two lncRNAs was still significantly associated with patient survival (P < 0.05, Table 1).

## Association of Deregulated Long Non-Coding RNAs With Biological Pathways and Processes

To better understand the biological function of the 17 DElncs, we constructed a co-expression network of deregulated protein-coding genes and lncRNAs. Using the Spearman's correlation coefficient above 0.8, a total of 32 deregulated mRNAs co-expressed with five lncRNAs (ZNF667-AS1, RP5-881L22.5, HOTAIR, HOTTIP, and HOXC-AS3) were acquired for GO enrichment analysis. Only nine protein-coding genes associated with HOTAIR, HOXC-AC3, and

TABLE 1 | Univariate and multivariate Cox regression analysis of variables associated with gastric cancer (GC) patient survival.

HR, hazard ratio; CI, confidence interval.

\*Using the mean expression of genes as the cut-off value.

† Using the median expression of genes as the cut-off value.

ZNF667-AS1 were assigned GO molecular function, involving oxidoreductase activity, nucleotide diphosphatase activity, and transcriptional activator activity (Figure 5A). In addition, GSEA analysis was performed to identify the associated biological processes and signaling pathways for these deregulated lncRNAs. As an example, we explored a functionally unknown lncRNA, RP5- 881L22.5, with significantly hypomethylated promoter regions that was upregulated in GC tissue and implicated with prognosis (Figure 5B). The expression of RP5-881L22.5 was positively correlated with "VECCHI\_GASTRIC\_CANCER\_EARLY\_UP" set, in which upregulated genes could differ early GC and normal tissue samples. The "NABA\_ECM\_GLYCOPROTEINS" set was enriched in the RP5-881L22.5 low expression group, which implied that this lncRNA could suppress tumor metastasis (Figure 5C).

## DISCUSSION

Cancer involves a complex regulatory network, and therefore integrating multiple omics data is required in the era of precision medicine (Olivier et al., 2019). The increasing applications of multi-omic profiling of GC have delivered new insight into the dynamics of this cancer type. In this study, we characterized DNA methylation in the promoters of lncRNA-encoding genes and inferred the potential lncRNAs regulated by aberrant DNA methylation in GC. Differential analyses were performed to compare DNA methylation and gene expression patterns between GC and normal tissues, and 548 DMCs and 2,399 DElncs were obtained. We thus identified lncRNAs (such as HOTAIR, HOTTIP, HOXA11-AS, HOXB-AS4, and HOXC-AS3) that could be regulated by aberrant DNA methylation via combination analysis. We further divided the potentially epigenetic regulated lncRNAs into different groups to explore their biological and clinical relationships with GC and found that the expressions of RP11-366F6.2 and RP5-881L22.5 were related to the prognosis of GC.

Methylation that interferes with transcription machinery binding to DNA has been reported to be highly associated with repression of gene transcription (Suzuki and Bird, 2008). Our research showed that a large number of lncRNAs were epigenetically deregulated by promoter methylation, and lncRNAs were globally hypomethylated, which was consistent with the previous observation that increased global DNA hypomethylation was a major event for the development and progression of cancer (Zeng et al., 2017). HOX genes are a subset of evolutionarily conserved homeobox genes that encode a class of important transcription factors that function in numerous developmental processes (Quinonez and Innis, 2014). By characterizing the transcriptional landscape of the four human HOX loci (A–D) at five base pair resolution in 11 anatomic sites, Rinn et al. identified 231 HOX non-coding RNAs (Rinn et al., 2007). HOTAIR is located in the HOXC cluster and serves as a scaffold protein by binding Polycomb repressive complex 2 (PRC2, including SUZ12, EED, and EZH2) via its 5′-domain and the LSD1/CoREST/REST complex via its 3′-domain, mediating gene silencing and reprogramming the overall chromatin dynamics in GC (Qu et al., 2019). Like HOTAIR, HOXA11-AS recruits EZH2 along with the histone demethylase LSD1 or DNMT1, which promotes proliferation and invasion of GC (Sun et al., 2016). HOTTIP enhances the expression of neighboring HOXA genes, particularly HOXA13 (Chang et al., 2016). HOXC-AS3, an antisense transcript of HOXC10, mediates gene transcriptional regulation in the tumorigenesis of GC by binding to YBX1 (Zhang et al., 2018). Genome-wide screening isolated HOXB-AS4 as specifically methylated in pancreatic cancer cells, which was useful to assess a cancer cell fraction in DNA samples (Ishihara et al., 2018).

Therapeutic targets and prognosis prediction from a comprehensive analysis of multi-omics data and clinical profiles is a critical for better understanding the biological complexity of GC. We identified two hypomethylated DElncs (RP11-366F6.2 and RP5-881L22.5) in GC that were significantly associated with overall survival. RP11-366F6.2, also called MAGEA4-AS1, is located in chrX and was reported to have significantly high expression in several tumor tissues, such as breast cancer and laryngeal squamous cell carcinoma (Yuan et al., 2017; Liu and Ye, 2019). Although pre-ranked GSEA analysis for RP11-366F6.2 returned no significantly gene sets, examining the functional roles of deregulated genes (such as MAGEA4, MAGEA10, HOXD10 and IGF2BP1) in the leading edge set indicated RP11-366F6.2 might be associated with tumor invasion and metastasis (Suzuki et al., 2008; Schultz-Thater et al., 2011; Xu et al., 2019). Regarding RP5-881L22.5, Zhu et al. developed an eleven-lncRNA signature, including this lncRNA, which could provide an effective individual mortality risk prediction and risk stratification in GC patients (Zhu et al., 2018). However, the biological functions of RP5-881L22.5 have

FIGURE 5 | Functional annotation for differentially methylated differentially expressed long noncoding RNAs (DElncs). (A) Circular plot of molecular function regulator Gene Ontology Term. (B) RP5-881L22.5 promoter hypermethylation in gastric cancer (GC) tumors compared with normal tissues. (C) Gene set enrichment analysis of RP5-881L22.5 in the The Cancer Genome Atlas (TCGA) dataset.

not been determined. GSEA results revealed that RP5-881L22.5 was likely to be involved in an extracellular matrix (ECM) interaction pathway. Glycoproteins make the ECM a cohesive network of molecules by linking cells together with structural components (Nallanthighal et al., 2019). Adhesive glycoproteins can bind to ECM components to activate downstream signaling pathways to regulate epithelial-mesenchymal transition, selfrenewal, migration, differentiation, and proliferation (Naba et al., 2016; Song et al., 2017). For example, the adhesion of cancer cells to fibronectin, a major adhesive ECM glycoprotein, remodels the tumor vasculature, enhances tumorigenicity, and facilitates metastasis. This mechanism could partly explain why decreased RP5-881L22.5 expression indicated a poor prognosis for GC patients.

To investigate the effect of epigenetically deregulated lncRNA in biological processes and pathways, an integrated analysis of DElncs and predicted mRNAs expression was performed. GO analysis revealed that these lncRNAs were involved in dysregulated transcriptional programs that invariably lead to cancer. We also found that predicted mRNAs HOXC10 and MSRB3 in GO analysis were significantly associated with overall survival in GC patients (Kim et al., 2019; Ma et al., 2019).

Several limitations in the present study should be pointed out. First, integrated analysis of genome-wide DNA methylation and lncRNA expression was based on the 27K Illumina array platform, which only contained 27,578 individual registered probes, and thus some possibly important methylation differences may be lacking from the current results. Second, the results of the present study are preliminary and primarily derived from bioinformatics analysis, and lack functional validation of the epigenetically deregulated lncRNAs. Third, due to limited availability of clinical data, it was not possible to obtain deeper insights into characterizing phenotype-genotype relationships.

In conclusion, the present results provide evidence for the changes of widespread DNA methylation of lncRNA-encoding genes in GC patients. The candidate factors identified in this study might function as potential molecular phenotypic biomarkers, especially RP11-366F6.2 and RP5-881L22.5, which were associated with prognosis. Our results help elucidate a more detailed explanation of epigenetic mechanisms for GC and deepen our understanding of the aberrantly methylated patterns in lncRNA-encoding genes.

## DATA AVAILABILITY STATEMENT

The following information was supplied regarding data availability:.

The level 3 TCGA data for DNA methylation arrays and lncRNA expression are available in Xena website.

HumanMethylation27: https://gdc.xenahubs.net/download/ TCGA-STAD.methylation27.tsv.gz

## REFERENCES

The microarray data that support this study are available through the NCBI database under accession GSE30601.

The data used to support the findings of this study are included within the article.

## AUTHOR CONTRIBUTIONS

PS and WG designed the study. PS and LW collected, analysed and interpreted the data. PS and LW wrote the draft. PS and WG edited the manuscript.

## ACKNOWLEDGMENTS

This research received no specific grant from any funding agency in the public, commercial, or not-for-profit sectors.

## SUPPLEMENTARY MATERIAL

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

Conflict of Interest: 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 © 2020 Song, Wu and Guan. Thisis 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.

# Identification of Familial Hodgkin Lymphoma Predisposing Genes Using Whole Genome Sequencing

Aayushi Srivastava1,2,3,4†, Sara Giangiobbe1,4†, Abhishek Kumar <sup>1</sup> , Nagarajan Paramasivam<sup>5</sup> , Dagmara Dymerska<sup>6</sup> , Wolfgang Behnisch<sup>7</sup> , Mathias Witzens-Harig<sup>4</sup> , Jan Lubinski <sup>6</sup> , Kari Hemminki 1,8, Asta Försti 1,2,3 and Obul Reddy Bandapalli 1,2,3,4 \*

<sup>1</sup> Division of Molecular Genetic Epidemiology, German Cancer Research Center (DKFZ), Heidelberg, Germany, <sup>2</sup> Hopp Children's Cancer Center (KiTZ), Heidelberg, Germany, <sup>3</sup> Division of Pediatric Neurooncology, German Cancer Research Center (DKFZ), German Cancer Consortium (DKTK), Heidelberg, Germany, <sup>4</sup> Medical Faculty, Heidelberg University, Heidelberg, Germany, <sup>5</sup> Computational Oncology, Molecular Diagnostics Program, National Center for Tumor Diseases (NCT), Heidelberg, Germany, <sup>6</sup> Department of Genetics and Pathology, International Hereditary Cancer Centre, Pomeranian Medical University, Szczecin, Poland, <sup>7</sup> Department of Pediatric Oncology, Hematology and Immunology, University of Heidelberg, Heidelberg, Germany, <sup>8</sup> Faculty of Medicine and Biomedical Center in Pilsen, Charles University in Prague, Pilsen, Czechia

### Edited by:

Lavanya Balakrishnan, Mazumdar Shaw Medical Centre, India

### Reviewed by:

Prashanth N. Suravajhala, Birla Institute of Scientific Research, India Raghu Metpally, Geisinger Health System, United States

### \*Correspondence:

Obul Reddy Bandapalli<EMAIL_ADDRESS>
†These authors have contributed equally to this work

### Specialty section:

This article was submitted to Bioinformatics and Computational Biology, a section of the journal Frontiers in Bioengineering and Biotechnology

Received: 03 December 2019 Accepted: 21 February 2020 Published: 06 March 2020

### Citation:

Srivastava A, Giangiobbe S, Kumar A, Paramasivam N, Dymerska D, Behnisch W, Witzens-Harig M, Lubinski J, Hemminki K, Försti A and Bandapalli OR (2020) Identification of Familial Hodgkin Lymphoma Predisposing Genes Using Whole Genome Sequencing. Front. Bioeng. Biotechnol. 8:179. doi: 10.3389/fbioe.2020.00179 Hodgkin lymphoma (HL) is a lymphoproliferative malignancy of B-cell origin that accounts for 10% of all lymphomas. Despite evidence suggesting strong familial clustering of HL, there is no clear understanding of the contribution of genes predisposing to HL. In this study, whole genome sequencing (WGS) was performed on 7 affected and 9 unaffected family members from three HL-prone families and variants were prioritized using our Familial Cancer Variant Prioritization Pipeline (FCVPPv2). WGS identified a total of 98,564, 170,550, and 113,654 variants which were reduced by pedigree-based filtering to 18,158, 465, and 26,465 in families I, II, and III, respectively. In addition to variants affecting amino acid sequences, variants in promoters, enhancers, transcription factors binding sites, and microRNA seed sequences were identified from upstream, downstream, 5′ and 3′ untranslated regions. A panel of 565 cancer predisposing and other cancer-related genes and of 2,383 potential candidate HL genes were also screened in these families to aid further prioritization. Pathway analysis of segregating genes with Combined Annotation Dependent Depletion Tool (CADD) scores >20 was performed using Ingenuity Pathway Analysis software which implicated several candidate genes in pathways involved in B-cell activation and proliferation and in the network of "Cancer, Hematological disease and Immunological Disease." We used the FCVPPv2 for further in silico analyses and prioritized 45 coding and 79 non-coding variants from the three families. Further literature-based analysis allowed us to constrict this list to one rare germline variant each in families I and II and two in family III. Functional studies were conducted on the candidate from family I in a previous study, resulting in the identification and functional validation of a novel heterozygous missense variant in the tumor suppressor gene DICER1 as potential HL predisposition factor. We aim to identify the individual genes responsible for predisposition in the remaining two families and will functionally validate these in further studies.

Keywords: familial Hodgkin lymphoma, whole genome sequencing, predisposing genes, germline variants, variant prioritization, next generation sequencing, genetic predisposition to disease

## INTRODUCTION

Hodgkin lymphoma (HL) is a lymphoproliferative malignancy originated in germinal center B-cells and is reported to account for about 10% of newly diagnosed lymphomas and 1% of all de novo neoplasms worldwide with an incidence of about 3 cases per 100,000 people in Western countries (Diehl et al., 2004). It is one of the most common tumors in young adults in economically developed countries, with one peak of incidence in the third decade of life and a second peak after 50 years of age.

Based on differences in the morphology and phenotype of the lymphoma cells and the composition of the cellular infiltrate, HL is subdivided into classical Hodgkin lymphoma (cHL) that accounts for about 95% of cases and nodular lymphocytepredominant Hodgkin lymphoma (NLPHL) that accounts for the remaining 5% of cases (Kuppers, 2009).

Although familial risk for HL is reported to be among the highest of all cancers (Kharazmi et al., 2015), not many genetic risk factors have been identified. An association between various HLA class I and class II alleles and increased risk of HL has been reported (Diepstra et al., 2005), while other non-HLA susceptibility loci have been detected through genome-wide association studies (Frampton et al., 2013; Cozen et al., 2014; Kushekhar et al., 2014). The identification of major predisposing genes is a more daunting task, however, rare germline variants in KLDHC8B, NPAT, ACAN, KDR, DICER1, and POT1 gene have been reported by different groups in high-risk HL families (Salipante et al., 2009; Saarinen et al., 2011; Ristolainen et al., 2015; Rotunno et al., 2016; Bandapalli et al., 2018; Mcmaster et al., 2018).

Here we report the results of whole genome sequencing (WGS) performed in three families with documented recurrence of HL. We used our Familial Cancer Variant Prioritization Pipeline (FCVPPv2) (Kumar et al., 2018) as well as two gene/variant panels based on cancer predisposing genes and variants prioritized in the largest familial HL cohort study to date in order to identify possible disease-causing high-penetrance germline variants in each family (Zhang et al., 2015; Rotunno et al., 2016). Pathway and network analyses using Ingenuity Pathway Analysis software also allowed us to gain insight into the molecular mechanisms of the pathogenesis of HL. We hope that these results can be used in the development of targeted therapy and in the screening of other individuals at risk of developing HL.

## MATERIALS AND METHODS

## Patient Samples

Three families with documented recurrence of HL were analyzed in this study, with a total number of 16 individuals (7 affected and 9 unaffected). HL family I and family III were recruited at the University Hospital of Heidelberg, Germany, while family II was recruited at the Pomeranian Medical University, Szczecin, Poland.

The study was approved by the Ethics Committee of the University of Heidelberg and Pomeranian Medical University, Poland. Collection of blood samples and clinical information from subjects was undertaken with a written informed consent in accordance with the tenets of the Declaration of Helsinki.

Germline DNA samples used for genome sequencing were isolated from peripheral blood using QIAamp <sup>R</sup> DNA Mini kit (Qiagen, Cat No. 51104) according to the manufacturer's instructions.

## Whole Genome Sequencing, Variant Calling, Annotation and Filtering

Whole genome sequencing (WGS) of available affected and unaffected members of the three HL families was performed using Illumina-based small read sequencing. Mapping to reference human genome (assembly version Hs37d5) was performed using BWA mem (version 0.7.8) (Li and Durbin, 2009) and duplicates were removed using biobambam (version 0.0.148). The SAMtools suite (Li, 2011) was used to detect single nucleotide variants (SNVs) and Platypus (Rimmer et al., 2014) to detect indels. Variants were annotated using ANNOVAR, 1000 Genomes, dbSNP, and ExAC (Smigielski et al., 2000; Wang et al., 2010; The Genomes Project Consortium et al., 2015; Lek et al., 2016). Variants with a quality score >20 and a coverage >5×, SNVs that passed the strand bias filter (a minimum one read support from both forward and reverse strand) and indels that passed all the Platypus internal filters were evaluated further for minor allele frequencies (MAFs) with respect to the 1,000 Genomes Phase 3 and non-TCGA ExAC data. Variants with a MAF <0.1% were deduced from these two datasets. A pairwise comparison of shared rare variants was performed to check for sample swaps and family relatedness.

## Data Analysis and Variant Prioritization Prioritization of Coding Variants

Variant evaluation was performed using the criteria of our in-house developed variant prioritization pipeline (FCVPPv2) (Kumar et al., 2018). Shortly, variants with MAF < 0.1% were first filtered based on the pedigree data considering cancer patients as cases and unaffected persons as controls. The probability of an individual being a Mendelian case or true control was considered.

Variants were then ranked using the CADD tool v1.3 (Kircher et al., 2014). Only variants with a scaled PHRED-like CADD score >10, i.e., variants belonging to the top 1% of probable deleterious variants in the human genome, were considered further. Genomic Evolutionary Rate Profiling (GERP) (Cooper et al., 2005), PhastCons (Siepel et al., 2005), and PhyloP (Pollard et al., 2010) were used to evaluate the evolutionary conservation of a particular variant. GERP scores > 2.0, PhastCons scores > 0.3, and PhyloP scores ≥ 3.0 were indicative of a good level of conservation and were therefore used as thresholds in the selection of potentially causative variants.

Next, all missense variants were assessed for deleteriousness using 10 tools accessed using dbNSFP (Liu et al., 2016), namely SIFT, PolyPhen V2-HDV, PolyPhen V2-HVAR, LRT, MutationTaster, Mutation Assessor, FATHMM, MetaSVM, MetLR, and PROVEAN. Variants predicted to be deleterious by at least 60% of these tools were analyzed further. Prediction scores for nonsense variants were attained via VarSome (Kopanos et al., 2018), the final verdict on pathogenicity offered by VarSome was based on the following tools: DANN, MutationTaster, FATHMM-MKL, FATHMM-XF, ALoFT, EIGEN, EIGEN PC, and PrimateAI.

Lastly, three different intolerance scores derived from NHLBI-ESP6500 (Petrovski et al., 2013), ExAC (Lek et al., 2016) and a local dataset, all of which were developed with allele frequency data, were included to evaluate the intolerance of genes to functional mutations. However, these scores were used merely to rank the variants and not as cut-offs for selection. The ExAC consortium has developed two additional scoring systems using large-scale exome sequencing data including intolerance scores (pLI) for loss-of-function variants and Z-scores for missense and synonymous variants. These were used for nonsense and missense variants, respectively.

Structural variants were analyzed using Canvas (version <IP_ADDRESS>3) (https://academic.oup.com/bioinformatics/article/ 32/15/2375/1743834) program's SmallPedigree-WGS separately to detect the larger copy number variants. The joint genotyped VCF for all the samples in a family generated via Platypus was used as the b-allele input file along with the BAM files, and the rest of the parameters were kept default. Variants with "PASS" filters and present in all the cases in a family were processed further and variants overlapping common structural variants (AF > 1%) from gnomAD (version 2.1) were marked as common and removed. The remaining rare structural variants that affects the known cancer predisposition genes were selected for the manual inspection in IGV.

### Analysis of Non-coding Variants

Variants located in the 3′ and 5′ untranslated regions (UTRs) were prioritized based on their location in regulatory regions. Putative miRNA targets at variant positions within 3′ UTRs and 1 kb downstream of transcription end sites were detected by scanning the entire dataset of the human miRNA target atlas from TargetScan 7.0 (Agarwal et al., 2015) using the intersect function of bedtools. Similarly, 5′ UTRs and regions 1 kb upstream of transcription start sites were scanned for putative enhancers and promoters using merged enhancer and promoter data from the FANTOM5 consortium as well as super-enhancer data from the super-enhancer archive (SEA) and dbSUPER. These regions were also scanned for transcription factor binding sites using SNPnexus (Dayem Ullah et al., 2018).

The regulatory nature and the possible functional effects of non-coding variants were evaluated using CADD v1.3, HaploReg V4 (Ward and Kellis, 2012), and RegulomeDB (Boyle et al., 2012), primarily based on ENCODE data (Birney et al., 2007). Epigenomic data and marks from 127 cell lines from the NIH Roadmap Epigenomics Mapping Consortium were accessed via CADD v1.3, which gave us information on chromatin states from ChromHmm and Segway. CADD also provided mirSVR scores to rank predicted microRNA target sites by a downregulation score. These scores are based on a new machine learning method based on sequence and contextual features extracted from miRanda-predicted target sites (Betel et al., 2010). Furthermore, SNPnexus was used to access non-coding scores for each variant and to identify regulatory variants located in CpG islands.

The final selection of 3′ UTR and downstream variants was based on their CADD scores > 10 and whether or not they had predicted miRNA target site matches. Similarly, upstream and 5′ UTR variants in enhancers, promoters, super-enhancers or transcription factor binding sites with CADD scores >10 were short-listed.

### Presence of Candidate Variants in 565 Cancer Predisposing and Other Cancer-Related Genes

In a study on cancer predisposing genes (CPGs) in pediatric cancers, Zhang et al. compiled 565 CPGs based on review of the American College of Medical Genetics and Genomics (ACMG) and medical literature (Zhang et al., 2015). The categories included genes associated with autosomal dominant cancer-predisposition syndromes (60), genes associated with autosomal recessive cancer-predisposition syndromes (29), tumor-suppressor genes (58), tyrosine kinase genes (23), and other cancer genes (395). We checked a list of genes corresponding to our shortlisted coding and non-coding variants for their presence in the list of genes in the aforementioned study.

## Presence of Candidate Variants in Prioritized HL Genes From a Large WES-Based Familial HL Study

In a study by Rotunno et al. (2016) 2,699 variants corresponding to 2,383 genes were identified in 17 HL discovery families after filtering and prioritization. We intersected our list of candidate genes with this list of 2,383 HL genes to identify coding and non-coding variants from our shortlist in potentially causative HL genes.

## Variant Validation

Specific variants of interest mentioned throughout the text (DICER1, HLTF, LPP, PLK3, RAD51D, RELB, SH3GL2, and SPTAN1) and highlighted as bold in the tables were validated using specific primers for polymerase chain reaction amplification designed with Primer3 (http://bioinfo.ut.ee/ primer3-0.4.0/) and Sanger sequencing on a 3,500 Dx Genetic Analyzer (Life Technologies, CA, USA), using ABI PRISM 3.1 Big Dye terminator chemistry, according to the manufacturer's instructions. The electrophoretic profiles were analyzed manually. Segregation analysis of the prioritized variants was performed in additional family members when DNA was available. Primer details are available on request.

## Ingenuity Pathway Analysis (IPA)

IPA (Qiagen; http://www.qiagen.com/ingenuity; analysis date 15/10/2019) was used to perform a core analysis to identify enriched canonical pathways, diseases, biological functions, and molecular networks among genes that passed the allele frequency cut-off, fulfilled family-based segregation criteria, met the CADD score cut-off and were not intergenic or intronic variants. Data were analyzed for all three families together. Top canonical pathways were identified from the IPA pathway library and ranked according to their significance to our input data. This significance was determined by p-values calculated using the right tailed Fisher's exact test.

IPA was also used to generate gene networks in which upstream regulators were connected to the input dataset genes while taking advantage of paths that involved more than one link (i.e., through intermediate regulators). These connections represent experimentally observed cause-effect relationships that relate to expression, transcription, activation, molecular modification and transport as well as binding events. The networks were ranked according to scores that were generated by considering the number of focus genes (input data) and the size of the network to approximate the relevance of the network to the original list of focus genes.

## RESULTS

## Whole Genome Sequencing Results

In our study, we analyzed three families with reported recurrence of Hodgkin lymphoma. Their respective pedigrees are shown in **Figure 1**.

In family I (**Figure 1A**), the proband (III-1) and her mother (II-2) were diagnosed with two different histological subtypes of classical Hodgkin lymphoma (cHL) at the ages of 7 and 34, respectively. The daughter was diagnosed with nodular sclerosis cHL and the mother with lymphocyte-rich cHL. The sample of the unaffected father (II-1) was also sequenced. Family II (**Figure 1B**) is characterized by a strong recurrence of HL. Five family members were diagnosed with HL (II-3, II-4, III-3, III-4, and III-5), of which three (III-3, III-4, and III-5) underwent WGS. In addition, the family member (II-6), who was considered as an obligatory carrier of the mutation, was sequenced as were samples and four healthy family members (III-1, III-2, III-6, and III-8) and one family member diagnosed with uterine cancer (II-1) as controls. In family 3 (**Figure 1C**), II-1 and II-2 were diagnosed with cHL, at the age of 27 and 24, respectively. Their parents (I-1, I-2) were not affected, however one of them is expected to be a carrier and analyzed accordingly.

WGS of 7 affected and 9 unaffected members from the three studied families identified a total number of 98,564, 170,550, and 113,654 variants which were reduced by pedigree-based filtering to 18,158, 465, and 26,465 in families I, II, and III, respectively.

## Prioritization of Candidates According to the FCVPPv2

After pedigree-based filtering, 130, 7, and 196 exonic variants were left in families I, II, and III, respectively, with a prevalence of non-synonymous and synonymous SNVs. The predominant type of substitution was the C>T transition. Among exonic variants fulfilling pedigree-based criteria, only variants with CADD scores >10 were taken into further consideration and prioritized according to deleteriousness, intolerance, and conservational scores, as detailed in the methods section. At the end of this process, 37 potential missense variants and 9 potential nonsense mutations were prioritized for families I—III and are shown in **Tables 1**, **2**.

Pedigree-based filtering also reduced the number of potentially interesting variants located in the untranslated regions to 523 for 5'UTR variants (130 in family I, 5 in family II, and 314 in family III) and 854 for 3'UTR variants (347 in family I, 10 in family II, and 497 in family III). These variants were further prioritized based on their CADD score >10 and their localization in known regulatory regions (**Supplementary Table 1**). 5′UTR variants were analyzed by the SNPNexus tool, which allowed us to identify 4 variants located in transcription factors binding sites. In addition, the intersect function of bedtools was used to identify further 15 variants located in promoter regions and 4 located in super-enhancer regions. Among variants located in the 3′UTR region, 56 variants located in miRNA seed sequences were selected.

Analysis of structural variants resulted in identification of a large deletion in exons 9 and 10 (del5395) of Chek2 kinase gene (CHEK2) in family 1 that segregates with the disease.

## Candidate Variants in 565 CPGs and 2383 Potentially Causative HL Genes

Intersecting our prioritized list of candidate genes with the list of 565 CPGs, we identified 11 variants in nine genes in coding and selected non-coding regions (upstream and downstream variants, 3′ and 5′ UTRs) of the known CPGs. These include FUBP1, SEPT6, DICER1, EZR, and NCOA1 from family 1 and BCL6, RAD51D, LPP, and PTCH1 from family 3 (**Table 3**). DICER1 and PTCH1 are known in autosomal dominant cancerpredisposition syndromes, whereas the rest are categorized as being "other cancer genes."

In addition to the identification of 11 variants in CPGs, we intersected our prioritized list of genes with a list of 2,383 genes with potentially causative variants from a large WES-based familial HL study. We found 25 variants in the coding and noncoding regions in 23 of the HL genes, with 7 coming from family I and 18 from family III (**Table 4**).

## Network and Pathway Analysis With IPA

Pathway analysis of the selected variants performed with IPA showed an enrichment of mutations in genes involved in pathways essential for B-cell proliferation and activation, specifically B-cell receptor signaling, and PI3K signaling in B lymphocytes and B cell activating factor signaling (**Supplementary Table 2A**, **Supplementary Figures 1A,B**).

Similarly, the IPA network analysis generated a comprehensive picture of possible gene interactions between our candidate genes (**Supplementary Table 2B**). The top network is related to cancer, hematological disease and immunological disease, which is in complete coherence with the pathogenesis of HL. Many genes from the prioritized list of top candidates are shown to play a role in the top networks (**Figure 2**).

## Literature Mining, Consolidation of Results, and Selection of Candidates

With the aim of identifying one highly penetrant dominant variant per family, we used our pipeline results and literaturebased mining to determine the genes' link to Hodgkin lymphoma or immune-related processes. For family 1, we have short-listed 5 potential candidates (DICER1, HLTF, NOTCH3, PLK3, and RELB). Based on segregation, confirmation and functional validation, we identified DICER1 as a candidate HL predisposing gene by showing significant down-regulation of

tumor suppressor miRNAs in DICER1-mutated family members (Bandapalli et al., 2018). The presence of DICER1 in the list of 565 CPGs also reinforces its status as the disease-causing variant in this family.

In family 2, three exonic variants made it to the final list (ALAD, CERCAM, and SPTAN1) of which SPTAN1 was shown to be among the genes in one of the top IPA networks (Network 3; **Figure 2C**). No coding or non-coding variants intersected with the panel of CPGs or HL candidate genes.

Two genes stand out in family 3, namely LPP and RAD51D. Both genes were found in the list of 565 CPGs and LPP was additionally found in the gene list from the large cohort of HL families. Three variants in LPP were prioritized by the FCVPPv2 and made it to the shortlist including one stopgain variant (3\_188123978\_G\_T), one 3′ UTR (3\_188608373\_A\_T) and one non-synonymous missense variant (3\_188123979\_A\_T). LPP (LIM domain containing preferred translocation partner in lipoma) is a member of the zyxin family of LIM proteins that is characterized as a promoter of mesenchymal/fibroblast cell migration. LPP has been shown to be a critical inducer of tumor cell migration, invasion and metastasis by virtue of its ability to localize to adhesions and to promote invadopodia formation (Ngan et al., 2018). A genome-wide association study of 253 Chinese individuals with B-cell NHL also identified a new susceptibility locus between BCL6 and LPP that was significantly associated with the increased risk of B-cell NHL (Tan et al., 2013). On the other hand, there are no reports of an association between RAD51D and lymphomas; however, it is a well-established susceptibility gene in Breast-Ovarian Cancer, Familial 4 and Hereditary Breast Ovarian Cancer Syndrome (Loveday et al., 2011; Chen et al., 2018). The final selection of a candidate in this family will be based on further functional studies.

## DISCUSSION

In summary, WGS data analysis of three families with reported recurrence of HL allowed us to prioritize 45 coding and 79 non-coding variants from which we subsequently selected and validated one for family I (DICER1), short-listed three in family II (ALAD, CERCAM, and SPTAN1) and two in family III (RAD51 and LPP), to investigate further with validation and functional studies. For family I we have already functionally validated DICER1 as the candidate predisposing gene in a previous study (Bandapalli et al., 2018). However, it was important to include the family in this paper, especially with regard to the integrity of the pathway and network analyses. We identified pathways related to

### TABLE 1 | Top missense variants prioritized using the FCVPPv2.

Chromosomal positions, classifications, PHRED-like CADD scores, protein changes and the number of positive intolerance (Int) and deleteriousness (Del) scores are shown for each variant. Variants of interest that were validated by Sanger sequencing in the provided family samples have been shown in bold.

B-cell proliferation and networks related to cancer, hematological disease, immunological disease, hereditary disorders, cell death and cell survival using IPA software, helping us to prioritize genes with functions in the pathogenesis of HL. Interestingly, several genes in our gene list were related to DNA repair (e.g., NOTCH3, RAD51, and SPTAN1).

In the current study, we also identified a deletion of exon 9 and 10 in CHEK2 in family 1. The same deletion has been reported in several unrelated patients with breast cancer of Polish origin. In that study the deletion of exon 9 and 10 in CHEK2 was shown to lead to a premature protein truncation at codon 381 and to evoke a 2-fold increase in the risk of prostate cancer and a 4 fold increase in the risk of familial prostate cancer (Cybulski et al., 2006). The detection of mRNA of abnormal length suggests that the deletion does not lead to complete transcript loss and therefore, the effect of this truncating mutation on cancer risk may differ or work in tandem with another genetic effect, may be with DICER1 in this family but warrants further experiments. Personalized medicine is an upcoming and promising field of medicine in which medical decisions, practices, interventions,

### TABLE 2 | Top non-sense variants prioritized using the FCVPPv2.

[I]VarSome Scores.

PM2, Pathogenic Moderat; PP3, Pathogenic Supporting (no. of scores predicting pathogenicity); Uncertain Significance: No scores could be found for the variant in question, PVS1, Pathogenic Very Strong.

Chromosomal positions, classifications, PHRED-like CADD scores, protein changes, the number of positive intolerance (Int) and VarSome prediction scores are included for each variant. Variants of interest that were validated by Sanger sequencing in the provided family samples have been shown in bold.

TABLE 3 | Variants corresponding to genes present in the panel of 565 known cancer predisposition genes from a study by Zhang et al. (2015).

and products are tailored to the individual patient based on their predicted response or risk of disease. The scope of this field has advanced rapidly with the advent of genomics and other omics and the possibility of implicating one gene or a set of genes in the pathogenesis of a particular disease. Thus, the identification of germline predisposing genes could be of great value in the screening of individuals at risk of developing HL, as well as in the development of personalized adjuvant therapies based on the affected pathways. In this aspect, delta-aminolevulinate dehydratase (ALAD) from family 2 is interesting, as it is

Variant details from both databases the present study and the study by Rotunno et al. (2016) are shown.

involved in the catalysis of the second step in the biosynthesis of heme and also acts as an endogenous inhibitor of the 26 S proteasome, a multi-catalytic ATP-dependent protease complex that functions as the degrading arm of the ubiquitin system, which is the major pathway for regulated degradation of proteins in all eukaryotes. Down regulation of ALAD is shown to be associated with poor prognosis in patients with breast cancer (Ge et al., 2017) whereas the existing data on non-erythroid spectrin αII (SPTAN1) suggest that overexpression of SPTAN1 in tumor cells reflects neoplastic and tumor promoting activity or tumor suppressing effects by enabling DNA repair through interaction with DNA repair proteins (Ackermann and Brieger, 2019). CERCAM is known as an unfavorable prognostic marker in urothelial, renal, and ovarian cancers implying the importance of the variants in these genes (Ma et al., 2016). RAD51D from family III is particularly interesting since it is involved in DNA repair through homologous recombination. Therefore, it is possible that carcinomas arising in patients carrying mutations in this gene will be sensitive to chemotherapeutic agents that target this pathway, such as cisplatin and the PARP (poly (ADP-ribose) polymerase) inhibitor olaparib. This has already been demonstrated in BRCA1/2 mutation-carrier cancer patients

cellular movement. Genes from our input-data are shown in gray, genes from our prioritized candidate list are highlighted in peach.

(Banerjee et al., 2010; Loveday et al., 2011). This approach can also be applied to target pathways affected by the mutated genes. Several candidate genes were identified by IPA pathway analysis in B cell receptor pathways, offering a valuable target for other pharmaceutical drugs. The B cell receptor (BCR) signaling pathway, when dysregulated, is a potent contributor to lympomagenesis and tumor survival (Valla et al., 2018). This pathway has been targeted in B-cell lymphomas and leukemias with several BCR-directed agents, such as inhibitors of Bruton's tyrosine kinase (BTK9), spleen tyrosine kinase (SYK) and phosphatidylinositol-3-kinase (PI3K) (Buggy and Elias, 2012; Dreyling et al., 2017; Liu and Mamorska-Dyga, 2017). In one study, excellent response rates could be demonstrated in certain non-Hodgkin lymphoma subtypes, however, issues related to the development of resistance to BTK inhibitors need to be addressed (Valla et al., 2018).

Advancements in the field of genomics have allowed WGS to become the state-of-the-art tool for the identification of novel cancer predisposing genes in Mendelian diseases. It is still a challenge to appropriately interpret the immense amount of data generated by WGS, especially with respect to non-coding variants. In our study, we have attempted to interpret a selection of non-coding variants using in silico and bioinformatic tools, however, the adequate analysis of intronic and intergenic variants remains a challenge. There are several reports of WGS being successfully implemented to implicate rare, high-penetrance germline variants in cancer, for example POT1 mutations in familial melanoma and Hodgkin lymphoma (Mcmaster et al., 2018; Wong et al., 2019) and POLE and POLD1 mutations in colorectal adenomas or carcinomas (Palles et al., 2013). In a previous study, we have used our pipeline (FCVPPv2) to prioritize novel variants in non-medullary thyroid cancer prone families (Srivastava et al., 2019). We have also successfully combined our pipeline with literature review and functional studies to identify DICER1 as a candidate predisposing gene in one Hodgkin lymphoma family (Bandapalli et al., 2018). We aim to apply these methods in the remaining Hodgkin lymphoma families and hope that these results will facilitate personalized therapy in the studied families and contribute to the screening of other individuals at risk of developing HL.

## DATA AVAILABILITY STATEMENT

Unfortunately, for reasons of ethics and patient confidentiality we are not able to provide the sequencing data into a public data base. The data underlying the results presented in the study are available from the corresponding author or from Dr. Asta Försti (Email: a.foersti@kitz-heidelberg.de).

## ETHICS STATEMENT

The studies involving human participants were reviewed and approved by Ethics Committee of the University of Heidelberg, Germany & Ethics Committee of the Pomeranian Medical University, Poland. Written informed consent to participate in this study was provided by the participants' legal guardian/next of kin.

## REFERENCES

## AUTHOR CONTRIBUTIONS

OB, AF, and KH conceived and designed the study. WB, MW-H, DD, and JL provided the HL family samples. NP ran WGS pipeline and CNVs analysis. AS, OB, SG, and AK analyzed the data. OB and SG performed the experiments. AS and OB wrote the first draft of the manuscript. All authors read, commented on, and approved the manuscript.

## FUNDING

This study was supported by the Harald Huppert Foundation and Transcan ERA-NET funding from the German Federal Ministry of Education and Research (BMBF).

## ACKNOWLEDGMENTS

The authors thank the Genomics and Proteomics Core Facility (GPCF) of the German Cancer Research Center (DKFZ) for providing excellent library preparation and sequencing services and the Omics IT and Data Management Core Facility (ODCF) of the DKFZ for the whole genome sequencing data management.

## SUPPLEMENTARY MATERIAL

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

and survival in BRCA1/2-negative breast cancers. Ann. Oncol. 29, 2046–2051. doi: 10.1093/annonc/mdy338

**Conflict of Interest:** 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 © 2020 Srivastava, Giangiobbe, Kumar, Paramasivam, Dymerska, Behnisch, Witzens-Harig, Lubinski, Hemminki, Försti and Bandapalli. 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.

# A Computational Approach for Mapping Heme Biology in the Context of Hemolytic Disorders

Farah Humayun1,2† , Daniel Domingo-Fernández<sup>2</sup> \* † , Ajay Abisheck Paul George<sup>1</sup> , Marie-Thérèse Hopp<sup>1</sup> , Benjamin F. Syllwasschy<sup>1</sup> , Milena S. Detzel<sup>1</sup> , Charles Tapley Hoyt<sup>2</sup> , Martin Hofmann-Apitius<sup>2</sup> \* and Diana Imhof<sup>1</sup> \*

<sup>1</sup> Pharmaceutical Biochemistry and Bioanalytics, Pharmaceutical Institute, University of Bonn, Bonn, Germany, <sup>2</sup> Department of Bioinformatics, Fraunhofer Institute for Algorithms and Scientific Computing (SCAI), Sankt Augustin, Germany

Heme is an iron ion-containing molecule found within hemoproteins such as hemoglobin and cytochromes that participates in diverse biological processes. Although excessive heme has been implicated in several diseases including malaria, sepsis, ischemiareperfusion, and disseminated intravascular coagulation, little is known about its regulatory and signaling functions. Furthermore, the limited understanding of heme's role in regulatory and signaling functions is in part due to the lack of curated pathway resources for heme cell biology. Here, we present two resources aimed to exploit this unexplored information to model heme biology. The first resource is a terminology covering heme-specific terms not yet included in standard controlled vocabularies. Using this terminology, we curated and modeled the second resource, a mechanistic knowledge graph representing the heme's interactome based on a corpus of 46 scientific articles. Finally, we demonstrated the utility of these resources by investigating the role of heme in the Toll-like receptor signaling pathway. Our analysis proposed a series of crosstalk events that could explain the role of heme in activating the TLR4 signaling pathway. In summary, the presented work opens the door to the scientific community for exploring the published knowledge on heme biology.

### Keywords: heme, hemolytic disorders, signaling pathways, knowledge graphs, biological expression language

## INTRODUCTION

Heme is an iron ion-coordinating porphyrin derivative essential to aerobic organisms (Zhang, 2011). It plays a crucial role as a prosthetic group in hemoproteins involved in several biological processes such as electron transport, oxygen transfer, and catalysis (Smith and Warren, 2009; Zhang, 2011; Kühl and Imhof, 2014; Poulos, 2014). Besides its indispensable role in hemoproteins, it can act as a damage-associated molecular pattern leading to oxidative injury, inflammation, and consequently, organ dysfunction (Jeney, 2002; Wagener et al., 2003; Dutra and Bozza, 2014). Plasma scavengers such as haptoglobin and hemopexin bind hemoglobin and heme, respectively, thus keeping the concentration of labile heme at low concentrations (Smith and McCulloh, 2015). However, at high concentrations of hemoglobin and, consequently heme, these scavenging proteins get saturated, resulting in the accumulation of biologically available heme (Soares and Bozza, 2016). With respect to hemolytic diseases, the formation of labile heme at harmful concentrations has been a subject of research for some years now (Roumenina et al., 2016; Soares and Bozza, 2016; Gouveia et al., 2017).

### Edited by:

Jyoti Sharma, Institute of Bioinformatics (IOB), India

fbioe-08-00074 March 6, 2020 Time: 14:16 # 1

### Reviewed by:

Yu Lin, Center for Devices and Radiological Health (CDRH), United States Obul Reddy Bandapalli, Hopp Children's Cancer Center Heidelberg (KiTZ), Germany

### \*Correspondence:

Daniel Domingo-Fernández daniel.domingo.fernandez@ scai.fraunhofer.de; danieldomingofernandez@ hotmail.com Martin Hofmann-Apitius martin.hofmann-apitius@ scai.fraunhofer.de Diana Imhof<EMAIL_ADDRESS>†These authors have contributed equally to this work

### Specialty section:

This article was submitted to Bioinformatics and Computational Biology, a section of the journal Frontiers in Bioengineering and Biotechnology

> Received: 14 October 2019 Accepted: 28 January 2020 Published: 06 March 2020

### Citation:

Humayun F, Domingo-Fernández D, Paul George AA, Hopp M-T, Syllwasschy BF, Detzel MS, Hoyt CT, Hofmann-Apitius M and Imhof D (2020) A Computational Approach for Mapping Heme Biology in the Context of Hemolytic Disorders. Front. Bioeng. Biotechnol. 8:74. doi: 10.3389/fbioe.2020.00074

**80**

Biomedical literature is an immense source of heterogeneous data that are dispersed throughout hundreds of journals. Furthermore, the majority of the results are scattered and published as unstructured free-text, or at best, presented in tables and cartoons representing the experimental study or biological processes and pathways. These shortcomings, combined with the exponential growth of biomedical literature, prevent the healthcare community and individual researchers from being aware of all the available information and knowledge in the literature. With the introduction of new technologies and experimental techniques, researchers have made significant advances in heme-related research and its role in the pathogenesis of numerous hemolytic diseases such as sepsis (Larsen et al., 2010; Effenberger-Neidnicht and Hartmann, 2018), malaria (Ferreira et al., 2008; Dey et al., 2012), and β-thalassemia (Vinchi et al., 2013; Conran, 2014; Garcia-Santos et al., 2017). In these diseases, large amounts of heme are released from ruptured erythrocytes and can potentially wreak havoc (Tolosano et al., 2010). Thus, it is crucial to develop new strategies that capture and exploit the vast amount of literature knowledge surrounding heme to better understand its mechanistic role in hemolytic disorders.

Biological knowledge formalized as a network can be used by clinicians as research and information retrieval tools, by biologists to propose in vitro and in vivo experiments, and by bioinformaticians to analyze high throughput -omics experiments (Catlett et al., 2013; Ali et al., 2019). Further, they can be readily semantically integrated with databases and other systems biology resources to improve their ability to accomplish each of these tasks (Hoyt et al., 2018). However, enabling this semantic integration requires organizing and formalizing the knowledge using specific vocabularies and ontologies. Although this endeavor involves significant curation efforts, it is key to the success of the subsequent modeling steps. Therefore, in practice, knowledge-based disease modeling approaches have been conducted only for major disorders such as cancer (Kuperstein et al., 2015) or neurodegenerative disorders (Mizuno et al., 2012; Fujita et al., 2014). In summary, while the scarcity of mechanistic information and the necessary amount of curation often impede launching the aforementioned approaches, modeling and mining literature knowledge provide a holistic picture of the field of interest. Furthermore, the underlying models derived from such approaches have a broad range of applications including hypothesis generation, predictive modeling and drug discovery.

Here, we present two resources aimed at assembling mechanistic knowledge surrounding the metabolism, biological functions, and pathology of heme in the context of selected hemolytic disorders. The first resource is a terminology formalizing heme-specific terms that have until now not been covered by other standard controlled vocabularies. The second resource is a heme knowledge graph (HemeKG), that is, a network comprising more than 700 nodes and more than 3,000 interactions. It was generated from 46 selected articles as the first attempt of modeling the knowledge, which is available from more than 20,000 heme-related publications. Finally, we demonstrate both resources by analyzing the crosstalk between heme biology and the TLR4 signaling pathway. The results of this analysis suggest that the activation profile for labile heme as an extracellular signaling molecule through TLR4 induces cytokine and chemokine production. However, the underlying molecular mechanism and individual pathway effectors are not fully understood and need further exploration.

## MATERIALS AND METHODS

This section describes the methodology used to generate the mechanistic knowledge graph and its supporting terminology. Subsequently, it outlines the approach followed to conduct the pathway crosstalk analysis. A schematic diagram of the methodology is presented in **Figure 1**.

## Knowledge Modeling

In order to identify recently published articles (i.e., published in the last 10 years) describing the role of heme in hemolytic disorders, PubMed was queried with the following: ("heme" AND "hemolysis") OR ("heme" AND "thrombosis") OR ("heme" AND "inflammation") AND ("2009"[Date – Publication]: "3000"[Date – Publication]). The resulting 3,108 articles were manually filtered by removing articles that were deemed too general or lacked a biochemical focus, as judged by expert opinion. After this filtering step, 6 reviews and 40 original research articles were selected for knowledge extraction and modeling. Knowledge was manually extracted and curated from this selected corpus using the official Biological Expression Language (BEL) curation guidelines from http://openbel.org/language/version\_2.0/bel\_specification\_ version\_2.0.html and http://language.bel.bio as well as additional guidelines from https://github.com/pharmacome/curation.

Next, evidence from this selected corpus is extracted and translated into BEL to generate a computable knowledge assembly model, HemeKG. In parallel to the modeling task, a terminology to support knowledge extraction of articles about the heme molecule was built. Finally, HemeKG can be used for numerous tasks such as hypothesis generation, predictive modeling and drug discovery.

Evidence from the selected corpus was manually translated into BEL statements together with their contextual information (e.g., cell type, tissue and dosage information). For instance, the evidence "Heme/iron-mediated oxidative modification of LDL can cause endothelial cytotoxicity and – at sublethal doses – the expression of stress-response genes" (Nagy et al., 2010) corresponds to the following BEL statement:

SET Cell = "endothelial cell" a(CHEBI:"oxidised LDL") pos bp(MESH:"Cytotoxicity, Immunologic").

## Generation of a Supporting Terminology