Employing single-cell RNA sequencing coupled with an array of bioinformatics approaches to ascertain the shared genetic characteristics between osteoporosis and obesity

Abstract

Aims: This study examined the relationship between obesity (OB) and osteoporosis (OP), aiming to identify shared genetic markers and molecular mechanisms to facilitate the development of therapies that target both conditions simultaneously.

Methods: Using weighted gene co-expression network analysis (WGCNA), we analyzed datasets from the Gene Expression Omnibus (GEO) database to identify co-expressed gene modules in OB and OP. These modules underwent Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment and protein-protein interaction analysis to discover Hub genes. Machine learning refined the gene selection, with further validation using additional datasets. Single-cell analysis emphasized specific cell subpopulations, and enzyme-linked immunosorbent assay (ELISA), protein blotting, and cellular staining were used to investigate key genes.

Results: WGCNA revealed critical gene modules for OB and OP, identifying the Toll-like receptor (TLR) signalling pathway as a common factor. TLR2 was the most significant gene, with a pronounced expression in macrophages. Elevated TLR2 expression correlated with increased adipose accumulation, inflammation, and osteoclast differentiation, linking it to OP development.

Conclusion: Our study underscores the pivotal role of TLR2 in connecting OP and OB. It highlights the influence of TLR2 in macrophages, driving both diseases through a pro-inflammatory mechanism. These insights propose TLR2 as a potential dual therapeutic target for treating OP and OB.

Fig. 1

Study flowchart. DEG, differentially expressed gene; GS1, gene set 1; GS2, gene set 2; GSE, Gene Expression Omnibus; KEGG, Kyoto Encyclopedia of Genes and Genomes; limma, linear models for microarray data; OB, obesity; OP, osteoporosis; PPI, protein-protein interaction; ROC, receiver operating characteristic; WGCNA, weighted gene co-expression network analysis.

Fig. 2

Weighted gene co-expression network analysis (WGCNA) identified modules of genes associated with osteoporosis and obesity in the test set. a) Both graphs elucidate the scale-free properties of the gene co-expression network and changes in the mean node connectivity, opting for 7 as the soft-threshold value. b) Both graphs collectively showcase the verification of the scale-free nature of the gene co-expression network, where R^2 = 0.87 and the slope = -2.1, aligning with the established criteria. c) Analogous to Figure 2a, albeit with different values; 20 is selected as the soft threshold. d) Mirroring Figure 2b, but with distinct values, R^2 = 0.81 and slope = -1.32, both meeting the prerequisites. e) Topological overlap matrix (TOM) heatmap delineates the TOM of all genes examined. Paler hues signify low overlap, transitioning to deeper reds indicating increased overlap. Darkened blocks along the diagonal represent modules. The dendrogram of genes and module allocation is also presented on the left and top. The TOM plot facilitates a refined comprehension of inter-gene relationships. f) and g) The module clinical trait correlation heatmap, a visualization method within WGCNA analysis, manifests the associations between modules and clinical features. Within the heatmap, the depth of colour directly corresponds to the magnitude of correlation: red depicts positive correlation, while blue signifies inverse associations. Each cell enumerates the correlation and its statistical significance. All p-values were calculated with Pearson correlation coefficient.

Fig. 3

Analysis of gene intersections within the test set module. a) to c) Bar charts are customarily employed to demonstrate the enrichment levels of each pathway, whereas bubble plots elucidate both the enrichment levels and the number of genes for each pathway. In these bubble plots, the colour of the bubble represents the p-value (or analogous statistical values), and the size signifies the number of overlapping genes (i.e. the genes provided that overlap with genes in a given pathway/Gene ontology (GO) term). Circular diagrams offer a vivid representation of the relationships among pathways. d) The protein-protein interaction (PPI) network diagram is depicted, where nodes symbolize proteins, and edges represent interactions between these proteins. Distinct colours and thicknesses of edges denote different sources of interaction. Sky-blue lines: known interactions from curated databases. Purple lines: known interactions from experimental studies. Blue lines: predicted interactions based on gene co-occurrence. Green lines: derived from the literature. e) The Hubgene network map, obtained using cytoHubba through diverse algorithms, assists in refining our comprehension of PPI dynamics and in unveiling pivotal biological processes and pathways. f) Random forest (RF) map of GSE151839 and GSE56814. g) RF graphs showing variable importance of GSE151839 and GSE56814. In Figure 3f, The RF map displays a model’s error rate as the tree count increases, with the x-axis labelled "trees" (0 to 500) and the y-axis labelled "Error" (0 to 0.5). Three lines, each a different colour and style, show varying error trends: the black line represents the overall error rate, while the red and green lines represent the out-of-bag (OOB) error rates for the two factor levels of the group variable. Error rates tend to decrease as more trees are added. The RF graphs of variable importance in Figure 3g show the importance of variables in a RF, listing genes or variables on the y-axis and "MeanDecreaseGini" on the x-axis. This measure reflects the mean reduction in the model’s Gini impurity caused by each variable, with larger values indicating greater importance for prediction in the model. All p-values were calculated with hypergeometric test. EPC, edge percolated component.

Fig. 4

Validation set: correlational analysis of differentially expressed genes (DEGs) associated with osteoporosis (OP) and obesity (OB). a) and b) Heatmap of DEGs, which visually displays the magnitude and clustering of gene expression levels. In the figure, a deeper colour denotes higher gene expression, while a lighter colour indicates lower expression. c) Volcano plot of DEGs combines statistical significance p-values with the fold change (logFC), facilitating a quick and intuitive identification of genes that are significantly altered and hold statistical relevance. d) and e) Bar graphs are commonly used to illustrate the enrichment level of each pathway, while bubble plots can showcase the enrichment levels alongside the number of genes. Within the bubble plot, the colour represents the p-value (or other metrics such as q-value), and the size signifies the number of genes, indicating the overlap between the submitted genes and those in a particular pathway/Gene Ontology (GO) term. Circle diagrams offer a more vivid depiction of the relationships between pathways. f) Box plot elucidates high expression levels of the Toll-like receptor 2 (TLR2) gene in both the OP and OB groups in the validation set, underscoring its statistical significance in gene expression. g) Receiver operating characteristic (ROC) curve of a single gene serves to assess the sensitivity and specificity of that gene as a biomarker for survival. The TLR2 gene demonstrates promising performance in the validation set, as portrayed by its ROC curve. All p-values were calculated with independent-samples t-test.

Fig. 5

Immune infiltration analysis and single-cell resolution studies. a) and b) Heatmaps delineating immune infiltration provide insights into the abundance of immunological cells. They provide a lucid visual representation of the density and clustering scenario of diverse immune cells within the tissues. In the heatmap, a darker hue signifies heightened immune cell abundance, whereas a lighter shade indicates reduced immune cell presence. c) t-distributed stochastic neighbor embedding (t-SNE) analysis of the osteoporosis (OP) single-cell RNA sequencing (scRNA-seq) data reveals 11 distinct cell subpopulations. Density plots highlight the pronounced expression of the Toll-like receptor 2 (TLR2) gene predominantly in macrophages. d) t-SNE examination of the OB scRNA-seq data uncovers nine unique cell subgroups. Density illustrations underscore the preeminent expression of the TLR2 gene mainly in macrophages. All p-values calculated with permutation test. RMSE, root mean square error.

Fig. 6

Elevated expression of Toll-like receptor 2 (TLR2) in macrophages precipitates the onset of obesity and osteoporosis. a) Following transfection, the secretion levels of interleukin (IL)-6 and tumour necrosis factor-α (TNF-α) in the supernatant were quantified using enzyme-linked immunosorbent assay (ELISA). b) Alterations in macrophage morphology were observed before and after treatment with palmitic acid (PA) and small interfering RNA (siRNA)-TLR2. Produced with crystal violet staining, magnification 200×. c) At the protein level, evaluate the silencing efficiency of TLR2 and the expression levels of ARG1. d) The protein expression of TLR2/ARG1 compared with the control group was examined. e) After treating macrophages with PA and siRNA-TLR2, we observed that the localization and expression of iNOS was in a more aggregated state. f) Fluorescence intensity of iNOS was detected using flow cytometry. g) Upon staining cells with tartrate-resistant acid phosphatase (TRAP), the number of multinucleated cells increased in the PA-induced group. This change was reversed in the siRNA-TLR2 group. h) We employed Western blot analysis to detect the expression of proteins associated with osteoclast differentiation. i) The protein expression of Cathepsin K (CTSK)/matrix metalloproteinase-9 (MMP9) compared with the control group was examined. Data are shown as mean and SD. *p < 0.05, **p < 0.005, ***p < 0.001 compared with control cells, calculated with one-way analysis of variance. All representative data from three independent experiments are shown.

Fig. 7

Competitive endogenous RNA (ceRNA) network diagram delineates the Toll-like receptor 2 (TLR2)-associated microRNAs (miRNAs) and long non-coding RNAs (lncRNAs).

Fig. 8

Illustration of the signalling pathway. Arg-1, Arginase 1; IL-6, interleukin-6; iNOS, inducible nitric oxide synthase; OPG, osteoprotegerin; RANK, receptor activator of nuclear factor kappa B; RANKL, RANK ligand; TNF-α, tumour necrosis factor alpha.