Network Analysis of Osteoarthritis Progression Using a Steiner Minimal Tree Algorithm

Introduction

Osteoarthritis (OA) is the leading cause of disability and source of societal cost in older adults.¹ Risk factors associated with OA include joint injuries, aging, and obesity.² It is projected that by 2040, OA will affect 78 million people.³ The impact of OA on individuals includes pain, loss of mobility, and loss of independence, with 25% of patients unable to perform normal daily activities. OA is characterized by various factors including subchondral bone remodeling,⁴ degeneration of the meniscus,^5,6 inflammation, and fibrosis of the infrapatellar fat pad and synovial membrane.⁷

OA will remain a large and growing global problem for clinical and public health systems. Although OA is a significant medical condition, there are currently no registered therapies available that can stop the structural damage caused by the disease. Symptom-modifying interventions only offer moderate long-term effects, at best.⁸ To advance in the treatment of OA, it is crucial to have a deep understanding of the cellular and molecular pathophysiology.

Several genes and pathways associated with mechanical, inflammatory, and synovial factors have been identified in OA.⁹ The identification of the possible essential genes such as MMP-3, MMP-13, VEGF, TIMP-1, TGFB1, FGF18 has important influence on molecular causes.^10–13

Although studies have identified important pathogenic genes in the development of OA, the underlying molecular mechanism remains unclear.¹⁴ Lots of genes has effects on OA progress rather than one or two major genes with large effects.¹⁵ Network and pathway analysis reveals crosstalk between genes. So, a comprehensive analysis of potentially causal genes can provide many important insights finding.¹⁶ We used a comprehensive analysis of genes associated with OA. Then, functional enrichment analyzes explore the significant biological themes within these genetic factors. We intended to analyze the topological characteristics of related genes in the context of a human protein-protein interaction network. To accomplish this, we first inferred the specific molecular network using the Steiner minimal tree algorithm. By isolating the network specific to Osteoarthritis from the larger human interaction group network, we gained insights into the potential pathological molecular network associated with the Osteoarthritis gene set. Additionally, we compared the significant RNA findings with previous results obtained from single-cell RNA-seq analysis.¹⁷ This comprehensive analysis provides systematic insights into the mechanisms underlying OA.

Materials and Methods

Identification of OA-Related Genes

Candidate genes associated with OA were curated by retrieving the human genetic association studies deposited in MEDLINE (http://www.ncbi.nlm.nih.gov/pubmed/). We searched OA with the term (osteoarthritis [MeSH]) and (polymorphism [MeSH] or genotype [MeSH] or alleles [MeSH]) not (neoplasms [MeSH]). Non-English articles studies were excluded. Two authors screened the titles and abstracts of the articles, and the full text of all eligible studies was reviewed for inclusion and exclusion criteria. Any discrepancies in literature screening were resolved by the third author. By 2020–12-01, we found 1203 publications for OA. Then, we collected the association studies after reviewing the abstracts, and narrowed the selection publication by focusing on OA. We only included publication that results support the conclusions. The negative or insignificant associations reports were excluded. The genes significantly associated with OA in these studies were selected for this study.

Functional Enrichment Analysis of OA-Related Genes

The functional features of the OA-related genes were examined by ClusterProfiler based on Gene Ontology (GO), BioCarta and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases. GO terms of biological processes with a p-adjust value smaller than 0.01 were kept as the significantly enriched ones. All the pathways with one or more genes overlapping the candidate genes were extracted. A p-value was assigned to each pathway by Fisher’s test to indicate the significant overlap between the pathway and the input gene. The pathways were significantly enriched with a p-adjust < 0.01.

Pathway Crosstalk Analysis

We further performed pathway crosstalk following the former way.¹⁸ To describe the overlap between any given pair of pathways, two measurements were computed. The Jaccard Coefficient (JC) \ (JC=|\frac {A\cap B {A\cup B} |\) and the Overlap Coefficient (OC) \ (OC=\frac {|A\cap B|} {min (|A|, |B|)} \). Where A and B are the lists of genes included in the two tested pathways. To construct the pathway crosstalk, we implemented the following procedure:

Select a set of pathways for crosstalk analysis. Only the pathways with a p-adjust < 0.01 were used. Meanwhile, the pathways containing less than six candidate genes were removed.

Count the number of shared candidate genes between any pair of pathways. Pathway pairs with less than seven overlapped genes were removed.

Calculate the overlap of all pathway pairs. All the pathway pairs were ranked according to their JC and OC values. The overlapping level between the two pathways was measured by the JC and OC average scores. The pathways could be grouped into one major module by their crosstalk analysis results.

Visualize the selected pathway crosstalk with the packages cytoscape.¹⁸

Construction of the Human Interactome

We constructed a comprehensive and reliable human interactome to investigate the correlation between the genes. First, we downloaded the human protein-protein interaction (PPI) data from the Protein Interaction Network Analysis (PINA) platform (release version: 21 May 2014). The database pooled and curated non-redundant physical interaction data from six databases: IntAct, BioGRID, DIP, HPRD, MINT, and MIPS/MPact. After excluding the redundant and self-interacting pairs and using org.Hs.ens.db R packages to map these interactome data onto gene SYMBOL, we constructed a comprehensive human physical interactome by merging the two data sets, which contained 17,457 nodes and 233,333 edges.

Construction of OA-Specific Network

We intended to explore the disease subnetwork extraction. We used the Steiner Net to extract specific potential pathological networks and take the collected OA genes as input seeds. Then, we generated 1000 random networks with the same number of interactions as the OA-specific network using the Erdos-Renyi model in the R igraph package to assess the non-randomness. Next, we calculated the average values of the shortest-path distance and clustering coefficient. The significance level of non-randomness needs two requirements; 1) the number of random networks with average shortest-path distance (ND) was smaller than the OA-specific networks, 2) the number of random networks with an average clustering coefficient (NC) was higher than the observed clustering coefficient. Finally, the empirical p-value was calculated by using ND/1000 and NC/1000.

Single-Cell RNA-Seq Analysis Results Identify

The important RNA was compared with former single-cell RNA-seq analysis results.¹⁷ Quanbo Ji has identified seven molecularly defined populations of chondrocytes in the human OA cartilage, including ECs, effector chondrocytes (RegCs), regulatory chondrocytes (ProCs), proliferative chondrocytes (preHTCs), pre-hypertrophic chondrocytes (FCs), fibrocartilage chondrocytes (HTCs), hypertrophic chondrocytes (HomCs) and homeostatic chondrocytes. They also presented gene expression profiles at different OA stages at single-cell resolution. We identified our crosstalk results from these expression profiles.

Results

Identification of Genes Associated with OA

(Supplementary Table 1) provided all genes related to OA. Altogether, a gene set with 177 members significantly associated with OA was collected from 187 studies.

Significant Enrichment Pathways for OA

Biological Functions and Pathway Enrichment Analysis Enriched in Gene Set Some GO terms significantly enriched in the candidate genes (Supplementary Table 2), including those associated with skeletal system development, cytokine-mediated signaling pathway, inflammatory response, cartilage development, extracellular matrix organization, negative regulation of cell population proliferation, ossification, positive regulation of tyrosine phosphorylation of STAT protein and positive regulation of epithelial to mesenchymal transition. We found 72 significant enrichment pathways for OA (Supplementary Table 3) through.

Crosstalk Among Significantly Enriched Pathways

To understand how these 72 pathways interact with each other, we performed a pathway crosstalk analysis. A total of 66 pathways were containing six or more members in the OA gene set. A number of 44 pathways met the criterion for crosstalk analysis. Each pathway shared at least seven genes with one or more other pathways (Figure 1).

Figure 1 Pathway crosstalk among significantly enriched pathways Nodes represent pathways and edges represent crosstalk between pathways. To understand how these 72 pathways interact with each other, we performed a pathway crosstalk analysis. A total of 66 pathways were containing six or more members in the OA gene set. A number of 44 pathways met the criterion for crosstalk analysis. Each pathway shared at least seven genes with one or more other pathways.

Gene Set Network Topological Characteristics

PPI analyzed the topological properties of nodes and interactions. The degree of the constructed network was analyzed in the context of the constructed human interactome. We intended to analyze the two other gene sets network features including the osteoporosis list with 323 members and the cancer genes with 506 genes (Supplementary Table 4). The nodes with specific degrees were scattered in a range from 1 to 836. The nodes were measured by the number of genes connecting with a given gene. The 151 of 177 genes mapped onto the human interactome network with the 52.51 mean degree, which each gene connected with 52.51 other genes on average. The 300 of 323 genes could be mapped onto the human interactome in osteoporosis, with 54.05 average degree. The 485 out of 506 had the corresponding nodes in the cancer genes with the 74.97 average degree. Further, 52.98% (80/151) and 51% (153/300) genes fell in the degree interval of 1–20, respectively for OA and Osteoporosis, while the cancer genes only 30.31%.

OA-Specific Molecular Network Inference

Steiner minimal tree algorithm linked the maximal and minimal members of genes nodes. The OA network contained 182 nodes and 181 edges (Figure 2). For these random subnetworks, the mean shortest-path distance was 6.51, significantly larger than that of the OA-related network (shortest-path distance, 4.06; empirical p=0). The average clustering coefficient of the random networks was 0.01, statistically significantly less than the OA distinctive network (clustering coefficient, 0; empirical p=0.72). As specified, 151 of the 177 genes in the OA gene set were included in the extracted OA-specific network. A total of 85.31% of the genes in the OA gene set and 82.96% of 182 genes in the OA-specific network, indicating a high coverage of the OA gene set in the subnetwork. The remaining 31 genes in the OA-specific network were not excluded, as they were non-input genes screened from the Steiner minimum network. These 31 genes are not part of the Osteoarthritis gene concentration, but they are involved in the biological process and have close interactions with each other. These genes provide potential candidate genes for further exploration in the OA biological processes. In addition, the expression of CLU, ENO1, SRRM1, UBC, HMGB1, NR3C1, NOTCH2NL, and CBX5 were significantly increased in seven molecularly defined populations of OA cartilage, which is consistent with the results of former single-cell RNA-seq analysis, making the consistency more convincing. Figure 3 presents violin plots depicting the expression levels of specific marker genes across seven different populations in human OA cartilage. These populations include effector chondrocytes, regulatory chondrocytes, proliferative chondrocytes, prehypertrophic chondrocytes, fibrocartilage chondrocytes, hypertrophic chondrocytes, and homeostatic chondrocytes.

Figure 2 We utilized the Steiner minimal tree algorithm to construct a specific network for osteoarthritis (OA) from the human interactome network. This algorithm effectively connected the maximal and minimal members of gene nodes. The resulting OA network consisted of 182 nodes and 181 edges.

Figure 3 Violin plots displaying the expression levels of candidate marker genes for the seven distinct molecularly defined populations within human OA cartilage are presented. The genes showcased in each figure are as follows: (a) CBX5 gene, (b) CLU gene, (c) CTNNB1 gene, (d) ENO1 gene, (e) HMGB1 gene, (f) NOTCH2NL gene, (g) NR3C1 gene, (h) SRRM1 gene, (i) UBC gene.

Discussion

OA is characterized by cartilage degeneration, synovial and infrapatellar fat pad inflammation, and subchondral bone remodeling.^19,20 In this study, we aimed to conduct a comprehensive collection of pathogenic associated with OA by including three databases. We found the intersection of all the signaling pathways through the GO and the KEGG functional enrichment analyses. Then, we further analyzed these novel genes in the context of the human OA protein-protein interaction network. Next, the Steiner minimal tree algorithm was used in OA-specific molecular network analysis.

Former research¹⁷ has identified seven molecularly defined populations of chondrocytes in the human OA cartilage, including three novel phenotypes with distinct functions. It presented gene expression profiles at different OA stages at single-cell resolution. These genes had close interaction with genes known to be related to the biological processes involved in OA.²¹ Thus, these genes provided a list of potential candidates for further exploration. After reviewing all these genes, we finally identified 31 important genes. Notably, CLU, ENO1, SRRM1, UBC, HMGB1, NR3C1, NOTCH2NL, and CBX5 showed significantly increased expression in OA chondrocytes.

Clusterin (CLU) has been found as a ubiquitous glycoprotein and its upregulation has been found in OA cartilage.²² Also known as apolipoprotein J, CLU acts as an ATP-independent holdase chaperone, preventing protein aggregation and proteotoxicity.²³ It is a disulfide-linked heterodimeric protein of approximately 60 kDa, involved in the clearance of cellular debris and the regulation of apoptosis. Recent studies suggest that the intracellular form of CLU may suppress stress-induced apoptosis and the secreted form of CLU functions as an extracellular chaperone that prevents protein aggregation.²⁴ Increased levels of CLU mRNA and protein in the joint environments of knee OA may reflect the severity of the condition. Csaba Matta’s findings demonstrate that the secretion of clusterin can be regulated by interleukin-1β and tumor necrosis factor-α, influencing cartilage degradation.²⁵ Furthermore, three neopeptides associated with Clusterin show potential as markers for OA disease.²⁶

Enolase 1(ENO1) has been implicated in autoimmune and inflammatory diseases.²⁷ ENO1 has been identified as a novel stimulatory receptor on monocytes.²⁷ ENO1 is a multifunctional glycolytic enzyme expressed in the cell-surface and has been found to be increased rapidly in response to inflammatory loops in rheumatoid arthritis.²⁸ Weijuan Ma investigated the protein changes in Kashin-Beck disease cartilage and founded the candidate proteins α-enolase.²⁹ Additionally, proteomic analysis of human osteoarthritic chondrocytes reveals protein changes, and Western blotting and immunohistochemistry in tissue cartilage were used to identify a decrease in ENO1.³⁰

Evidence supports increased expression of high mobility group box 1 (HMGB1) in OA cartilage.³¹ HMGB1 protein is associated with chemoattractant properties.³² Chondrocytes release HMGB1, which has migratory effects on chondrogenic progenitor cells.³³ Additionally, increased expression of HMGB1 has been observed in the synovium of rat models with anterior cruciate ligament transection-induced KOA.³⁴ HMGB1 is a downstream product of pyroptosis and may have potential inflammatory regulating effects. Chromobox 4 (CBX4), a component of polycomb repressive complex 1, is involved in maintaining cell identity and organ development through gene silencing. The expression of CBX4 has been shown to attenuate the development of OA in mice.³⁵

Ubiquitin C(UBC) is a gene encoding ubiquitin precursor in mammals.³² Conjugation of ubiquitin monomers or polymers can lead to diverse effects within a cell, depending on the residues to which ubiquitin is conjugated.³⁶ A study suggests that UBC could be a promising therapeutic target for ovarian cancer patients with recurrent ubiquitin B(UBB) silencing.³⁷ Hypoxia-inducible factor 1α (HIF-1α) plays a crucial role in preventing apoptosis in chondrocytes. The protein UCHL1 can inhibit apoptosis in chondrocytes by increasing the levels of HIF-1α, which in turn promotes mitophagy and preserves mitochondrial function.³⁸ Bobin Mi reported one of the top 10 central genes of female patients with OA was UBC, and the upregulation of UBC is a therapeutic strategy for maintaining the normal function of cells in the synovium.³⁹

Notch is a single-pass transmembrane cell surface receptor, it regulates differentiation and apoptosis during embryogenesis, as well as in various developmental systems such as neurogenesis and hematopoiesis.⁴⁰ The proliferation and differentiation of disc cells are dependent on the Notch signaling pathway.⁴¹ In a study, miR-485-3p was found to promote proliferation and prevent apoptosis, ECM degradation, inflammation, and oxidative stress in OA chondrocytes by inhibiting Notch2 and the NF-jB pathway.⁴² The inhibition of Notch2 expression resulted in increased proliferation and reduced apoptosis rate.⁴² It has been reported that the activation of Notch2 in specific chondrocytes can stimulate the progression of OA.⁴³ Notch2 serves as the receptor for the canonical NF-jB pathway, which is essential for chondrocytes to express MMPs, release ECM components, trigger inflammatory responses, and further contribute to cartilage destruction.⁴⁴

Serine and arginine repetitive matrix 1(SRRM1), also known as Srm160,⁴⁵ is an SR-related protein.⁴⁶ SRRM1 was initially identified as a nuclear matrix antigen highly concentrated in interphase nuclear “speckle” domains, which are enriched in splicing components.⁴⁷ It functions as a coactivator of both constitutive and exon enhancer-dependent splicing by forming cross-intron interactions with multiple splicing factors bound directly to pre-mRNA.⁴⁶ SRRM1 has been found to be correlated with tumor malignancy in patients with poor clinical prognosis.⁴⁸

We first made these systems biology-based analyses to explore complicated OA phenotypes. Steiner minimal tree algorithm results might provide m important insights beyond the conventional single-gene analyses and evaluated. This analysis provided important findings to understand OA progress in a systems biological view.

Although no OA genetic study has yet led to a diagnostic tool or treatment, it has been observed that these genes show a significant increase in expression in seven distinct molecularly defined populations of OA cartilage. Therefore, the current gap between discovery and utility is expected to narrow.

Conclusion

The Steiner tree-based approach identified new biological molecules that are associated with OA genes. Specifically, the CLU, ENO1, SRRM1, UBC, HMGB1, NR3C1, NOTCH2NL, and CBX5 genes were found to be associated with the progression of OA cartilage. In this study, it was observed that these genes exhibited significantly increased expression in seven distinct molecularly defined populations of OA cartilage.