The Biological Interaction of SARS-CoV-2 Infection and Osteoporosis: A Preliminary Study

The COVID-19 pandemic caused by the severe acute coronavirus disease 2 (SARS-CoV-2) virus represents an ongoing threat to human health and well-being. Notably, many COVID-19 patients suffer from complications consistent with osteoporosis (OP) following disease resolution yet the mechanistic links between SARS-CoV-2 infection and OP remain to be clarified. The present study was thus developed to explore the potential basis for this link by employing transcriptomic analyses to identify signaling pathways and biomarkers associated with OP and SARS-CoV-2. Specifically, a previously published RNA-sequencing dataset (GSE152418) from Gene Expression Omnibus (GEO) was used to identify the differentially expressed genes (DEGs) in OP patients and individuals infected with SARS-CoV-2 as a means of exploring the underlying molecular mechanisms linking these two conditions. In total, 2,885 DEGs were identified by analyzing the COVID-19 patient dataset, with shared DEGs then being identified by comparison of these DEGs with those derived from an OP patient dataset. Hub genes were identified through a series of bioinformatics approaches and protein-protein interaction analyses. Predictive analyses of transcription factor/gene interactions, protein/drug interactions, and DEG/miRNA networks associated with these DEGs were also conducted. Together, these data highlight promising candidate drugs with the potential to treat both COVID-19 and OP.


INTRODUCTION
The COVID-19 (coronavirus disease 2019) pandemic, caused by the single-stranded RNA severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) virus, has caused extensive death and suffering throughout the world Shapira et al., 2022;Swets et al., 2022). Notably, many older adults that recover from SARS-CoV-2 infection have been reported to develop new-onset or aggravated osteoporosis (OP). As a systemic bone disease, OP is associated with reductions in bone density, increased bone fragility, and the degradation of the bone microstructure . Several potential mechanisms may explain the observed link between COVID-19 and OP development. For one, a range of antiviral drugs are employed in the treatment of COVID-19 patients, including corticosteroids, which can contribute to OP onset or aggravation (Xiong et al., 2020). Infected patients also exhibit abnormal increases in miR-4485 expression in their bone marrow, and this microRNA (miRNA) has been reported to target TLR for and to interfere with appropriate osteogenic remodeling (Mi et al., 2021). COVID-19 can also interfere with patient exercise and medication use, both of which can contribute to higher OP rates in elderly individuals (Salvio et al., 2022). There is thus a clear need to better clarify the mechanisms linking SARS-CoV-2 infection and OP development.
The angiotensin-converting enzyme 2 (ACE2) receptor has been confirmed to play a key role as a receptor that SARS-CoV-2 utilizes to enter into and infect target cells (Casciola-Rosen et al., 2022;Popescu et al., 2022). Bone tissue, however, exhibits only low levels of ACE2 expression, and whether COVID-19-related bone damage is the result of direct infection or is secondary to pathological changes in other organs thus remains to be conclusively established. While certain bone marrow cells do exhibit ACE2 receptor expression, a biological link between SARS-CoV-2 and bone remodeling has yet to be demonstrated.
The present study was developed to explore the potential mechanisms whereby SARS-CoV-2 infection may influence OP development or severity. To that end, extant RNA-sequencing (RNA-seq) data were utilized to conduct transcriptomic analyses of signaling pathways and biomarkers associated with SARS-CoV-2 infection and OP. Differentially expressed genes (DEGs) expressed in OP patients infected with SARS-CoV-2 were identified using the GSE152418 dataset from the GEO database, and the resultant gene list was used to explore relevant signaling pathways as well as candidate drugs with the potential to aid in COVID-19 treatment. Hub regulatory genes in this pathological context were further identified through a series of bioinformatics analyses and through the construction of a protein-protein interaction (PPI) network. Moreover, efforts to clarify the biological link between SARS-CoV-2 infection and OP incidence were made by using the identified DEGs to conduct predictive transcription factor/ gene interaction, protein/drug interaction, and DEG/miRNA network analyses with the goal of elucidating the underlying molecular mechanism.

Data Search and DEG Identification
COVID-19) ( Figure 1A). The GSE100609 dataset consisting of 4 OP patients and 4 healthy controls was subsequently retrieved from the GEO database and analyzed, leading to the identification of 494 DEGs, of which 296 and 198 were respectively upregulated and downregulated in OP patients ( Figure 1B).

Identification of DEGs Shared Between COVID-19 and OP
The JVenn program was next used to identify 73 DEGs shared between these COVID-19 and OP datasets, suggesting that they may be associated with the biological link between these two conditions ( Figure 2). Given this overlap, it is possible that drugs capable of treating OP may offer value in the treatment of COVID-19.

Gene Enrichment Analyses of Shared COVID-19-and OP-Related DEGs
Next, GO and KEGG enrichment analyses of these shared DEGs were conducted using the Enrichr application. These DEGs were significantly enriched in biological process GO terms including regulated exocytosis, platelet degranulation, cytokine-mediated signaling pathway, regulation of macrophage-derived foam cell differentiation, and positive regulation of tau-protein kinase activity. They were also enriched for cellular component GO terms including platelet alpha granule, platelet alpha granule lumen, muscle myosin complex, platelet alpha granule membrane, high-density, and lipoprotein particle. With respect to molecular function GO terms, these shared DEGs were enriched in the calcium-dependent phospholipid binding, CXCR chemokine receptor binding, glutamate receptor activity, ligand-gated channel activity, and ligand-gated ion channel activity terms ( Figure 3; Table 1). Common DEGenriched pathways were further explored using the KEGG, Reactome, WikiPathway, and BioCarta databases, with the top 10 enriched terms being compiled in Figure 4 and Table 2.

PPI Network Construction and Hub Gene Identification
Following the importation of the 73 shared DEGs into the STRING database, a PPI network consisting of 73 nodes and 235 edges was generated ( Figure 5A), with the Cytoscape 3.7.2 tool being used to reorganize this network ( Figure 5B) The Cytohubba plugin was then used to identify the 10 genes with the highest degree value in this network, with these genes (CXCL5, MMP1, PF4, VWF, IL11, ITGB3, HP, IL1R2, FASLG, MMP10) being identified as hub genes ( Figure 5C).

Identification of Transcription Factors and miRNAs Associated With Shared DEGs
To more fully explore the potential regulation of the identified hub DEGs at the transcriptional level, the NetworkAnalyst tool was used to identify shared transcription factors and miRNAs associated with these DEGs that may regulate their expression ( Figure 6).

Exploration of Potential Therapeutic Drugs and Gene-Disease Associations
Lastly, the DSigDB database was used to identify candidate drugs associated with these shared DEGs, working under the assumption that these genes may represent promising therapeutic targets associated with both COVID-19 and OP. The top 10 retrieved drug compounds were ARSENIC CTD 00005442, MS-275 PC3 UP, camptothecin PC3 UP, fluoride CTD 00005982, 1,4-chrysenequinone PC3 UP, azacitidine PC3 UP, sanguinarine HL60 UP, benzo [a]pyrene CTD 00005488, ellipticine PC3 UP, and CP-690334-01 PC3 UP ( Figure 7; Table 3). Gene-disease association analyses conducted with Network-Analyst further revealed the Hypersensitivity, Prostatic Neoplasms, Atherosclerosis, Diabetic Angiopathies, Anemia, Myocardial Infarction, Schizophrenia, Hemorrhage, Unipolar Depression, and Major Depressive Disorder disease states to be most associated with the identified COVID-19/ OP-related hub genes ( Figure 8).

DISCUSSION
OP is a systemic bone disease that results in reductions in bone density, increased bone fragility, and the degradation of the bone microstructure. Recent research suggests that there may be an important link between COVID-19 and the onset or aggravation of OP among older adults. The present study was thus conducted to identify potential molecular biomarkers that may be shared between COVID-19 and OP in an effort to clarify the link between these two conditions and to guide future therapeutic efforts. High-throughput gene expression profiling is commonly used to identify molecular biomarkers associated with a range of disease states (Mahmud et al., 2021). In the present study, 73 DEGs shared between SARS-CoV-2-infected and OP patients were identified through transcriptomic analyses, and these genes were then subject to further enrichment analyses aimed at better exploring the link between OP and COVID-19. GO analyses enable the assessment of the regulatory relationships for particular genes based on theoretical models of associated genes and internal relationships (Nashiry et al., 2021), with these predictive enrichment analyses being based on a growing base of knowledge regarding gene function and associated ontological classes (Wang et al., 2022). GO analyses broadly classify gene characteristics based on associated biological process (BP), cellular component (CC), and molecular function (MF) terms (Cao et al., 2022). The shared DEGs identified in this study were enriched in BP terms including regulated exocytosis (8 genes) and the platelet degranulation signaling pathway (6 genes) are among the top GO terms. Exocytosis is closely tied to the survival and function of osteoblasts (Oppert et al., 2022), and the regulation of osteoblast exocytosis by vitamin D3 is reportedly significant in the context of bone remodeling . The shared DEGs were also highly enriched in the platelet alpha granule (6 genes) and platelet alpha granule lumen (4 genes) CC terms. Platelet alpha granules impact endothelial exocytosis and OP development, and platelets play a well-documented role in regulating bone remodeling and vascular function (Li et al., 2021;Salamanna et al., 2021). The top MF terms in which these shared DEGs were enriched included calcium-dependent phospholipid binding (3 genes) and CXCR chemokine receptor binding (2 genes). Notably, calcium-dependent phospholipid binding plays a key role in viral replication (Schleiss et al., 2021), while the CXCR chemokine receptor family can recruit stem cells to facilitate subsequent osteogenic differentiation (Khokhar et al., 2022).
KEGG enrichment analyses were additionally used to explore shared pathways associated with these 73 common DEGs evident in the analyzed COVID-19 and OP datasets. The top 10 KEGG pathways associated with these DEGs included the hematopoietic cell lineage, cytokine-cytokine receptor interaction, platelet activation, complement and coagulation cascades, ECMreceptor interaction, rheumatoid arthritis, human papillomavirus infection, tyrosine metabolism, long-term depression, and phenylalanine, tyrosine and tryptophan biosynthesis pathways. A PPI network was further constructed based on the identified DEGs, with the hub proteins in this network (CXCL5, MMP1, PF4, VWF, IL11, ITGB3, HP, IL1R2, FASLG, MMP10) being regarded as the most critical shared regulators of OP and SARS-CoV-2 infection. CXCL5 signaling is associated with cellular proliferation and differentiation (Kawagoe et al., 2020), and the aberrant expression of this chemokine can lead to a range of dysfunctional outcomes, potentially contributing to the incidence of infection and/or OP development (Zhang et al., 2020).
Transcription factors and miRNAs that may function as upstream regulators of these DEGs were additionally identified in an effort to gain more insight into the pathological basis of these disease states. Many drugs have   (Sing et al., 2020). Many older adults suffering from COVID-19 also exhibit comorbid OP. The identification of pathogenic factors shared by these two diseases may enable the simultaneous administration of antiviral and anti-OP agents with the potential to provide immense clinical benefit to this patient population.
There are some limitations of this study. Firstly, we didn't perform any in vivo or in vitro experiments to validate our bioinformatics result, because it is a preliminary study. Validations will be carried out in our further research. Secondly, the number of cases in GSE100609 is relatively small. In addition, for clinical traslation of our findings, it is really a long way from the bench side to clinical use. Osteoporosis is one of the most common disorders in the elderly, who are also relatively vulnerable to COVID-19. The drugs identified in this

Data Search and DEG Identification
The GEO (Gene Expression Omnibus, https://www.ncbi.nlm.nih.gov/ geo/) database repository, which compiles MIAME-compliant data submissions, was searched to identify RNA-seq data corresponding to COVID-19 and OP patients. Genes differentially expressed between COVID-19 or OP patients and healthy controls were identified with the R limma package and DESeq2 using the following criteria: p < 0.05 and |logFC| > 1.0. A heatmap of the identified DEGs was constructed with the R pheatmap package applied to detect significant DEGs. The heatmap was drawn with the pheatmap package in R.

Identification of Shared DEGs Between COVID-19 and OP Patients
All DEGs identified in the initial comparison of COVID-19 and OP patient datasets were imported into the online JVenn  Frontiers in Cell and Developmental Biology | www.frontiersin.org May 2022 | Volume 10 | Article 917907 7 tool (http://jvenn.toulouse.inra.fr/app/example.html), which was used to identify mutual DEGs shared between these two datasets. The shared DEGs were then represented in a pie chart generated by JVenn.

Functional Enrichment Analyses of Shared DEGs
GO (Gene ontology) and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway analyses enable researchers to gauge the potential functional roles of particular genes of interest. The Enricher tool is a web-based platform that enables these functional enrichment analyses following the importation of gene lists of interest. For the present study, Enrichr analyses of mutual DEGs were conducted using four databases (KEGG,

Protein-Protein Interaction Network Analysis
The STRING (https://string-db.org/) database was used to construct a PPI network incorporating the identified mutual DEGs, with interactions exhibiting a combined score >0.5 being incorporated into the generated network. The open-source Cytoscape (v 3.7.2) platform was used to construct and visualize the resultant network.

Hub Gene Identification
The Cytohubba plugin for Cytoscape was used to extract network features using an MMC (Maximal Clique Centrality) approach, with the top 10 genes identified within the PPI network being selected as hub genes for further analysis.

Retrieval of Transcription Factors and miRNAs That Interact With Mutual DEGs
Transcription factors regulate the expression of specific target genes, and identifying shared transcriptional regulators of particular genes of interest can thus offer valuable molecular insight. The web-based NetworkAnalyst (http://www. networkanalyst.ca) tool was used with the JASPAR database to identify transcription factors with the potential to bind to the shared DEGs associated with OP and COVID-19 in the present study. In addition, miRNAs with the potential to regulate these DEGs were identified through the use of the Tarbase and mirTarbase databases. Cytoscape 3.7.2 was then used to visualize the identified transcription factor/gene and miRNA/ gene interaction networks produced through these analyses.

Exploration of Potential Therapeutic Drugs and Gene-Disease Associations
The gene set enrichment analysis-based DSigDB database was used with EnrichR to conduct a protein/drug interaction analysis with the goal of identifying pharmacological compounds with the potential to regulate these target DEGs. In total, DSigDB incorporates 22,527 gene sets and 17,389 unique compounds associated with 19,531 genes. In addition, DEG/disease relationships were assessed with DisGeNTET via NetworkAnalyst.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. This data can be found here: GEO.

AUTHOR CONTRIBUTIONS
HZ conceived and designed the study; XK and XW performed bioinformatics analysis; JL and LL analyzed the data; YZ and QW provided advice and technical assistance; XK and XW wrote the manuscript. All authors approved the final manuscript.