Single-Cell Analyses of Heterotopic Ossification: Characteristics of Injury-Related Senescent Fibroblasts.

Introduction: Injury-related cellular senescence may be involved in heterotopic ossification, and no research has been performed about this before. Methods and
Results: The study utilized integrated single-cell RNA-sequencing (scRNA-seq) data from heterotopic ossification samples. The number of senescent cells increased from day 3 and reached the highest level at day 21. However, the expression level of Cyclin Dependent Kinase Inhibitor 2A (Cdkn2a) has no such tendency as the change of cell amount, indicating that the expression level of Cdkn2a may be different in different types of senescent cells or the same time of senescent cell at different time points. The expression level of SASPs (senescence associated secret phenotypes) was also different in different types of senescent cells or at different time points. The GO (gene ontology) analysis revealed that the senescent cells were significantly correlated with the ossification processes, like ECM organization, cell adhesion, ossification, cartilage development, etc. Trajectory analysis showed that injury-related senescent fibroblasts (day 7 and 21) and age-related senescent fibroblasts (day 0 and 42) were in different branches. GO analysis demonstrated that injury-related senescent fibroblasts were mainly related to ossification and ECM remodeling. The KEGG (Kyoto Encyclopedia of Genes and Genomes) results revealed that the ossification was significantly corrected with protein processing in PI3K-Akt signaling, MAPK signaling, focal adhesion, etc.
Conclusion: Consequently, we demonstrated that, unlike age-related senescence, the injury-related senescence demonstrated significantly different SASP phenotypes. The injury-related senescence of fibroblasts is associated with heterotopic ossification formation and may act through PI3K/Akt-induced SASPs.

Introduction

Heterotopic ossification (HO) is common in clinic, with the pathogenesis remains unknown and the treatments limited. HO can occur in any soft tissue of the body, no matter muscle, tendon, ligament, vascular wall, etc. The etiology of HO can be divided into hereditary, neurogenic, traumatic, degenerative, and so on. The pathogenesis can be divided into two categories: 1) Non cell mediated heterotopic ossification, characterized by direct calcium deposition, without osteoblast participation, such as dystrophic ossification. 2) Cell mediated heterotopic ossification, which can also be divided into two kinds, heterotopic intramembranous ossification, and heterotopic endochondral ossification. Heterotopic endochondral ossification is the main type that happens in clinic, such as neurogenic, traumatic, and degenerative ectopic ossification, which are all formed by endochondral ossification. The pathologies of heterotopic endochondral ossification remain unclear; however, it is widely believed that heterotopic endochondral ossification shares much common points with skeletal development and fracture healing. The pathology of heterotopic endochondral ossification is a precise and programmed pathological process, composed of four consecutive processes: injury/inflammation, stem cell recruitment, chondrogenic differentiation, and finally ossification. The injury-inflammation-stem cell recruitment-regeneration process is common in wound healing, while how the stem cells differentiated into chondro/osteogenic systems but not tissue regeneration on HO is far from understood.

Senescence widely exists in human, primate, and rodent tissues, and can appear under a variety of diseases or inducements, such as DNA damage, telomere dysfunction, oncogene activation, organelle damage and so on. In the face of these pressures, especially DNA damage, cells will start self-protection and stop self-replication or self-division. However, these cells do not die obediently, but exist like “zombies”, lose the function of original cells, and become the so-called senescent cells. It was previously thought that senescent cells would continuously secrete some proteins (senescence associated secret phenotype, SASP) to resist the arrival of death. It is usually a mixture of various cytokines, chemokines, growth factors and matrix basic proteases. By cooperating with the immune system, these factors not only protect the senescent cell from death but also affect the proliferation and differentiation of adjacent cells, resulting in negative effects such as organ aging, tissue inflammation, tumorigenesis, etc. This negative effect has been reported in a lot of diseases, such as osteoarthritis, Alzheimer’s disease, atherosclerosis and so on. Targeted removal of these aging cells will also significantly delay the occurrence of aging diseases or reduce the symptoms of aging diseases. For example, the disease-free life span of mice can be prolonged by scavenging aging cells (senolytic drugs), and the average total life span is prolonged by 17%–35%. For mice in the early stage of atherosclerosis, clearing the aging cells in the body can eliminate its early lesions, and even for the late lesions, the senolytic drugs can greatly reduce the formation of atherosclerotic plaque. However, in addition to these negative reactions, senescent cells may also play a two-way role in promoting tissue regeneration. More and more evidence shows that the accumulation of senescent cells after tissue injury may induce cell reprogramming and tissue remodeling. In muscle tissue, both acute and chronic injury can significantly induce transcription mediated cell reprogramming, which is closely related to SASP component IL-6.

The role of cellular senescence in HO has not been researched before, but several traces could be found, indicating its potential role in inducing HO. Firstly, as we all know, adult terminal differentiated cells, whether in vivo or in vitro, can undergo cell reprogramming, dedifferentiation, or trans-differentiation into undifferentiated cells (cell plasticity) under some special conditions, thus participating in some physiological and pathological processes.1 Cell reprogramming is important in tissue regeneration, and senescence-related reprogramming has been proven to be an important part of tissue reconstruction.2– SASP, including IL-6, activin A, and some other unknown content, can dedifferentiate surrounding cells into pluripotent cells by regulating the expression of Oct4, Sox2, KLF4, and c-Myc (OSKM), so as to establish the basis of tissue repair.5,6 In p16INK4a/ARF knockout mice (lacking senescent cells), the ability of cell reprogramming was seriously reduced, and the inhibition of NF-kB activity (an important upstream of SASP) will also greatly reduce the occurrence of cell reprogramming.7,8 Secondly, senescence has also been directly linked with chondrogenic differentiation. Activin A is a member of the TGF-β superfamily, acts as an activator of the BMP pathway, which is a known key pathway related to chondrogenic differentiation. The role of activin A in HO has also been widely researched. In fibrodysplasia ossificans progressiva (FOP), activin A can directly promote proximal chondrogenic differentiation and thus inducing FOP. IL-6/STAT3 signaling is also able to activate HIF1α (known strong chondrogenic differentiation factor), thus participating in the chondrogenesis process. Similarly, other components of SASP, such as PDGF and MMPs, have also been widely proven in previous studies to promote the formation of fibrocartilage matrix.9 Most importantly, our previous research demonstrated that injury-induced senescent cell burden and the SASP contribute to FOP lesion formation and that tissue reprogramming in FOP is mediated by cellular senescence, altering myogenic cell fate toward a chondrogenic cell fate. Furthermore, pharmacological removal of senescent cells abrogates tissue reprogramming and HO formation.10

Consequently, we hypothesized that injury-related cellular senescence was involved in heterotopic ossification, and we utilized bioinformatics analysis of transcriptomic data of HO model to verify our hypothesis.

Materials and Methods

Data Acquisition and Preprocessing

The raw scRNA-seq data (GSE150995 and GSE126060) was downloaded from GEO database (). Both datasets were based on mouse burn and tenotomy model. Briefly, a complete transection of the Achilles tendon together with a dorsum burn injury were performed. The samples of the Achilles tendon were collected at d0, d3, d7, d21, and d42 post-surgery. A total of 20 samples and 42,826 cells were included in the final analysis. The present study involves only mouse data and is a study based on GEO data; thus, there is no need of an Institutional Review Board (IRB) or local ethics review board.

scRNA-Seq Data Analysis

The raw scRNA-seq data was processed using the Seurat v4.0, SingleR 1.6.1, and Monocle 2.20.0 packages in R 4.0.5. First, 3877 low-quality cells were excluded based on the following quality control standards: 1) genes detected in <3 cells were excluded; 2) cells with <200 or >5000 total detected genes were excluded; and 3) cells with ≥25% of mitochondria-expressed genes were excluded. Then, the gene expression of the remaining cells was normalized using a linear regression model. Highly variable genes were detected by the “FindVariableFeatures” function, and then the principal component analysis (PCA) was performed to identify significantly available dimensions with a P value <0.05.

Then, the Uniform Manifold Approximation and Projection for Dimension Reduction (UMAP) algorithm was applied for dimensionality reduction with 20 initial PCs (as demonstrated by Jackstraw and Elbow) and for performing cluster classification analysis across all cells. The differential expression analysis among all genes within cell clusters was performed using the FindAllMarkers algorithm in R to identify the marker genes of each cluster. An adjusted P value <0.05, expression percent >0.25, and | log2[fold change (FC)] | >0.25 were set as the cutoff criteria for identifying marker genes.

Afterwards, different cell clusters were determined and annotated by the singleR package according to the composition patterns of the marker genes and were then manually verified and corrected according to the markers reported in published papers.11,12 The Seurat functions DotPlot, Vlnplot, FeaturePlot, and Heatmap were used to visualize the gene expression with dot plot, violin plot, feature plot, and heatmap, respectively. The ggpubr package was used for statistical analysis of the target gene expression between groups. The Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) enrichments were also performed with the ClusterProfiler 4.0.2 package to detect the pathways in senescence cells.

Trajectory Analysis

Single-cell pseudotime trajectories were constructed using the Monocle 2.20.0 algorithm.14 This algorithm adopts a machine learning technique, learning a parsimonious principal graph to reduce the given high-dimensional expression profiles to a low-dimensional space. Single cells were projected onto this space and ordered into a trajectory with branch points. For data interpretation, the cells in the same branch were generally considered to be in the same differentiation state, while cells located in different branches were considered to have different cell differentiation characteristics. In addition, differential expression analysis was performed between branches and genes that showed differential expression levels were defined as branch-dependent or state-specific genes or marker genes.

Results

Following the quality control standard and the normalization of scRNA-seq data, 4137 low-quality cells were excluded, and 38,689 cells from different time points were included in the final analysis (Figure 1A and B). Unsupervised clustering was applied, and 41 transcriptionally unique cell clusters were identified at the injury site, and then annotated by singleR, and manually corrected according to the reported markers (). Finally, a total of 16 clusters were found (Figure 1C), including macrophage cells, DC cells, granulocyte cells, NK cells, skeletal muscle cells, pericyte/smooth muscle cells, endothelial cells, tenocytes, fibroblasts, MSC, neuromuscular cells, nerve cells, skin fibroblasts, and a remaining small uncharacterized cluster.

Figure 1

Identification of 13 cell clusters in HO samples based on single-cell RNA-seq data. (A and B) After quality control of the 42,826 cells from 20 Achilles tendon samples, 38,689 cells were finally included in the analysis. (C) The uMap algorithm was applied and 16 cell clusters were successfully classified.

As the expression of Cdkn2a and its protein product p16 is one of the most important and frequently used markers of senescence and aging, the proportion and expression level of Cdkn2a were analyzed in all cells of different time points. The proportion of cells expressing Cdkn2a elevated from day 3, reached the highest level at day 7 and day 21, and then dropped significantly at day 42 (Figure 2A). However, there was no significant difference in the expression levels of Cdkn2a at different time points (Figure 2B). The expression level is interesting and confusing as it does not increase significantly along with time as we suppose. On the contrary, the expression level decreased significantly after surgery (d3, d7, and d21 vs d0), though the cdkn2a positive cells increased significantly. Another fact is that the expression level of cdkn2a is similar to d0 and d42, and as we all know, most of the healing process completes before d42. An explanation is that senescent cell types at day 42 are similar with day 0 (normally aging related senescent cells) but much different from injury-related senescence (d3, d7, and d21). Therefore, the difference of the senescent cell types or the different reasons of senescence may influence the expression of cdkn2a.

Figure 2

The proportion and expression level of Cdkn2a in all cells of different time points. (A) The proportion of cells expressing Cdkn2a. (B) The expression levels of Cdkn2a at different time points. *p < 0.05; **p < 0.01; ***p < 0.001).

Then, we divided all the cells into two subgroups according to the expression level of cdkn2a, the senescence group (cdkn2a count >0, 1325 cells) and the non-senescence group (cdkn2a count = 0.37364 cells). Using a list of previously characterized senescence markers, we also plotted the expression level of these markers in cdkn2a+ cells against cdkn2a− cells (Table 1). Unlike those previously published data, a lot of the reported markers found no difference between the groups, like Parp6, Il6, etc. As we all know, the senescence marker differs significantly across different pathologies and tissues, and no consensus has been reached about markers of the injury-related senescence.

Table 1

The Expression Level of Senescence Markers in cdkn2a+ and cdkn2a− Cells (Ns, Not Statistically Significant; *p < 0.05;**p < 0.01; ***p < 0.001; ****p < 0.0001)

Gene Symbol	p	p. Adj	p. Format	p. Signif
Sirt6	4.14E-07	4.1E-07	4.1E-07	****
Parp6	0.366906	0.37	0.37	ns
Serpine1	3.02E-55	3E-55	<2e-16	****
Csf1	0.03657	0.037	0.037	*
Il12b	0.002411	0.0024	0.0024	**
Parp11	0.013418	0.013	0.013	*
Trp53	2.08E-31	2.1E-31	<2e-16	****
Wbscr27	1.6E-10	1.6E-10	1.6E-10	****
Cdkn2a	0	0	<2e-16	****
Dao	0.001294	0.0013	0.0013	**
Zc3hav1	0.006637	0.0066	0.0066	**
Parp1	0.899342	0.9	0.9	ns
Nampt	2.47E-52	2.5E-52	<2e-16	****
Nmnat1	0.39344	0.39	0.39	ns
Sirt5	0.000126	0.00013	0.00013	***
Pvrl4	6.42E-15	6.4E-15	6.4E-15	****
Angptl2	4.5E-206	4.5E-206	<2e-16	****
Bcl2l2	9.73E-07	9.7E-07	9.7E-07	****
Sirt2	0.20237	0.2	0.2	ns
Ppm1d	0.000552	0.00055	0.00055	***
Evl	4.77E-07	4.8E-07	4.8E-07	****
Ccdc74a	1.22E-07	1.2E-07	1.2E-07	****
Cd38	1.5E-131	1.5E-131	<2e-16	****
Nxph4	3.47E-05	0.000035	0.000035	****
Retnla	4.4E-37	4.4E-37	<2e-16	****
Sirt1	4.55E-18	4.5E-18	<2e-16	****
Fabp4	1.85E-71	1.8E-71	<2e-16	****
Nmnat3	1.76E-07	1.8E-07	1.8E-07	****
Prodh	0.000356	0.00036	0.00036	***
Ier5	1.71E-52	1.7E-52	<2e-16	****
Itgax	2.31E-47	2.3E-47	<2e-16	****
Il1b	1.1E-192	1.1E-192	<2e-16	****
Parp9	7.32E-05	0.000073	0.000073	****
Nfkbiz	1E-46	1E-46	<2e-16	****
Tnks	0.000174	0.00017	0.00017	***
Il18	1.31E-10	1.3E-10	1.3E-10	****
Arg1	1.76E-63	1.8E-63	<2e-16	****
Epn3	0.000152	0.00015	0.00015	***
Btg2	2.19E-59	2.2E-59	<2e-16	****
Il10	1.81E-29	1.8E-29	<2e-16	****
Crabp2	1.61E-89	1.6E-89	<2e-16	****
Itgam	1.3E-270	1.3E-270	<2e-16	****
Fcgr1	2.85E-235	2.9E-235	<2e-16	****
Adgre1	4.39 E-214	4.40 E-214	<2e-16	****
Cdkn1a	2E-157	2E-157	<2e-16	****
Parp14	1.7E-64	1.7E-64	<2e-16	****
Parp16	0.00251	0.0025	0.0025	**
Parp4	0.076721	0.077	0.077	ns
Il6	0.412813	0.41	0.41	ns
Parp10	1.97E-20	2E-20	<2e-16	****
Ccl2	5.3E-08	5.3E-08	5.3E-08	****
Cdkn2b	1.2E-108	1.2E-108	<2e-16	****
Sirt7	9.55E-10	9.6E-10	9.6E-10	****
Nadk	3.19E-10	3.2E-10	3.2E-10	****
Ly6d	1.55E-64	1.5E-64	<2e-16	****
Lyve1	1.76E-27	1.8E-27	<2e-16	****
Bmi1	1.81E-15	1.8E-15	1.8E-15	****
Nfkb1	3.4E-34	3.4E-34	<2e-16	****
Il12a	0.021151	0.021	0.021	*
Fcna	1.31E-43	1.3E-43	<2e-16	****
Parp8	9.43E-14	9.4E-14	9.4E-14	****
Igfbp2	3.86E-08	3.9E-08	3.9E-08	****
Slc52a2	2.2E-15	2.2E-15	2.2E-15	****
Parp12	0.105807	0.11	0.11	ns
Irg1	1.99E-25	2E-25	<2e-16	****
Parp3	2.02E-25	2E-25	<2e-16	****
E2f2	0.018451	0.018	0.018	*
Mrc1	2.4E-224	2.4E-224	<2e-16	****
E2f7	0.076711	0.077	0.077	ns
Sarm1	0.232234	0.23	0.23	ns
Sirt4	0.001185	0.0012	0.0012	**
Lmnb1	6.39E-09	6.4E-09	6.4E-09	****
Tiparp	3.93E-64	3.9E-64	<2e-16	****
Nnmt	1.3E-168	1.3E-168	<2e-16	****
Bst1	3.1E-185	3.1E-185	<2e-16	****
Slc48a1	0.747632	0.75	0.75	ns
Tnf	2.1E-175	2.1E-175	<2e-16	****
Hmgb1	2.5E-121	2.5E-121	<2e-16	****
Nt5e	3.47E-15	3.5E-15	3.5E-15	****
Sulf2	6.2E-147	6.2E-147	<2e-16	****
Cd68	1E-271	1E-271	<2e-16	****
Cd36	8.7E-225	8.7E-225	<2e-16	****
Sirt3	0.002598	0.0026	0.0026	**
Mgl2	1.44E-42	1.4E-42	<2e-16	****
Nmnat2	0.002517	0.0025	0.0025	**
Parp2	0.79456	0.79	0.79	ns
Loxl4	0.085264	0.085	0.085	ns

Then, these senescent cells were further analyzed to explore the senescent cell types. A total of 2000 highly variable genes were detected with the variance analysis (Figure 3A). PCA was performed to identify available dimensions and screen correlated genes. The top 20 significantly correlated genes are displayed as dot plots and heatmaps (Figure 3B and C). However, the PCA results did not demonstrate clear separations among cells (Figure 3D). We selected 20 principal components (PCs) with an estimated P value <0.05 for subsequent analysis (Figure 3E and F).

Figure 3

Identification of 6 cell clusters with diverse annotations based on single-cell RNA-seq data. (A) The variance diagram shows 21,289 corresponding genes throughout all cells. The red dots represent highly variable genes, and the black dots represent nonvariable genes. (B and C) The top 20 significantly correlated genes are displayed as dot plots and heatmaps. (D) PCA did not demonstrate clear separations of the senescent cells. (E and F) PCA identified the 20 PCs with an estimated P value < 0.05. (G and H) The uMAP algorithm was applied for dimensionality reduction with the 20 PCs, 14 cell clusters were successfully classified, and finally 6 cell clusters were identified.

Afterwards, the UMAP algorithm was applied, and cells were successfully classified into 14 separate clusters (Figure 3G). Differential expression analysis was performed, and a total of 7449 marker genes from all clusters were identified ( and ). According to the expression patterns of the marker genes, these clusters were annotated by singleR, and then manually corrected according to the reported markers (Figure 3H, ). Clusters 0–7, containing 1098 cells, were annotated as fibroblasts; cluster 8, containing 65 cells, was annotated as macrophages, gran, or DC cells; cluster 9, containing 52 cells, was annotated as skMuse/Peri/SMC; clusters 10 and 11, containing 71 cells, were annotated as MSCs; cluster 12, containing 21 cells, was annotated as T cells; and cluster 13, containing 18 cells, was annotated as endothelial cells.

As we all know, different types of senescent cells may have different SASP components. It is interesting to figure out the difference of the SASP among the different cell types. The ANOVA analysis was used to compare the mean expression level of the SASPs of different cell types (Table 2). The expression of some of the most discussed SASPs (TNF-α, IL-6, IL-1β, CXCL10, CCL5, CCL2, MMP3, and PAI-1/Serpine1) are shown in Figure 4A–H.

Table 2

The Mean Expression Level of the SASPs of Different Senescent Cell Types (Ns, Not Statistically Significant; *p < 0.05;**p < 0.01; ***p < 0.001; ****p < 0.0001)

Adipoq	0.981954	0.98	0.98	ns	Anova
Ang	0.001664	0.0017	0.0017	**	Method
Areg	3.8E-30	3.8E-30	<2e-16	****	Method
Areg	3.8E-30	3.8E-30	<2e-16	****	Method
Bcl2	0.0316	0.032	0.032	*	#N/A
Btc	0.735402	0.74	0.74	ns	Anova
Ccl1	0.997819	1	1	ns	Anova
Ccl2	0.004025	0.004	0.004	**	Method
Ccl20	0.910174	0.91	0.91	ns	Anova
Ccl27a	0.015168	0.015	0.015	*	Method
Ccl3	1.44E-50	1.4E-50	<2e-16	****	Method
Ccl5	1.28E-44	1.3E-44	<2e-16	****	Method
Ccl7	1.75E-06	1.7E-06	1.7E-06	****	Method
Ccl8	1.86E-15	1.9E-15	1.9E-15	****	Method
Cd44	9.37E-10	9.4E-10	9.4E-10	****	Method
Csf1	0.157989	0.16	0.16	ns	Anova
Csf2	0.991921	0.99	0.99	ns	Anova
Csf3	0.338366	0.34	0.34	ns	Anova
Ctgf	1.74E-07	1.7E-07	1.7E-07	****	Method
Cxcl12	1.11E-32	1.1E-32	<2e-16	****	Method
cxcl2	5.42E-26	5.4E-26	<2e-16	****	Method
Cxcl5	5.18E-06	5.2E-06	5.2E-06	****	Method
Cxcr2	0.936235	0.94	0.94	ns	Anova
Cyr61	2.75E-22	2.7E-22	<2e-16	****	Method
Egfr	2.12E-17	2.1E-17	<2e-16	****	Method
Fas	0.668002	0.67	0.67	ns	Anova
Fgf1	0.276775	0.28	0.28	ns	Anova
Fgf2	0.003874	0.0039	0.0039	**	Method
Fgf7	4.44E-08	4.4E-08	4.4E-08	****	Method
Fn1	3.9E-101	3.9E-101	<2e-16	****	Method
Gdnf	0.718293	0.72	0.72	ns	Anova
Hgf	0.242066	0.24	0.24	ns	Anova
Icam1	1.95E-17	2E-17	<2e-16	****	Method
Ifng	1.88E-28	1.9E-28	<2e-16	****	Method
Igfbp2	0.441042	0.44	0.44	ns	Anova
Igfbp3	1.62E-15	1.6E-15	1.6E-15	****	Method
Igfbp4	6.31E-60	6.3E-60	<2e-16	****	Method
Igfbp5	1.08E-06	1.1E-06	1.1E-06	****	Method
Igfbp6	4.93E-33	4.9E-33	<2e-16	****	Method
Igfbp7	1.1E-141	1.1E-141	<2e-16	****	Method
Il11	0.161522	0.16	0.16	ns	Anova
Il15	0.002351	0.0024	0.0024	**	Method
Il1A	0.0151	0.015	0.015	*	Method
Il1B	4.2E-59	4.2E-59	<2e-16	****	Method
Il1R1	3.23E-07	3.2E-07	3.2E-07	****	Method
Il2Ra	2.7E-158	2.7E-158	<2e-16	****	Method
Il6	5.24E-22	5.2E-22	<2e-16	****	Method
Il6st	1.91E-14	1.9E-14	1.9E-14	****	Method
Il7	0.958777	0.96	0.96	ns	Anova
Lif	0.082219	0.082	0.082	ns	Anova
Mcpt4	0.553296	0.55	0.55	ns	Anova
Mif	6.59E-14	6.6E-14	6.6E-14	****	Method
Mmp10	0.835597	0.84	0.84	ns	Anova
Mmp12	0.000179	0.00018	0.00018	***	Method
Mmp13	0.143169	0.14	0.14	ns	Anova
Mmp14	3.26E-45	3.3E-45	<2e-16	****	Method
Mmp2	2.2E-121	2.2E-121	<2e-16	****	Method
Mmp3	0.000762	0.00076	0.00076	***	Method
Nap1l4	6.71E-05	0.000067	0.000067	****	Method
Osm	9.14E-46	9.1E-46	<2e-16	****	Method
Pdgfa	9.58E-45	9.6E-45	<2e-16	****	Method
Pdgfb	6.6E-100	6.6E-100	<2e-16	****	Method
Pigf	0.037762	0.038	0.038	*	Method
Plat	8.33E-16	8.3E-16	8.3E-16	****	Method
Plau	0.006943	0.0069	0.0069	**	Method
Plaur	3.99E-09	4E-09	4E-09	****	Method
Ptger2	3.12E-09	3.1E-09	3.1E-09	****	Method
Saa1	0.801068	0.8	0.8	ns	Anova
Serpinb2	0.973808	0.97	0.97	ns	Anova
Serpine1	2.12E-05	0.000021	0.000021	****	Method
Spink1	0.982182	0.98	0.98	ns	Anova
Tgfb1	1.82E-36	1.8E-36	<2e-16	****	Method
Thpo	0.859743	0.86	0.86	ns	Anova
Timp1	8.25E-50	8.3E-50	<2e-16	****	Method
Timp2	3.27E-83	3.3E-83	<2e-16	****	Method
Tnf	2E-48	2E-48	<2e-16	****	Method
Tnfrsf18	1.3E-221	1.3E-221	<2e-16	****	Method
Tnfrsf1a	2.91E-13	2.9E-13	2.9E-13	****	Method
Tnfrsf1b	8.41E-65	8.4E-65	<2e-16	****	Method
Tnfsf10	2.99E-20	3E-20	<2e-16	****	Method
Tnfsf11	0.000736	0.00074	0.00074	***	Method
Tsc1	0.526758	0.53	0.53	ns	Anova
Vegfa	0.003018	0.003	0.003	**	Method
Wnt16	0.024772	0.025	0.025	*	Method
Wnt2	0.462268	0.46	0.46	ns	Anova
Wnt5A	9.13E-05	0.000091	0.000091	****	Method

Figure 4

The expression level of SASPs in different senescent cell types and in fibroblasts at different time points. (A–H) The expression level of SASPs in different senescent cell types. (I) The expression level of SASPs in fibroblasts at different time points.

Also, as discussed above, the same type of cell may express different types of SASP under different inducements of senescence. The expression level of cdkn2a differs significantly at day 0/42 (normal senescence or age-related senescence) and day 3/7/21 (injury-related senescence). In this study, we compared the expression of SASPs in fibroblasts at different time points (Figure 4I). The expression of SASPs in d0 was much different from that in other time points, which means the age-related and injury-induced senescence were totally different molecular processes. Though the number of senescent cells in day 3 is not high, the expression level of SASPs is the highest in all time points, except Ccl5 and Pai1.

KEGG and GO analyses were performed to identify the potential molecular mechanisms and pathways of senescence. The GO results revealed that the senescent cells were significantly correlated with the ossification processes, like ECM organization, cell adhesion, ossification, cartilage development, etc. (Figure 5, ). The KEGG results revealed that the senescence was significantly corrected with protein processing in endoplasmic reticulum, PI3K-Akt signaling, ribosome, focal adhesion, etc. Endoplasmic reticulum stress has been widely associated with the process of senescence in a number of diseases.13,14 PI3K-Akt signaling has also been associated with senescence regarding the p21 activity.15,16

Figure 5

The KEGG and GO enrichments of senescent fibroblasts to identify related molecular mechanisms and pathways.

Trajectory analysis was performed to project all senescent fibroblasts onto one root and three branches, termed branches I, II, III, and IV (Figure 6A and B). The results showed that the senescent fibroblasts of day 7 were mainly located in the root and branch I, the senescent fibroblasts of day 21 were mainly located in branch I and II. The senescent fibroblasts of day 0 and d42 were mainly located in branch III, which means the kind or degree of senescent cells were different at different time points. At day 7, the cells were mainly injury-related senescent cells, but at d42, the cells were mainly age-related senescent cells (similar to day 0). Differential analysis was performed, and the significant marker genes were identified as shown in figure (Figure 6C, ). KEGG and GO analysis were also performed to identify related molecular mechanisms and pathways of the different branches (Figure 7, ). It is so interesting that cluster 2 is related to ossification and cluster 3 is related to ECM remodeling, as shown in GO analysis. Especially the cluster 2, which is mainly composed of day 21 senescent cells, is consistent with the pathological findings (ossification-related processes). The KEGG results revealed that the cluster 2 was significantly corrected with protein processing in PI3K-Akt signaling, MAPK signaling, focal adhesion, etc. These findings indicate that senescent cells in distinct differentiation states demonstrate distinct biological characteristics, which might provide new evidence for the molecular signatures and biology of senescence, including both intrinsic properties and the regulation of related pathways. The most important finding is the relationship of branch II with ossification.

Figure 6

Cell trajectory analysis of senescent fibroblasts. (A–C) The senescent fibroblasts of day 7 were mainly located in the root and branch I, the senescent fibroblasts of day 21 were mainly located in branch I and II. The senescent fibroblasts of day 0 and d42 were mainly located in branch III, which means the kind or degree of senescent cells were different at different time points. At day 7, the cells were mainly injury-related senescent cells, but at day 42, the cells were mainly age-related senescent cells (similar with day 0). (D) Differential expression analysis of the five different subsets.

Figure 7

GSEA of the senescent fibroblast subsets found in cell trajectory analysis was performed to identify related molecular mechanisms and pathways. An FDR ≤ 0.05 was considered statistically significant (cluster 1, root; cluster 2, branch I; cluster 3, branch IV, cluster 4, branch III; cluster 5, branch II).

Discussion

The present study found that injury-related senescence is associated with heterotopic ossification formation. Cell aging-related secretory phenotype SASP is a complex mixture composed of a variety of cytokines, chemokines, growth factors and matrix basic proteases, including IL-6, MMP-3, MMP-13, COL1A1, IL-8, TNFa, HMGB1, activin A, etc. The exact composition of SASP varies significantly with cell type, tissue environment and stimulation of aging induction. These secretory factors interact with adjacent cells and immune system and finally affect the fate of aging cells and surrounding normal cells. The present study found that injury-related senescent fibroblasts may be related to heterotopic ossification, and the mechanism how senescence participates in ossification remains unknown. The senescent cells have no ability to self-replication or self-division, so there is no possibility for them to directly differentiate into chondro/osteogenic cells. The only way they can go is secreting SASP and then influencing the adjacent cells.

The PI3K/Akt pathway is one of the most important signaling pathways for regulating cellular senescence. It regulates cell growth by acting on TSC1/TSC2 complex and mTORC signaling. PI3K/AKT pathway triggers senescence in human cells by promoting enhanced p53 protein synthesis via mTORC1.17 Akt is a major regulator of cell survival, which is regulated by directly inhibiting pro-apoptotic proteins (such as Bad) or inhibiting pro-apoptotic signals produced by transcription factors (such as FoxO). Akt affects cell proliferation by phosphorylating CDK inhibitors p21 and p27.18 And the PI3K/AKT/mTOR signaling pathway is also important for regulating cellular senescence.19 PI3K/Akt pathway is also important in inducing SASP secretion. High-throughput genetic screening studies combined with gene expression analysis and cytokine arrays indicate that the transcription of SASP inflammatory genes is regulated by NF-κB. NF-κB suppression causes bypass of Akt-induced senescence in fibroblasts and a diminished SASP and NF-κB inhibition contributes to chemoresistance in murine lymphomas.20 If we treat the senescent cells with the PI3K/Akt downstream mTORC1 inhibitor rapamycin, the SASP cytokine secretion will be dramatically impaired.21 Rapamycin could decrease the secretion of 35% of the evaluated SASP factors as well as NF-κB transcriptional activity. Akt-induced senescence also results in an induction of the SASP components IL-1α, IL-1β, IL-6 and IL-8,17,20 which could be inhibited by mTORC1 inhibitor rapamycin. However, as discussed above, the effect of senescence on ossification is mainly through SASP; it is so interesting and important if PI3K/Akt/mTOR signaling is associated with secreting ossification-related SASPs. The PI3K/Akt/mTOR has been linked with multiple SASPs like IL-6, IL-8, ANGPTL4, SPINK1, SERP2, MMP3, AREG, etc. Lots of them have been proven effective in inducing osteogenesis.

However, the present study is only a bioinformatics analysis of transcriptomic data of HO model, and no corresponding verification experiment was performed. In our recently published paper, injury-induced senescent cell burden and the SASPs may contribute to fibrodysplasia ossifican progressiva (FOP) lesion formation. Furthermore, pharmacological removal of senescent cells abrogates tissue reprogramming and HO formation.10 Thus, the relationship between senescence and HO may be exact, but the mechanism or signaling pathway should be further explored by molecular, cell, and animal experiments. Furthermore, though the senescent fibroblast is the main cell type, will the other cell types also contribute to HO formation is still worth exploring. Most importantly, what is the target cell, or in other word, what cell is affected by SASPs and thus reprogrammed to HO progenitors? Further studies should be performed upon these topics.

Conclusion

Taken together, the present study found that injury-related senescence of fibroblasts is associated with heterotopic ossification formation and may act through PI3K/Akt-induced SASPs. Unlike age-related senescence, the injury-related senescence demonstrated significantly different SASP phenotypes.