Identification of key genes and pathways associated with sex difference in osteoarthritis based on bioinformatics analysis.

OBJECTIVES: The present study aimed to identify different key genes and pathways between postmenopausal females and males by studying differentially expressed genes (DEGs).
METHODS: GSE32317 and GSE55457 gene expression data were downloaded from the GEO database, and DEGs were discovered using R software to obtain overlapping DEGs. The interaction between overlapping DEGs was further analyzed by establishing the protein-protein interaction (PPI) network. Finally, GO and KEGG were used for enrichment analysis.
RESULTS: 924 overlapping DEGs between postmenopausal women and men with osteoarthritis (OA) were identified, including 674 up-regulated genes and 249 down-regulated ones. And 10 hub genes were identified in the PPI network, including BMP4, KDM6A, JMJD1C, NFATC1, PRKX, SRF, ZFX, LAMTOR5, UFD1L and AMBN. The findings of the functional enrichment analysis suggested that these genes were predominantly expressed in MAPK signaling pathway as well as the Thyroid hormone signaling pathway, indicating that those two pathways may be involved in onset and disease progression of OA in postmenopausal patients.
CONCLUSION: BMP4, KDM6A, JMJD1C, PRKX, ZFX and LAMTOR5 are expected to play crucial roles in disease development in postmenopausal patients and may be ideal targets or prognostic markers for the treatment of OA.

Introduction

Osteoarthritis (OA), the most common joint disorder, is characterized by cartilage deterioration, synovitis, subchondral bone sclerosis and persistent pain[1-4]. Limited movement of joints often occurs in severe cases. Other than relieving patient’s pain, there is no effective treatment that can delay or stop the disease from progressing. Joint replacement surgeries are only offered to patients with advanced stage of OA.

It is estimated that about 18% of females and 6% of males aged 60 or older have OA worldwide[5]. The incidence of OA is about 3 times higher in women than in men. The fact that the disease is more prevalent in female patients than in male patients suggests that sex hormones have a role in OA’s pathogenesis; hence, estrogen is thought to be implicated in the development of the disease. Endogenous estrogen was linked to radiological OA and cartilage turnover in a comprehensive retrospective analysis of 27 trials, while testosterone was linked to cartilage volume. Furthermore, there has been evidence of a link between exogenous estrogen and cartilage and bone turnover. However, it’s still unknown how estrogen could alter radiological and MRI structure, and joint replacement as well. In addition, an association has also been found between gene polymorphisms of estrogen receptors α and β, and OA[6].

High-throughput cellular microarray platforms have been increasingly employed for profiling gene expression. Studies have been conducted extensively to identify or evaluate specific genes and pathways involved in OA progression, but most have focused on comparing OA patients’ gene expression to healthy individuals. Literature on patients’ gene expression related to sex differences is limited.

We aimed to identify the genes with differential expression (DE) in postmenopausal females and males with OA through comprehensive bioinformatics analysis. The Gene Expression Omnibus (GEO) database was used to get gene expression profiles from postmenopausal females and males with OA. Then protein interaction analysis was utilized to determine hub genes of DEGs. This study contributes to a better knowledge of sex-specific mechanisms of OA, as well as prospective biomarkers and therapeutic targets to treat OA.

Materials and Methods

Data source

In this study, we searched two independent OA data sets, GSE55457 and GSE32317, using keywords “osteoarthritis” [MeSH Terms] OR “osteoarthritis” [All Fields] AND “Homo sapiens”[porgn] AND “Expression profiling by array” [Filter]. The two datasets were obtained from two independent cohorts, including synovial samples from females and males with OA. Samples from female and male OA patients were screened for different genes.

Data preprocessing and identification of DEGs

R software (v. 4.1.1) and associated R package was used to compare DEGs between female and male OA patient groups in GSE55457 and GSE32317, respectively. Firstly, background-adjustment was performed for the two datasets, followed by normalization using log2 transformation, as shown in Figure 1. Background correction, quantile normalization and raw microarray data were summarized by the RMA algorithm[7] in the microarray data linear model (Limma) package[8]. The false detection rate (FDR) adjusted P values using the Benjamini and Hochberg (BH) Procedure[9] in Limma package[8] were used for comparison. DEGs with P values less than 0.05 and |log2 fold change (FC)| values greater than 0 were selected as common thresholds[10,11]. The ggplot2 software program in the R software was used to create volcano plots of DEGs. The overlapping DEGs from both two datasets were selected for further study.

Figure 1

Two datasets were background adjusted and normalized by the log2 transformation, (A) GSE55457 before adjustment, (B) GSE55457 after adjustment; (C) GSE32317 before adjustment, and (D) GSE32317 after adjustment.

Analysis on GO enrichment and KEGG pathway

The KEGG pathway and GO enrichment of DEGs were analyzed with bioinformatics resources (Database for Annotation, Visualization and Integrated Discovery, DAVID), to explain the biological processes of DEGs between females and males[12,13].

Construction of protein-protein interaction (PPI) network and identification of hub genes

DEGs were imported into the interaction gene retrieval tool (STRING, https://www.string-db.org/) to create a PPI network, with a default threshold of a comprehensive score greater than 0.4, in order to understand the mutual functions of DEGs. To visualize the network, the PPI was then imported into the Cytoscape software (v.3.8.2). Hub genes (methods) were screened and visualized with the plug-in software cytohubba[14].

Results

Identification of DE-miRNAs

GSE55457 and GSE32317 datasets yielded 946 and 984 DEGs, respectively, by comparing the genes expressed in synovium of postmenopausal female OA patients and that of male patients. These included 694 up-regulated genes and 252 down-regulated genes from GSE55457, whereas GSE32317 had 648 up-regulated genes and 336 down-regulated genes. The volcano plots are shown in Figure 2, with up-regulated genes in red and down-regulated ones in blue. After crossover, there were 64 overlapping genes within the two datasets, with 10 genes responsible for up-regulation and 54 for down-regulation. As shown in Figure 3, the Venn diagram shows overlapping genes from both two data sets.

Figure 2

Volcano plots of differential genes between females and males with OA in (A) GSE32317 and (B) GSE55457.

Figure 3

GO enrichment analysis of overlapping DEGs with postmenopausal females and males in OA.

Enrichment analysis of target genes

GO enrichment analysis indicated that, DEGs were enriched in the cell adhesion, cell migration, ureteric bud development, BMP signaling pathway and translational initiation in BP. In CC, DEGs were enriched in both nucleus and proteinaceous extracellular matrix. And in MF, DEGs were shown to have higher level of DNA binding, chromatin DNA binding, heparin binding, complement receptor activity and proteoglycan binding (Figure 4, Table 1).

Figure 4

KEGG pathway with overlapping DEGs in postmenopausal females and males with OA.

Table 1

GO enrichment analysis of sex-related overlapping DEGs in OA.

Term	Biological function	Count	PValue
GO:0005634	Nucleus	23	0.066427
GO:0003677	Dna binding	10	0.061307
GO:0007155	Cell adhesion	5	0.056515
GO:0016477	Cell migration	4	0.017014
GO:0005578	Proteinaceous extracellular matrix	4	0.047148
GO:0001657	Ureteric bud development	3	0.00639
GO:0030509	Bmp signaling pathway	3	0.024005
GO:0006413	Translational initiation	3	0.069535
GO:0031490	Chromatin dna binding	3	0.01328
GO:0008201	Heparin binding	3	0.084231
GO:1901896	Positive regulation of calcium-transporting atpase activity	2	0.012567
GO:0051150	Regulation of smooth muscle cell differentiation	2	0.018792
GO:0072358	Cardiovascular system development	2	0.024978
GO:0002430	Complement receptor mediated signaling pathway	2	0.031126
GO:0001759	Organ induction	2	0.037237
GO:0060425	Lung morphogenesis	2	0.040277
GO:0002042	Cell migration involved in sprouting angiogenesis	2	0.046331
GO:0070584	Mitochondrion morphogenesis	2	0.058326
GO:1902895	Positive regulation of pri-mirna transcription from rna polymerase ii promoter	2	0.061302
GO:0048754	Branching morphogenesis of an epithelial tube	2	0.070173
GO:0051491	Positive regulation of filopodium assembly	2	0.078962
GO:0004875	Complement receptor activity	2	0.017993
GO:0043394	Proteoglycan binding	2	0.032745
GO:0051213	Dioxygenase activity	2	0.061599
GO:0032452	Histone demethylase activity	2	0.067268

According to the analysis of the KEGG pathway, MAPK and Thyroid hormone signaling pathways were the pathways associated with the sex differences in OA (Figure 5, Table 2).

Figure 5

Protein-protein interaction (PPI) network of top 10 hub genes.

PPI network

PPI network was created using the STRING database and then loaded into the Cytoscape program (v3.8.2) for visualization to assist the study of the relationships between the overlapping DEGs. Degree value was calculated for DEGs with the plug-in software Cytohubba and sorted, and those with higher degree values were hub genes, which were more likely to be associated with the disease. As shown in Figure 6, 10 hub genes with the largest degree values were BMP4, KDM6A, JMJD1C, NFATC1, PRKX, SRF, ZFX, LAMTOR5, UFD1L and AMBN. Among them, the levels of BMP4, KDM6A, JMJD1C, PRKX and ZFX were higher in postmenopausal female patients than those in males patients. LAMTOR5 was down-regulated in postmenopausal female patients compared with those in males with OA. NFATC1, AMBN and SRF were co-expressed genes for OA in postmenopausal females and males.

Discussion

OA is a degenerative joint disease that affects men and women differently. OA tends to be more common in female than male, and the incidence is higher in postmenopausal females. Sex hormones help to explain the sex difference in OA. After menopause, estrogen levels continue to decline, and the prevalence of OA increases sharply, indicating that estrogen may play a preventive role against OA. Prolonged estrogen usage, on the other hand, does not prevent or reverse tumor growth, but it does raise the risk of endometrial and breast cancer[15,16]. The fundamental mechanisms that result in different OA prevalence in males and females are still unclear. Therefore, it is critical to investigate the biochemical processes behind the gender difference in OA.

In this study, comprehensive bioinformatics analysis was conducted with two GEO data gene expression profile datasets, to identify sex-differential biological mechanisms associated with the OA pathogenesis. 64 overlapping genes were identified in postmenopausal females and males with OA, including 10 up-regulated and 54 down-regulated genes. Among them, BMP4, KDM6A, JMJD1C, NFATC1, PRKX, SRF, ZFX, LAMTOR5, UFD1L and AMBN are hub genes. GO enrichment analysis on overlapping genes revealed that hub genes were mainly enriched in cell adhesion, cell migration, ureteric bud development, BMP signaling pathway, translational initiation, nucleus, proteinaceous extracellular matrix DNA binding, chromatin DNA binding, and heparin binding, complement receptor activity and proteoglycan binding. According to these functional enrichment analyses, DEGs were enriched in chondrogenesis, mesenchymal stem cells differentiation into chondrocytes and angiogenesis in OA[17-19], contributing to the pathogenesis of the disease. Furthermore, KEGG pathways with overlapping DEGs showed that these hub genes were primarily found in the thyroid hormone signaling pathway. Previous studies on OA also suggested that both MAPK[20] and Thyroid hormone signaling pathways[21,22] contribute to OA pathogenesis.

In the constructed PPI network, 10 key genes were identified, including BMP4, KDM6A, JMJD1C, NFATC1, PRKX, SRF, ZFX, LAMTOR5, UFD1L and AMBN, which may be involved in the disease’ development. More evidence in previous studies proves that NFATC1, AMBN, SRF, BMP4, KDM6A and PRKX are also associated with OA pathogenesis. Angiogenesis is crucial in pathophysiology of inflammation of joints, such as OA. By promoting inflammatory cells invasion and the increase of local pain receptors, angiogenesis can result in pain as well as structural damage. Vascular endothelial growth factor (VEGF) is one of the mediators involved in angiogenesis[23]. It induces angiogenesis by activating the migration and proliferation of endothelial cells, which promotes macrophage recruitment and angiogenic response during inflammation. Angiogenesis and osteogenic coupling require optimal VEGF levels in regions where intramembrane ossification repair occurs. VEGF may operate as a paracrine factor during this process, since the loss of Vegfr2 in osteoblasts enhances osteoblast maturation and mineralization. Furthermore, during the endochondral osteogenic stage, the recruitment of osteoclasts and blood vessels are stimulated by mast chondrocyte-derived VEGF and osteoblast, promoting cartilage resorption at the repair site[24]. Serum response factor (SRF) is an essential transcription factor that plays a role in VEGF-induced angiogenesis in endothelial cells as a downstream mediator of VEGF signaling. PRKX is a vital protein kinase responsible for regulating angiogenesis, and it is involved in various pathological and physiological conditions associated with angiogenesis[25]. Ambn inhibits rankl-induced osteoclastogenesis by regulating the NFATC1 axis[26]. It has been found that bone morphogenetic proteins (BMPs) act as critical morphogenic factors because of their pleiotropic functions, which contribute to the regulation of tissue development, homeostasis as well as tissue repair[27]. BMP2 is a BMP4 homologous protein, which causes hypertrophy of chondrocytes and degradation of cartilage by activating LRP-5-induced Wnt/β-catenin signaling[28]. Cartilage specific ablation activates the T-nuclear factor C1 (Nfatc1), leading to early, invasive OA affecting multiple joints[29]. KDM6A, SOX9 and Aggrecan can significantly increase in ACCs co-cultured with bone marrow mesenchymal stem cells, collagen type 2 and BMSCs[30]. What role JMJD1C plays in the development of OA needs to be further explored, but JMJD1C is believed to help people understand how circulating androgens are regulated and provide potential therapeutic targets for androgen therapy[31], which may produce different pathogenesis in postmenopausal females with OA versus males with OA. In addition to the genes discussed, there is a lack of data on whether ZFX, LAMTOR5 and UFD1L are strongly linked to the onset and development of OA. Given the specific pathophysiological roles of NFATC1, AMBN, SRF, BMP4, KDM6A and PRKX in joint development, it is also valuable to investigate the targeted inhibition of receptors, including ZFX, LAMTOR5, UFD1L and JMJD1C in OA.

The KEGG pathway showed that the majority of hub genes were found to be abundant in the MAPK and the Thyroid hormone signaling pathways. The latter pathway mainly participates in the onset of OA by regulating microvessels in synovial, osteophytes and meniscus[32]. Scholars have also found decreased expression of autophagy markers Beclin-1 and lc3, decreased autophagosomes and p62 protein accumulation in chondrocytes stimulated with high TSH, indicating impaired autophagy flux. More interestingly, the mTOR was up-regulated with decreased AMPK activity in TSH-stimulated PMCs, indicating that mTOR/AMPK pathway is associated with TSH autophagy regulation in PMCs. It was also found that primary TSH-stimulated chondrocytes have increased apoptosis, autophagy is inhibited and involved in the OA progression[21]. The ability of the MAPK signaling pathway also includes the regulation of master transcription factors for the occurrence and function of chondrocytes, osteoblasts, and osteoclasts, demonstrating that this pathway is essential in terms of physiological bone development as well as homeostasis[20].

In conclusion, the findings of this study show that pathways such as the MAPK signaling pathway and the Thyroid hormone signaling pathway may be necessary in occurrence and development of postmenopausal OA. BMP4, KDM6A, JMJD1C, PRKX and ZFX may be key genes associated with the progression of OA in postmenopausal females and may be ideal targets for the treatment of OA in postmenopausal females. In addition, BMP4, KDM6A, JMJD1C, PRKX, ZFX, LAMTOR5 and UFD1L may be associated with therapeutic OA.

Ethics approval

The study was approved by the Ethics Committee of Xiangyang No.1 People’s Hospital, Hubei University of Medicine.

Authors’ contributions

JX and ZY conceived and designed the study, and drafted the manuscript. JX, ZY, GW, YZ, XL and FZ collected, analyzed and interpreted the experimental data. ZY and GW revised the manuscript for important intellectual content. All authors read and approved the final manuscript.