Integrative Analyses of Biomarkers and Pathways for Osteosarcopenia

Osteosarcopenia is a geriatric syndrome coexistence of osteoporosis and sarcopenia. However, the molecular mechanism underlying osteosarcopenia have not been fully elucidated. Differentially expressed genes (DEGs) for osteoporosis and sarcopenia were respectively identied by analyzing four expression datasets from the GEO. We extracted the gene expression datasets GSE56814 and GSE56815 for osteoporosis, GSE1428 and GSE8479 for sarcopenia. 133 co-expressed DEGs were included in osteoporosis and sarcopenia datasets. Furthermore, functional enrichment analyses and PPI network construction were performed to explore the potential biological function of the DEGs and identify hub genes. S100 protein binding (GO:0044548; p-value = 1.83E-06) and regulation of mRNA metabolic process (GO:1903311; p-value = 2.30E-05) were signicantly enriched in gene ontology(GO) analysis. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis showed that salmonella infection (hsa05132; p-value = 1.05E-04) and AMPK signaling pathway (hsa04152; p-value = 2.18E-03) were signicantly enriched. According to the results of PPI, we nally identied ve most critical genes as the hub genes, including AKT1, ANXA2, VIM, S100A6 and S100A11. The integrated analysis will contribute to the understanding of comprehensive molecular changes in osteosarcopenia and the development of new target therapies.


Introduction
Population aging is accelerating rapidly worldwide. Sarcopenia and osteoporosis are two conditions that are associated with aging, with similar risk factors that include genetics, endocrine function, and mechanical factors 1 . Osteoporosis is a metabolic bone disorder that is characterized by low bone mass and micro-architectural deterioration of bone tissue which can lead to an increased risk of fracture 2 . Sarcopenia means loss of muscle either mass or function in the elderly, and can reduce mobility, diminish quality of life, and lead to fall-related injuries 3,4 . Usually, osteoporosis and sarcopenia occur in consort, it is recently referred as 'osteosarcopenia' 1,5 .
Bone and muscle closely interact with each other, not only anatomically, but also chemically and metabolically 6,7 . Additionally, speci c pathophysiological nding, such as fat in ltration, is common to both diseases. It suggests that sarcopenia and osteoporosis are closely linked 8 .
Despite sharing common risk factors and biological pathways, the inherent molecular mechanisms underlying the pathogenesis of osteosarcopenia have not been fully clari ed. Further researches are needed to provide deeper understanding to explore more advantageous therapeutic targets.
In the past few years, microarray technology had been widely used for gene expression pro ling in osteoporosis and sarcopenia patients 9,10 . However, there was no microarray study in osteosarcopenia patients. Thus, we conducted a gene expression meta-analysis between osteoporosis and sarcopenia using integrated bioinformatics methods. Additionally, based on the result of this analysis, gene enrichment and pathway annotation analysis were performed and hub genes were also identi ed in the DEGs.

Identi cation of DEGs
After standardization of the datasets (Supplementary Figure 1

PPI network analysis and hub gene selection
Based on the STRING database, the PPI network of DEGs were gathered as a cluster consisting of 133 nodes and 116 edges ( Figure 5). MCODE plug-in of Cytoscape software was applied to identify the most signi cant module that was comprised of 5 nodes (MCODE score ≥ 2). We nally de ned the ve hub genes, including AKT1, ANXA2, VIM, S100A6 and S100A11.

Discussion
Osteosarcopenia is a recent terminology, which means a person under osteoporosis and sarcopenia at the same time. Some professors also called it as sarco-osteoporosis 11 . Previous studies demonstrated that sarcopenia patients had a higher risk of having coexisting osteoporosis compared with nonsarcopenia individuals 12 . Some other studies showed that one with both osteoporosis and sarcopenia has a higher risk of fall and fracture than those with osteoporosis or sarcopenia alone 11,13 . Osteosarcopenia is recognized as a complex disease, not only associated with age, but also with endocrine disorders, malnutrition, obesity and use of corticosteroids 14 . Identifying the susceptibility gene of osteosarcopenia could help us to deeper understand this disease and to nd potential treatments.
Our study was the rst to systematically search and incorporate the microarrays on osteosarcopenia in GEO. In the present study, we included two microarrays in osteoporosis and sarcopenia respectively, compared gene expression pro les between old and young person, identi ed the co-expressed DEGs in the two diseases. Furthermore, functional and pathway enrichment analyses were performed to understand the potential biological function of the DEGs, and PPI network construction were performed to identify interactions between DEGs.
133 DEGs were ltered out across multiple datasets. Among these DEGs, DUSP1 and SARAF got signi cant difference in all the four datasets. AKT1, ANXA2, VIM, S100A6 and S100A11 were identi ed as the hub genes under PPI network analysis. The AKT1 gene could be transcribed to produce a protein called AKT1 kinase, which plays a critical role in many signaling pathways, helps regulate cell growth and division, and also helps control apoptosis. Several researches have reported that AKT1 could control osteoblast differentiation and development 15 , and also could in uence metabolic and differentiation of muscle 16 . It was also be validated in our study, that AKT1 was high expressed in osteoporosis patients while low expressed in sarcopenia patients.
ANXA2 encodes a member of a calcium-dependent phospholipid-binding protein family, which play a role in the regulation of cellular growth and signal transduction pathways. Previous studies have shown that elevated ANXA2 protein expression level could stimulate peripheral blood monocytes (PBM) differentiate into higher number of osteoclasts and resorb bone at higher rates to decreasing BMD 17 . In our study, ANXA2 was low expression in osteoporosis datasets. In sarcopenia datasets, the ANXA2 was low expression in GSE1428, while high expression in GSE8479. It is consistent with previous studies that ANXA2 is helpful to repair muscle cell membrane and also could mediate the occurrence of muscle in ammation [18][19][20] .
VIM gene could encode a type III intermediate lament protein, which maintaining the shape of cell and integrity of the cytoplasm, and stabilizing interaction of cytoskeleton. Previous study has report that VIM has a key role in regulating in ammation in macrophages 21 . And increased expression of VIM could be considered as a reliable marker of muscle regeneration 22,23 . The expression of VIM in sarcopenia datasets was difference in our study, that maybe the result of different micro-environment. VIM gene was up-regulated in ostoeporosis patients compared with contral group. The regulatory mechanism of VIM in osteoporosis still need to be further studied. S100A6 and S100A11 are members of the S100 family. Several studies demonstrated S100 proteins may participate in functional inhibition of osteobalst progenitor cells, in ammatory, and combined with ANXA1 to participate in membrane repair 18,24−26 . In our study, S100A6 and S100A11 were up-regulared in osteoporosis patients, but different expressed in sarcopenia datasets. This result showed that S100 protein has various regulatory functions in muscle tissue under different situation.
Further functional enrichment analyses demonstrated that DEGs in osteosarcopenia were signi cantly enriched in regulation of mRNA metabolic process and stability. Previous studies showed that some hormone, such as growth hormone/insulin-like growth factor-1 (GH/IGF1), play a key role in the development of osteosarcopenia 27 . Besides the AMPK signaling pathway, a pathway that affects the metabolism of various substances, the DEGs were signi cantly enriched in the salmonella infection with KEGG analyses. Several studies have con rmed that some pro-in ammatory cytokines could promote bone resorption, and in uence the development of sarcopenia through the activation of ubiquitinprotease pathway 28-30 . However, there was a limitation in our study that the four datasets were included from two diseases. Until now, there have been no genetic sequencing studies in osteosarcopenia patients. Further researches are needed to verify the changes of gene expression in osteosarcopenia patients in the future.

Conclusions
In summary, the study provided deeper insight to the comprehensive molecular changes in the pathogenesis of osteosarcopenia, and identi ed several potential therapeutic targets, including AKT1, ANXA2, VIM, S100A6 and S100A11. Through GO and KEGG pathway analysis, we found that the DEGs were mainly enriched in regulation of mRNA metabolic process, regulation of mRNA stability and S100 protein binding, and were involved in salmonella infection and AMPK signaling pathway. The results could help to anticipate adverse health outcomes and development of new target therapies.

Microarray Data
We selected GEO (http://www.ncbi.nlm.nih.gov/geo), which is a publicly available database of gene/microarray pro les for our study 31 . The search strategy ("osteoporosis, postmenopausal" [MeSH Terms] OR "osteoporosis" [MeSH Terms] OR osteoporosis [All Fields]) AND "Homo sapiens" [porgn] AND "gse" [Filter] was adopted for osteoporosis patients. And the search strategy (sarcopenia [All Fields] AND "Homo sapiens" [porgn] AND "gse" [Filter]) was adopted for sarcopenia patients. All the included datasets contain normal tissue as controls. We extracted the gene expression datasets GSE56814 and GSE56815 for osteoporosis, GSE1428 and GSE8479 for sarcopenia ( Table 1). The tissues of osteoporosis datasets are blood monocytes, which are progenitors of osteoclasts and produce important factors to bone metabolism. The tissues of sarcopenia datasets are vastus lateralis muscle, which can directly demonstrate the gene expression of skeletal muscles. with high BMD. Additionally, the platform for GSE1428 is also GPL96, which include older (n = 12) and younger (n = 10) male. And the platform for GSE8479 is GPL2700, Sentrix HumanRef-8 Expression BeadChip, which included older (n = 25) and younger (n = 26) adults. Series matrix les and related annotation document of GSE56814, GSE56815, GSE1428 and GSE8479 were downloaded.

Identi cation of differentially expressed genes
The corresponding annotation document was used to map the microarray probes to gene symbols. The mean value was adopted when multiple probes mapped to the same symbol. All the expression microarray datasets were rst standardized. And the DEGs were determined between osteoporosis or sarcopenia tissues and normal tissues in each microarray by using the limma V3.46.0 (linear models for microarray data) package of the R software program (version 4.0.3) 32 . The |log2 fold change (FC)| > 0.1 and p-value < 0.05 were regarded as the cut-off criteria to determine DEGs. Online tool Calculate and draw custom Venn diagram (http://bioinformatics. psb.ugent.be/webtools/Venn/) was applied to detect common DEGs among the 4 datasets. DEGs in any three of the four datasets were considered to the common DEGs.

Functional and Pathway Enrichment Analyses
Gene ontology (GO) analysis is used extensively to identify the characteristic biological attributes of genes, gene products, and sequences, including biological process (BP), cell components (CC) and molecular function (MF) 33 . Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis provides a comprehensive set of bio-interpretation of genomic sequences and protein interaction network information 34 .
In this study, GO analysis and KEGG pathway enrichment analysis of DEGs were automatically completed and visualized by the clusterPro ler V3.18.0 35 , org.Hs.eg.db V3.10.0, and Goplot V1.0.2 package 36 in the R software statistical analysis platform. P.adjust < 0.1 and q-value < 0.2 were regarded as the cut-off criteria in GO analysis, while p-value < 0.05 was regarded as the cut-off criteria in KEGG analysis.

Protein-Protein Interaction (PPI) Network Analysis
We searched the STRING database (http://www.string-db.org/) 37 to identify and predict interactions between DEGs, and construct the PPI network with the cut-off standard as a con dence > 0. 4

Con icts of Interest
No potential con icts of interest were reported by any of the authors in the conduct of this work.

Funding
This work was supported by the clinical medicine science and technology innovation plan of Jinan, Grant No. 202019202. Figure 1 Volcano plots of the four microarray datasets. Red points represented up-regulated genes, while blue points represented down-regulated genes. Gray points represented genes with no signi cant difference.

Figure 2
Venn diagram of common differentially expressed genes from the four datasets.   Supplementary Files