1.
Introduction
Osteoarthritis is chronic degenerative joint disease, which may affect the cartilage, synovium, joint ligaments and subchondral bone [1]. In 2005, approximately 27 million people in the United States had OA [2]. The traits of joint lesions are joint space narrowing, osteophytosis, subchondral sclerosis and osteophyte formation. The traditional treatment in the later stages of OA is joint replacement surgery, obviously, its drawbacks include the inability to effectively treat and prevent early-stage OA and the need for invasive the surgery. The exploration of the pathogenesis of early OA, and timely prevention and treatment of OA will be beneficial to improve patient survival condition. Various biochemical compounds involved in OA development have been regarded as disease pathogenesis consisting of several inflammatory response proteins (e.g., interleukin 1β, interleukin 6 and tumor necrosis factor), matrix-degrading enzymes and toll-like receptors secreted by chondrocytes in cartilage [3]. However, some researchers used isotopes and radiography of OA joints to show that osteophyte and sclerosis formation in subchondral bone precede cartilage change [4]. These changes may be mediated by osteoblasts, osteocytes and osteoclast in subchondral bone. Abnormal bone remodeling, bone matrix mineralization and skeletal architecture construction are dysfunctions in osteoarthritic subchondral bone and directly effect on cartilage [5,6,7]. The irregular expression of RANKL had been found in OA osteoblasts [8,9]. osteoblasts may be acted in OA pathogenesis, by producing abnormal levels of A disintegrin and metalloproteinase with thrombospondin motifs (ADAMTS) and MMPs [10]. Therefore, the molecular mechanisms underlying osteoblasts in OA bone should be investigated.
After extensive research, various genes associated with disease treatment were established, including non-coding RNA, DNA methylation, histone modifications and regulatory RNAs [11,12]. Earlier studies have searched for potential target genes involved in OA in bulk transcriptome profiles [13,14,15]. The identification and exploration of differential expression genes in osteoblasts between normal and OA could expand our overall understanding of the role of osteoblasts in OA development. For samples with high heterogeneity such as the joint, traditional high-throughput sequencing technologies can only provide the average of transcript levels of all cells in the sample, but single-cell transcriptome sequencing can accurately characterize the transcriptome of each cell in the sample. Chenlu Li [16] successfully utilized single-cell RNA sequencing and found that COL6A3 and ACTG were validated as key marker genes of fibroblastic-like chondrocytes and fibroblasts in OA. Zhen Wu [17] found that fibronectin1 as key gene of synovial fibroblasts and the pathway of marker genes activated in differential state cells through single-cell expression profiling. Subsequently, we explored significant hub genes between subpopulations osteoblasts of normal and OA sample from single-cell expression profiling.
miRNA consists of 20–24 nucleotides, which binds to the target messenger RNA binding to the 3' untranslated region of the target RNA, and it can inhibit the translation of the target messenger RNA and affect protein expression [18]. miRNAs are recommended as circulating markers for osteoarthritis and are considered as novel therapeutic materials for in clinical treatment [19]. However, limited studies have reported the use of miRNA in OA. Coutinho de Almeida R found that 238 target mRNAs of 62miRNAs were mainly enriched in the 'nervous system development' mediated by miRNA regulatory mechanisms. TAO [20] found that miR34a is overexpressed in cartilage of OA patients and regulates the PI3K/AKT pathway. Therefore, an in-depth study of OA mRNA-miRNA regulatory mechanisms can provide a comprehensive understanding of the pathogenesis of OA and identify potential therapeutic loci.
In our research, we explored the genomic signatures and marker genes in osteoblasts between normal and OA osteoblasts by using scRNA-seq data from osteoblasts cluster derived from OA patient and healthy person. After differential expression analysis, we explored the biological functions of the differential expression genes. We constructed the PPI network and core module to accurately investigate differential gene interactions and screen for hub genes. Then, based on several credible databases, we rigorously screened several potential miRNAs and constructed mRNA-miRNA regulatory networks, in which partial miRNAs were experimentally validated by GSE175961.
2.
Materials and methods
2.1. Processing of scRNA-seq data
The scRNA-seq data of human hip joint tissue of OA and normal were accessed from GSE169396 and GSE147390 with a reading depth of 10× genomics based on Illumina NextSeq 500 via the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) [21,22]. Four normal samples and one OA sample were screened using Seurat package (version 4.0.2) [23]. The cells were subjected to 1000 highly variable genes by using the Normalization function and FindVariableFeatures function in the Seurat package. After using the FindNeighbors function, FindClusters function and RunTSNE function, cells were divided into different clusters according to highly variable genes. Finally, normal and osteoarthritic osteoblasts which remarked by marker genes (ALPL and RUNX2) were merged into one matrix by using the cbind function of R.
2.2. Quality control
We utilized the Seurat package in R 4.1.1 to process osteoblasts matrix with subsequent analysis [23]. First, the osteoblasts matrix created the object and excluded low-quality cells by using the CreateSeuratObject and subset function. Low-quality osteoblasts were removed according to the following criteria: 1) genes identified < 3 cells; 2) cells of total cell expression genes < 200; 3) volume of total detected RNA < 500 or > 2500; and 4) The proportion of mitochondria and ribosomal expressed genes > 25%. The osteoblasts object was normalized using the NormalizeData function, and 1000 highly variable genes were screened out using the FindVariableFeatures function. To analyze multidmensinoal and complexity scRNA-seq data, we selected the ScaleData function and RunPCA function to linearly scale the data and reduction analysis [24].
2.3. Cell cluster analysis
To visualize and statistically analyze scRNA-seq data, we tested 20 principal component and were chosen for t-SNE. The function of FindNeighbors and FindClusters with the 0.5 clustering resolution were used for classification analysis of cell clusters. Then, we performed cluster identification and visualization in osteoblasts through RunTSNE function [25]. The osteoblast subpopulations, including pre-osteoblasts, mature osteoblasts and undetermined osteoblasts, were remarked through experimentally validated marker genes by manual annotation [22,26,27].
2.4. Differentially expressed genes trajectory analysis
The differentially expressed genes (DEGS) between normal and OA osteoblasts were evaluated using the FindMarkers function as implemented to identify the DEGS of osteoblasts samples. We filtered the marker genes in accordance with cutoff of p_value_adj < 0.05 and |log2[fold change (FC)]| > 1. To further demonstrate the differences among different cell clusters, we used the Monocle package in R 4.1.1 to demonstrate the pseudotime trajectories analysis [28]. Using the differentialGeneTest function to selected test genes. Then, the prerequisite object was constructed as Monocle package, and the multi-dimensional expression profile of cell clusters was reduced to a low-dimensional data and presented in the form of different branches by using the method of DDRTree. By using the orderCells function, each cell was projected onto a dimension and tracked with branching points according to test genes.
2.5. Functional enrichement analysis of marker genes
ClusterProfiler package (version 4.0.5) in R 4.1.1 was utilized to explore the comprehensive molecular biology mechanisms and pathway of osteoblasts [29]. GO and KEGG analysis of marker genes were identified by enrichGO and enrichKEGG function in ClusterProfiler based on the cutoffvalue < 0.01 and cutoffvalue < 0.05. However, GSEA used for the analysis of the enrichment of total DEGS arranged by logFC from largest to smallest in ClusterProfiler package with the P-value < 0.05.
2.6. Construction of the PPI network and identification of core module
The Search Tool for the Retrieval of Interacting Genes (STRING V11.5, https://string-db.org/) database was used to identify the hub gene and visualize the PPI among the marker genes. The confidence PPI network was set minimum required interaction score > 0.4 and exported the result as a tsv.file. The original result was imported to the Cytoscape software (version 3.8.2), and a PPI network was constructed according to the co-expression active interaction source [30]. The core module was selected by Molecular Complex Detection (MCODE) plugin of Cytoscape with the degree Cutoff = 2, node score Cutoff = 0.2, K-core = 7 and max depth = 100. All nodes in the core module were selected as hub genes. The KEGG enrichment analysis of hub genes from top 1 module were applied using CluePedia plugin of Cytoscape. The KEGG enrichment analysis terms, it should have contained at least three genes from top 1 module and at least 3% of all gene characteristic to each term. Regarding the statistical method of the enrichment analyses, a cutoff of P value < 0.05 and kappa score of 0.4 were identified as screening criteria.
2.7. Construction of potential miRNA-mRNA regulatory network
Hub genes were respectively predicted their potential miRNA depending on the three online miRNA databases including miRWalk (http://mirwalk.umm.uni-heidelberg.de/), TargetScan (http://www.targetscan.org/vert_80/) and miRDB (http://mirdb.org/). Not all miRNAs were selected, and the miRNA was found in at three databases to be considered as potential miRNA. The network was established according to bindingp score in Cytoscape software [30].
2.8. Expression levels of miRNA by GSE175961 Dataset
GSE175961 was downloaded from GEO database and used to analysis the miRNA expression levels of miRNA in human knee based on the platform of GPL20712 (Agilent-070156 Human miRNA). Significant differential expression levels of miRNAs between healthy and OA knee were analyzed by R language limma package [31]. The valid miRNA needs to meet two requirements simultaneously: 1) P value < 0.05; 2) The mRNA and miRNA have the opposite expression levels between healthy and OA knee.
3.
Results
3.1. Single-cell RNA-seq profiling of osteoblasts
After overlap the osteoblasts marker genes (ALPL, RUNX2) and each cluster marker genes with normal and osteoarthritis samples, we acquired 885 normal osteoblasts and 4221 osteoarthritic osteoblasts derived from four normal individuals and one patient (Figure S1, Table S1). After excluding low-quality cells with percentages of both ribosomal and mitochondrial sequencing count > 25%, we retained 577 normal osteoblasts and 3810 osteoarthritic osteoblasts for subsequent analysis. We combined the sequencing data of 5106 files into one matrix. The results of the quality control were presented in the form of violin diagrams, showing the detected gene numbers, the sequencing count and proportion of ribosomal and mitochondrial genes in each cell (Figure S2A). A sum of 19,535 corresponding genes were involved, from which 1000 highly variable genes were selected for cell subpopulation analysis, and the top 10 genes were displayed in Figure S2B. Then, we used principal component analysis (PCA) to selected optimal principal component and screened remarkably correlated genes in each component. According to the distribution of cell clusters in PC_1 and PC_2, scRNA-seq data from different sample have no significant batch effects and presented good cell clustering effects shown in Figure S2C. Furthermore, we selected 20 PCs with the P value < 0.05 for cell subpopulation visualization (Figure S2D). To directly visualize the multidimensional scRNA-seq data, we used TSNE projection to precisely classify osteoblasts into seven clusters (Figure 1A). And We used the split TSNE projection to show osteoblasts cells from different samples (Figure 1B). The manually annotated marker genes of three osteoblasts subpopulations are listed in Table S2. The marker genes of three osteoblasts subpopulations expression were shown in Figure 1C. We calculated the marker genes of seven clusters with criteria of p_val_adj < 0.05, in which LEPR enriched in clusters 0, 1, 2, VCAM1 enriched in clusters 0, 1; clusters; COL1A1, SPP1 and IFITM5 enriched in cluster 5; BGLAP and IBSP enriched in clusters 3, 4, 5 and NR4A1 and NR4A2 enriched in cluster 6 (Table S3). Taking into account the subpopulations marker genes at the cell cluster marker genes and cell expression levels, we annotated cell clusters 0, 1, 2 as pre-osteoblasts clusters 3, 4, 5 as mature osteoblast clusters, and 6 as undetermined osteoblast clusters. When the 7-cluster was analyzed for differences with pre-osteoblasts cells, most of the differential genes were ribosomal genes and therefore the 7-cluster was annotated as pre-osteoblasts clusters (Table S4). Pre-osteoblasts included 3009 cells; mature osteoblasts included 1208 cells and undetermined Osteoblasts included 170 cells (Table S5). Then, we exhibited the trajectory analysis to explore the osteoblasts subpopulations differentiation trends. A relatively high proportion of OA osteoblasts consisted of pro-osteoblasts compared with normal osteoblasts, but a relatively low proportion consisted of mature osteoblasts and undetermined osteoblasts (Figure 1E). Exactly as described by the cell annotation results, osteoblasts were correspondingly divided into three directions in the trajectory analysis of scRNA-seq results. From the left to the right of the diagram shows the differentiation of pre-osteoblasts to mature osteoblasts, with undetermined osteoblasts in the middle (Figure 1F).
3.2. Global gene expression characteristics of differently expressed genes
A total of 117 marker genes were included through differential expression analysis between normal and osteoarthritic osteoblasts with |avg_log2FC| > 1 and p_value_adj < 0.05 (Table S6).
GO enrichment analysis was explored the molecular biology functions of marker genes in biological processes, cellular component and molecular function category (Figure 2A, B). In the biological processes category, up-regulated marker genes were remarkably related with post-translational protein modification, regulation of leukocyte migration and positive regulation of leukocytes migration, while the enrichment items of down-regulated marker genes included humoral immune response, defense response to bacterium and neutrophil activation. In the cellular component category, most up-regulated marker genes were associated with collagen-containing extracellular matrix and down-regulated marker genes were associated with secretory granule lumen, cytoplasmic vesicle lumen and vesicle lumen. In the molecular function category, significant enrichment items of up-regulated marker genes included enzyme inhibitor activity, receptor ligand activity and signaling receptor activator activity, while the enrichment items of down-regulated marker genes included glycosaminoglycan binding. The distribution of top-one GO analysis results was shown in Figure 3A, B.
KEGG enrichment analysis was explored the potential pathways of marker genes. The significantly enrichment item of up-regulated marker genes was cholesterol metabolism (P < 0.05, Figure 3C). Down-regulated marker genes were significantly activated in the IL-17 signaling pathway, malaria, NOD-like receptor signaling pathway, AGE-RAGE signaling pathway, TNF signaling pathway, transcriptional mis-regulation in cancer and osteoclast differentiation (P < 0.05, Figure 3D).
3.3. GSEA analysis
GSEA is a genetic statistical model that calculates significant biological pathways after ranking complete DEGS ordered avg_log2FC. The DEGS significantly enriched in IL-17 signaling pathway, NOD-like receptor signaling pathway and ribosome (P < 0.05, Figure 4). As shown in the figure, ribosome mostly enriched in down-regulated marker genes in DEGS. IL-17 and NOD-like receptor signaling pathway share the same linear trend and peak and they are mostly enriched in down-regulated marker genes in DEGS.
3.4. The PPI network and identification of core module
For a good understanding of the interaction correlation of marker genes, we performed STRING online database to construct PPI network. Based on the minimum required interaction score = 0.04, the network consists of 756 edges and 95 nodes. All the genes were imported into degree algorithm method on CytoHubba plugin, and the weight of marker genes with more interactions with other maker genes was visualized (Figure 5A). The most significant cluster with highly interconnected regions in the PPI network was obtained from MCODE plugin. The module involved 18 nodes and 134 edges (Figure 5B). All nodes in the module such as CXCL12, CXCL8, NR4A1, LCN2, JUN, S100A8, IER2, ZFP36, HP, SOCS3, CAMP, JUNB, BTG2, FOSB, RETN, IRF1, EGR1, and ATF3 were regarded as hub genes for the next functional enrichment analysis by using the ClueGO plugin in Cytoscape. Hub genes correspondingly involved in the IL-17 signaling pathway, osteoclast differentiation, TNF signaling pathway, AGE-RAGE signaling pathway in diabetic complications and Pertussis and rheumatoid arthritis in KEGG enrichment analysis (Figure 5C, Table S7). The connection genes between the different KEGG terms include CXCL8, JUN, JUNB, SOCS3, and FOSB.
3.5. The potential miRNA-mRNA regulatory network
To explore the potential regulatory pathways of miRNA, we searched three online databases, namely, TargetScan, miRWalk and miRDB, to predict the corresponding miRNA of hub genes. A miRNA can only be a potential miRNA if it is predicted by all three databases at the same time. Furthermore, we identified 189 potential miRNAs including SOCS3 intersects 2 potential miRNAs, BTG2 intersects 30 potential miRNAs, CXCL12 intersects 64 potential miRNAs, EGR1 intersects 4 potential miRNAs, FOSB intersects 10 potential miRNAs, IER2 intersects 4 potential miRNAs and JUN intersects 1 potential miRNA. To accurately visualize the regulatory network, we imported the intersections between mRNA and miRNA into Cytoscape software (Figure 6A). has-miR-449a and has-miR-218-5p were not only predicted by the database, but were also experimentally validated to be potentially regulated in relation to FOSB and SOCS3 based on the principle that miRNA inhibited the translation of the target messenger RNA. FOSB expressed high levels in normal sample than in OA sample (P < 0.01), but hsa-miR-449a has the opposite profile (P < 0.01, Figures 6B, C). SOCS3 expressed high levels in normal samples (P < 0.01), but hsa-miR-218-5p is highly expressed in OA samples (P < 0.05, Figures 6D, E).
4.
Discussion
Osteoarthritis is a chronic disease in which lesions in the subchondral bone area are characterized by osteophytes formation and sclerosis. OA lesions in the subchondral bone region are mediated by osteocytes and bone matrix. The osteoblasts act in bone production and remodeling, regulate skeletal architecture and bone matrix mineralization by producing extracellular matrix proteins and induce osteoclastogenesis by producing cytokines or direct cell contact [32]. However, with the onset of OA, their function changes dramatically in the subchondral bone region. Traditional RNA bulk-sequencing mainly obtains the average expression of gene transcripts of all cells in a tissue, but not the expression profile of individual cells or a certain type or state of cells, thus obscuring the role of key cell subpopulations or intermediate states. scRNA-seq allows us to obtain the gene expression profiles of individual cells and analyze them to identify the transcriptional profiles of cell types and subtypes, and thus understand the occurrence of OA and explore novel therapeutic loci. In our research, we combined osteoblasts clusters from single-cell datasets GSE169396 and GSE147390 into the same object to explore gene expression differences, molecular mechanisms and miRNA regulatory networks in OA.
In the present study, we identified seven clusters and manually annotated them as pro-osteoblasts cluster, mature osteoblasts cluster and undetermined osteoblasts cluster according to marker gene. After the proposed time series analysis, the osteoblasts subpopulations were divided into three directions, in which the pro-osteoblasts cluster was distributed in the lower left, the mature osteoblasts cluster was distributed in the lower right and the undetermined osteoblasts cluster was distributed in the upper part, which further validated the manual annotation results of the three osteoblasts subpopulations. Undetermined osteoblasts may be in an intermediate stage between pre-osteoblasts and mature osteoblasts. In the present study, osteoblasts in normal and OA samples shared the same subpopulation typology and no subpopulation appeared in separate samples. However, the proportion of subpopulations in the samples changed remarkably, and a higher proportion of mature osteoblasts and undetermined osteoblasts in normal samples and a high proportion of pre-osteoblasts in osteoarthritic samples. Upon the activation of IGF1 and PGE2 hormones, normal pre-osteoblasts further differentiates into mature osteoblasts, while the mineralization of the extracellular matrix occurred [33]. Indira Prasadam et al. [7] found that OA osteoblasts were unable to produce normal minerals and extracellular mesenchyme resulting in the inability to express SOST and DMP late markers, demonstrating that osteoblasts were in an immature state. As an essential survival environment, the extracellular mesenchyme regulates cell proliferation, differentiation and migration through multiple signaling pathways. Significant abnormalities in the extracellular mesenchyme have been observed in OA [34,35,36]. In the present study, the result of GO analysis show that marker genes were remarkably enriched in collagen-containing extracellular matrix. Subsequently, we found that collagen-containing extracellular matrix was highly expressed in the pre-osteoblasts and undetermined osteoblasts 2 clusters, which further validated the above results. Therefore, the undetermined osteoblasts 2 clusters may be a state between pre-osteoblasts and mature osteoblasts. Furthermore, alterations in osteoblasts status and function triggered by the abnormal mass of the extracellular mesenchyme may be involved in the occurrence of OA.
In the bone matrix, overloading of the joint causes microdamage and osteoblasts detecting microdamage can initiate repair by osteoblasts and osteoclasts [6]. IL-17 and TNF-a can mediate the RNAKL-RANK-OPG system to control the communication between osteoblasts and osteoclasts to achieve bone homeostasis [37,38,39,40]. Osteoblasts are important cells in the repair of inflammatory bone tissue damage. Michiel Croes [41] found that IL-17 greatly promotes osteogenic differentiation of human bone marrow mesenchymal stem cells. Zhongxiu Wang [38] treated primary osteoblasts with IL-17 with high expression of mRNA levels of ALP and RUNX2 (markers of early osteoblasts differentiation) and osteocalcin (marker of late differentiation) in a mouse model. Both GSEA and KEGG analyses results, IL-17 signaling pathway was highly activated in down-regulated marker genes. Obtained KEGG analysis of the core module and marker genes, IL-17 and TNF signaling pathway were are jointly regulated by JUN. Osteoclast differentiation was jointly regulated by IL-17 and TNF signaling pathway by JUN. Therefore, we speculate that the abnormal synthesis of IL-17A receptor and TNF-receptor on the surface of osteoblasts results in the inability of IL-17 and TNF to affect osteoblasts effectively, resulting in impaired bone remodeling and hypofractionation of primary osteoblasts, which in turn causes the clinical features of osteoarthritis (e.g., bone redundancy and subchondral osteosclerosis). We hypothesize that JUN may be potential therapeutic loci.
The expression of mRNAs was negatively regulated and repressed by miRNAs. Hence, mRNAs and miRNAs have opposite expression trends. miRNA expression can influence the development of OA [42,43,44]. In recent years, Numerous miRNAs have been identified in osteoarthritic cartilage cells [45]. Soyocak [46] found that miR-146a and miR-155 expressed higher levels in the osteoarthritis human than that of normal persons subjects and increases with the disease duration. Wang [47] found that miR-411 and matrix metalloproteinase 13 were negatively correlated in cartilage tissue in osteoarthritis, and that matrix metalloproteinase 13 was the target gene of miR-411 which overexpression inhibited the expression of matrix metalloproteinase 13, thereby delaying the development of OA. In our research, we searched for 130 potential miRNAs, among which hsa-miR-449a and hsa-miR-218-5p were obtained from the GSE175961 dataset for experimental validation. We found that hsa-miR-449a and FOSB, hsa-miR-218-5p with SOCS3 all exhibited opposite expression trends in OA samples. Thus, hsa-miR-449a and hsa-miR-218-5p may accompany the development and progression of OA. More experiments are still needed to validate and define it.
In addition, cementoblasts, a manner comparable to osteoblasts, expressed various functional TLRs and nucleotide-binding oligomerization domain (NOD) proteins [48]. The expressions of TLRs and NODs were upregulated upon differentiation in both cementoblasts and osteoblasts in a similar manner. Both GSEA and KEGG analysis, showed the NOD-like receptor signaling pathway located in down-regulated area. Therefore, low NOD-receptor expression may delay osteoblasts differentiation, which still needs further study.
Our current research still has limitations. Patient details are incomplete, and some clinical parameters such as laboratory tests, medical and surgical records cannot be added to the discussion, because they cannot be downloaded. Although the genetic profile of osteoblasts and miRNAs are well validated, complementary basic experiments are still necessary to reveal the specific mechanisms underlying the development of OA. We have shared the code on the GitHub website and expect more researchers to communicate with us via email (GitHub: https://github.com/huanchangxiang/Code-of-scRNA-in-osteoarthritis.git).
In conclusion, our study is the first to screen marker gene and potential miRNA for OA osteoblasts based on single-cell transcriptomics, and in addition to this, we have identified OA osteoblasts subpopulations and genetic signatures to explore the role of osteoblasts in the pathogenesis of OA at the genetic level.
Data availability
All data in this study are already available from the Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/).
Conflicts of interest
The authors declare that there is no conflict of interest in this work.
DownLoad:
