Pathway-based network analyses and candidate genes associated with Kashin-Beck disease.

To perform a comprehensive analysis focusing on the biological functions and interactions of Kashin-Beck disease (KBD)-related genes to provide information towards understanding the pathogenesis of KBD.A retrospective, integrated bioinformatics analysis was designed and conducted. First, by reviewing the literature deposited in PubMed, we identified 922 genes genetically associated with KBD. Then, biological function and network analyses were conducted with Cytoscape software. Moreover, KBD specific molecular network analysis was conducted by Cytocluster using the Molecular Complex Detection Algorithm (MCODE).The biological function enrichment analysis suggested that collagen catabolic process, protein activation cascade, cellular response to growth factor stimulus, skeletal system development, and extrinsic apoptosis played important roles in KBD development. The apoptosis pathway, NF-kappa B signaling pathway, and the glutathione metabolism pathway were significantly enriched in the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway network, suggesting that these pathways may play key roles in KBD occurrence and development. MCODE clusters showed that in top 3 clusters, 54 of KBD-related genes were included in the network and 110 candidate genes were discovered might be potentially related to KBD.The 110 candidate genes discovered in the current study may be related to the development of KBD. The expression changes of apoptosis and oxidative stress-related genes might serve as biomarkers for early diagnosis and treatment of KBD.

Introduction

Kashin-Beck Disease (KBD) is a type of endemic, chronic, and disabling osteoarthropathy with high disability. KBD is also named as endemic osteochondropathy, as it is mainly prevalent at defined areas of North Korea, Russia, and China (mainly in northeastern and southwestern China).[ By 2015, there were 0.57 million KBD patients in China. Of these patients, nearly 13 thousand were children under 13 years old.[ KBD patients usually exhibit the clinical symptoms, such as enlargement of phalanges or other joints, joint pain, limited movement, deformed limb joints.[ Besides the detrimental effect on all aspects of quality of life, KBD places a heavy economic burden on society as a whole.[

Previous studies reported that environmental factors played important roles in KBD development, such as selenium deficiency in the environment of endemic KBD areas and cereal contamination by fungal mycotoxins.[ However, the etiology and pathogenesis have not been fully unraveled. In recent years, growing evidence has implicated that KBD is a complex disease due to a complicated interplay of environmental and genetic factors. Some studies suggested that some genes might contribute to the development of KBD, such as GPX1,[GPX4,[SELS,[SEPP,[ADAM12,[Col2A1,[Bcl-2, Bax, and Fas.[ But these studies only provided a limited explanation for the KBD occurrence and insufficiently revealed their roles and functions in the development of KBD. In the whole genome-wide, these genes do not play independent roles; instead, these genes form biological function and pathway networks through their intricate interactions, which can give us a more complete picture of their genetic function contributing to the development of KBD.[ Furthermore, genetic studies have indicated that the pathological molecular mechanisms may be attributed to hundreds of genes and their variants for KBD.[ With the rapid development of bioinformatics, such as biological function enrichment analysis, protein–protein interaction (PPI) and the Molecular Complex Detection Algorithm (MCODE), an increasing number of differentially expressed genes (DEGs) were identified between KBD cases and healthy controls.[ Therefore, a comprehensive analysis of potential genes, which uncovers more detailed information involved in pathways and networks as compared to conventional single-gene analysis, is necessary.

In the present study, we first formed a comprehensive collection of genes that were reported to be associated with KBD. Biological function enrichment analysis, based on the KBD-related genes, was then conducted to identify their significant biological functions. We also constructed a pathway network to explore possible crosstalk among the significant pathways. Moreover, we made a PPI network of these KBD-related genes using MCODE.[ This study may provide further information for understanding the molecular mechanisms of KBD based on a systematic bioinformatics analysis.

Materials and methods

Collection of KBD-related genes

KBD-related genes identified in human genetic association studies in PubMed (http://www.ncbi.nlm.nih.gov/pubmed/) were collected using the term (Kashin-Beck Disease or KBD) AND (polymorphism [MeSH] or genotype [MeSH] or alleles [MeSH]).[ KBD-related genes were retrieved and screened under the following criteria:

they must be reported in publications which matched the above retrieve strategy;

the KBD-related genes should be provided in official symbols;

genes that exhibited significant associations (P < .05) with KBD were included and that showed insignificant associations (P > .05) with KBD were excluded.

By July 19, 2017, a total of 58 publications that were matched the criteria were retrieved. We further included genes with a significant association with KBD (P < .05) after reading the full text or abstracts of the original articles. A KBD-related gene set with 922 members that was reported significantly associated with KBD was obtained from 24 studies. For microarray analyses or genome-wide association studies (GWAS) at a genome-wide significance level, the comprehensive and positive genes were obtained from the supplementary data of these studies.

Enrichment analysis of KBD-related genes

Gene Ontology Term (GO; www.geneontology.org) and Kyoto Encyclopedia of Genes and Genomes (KEGG; www.genome.ad.jp/kegg) pathway enrichment analyses were performed using Cytoscape 3.5.1. GO is usually used to explore the molecular function (MF), biological process (BP) and cellular component (CC) of candidate genes. In our present study, we mainly focused on BP of the candidate genes. KEGG is widely used for extracting the pathway information based on the molecular interactions. Fisher's exact test (two-side) or x2 test were performed to classify the pathway category. The false discovery rate (FDR) was used for the P-value correction based on the Benjamini and Hochberg method.[ Significantly enriched pathways and GO terms (P ≤ .05, number of enrichment genes ≥2) were identified through Cytoscape 3.5.1.

Establishment of pathway networks among KBD-related genes

The pathway network among KBD-related genes was established to understand their possible interactions and crosstalk. Due to differences in the sensitivity and specificity of detection methods or sources, we constructed the pathway networks based on the gene lists from GWAS, microarray analysis, and exome sequencing using ClueGO, which is a plugin for Cytoscape 3.5.1. The ClueGO network is analyzed with kappa statistics based on the similarity of their associated genes, which reveals the interactions among the pathways. The significance of the terms and groups is automatically calculated from Gene Ontology, KEGG, Reactome, and CluePedia by ClueGO based on Cytoscape 3.5.1.[

Construction of KBD-related PPI network via MCODE

In the present study, we first identified the PPI network relationships of the 922 KBD-related genes in the whole network based on the biological network from GeneMANIA, which used many large, publicly available biological datasets to find related genes. MCODE cluster analysis was then performed using Cytocluster software (degree cutoff = 2, node score cutoff = 0.2, K-core = 2, MAX depth = 100) to identify the most important MCODE clusters according to clustering scores. Finally, the screened cluster networks were visualized with Cytoscape 3.5.1 software. The importance of clustering networks was assessed by clustering scores.

Results

Identification of KBD-related genes from the literature

By carefully and comprehensively retrieving articles matching our criteria from the PubMed database, we screened genetic association studies related to KBD and selected only the publications that found gene(s) significantly associated with KBD. Original manuscripts reporting a negative or insignificant association were excluded. Finally, a KBD-related gene set consisting of 922 members reported to be significantly associated with KBD was retrieved from 24 studies (Table 1). The list of 922 genes is provided in Supplemental File 1.

Table 1

Sources of KBD-related genes.

Biological function and pathway enrichment analysis towards KBD-related genes

Results showed that the KBD-related genes were involved in 672 biological function terms in total. The top five biological function categories enriched in KBD-related genes detected by exome sequencing, GWAS and microarray analysis are listed in Table 2. The results showed that biological functions (collagen catabolic process, protein activation cascade, cellular response to growth factor stimulus, skeletal system development, extrinsic apoptosis, et al) were enriched among KBD-related genes, which indicates that collagen formation, skeletal development, and cellular programs (growth, differentiation, migration, apoptosis, et al) may affect KBD.

Table 2

Top five biological functions enriched in KBD-related genes detected by exome sequencing, GWAS and microarray analysis.

Among the significantly enriched pathways for KBD, several pathways related to cell death and apoptosis, for example, Sa_programmed cell death, Reactome_intrinsic pathway for apoptosis, Biocarta_WNT pathway, extracellular matrix (ECM) receptor interaction, were enriched in the KBD-related genes. Also, pathways related to immune or inflammation-related signaling were identified, such as, Adherens junction, Regulation of actin cytoskeleton, and Reactome signaling in immune system, suggesting that the immune system is also involved in the development of KBD. In addition, oxidative stress-related signaling pathways, such as Glutathione metabolism, were also significantly enriched, which indicates that KBD development is related to the body's oxidative stress response. Top five pathways enriched in KBD-related genes detected by exome sequencing, GWAS and microarray analysis were showed in Table 3.

Table 3

Top five pathways enriched in KBD-related genes detected by exome sequencing, GWAS and microarray analysis.

KEGG pathway network of KBD-related genes

The pathway network showed that significantly enriched KEGG pathways related to KBD were divided into 15 dependent clusters. Pathways in the same cluster exhibit the same or similar functions; however, the different clusters were relatively independent, suggesting they may exhibit different functions. Top three clusters with most related genes were exhibited in Fig. 1. Apoptosis pathway, NF-kappa B signaling pathway, and glutathione metabolism pathway were enriched in the predominant three clusters.

Figure 1

Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway network for Kashin-Beck disease (KBD)-related genes. A: From genes detected by genome-wide association studies (GWAS); B: From genes detected by microarray analysis; C: From genes detected by exome sequencing; pathways in the same cluster exhibit the same or similar functions; apoptosis pathway, NF-kappa B signaling pathway, and glutathione metabolism pathway were enriched in the predominant three clusters.

Construction of KBD-specific protein network via MCODE

The topological properties of PPI via pathway-based MCODE cluster analysis were performed to help understand the potential biological mechanisms associated with the network. The higher the clustering scores, the more important the biological function of this clustering in the development of KBD. Therefore, we analyzed the three clusters with the highest clustering scores in detail. The results are shown in Fig. 2. In Fig. 2A, the PPI cluster (MCODE cluster score = 6.114) had 36 genes in total, including 24 (66.67%) gene members in the KBD-related gene list; while 12 (33.33%) genes were not included in the list. In Fig. 2B, the PPI cluster (MCODE cluster score = 5.077) had 40 genes in total, including 25 (62.50%) gene members in the KBD-related gene list, but 15 (37.50%) genes were not in the list. In Fig. 2C, the PPI cluster (MCODE cluster score = 4.184) had 88 genes in total, in which 14 (15.91%) gene members were KBD-related but 74 (84.09%) genes have not been reported before. In these 3 PPI clusters, 54 of the KBD-related genes were included in the human interactome network, among which 110 candidates genes that are likely to be highly associated with KBD based on MCODE clusters. These 110 genes were listed in Table 4 , which provided a list of new potential candidates for KBD.

Figure 2

Kashin-Beck disease (KBD) top three specific PPI networks. A: From genes detected by genome-wide association studies (GWAS); B: From genes detected by microarray analysis; C: From genes detected by exome sequencing; the red nodes mean genes in the KBD-related gene list, the green nodes are genes not in the KBD-related gene list.

Table 4

Genes included in KBD top three specific PPI networks but not in the KBD-related gene list.

Discussion

By July 19, 2017, there were 24 articles on gene-related studies about KBD, and we assembled a list of 922 genes that were reported statistically associated with KBD. The detection methods used in these studies were mainly exome sequencing, microarray analysis, and GWAS. However, most of the biological effects of the 922 genes potentially involved in KBD development have not yet been evaluated. Therefore, a thorough understanding of the BPs related to the molecular pathogenesis of KBD is still far from complete.[ Fortunately, the developments of high-throughput microarray analysis and sequencing techniques in recent years provide an opportunity to understand the pathogenesis of KBD at a systems biology level based on the biomarkers that have been discovered.

Pathways related to apoptosis play an important role in KBD occurrence and development

Our biological function enrichment analysis identified specific BPs and pathways in which KBD-related genes are involved. For instance, collagen catabolic process, protein activation cascade, cellular response to growth factor stimulus, skeletal system development, and extrinsic apoptosis were significantly enriched among the KBD-related genes, indicating that biological events such as the maintenance of apoptosis, the degradation of ECM, and biology process related to oxidative stress play key roles in the development of KBD. Studies had reported that KBD is characterized as specific pathological changes in cartilage, causing apoptosis and necrosis of chondrocytes in the epiphyseal plate,[ which further demonstrated by our results. The current view about what factors lead to articular chondrocyte apoptosis in KBD focuses on varying degrees of oxidative stress occur in the articular chondrocytes, which is primarily caused by articular chondrocyte apoptosis.[ Our results also provided evidence for this view. Our enrichment results showed that oxidative stress-related signaling pathways such as glutathione metabolism, were also significantly enriched in KBD-related genes, which means that KBD development is also associated with the body's oxidative stress response. Due to the above results, we have more evidence to infer that a reduction in oxidative stress may reduce chondrocyte apoptosis, and thus prevent the occurrence of KBD.[ In addition, the agreements between our results and previous work indicated that our results are reliable.

Inflammation and immune-related pathways may be also associated with KBD occurrence and development

The KEGG network showed that the key pathways enriched in genes identified by GWAS were complement and coagulation cascades pathways, which are central to host defense. The complement pathway is part of innate immune system and the coagulation pathway is essential for clot formation and the prevention of excessive bleeding.[ The complement and coagulation cascades generate crucial enzymatic activities that, in turn, generate the complement effector molecules, with the opsonization of pathogens, the recruitment of inflammatory and immunocompetent cells, and the direct killing of pathogens as the main consequences.[ The results further consolidate ties between the pathology of KBD and immune-specific activities.[ In addition, the KEGG network also showed that dilated cardiomyopathy pathway was enriched in genes from exome sequencing. Dilated cardiomyopathy pathway represents the final common morphofunctional pathway of various pathological conditions, in which a combination of myocyte injury and necrosis associated with tissue fibrosis results in impaired mechanical function. To our knowledge, this is the first time that KBD has been suggested to have a close biological association with heart disease, which still needs further evaluation to confirm.

Apoptosis, cell structure and oxidative stress-related events might affect KBD occurrence at protein level

Top three important PPI networks of the KBD-related genes were extracted from reference interactome network (GeneMANIA network) by MCODE clustering, which is a density-based non-overlapping clustering algorithm and based on the seed nodes as the center to expand to its neighbor nodes. Then the nodes that may interact with the seed nodes were screened out to construct complexes in a PPI network.[ MCODE can also identify clusters (highly interconnected regions) within a network.[ It is worth noting that 110 extended genes were not included in the KBD-related genes, but appeared in the top three PPI networks, which were plausible genes that previously not been reported to be associated with KBD. For example, BAD (BCL-2 associated agonist of cell death) is a member of the BCL-2 family,[ which is known as a regulator of programmed cell death. Previous studies have reported that BAD can positively regulate cell apoptosis by forming heterodimers with BCL-xL and BCL-2.[COL1A2, COL4A2, COL4A3, COL4A5, and COL4A6 were identified as potential target genes in KBD, and these are related to the synthesis of type I collagen and strengthening as well as supporting many tissues, including cartilage and bone.[MAPK14 is a member of the MAP kinase family, which acts as integration points for multiple biochemical signals, and are involved in cell proliferation, differentiation, transcription regulation, and development. MAPK14 is activated by oxidative stress and proinflammatory cytokines, but this activation requires phosphorylation of MAPK14 by MAP kinase kinases (MKKs).[ Genes play their crucial biological roles in the form of proteins. Therefore, the above events or functions, which were reflected by the PPI network, may be closer to the practical biology processes of the KBD-related genes. Our results extend the findings of the previous studies, and provide further information about the pathogenesis of KBD. Meanwhile, our MCODE cluster analysis was based on a reference interactome network, which not only provided a meaningful inferred network of KBD-related genes but also can be used to identify potential candidate biomarkers.

According to our knowledge, our present study is the first report to integrate all the candidate genes from the previous studies together and perform a comprehensively analysis in order to answer the above key questions (which candidate gene is more important in KBD development, what are the possible signal pathways involved when these candidate genes play their roles, et al) as well as to provide more valuable information. Exploring the environmental susceptibility genes of KBD has become a hot topic in this field. The conclusions of this study may help to update the new understanding of the pathogenesis of KBD and to expand the biomarkers of KBD, which helps early identification and treatment of KBD.

In conclusion, 110 candidate genes discovered in the present study may be related to the development of KBD. The expression changes of apoptosis and oxidative stress related genes might serve as biomarkers for early diagnosis and treatment of KBD. We will further verify our findings through rigorous and independent experiments in big KBD populations.

Acknowledgments

We would like to thank Dr Zhenzhong Li from Beijing Compass Biotechnology Co., Ltd (Beijing, China; http://www.kangpusen.com/) for his kind help with the data analysis. Our study was supported by grants from the National Natural Science Foundation of China (grant nos. 81773372 and 81573104).

Author contributions

RQZ and YMX designed the study. RQZ, HG, and XLY conducted the analyses. RQZ, BRL, DDZ, ZFL, and HLJ performed the analyses and prepared the manuscript. AMY and JZ revised the manuscript. RQZ and YMX submitted the study.

Conceptualization: Yongmin Xiong.

Data curation: Rongqiang Zhang, Hao Guo.

Formal analysis: Rongqiang Zhang, Hao Guo, Xiaoli Yang, Dandan Zhang, Baorong Li, Zhaofang Li.

Funding acquisition: Yongmin Xiong.

Methodology: Rongqiang Zhang.

Software: Rongqiang Zhang.

Supervision: Yongmin Xiong.

Writing – original draft: Rongqiang Zhang, Yongmin Xiong.

Writing – review & editing: Rongqiang Zhang, Hao Guo, Xiaoli Yang, Dandan Zhang, Baorong Li, Zhaofang Li, Yongmin Xiong.

Table 4 (Continued)

Genes included in KBD top three specific PPI networks but not in the KBD-related gene list.