{'fraction_non_alphanumeric': 0.02787433155080214, 'fraction_numerical': 0.017179144385026737, 'mean_word_length': 5.257214554579674, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Identification of key genes and pathways associated with sex difference in osteoarthritis based on bioinformatics analysis

Abstract 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][2][3][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 sexspecific mechanisms of OA, as well as prospective biomarkers and therapeutic targets to treat OA.

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, backgroundadjustment 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.

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 .

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 downregulated 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 upregulation and 54 for down-regulation. As shown in Figure   3, the Venn diagram shows overlapping genes from both two data sets.

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).
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).

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][18][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 VEGFinduced 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.