Stochastic investigation of HIV infection and the emergence of drug resistance

1.   Introduction

Osteoporosis (OP) is the most common bone disease characterized by a loss of bone mass and quality that results in fragility fractures. Osteoporosis is responsible for over 8.9 million fractures each year worldwide, with most cases occurring in postmenopausal women ^[1]. The most frequent osteoporotic fractures are fractures of the hip, wrist and spine, although most fractures in the elderly are probably at least partly related to osteoporosis. To prevent the osteoporosis, clinical risk factors ^{[2,3,4,5,6,7]}, including age, female gender, rheumatoid arthritis, previous fragility fracture, genetic predisposition, smoking, > 14 units of alcohol per week, early-onset menopause (<45 years of age) and low BMI (<19 kg per m²), should be evaluated. To diagnose the osteoporosis, bone mineral density (BMD) measurement is widely used in clinical medicine as an indirect indicator of osteoporosis and fracture risk ^[8]. Specifically, the BMD measurement can category the subjects into four statuses, including normal, osteopenia, osteoporosis and severe osteoporosis.

Bone remodeling is accomplished by two specialized cells: bone-resorbing osteoclasts and bone-forming osteoblasts ^[9]. Osteoclasts can break down bone, while osteoblasts create new bone. Osteoclasts and osteoblasts can coordinate well for most of the human life ^[10]. In a pathological state of osteoporosis, this coordination can break down, and the osteoclasts begin to remove more bone than the osteoblasts can create ^[11]. To understand the molecular mechanism of osteoporosis, several studies introduced molecules, including SNP ^{[12,13,14,15]}, DNA methylation ^[16], RNA ^[17] and microRNA ^[18], to investigate the genetic, epigenetic and transcriptomic regulators involved in osteoclastogenesis or osteoblastogenesis.

Osteoclastogenesis has been reported to cause osteoporosis by previous studies ^[19,20,21], however, the critical functional modules and signaling pathways still need to be explicitly characterized, which is useful for the discovery of the therapeutic targets in osteoclasts for the treatment of skeletal diseases. In this study, we collected gene expression datasets of peripheral blood monocytes (PBMs) from low or high BMD subjects. Notably, PBMs, which are osteoclast progenitor cells and produce cytokines for osteoclastogenesis and bone resorption, could be used for osteoporosis study. Furthermore, network based analysis and overrepresentation enrichment analysis highlighted critical functional modules and signaling pathways, which improved our understanding of the molecular mechanism, and provided the potential drug targets for osteoporosis.

2.   Materials and method

2.1.   Gene expression data pre-processing

Gene expression datasets were obtained from the NCBI Gene Expression Omnibus (GEO) (http://www.ncbi.nlm.nih.gov/geo) with accession numbers, GSE56814 ^[22]. Prior to downstream analysis, we firstly mapped the array probes to the GENCODE gene annotation v19. We calculated the average expression values of genes matching multiple probes.

2.2.   Differential expression analysis

The differential expression analysis was conducted in R programming language (https://www.r-project.org) with student-t test. The differentially expressed genes were identified at the threshold P-value < 0.05. The up-or down-regulation status was determined based on the fold change for each gene.

2.3.   Overrepresentation enrichment analysis (OEA)

Overrepresentation enrichment analysis was also implemented at WEB-based Gene Set Analysis Toolkit (WebGestalt) ^[23]. The Gene Ontology biological processes and KEGG pathways were selected as the functional database.

2.4.   Protein-protein interaction network mapping

The Search Tool for the Retrieval of Interacting Genes/ Proteins (STRING) ^[24] online software (https://string-db.org) was used to assess the interactions. The interactions of the proteins encoded by the differently expressed genes were searched using STRING online software, and the combined score of > 0.4 was used as the cut-off criterion. The PPI network was visualized using Cytoscape software (http://www.cytoscape.org). The Cytoscape MCODE plug-in (version 3.4.0) was applied to search for clustered sub-networks of highly connected nodes from the PPI network. The resulting network was subjected to module analyses with the Plugin MCODE with the following default parameters: Degree cut-off ≥3.

2.5.   Construction of co-expression network between differentially expressed mRNAs and lncRNAs

The Pearson correlation coefficient of DEG-lncRNA pairs was calculated according to their expression values. The co-expressed DEG-lncRNA pairs with an absolute value of the Pearson correlation coefficient of ≥0.8 were selected and the co-expression network was visualized by using Cytoscape software.

3.   Results

3.1.   Identification of long non-coding RNAs by microarray probe reannotation

To identify the probes on the array corresponding to long non-coding RNA (lncRNA) sequences, we re-annotated all probes using GENCODE gene annotation with version 19. Finally, we identified a total of 830 long non-coding RNAs, with gene types such as processed transcript, antisense, long intergenic noncoding RNA (lincRNA) and pseudogene. We found that pseudogene accounted for the most of the identified lncRNAs (45%), followed by the long intergenic noncoding RNAs (29%, Figure 1A). In addition, the reannotation of the probes also mapped 17,051 probes to 17,743 protein-coding genes. We then compared the expression levels of lncRNAs with those of protein-coding genes, and found that expression levels were lower in lncRNAs than protein-coding genes (Figure 1B).

3.2.   Differentially expressed mRNAs and lncRNAs

To uncover genes that may be responsible for osteoporosis, we compared the expression profiles of low BMD subjects with that of high BMD subjects. We identified a total of 496 differentially expressed genes (DE-genes), including 120 up-regulated and 376 down-regulated genes in low BMD subjects (t-test, P-value < 0.05, Supplementary Table 1). We observed that four zinc finger genes, ZNF528, ZNF79, ZNF223 and ZNF765, were included in the top-ten up-regulated genes, indicating that zinc finger genes may play key roles in osteoporosis. Notably, 24 lncRNAs were differentially expressed (DE) between samples with low and high BMD, which were presented in Supplementary Table 1. Visualization of the expression levels of DE-genes revealed that a fraction of normal controls (n = 15) exhibited similar expression patterns with the cases with low BMD (Figure 2). Furthermore, a higher proportion of postmenopausal normal controls were observed in the 15 samples than other normal controls (proportion test, P < 0.05, 66.7% vs. 22.2%), suggesting that postmenopausal normal controls had a higher risk of osteoporosis.

3.3.   GO and KEGG enrichment analysis

Next, the differentially expressed mRNAs were subjected to GO and KEGG analyses (Figure 2). GO analysis indicated that the upregulated genes were mainly involved in regulating nucleobase metabolic process, regulation of organelle assembly, regulation of plasma lipoprotein particle levels, response to interferon-gamma, snRNA metabolic process, ncRNA transcription and ncRNA processing (Figure 3A). Furthermore, the downregulated genes were mainly enriched in categories associated with regulation of sequence-specific DNA binding transcription factor activity, neutral lipid metabolic process, cytokinesis, nucleoside triphosphate metabolic process, regulation of cell-cell adhesion, nucleoside monophosphate metabolic process, and glycoprotein metabolic process (Figure 3B). The above pathways may therefore participate in regulating the bone miner density.

KEGG pathway analysis revealed that upregulated genes were primarily enriched in pathways associated with the glycine, serine and threonine metabolism, ABC transporters, arginine and proline metabolism, Jak-STAT signaling pathway, RNA degradation, and herpes simplex infection (Figure 3C). In accordance with the enriched GO terms, the up-regulated genes were significantly enriched in inflammation-related pathways. Downregulated genes were mainly associated with adipocytokine signaling pathway, metabolic pathways, oxidative phosphorylation, cytosolic DNA-sensing pathway, and glycerophospholipid metabolism (Figure 3D). The GO and KEGG enrichment analysis of the down-regulated genes suggested that metabolic pathways above may function by maintaining the bone miner density.

3.4.   PPI network construction for DE-mRNAs

To investigate the interactions between the proteins encoded by the DEGs, we mapped the DEGs to protein-protein interaction network (Figure 4). The PPI network for the upregulated genes contained 38 nodes and 33 edges, and the hub nodes with the higher connectivity degree were SHMT1, MYLIP, MKRN1 and EXOSC8 (connectivity degree = 4, Figure 4A). The PPI network for the downregulated genes contained 198 nodes and 348 edges, and the hub nodes with the higher connectivity degree were IL8 (connectivity degree = 28), NFKBIA (connectivity degree = 13), HADC5 (connectivity degree = 12) and HDAC4 (connectivity degree = 12) (Figure 4B).

Furthermore, we also performed a module analysis of the network using MCODE plugin. With degree cut-off ≥3, four up-regulated proteins, including MEX3C, MKRN1, TRIM4 and MYLIP, were identified as key proteins/genes of a module in the up-regulated PPI network (Figure 5A). The high connectivity and components of a module suggested that MKRN1 and MYLIP may play important roles in regulating bone miner density. Similarly, the down-regulated proteins, including JDP2, HADC4, HDAC5, CDYL2, ACADVL, ACSL1 and BRD4, were identified as key proteins/genes in the down-regulated PPI network (Figure 5B). Moreover, IL8, CCL25, CNR1, HCAR2, HCAR1, GNG13 and GNG8 were also identified as a module in the down-regulated PPI network (Figure 5C). All these proteins except IL8 participated in G-alpha (ⅰ) signaling events, indicating that this pathway may be important for bone mineral homeostasis.

3.5.   Co-expression network analysis of differently expressed mRNAs and lncRNAs in osteoporosis

To predict the potential functional roles of the differentially expressed lncRNAs, the Pearson correlation coefficient for lncRNA-DE-mRNA pairs was first calculated according to their expression value. The co-expressed mRNA-lncRNA pairs with an absolute value of their Pearson correlation coefficient of ≥0.8 were selected. As presented in Figure 6A, the network included 7 differentially expressed lncRNAs and 70 DE-mRNAs.

To characterize the biological function of the 7 differentially expressed lncRNAs, we performed gene set enrichment analysis on the highly correlated DE-mRNAs. Based to the GO terms by enrichment analysis, we successfully annotated 6 of the 7 DE-lncRNAs (P < 0.05, Figure 6B). Specifically, AC009041.1 and SSTR5-AS1 were primarily enriched in DNA conformation change and DNA recombination. Particularly, RP11-498C9.17 was enriched in pathways involved in epigenetic regulation, such as macromolecule deacylation, epigenetic regulation of gene expression, chromatin remodeling, chromatin assembly or disassembly and protein-DNA complex subunit organization. Other lincRNAs may participate in pathways such as cell-cell signaling by Wnt, response to extracellular stimulus, autophagy, and cell cycle checkpoint.

4.   Discussion

Osteoclastogenesis and bone resorption play a fundamental role in osteoporosis pathogenesis. Better understanding the regulation of osteoclastogenesis is very important for the discovery of therapeutic targets for the treatment of osteoporosis.

In this study, we aimed to analyze the key mRNAs and lncRNAs in osteoporosis. To explore the molecular mechanism of regulating bone miner density, we collected gene expression datasets of peripheral blood monocytes (PBMs) from low or high BMD subjects. To identify the probes on the array corresponding to long non-coding RNA (lncRNA) sequences, we re-annotated all probes using GENCODE gene annotation, and identified a total of 830 long non-coding RNAs.

To investigate the biological differences between samples with low and high BMD, we performed differential expression analysis of the gene expression data from the discovery datasets, and identified a total of 496 differentially expressed genes (DEGs), including 120 up-regulated and 376 down-regulated in low BMD subjects, of which, 24 were long non-coding RNAs differentially expressed between samples with low and high BMD. Subsequently, the differentially expressed mRNAs were subjected to GO and KEGG analyses. The up-regulated mRNAs were significantly enriched in inflammation-related pathways, which was consistent with the conclusion that inflammatory condition could contribute to the differentiation from monocyte to osteoclast by previous study ^[25]. In contrast, the downregulated genes were mainly associated with adipocytokine signaling pathway, metabolic pathways, oxidative phosphorylation, cytosolic DNA-sensing pathway, and glycerophospholipid metabolism, suggesting that metabolic pathways may function by maintaining the bone miner density.

To investigate the interactions between the proteins encoded by the DEGs, we mapped the DEGs to protein-protein interaction network. We observed that the down-regulated proteins, including JDP2, HADC4, HDAC5, CDYL2, ACADVL, ACSL1 and BRD4, were key components in the down-regulated PPI network (Figure 5B). It is well established that the downregulation of histone deacetylases can promote osteoclastogenesis ^[26,27], indicating that the downregulation of histone deacetylases and cofactors, such as, HDAC4, HDAC5 and JDP2 may be critical regulators in osteoclastogenesis. Moreover, IL8, CCL25, CNR1, HCAR2, HCAR1, GNG13 and GNG8 were also identified as a module in the down-regulated PPI network (Figure 5C). All these proteins except IL8 participated in G-alpha (ⅰ) signaling events, which could inhibit cAMP dependent pathway ^[28]. Previous study ^[29] also reported that cAMP-PKA could regulate osteoclast differentiation, indicating that this pathway may be important for bone mineral homeostasis.

To predict the potential functional roles of the differentially expressed lncRNAs, the Pearson correlation coefficient for lncRNA-DE-mRNA pairs was calculated to construct the co-expression network for mRNAs and lncRNAs. Particularly, RP11-498C9.17 was highly correlated with epigenetic regulators, such as HDAC4, MORF4L1, HMGA1 and DND1. As we described above, the downregulation of histone deacetylases can promote osteoclastogenesis, indicating that the lncRNA RP11-498C9.17 may play a key role in bone mineral homeostasis via controlling osteoclastogenesis.

In conclusion, our integrative analysis revealed key mRNAs and lncRNAs, involved in bone mineral homeostasis and osteoclastogenesis. The results not only broaden our insights into lncRNAs in bone mineral homeostasis and osteoclastogenesis, but also improved our understanding of molecular mechanism.

Acknowledgements

This work were financially supported by the Doctoral Fund of Inner Mongolia Normal University under Grant No. 100900091719, the Scientific research projects of the Inner Mongolian higher educational system(NJZY19025)and Provincial Nature Science Research Project of Anhui Colleges (KJ2018A0331). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Conflict of interest

The authors declare that they have no competing interests.