Comparative proteomic analysis identifies differentially expressed proteins and reveals potential mechanisms of traumatic heterotopic ossification progression.

Background: Traumatic Heterotopic Ossification (tHO) is one of complications of elbow fractures to the detriment of patients' rehabilitation, and the severity of tHO corresponds to the size of ectopic bone. It has yet to be elucidated which proteins and pathways underlying the progression of tHO, and biomarkers to predict the severity of tHO at early stage of the disease also need further investigation.
Methods: In this study, a new rat model with distinct volume of ectopic bone was established first. Then a data-independent acquisition proteomics approach was used to investigate injured site tissues sequentially obtained from these rats (2, 7, 14, and 28 days post-injury). Differentially expressed analysis, functional annotation and co-expression analysis and protein-protein interaction network were performed to explore the pathways and hub proteins in the tHO progression. Clinical samples from a nest case-control study were used to validate the selected proteins for predicting the severity of tHO.
Results: The Achilles Tenotomy (AT) induced significantly larger sizes of ectopic bone compared to Partial Achilles Tenotomy (PAT) in rat models. A total of 3547 quantifiable proteins were screened for differential expression analysis among the AT, PAT and control groups. The hierarchical clustering and expression pattern analysis revealed more apparent difference in the pathways such as oxidative phosphorylation, mitochondrial function, and sirtuin signaling between AT and PAT group at the early stage (2 dpi) of tHO. The co-expression analysis identified five hub proteins, UBA1, EIF3E, RPL17, RPL27, and RPS28. qPCR assay, immunoblot assay and immunohistochemistry assay verified that these proteins had higher expression level in the tissue samples of clinically relevant HO patients and clinically irrelevant HO patients than HO negative patients.
Conclusion: The new established animal model and proteome profile could serve as a solid foundation for the comprehensive investigation of the progression of traumatic heterotopic ossification. And the identified 5 proteins (UBA1, EIF3E, RPL17, RPL27, and RPS28) may serve as potential biomarkers to predict the severity of tHO. The translational potential of this article: The proteins identified in this study may be the potential biomarkers and therapeutic targets for predicting and treating the tHO at early stage.

Introduction

Traumatic heterotopic ossification (tHO) refers to an orthopaedic condition defined by endochondral ossification in the extraskeletal tissues following orthopedic trauma, severe burns, brain and spinal cord injury, or combat-related damage. When ectopic osteogenesis develops around the joint, it may cause swelling, chronic pain, nerve compression, and joint contractures [[1], [2], [3]]. However, HO can be divided into two categories: clinically relevant HO and clinically irrelevant HO according to its severity. While clinically relevant HO may greatly impact the extremity function of patients or may even cause disability, clinically irrelevant HO lesions rarely disturb patients [[4], [5], [6], [7]]. The severity of tHO is mainly depend on the volume of the ectopic bone. According to the experience of Shanghai Sixth people's hospital, larger ectopic bone had higher risk resulting in elbow stiffness [8,9]. However, lacking of method to predict the severity of tHO at early stage, the physicians always have to indiscriminately administrate prophylactic treatment to elbow fracture patients. Up to now, the recognized prophylaxis approach of HO, including the non-steroidal anti-inflammatory drugs (NSAIDs) and low-dose radiation therapy, have side effects such as gastrointestinal toxicity and may result in nonunion of bone fracture [10,11]. For patients who only developed the clinically irrelevant HO, the side effects of prophylactic drugs may outweigh its benefits.

Previous researches on the mechanism of tHO progression had proven that inflammatory niches, biomechanical stimuli, and electromagnetic stimulation at the trauma site activate signaling pathways of the progenitor cells around the trauma site, including the members of the tumor growth factor (TGF)-β superfamily signaling cascades and mitogen-activated protein kinase (MAPK) pathway, which induce chondrogenesis/osteogenesis [12,13]. Other pathways, such as the PI3K-mTOR-AKT axis, retinoic acid receptor (RAR) family, and hypoxia inducible factor-1 alpha (HIF-1α) have also been shown to contribute to osteogenesis and angiogenesis that occur during bone development [[14], [15], [16]]. The terminal mass of ectopic bone is variable, and not always proportional to the severity of the original trauma [17,18]. However, little is known about the mechanism(s) regulating this process and the volume of ectopic bone produced. Considering the correlation between the volume of heterotopic ossification and the severity of tHO, it is vital to obtain a better understanding of the mechanism(s) that regulate the progression of tHO in order to find biomarkers and potential targets for diagnosing and treating tHO at early stage. Nevertheless, we still lack of suitable animal model to imitate the clinically relevant HO and clinically irrelevant HO for mechanism research yet.

Achilles tenotomy (AT) rat model is a classical animal traumatic heterotopic ossification model. In previous experiments, we unexpectedly found that partial Achilles tenotomy (PAT) could be used as a valid method to stably induce less ectopic bone than AT in rat models, without any other pharmacological or transgenic approaches. In current study, we performed Achilles tenotomy (AT) and partial Achilles tenotomy (PAT) in rat models to simulate the distinct size of ectopic bone observed in clinical courses and sequentially collected tissue samples from injury site with the progression of tHO. We then applied a data-independent acquisition (DIA) proteomic strategy, a newly developed proteomic technique with higher reproducibility and throughput compared to the canonical data-dependent acquisition (DDA) approach [19], to identify the proteins and signaling pathways associated with the development of heterotopic ossification. Meanwhile, clinical samples at early stage of tHO formation from a nest case–control study were collected to validate the protein that indicate the severity of tHO. This study provide insight into the mechanism of HO progression as well as the potential biomarkers and targets for monitoring and treating the disease at early stage.

Methods

Animal samples

The study was conducted using 200 ​± ​10 ​g male Sprague Dawley rats (Rattus norvegicus). The rats were provided and maintained by the Shanghai Jiao Tong University Animal Department (Shang Hai, China). All animal experiments were approved by the Animal Ethics and Welfare Committee of Shanghai Jiao Tong University and abided by the U.K. Animals Act and associated guidelines. Achilles' tenotomy (AT) and partial Achilles’ tenotomy (PAT) were inflicted on rats in the AT and PAT group respectively as previously described [20]. The rats were provided free access to food and water, and observed carefully after the operation. After 10 weeks, the rats in the AT and PAT group were euthanized by CO2 asphyxiation, both hind limbs were collected to measure the volume of ectopic bone by micro-CT, and the tissue of the injured site from the hind limbs was harvested for histological examination. For proteomic analysis, injured tissue samples from 2 dpi (days post injury), 7 dpi, 14 dpi, 28 dpi of AT and PAT rat models, and tissue from the sham group (control) were collected. The samples were cleaned in phosphate-buffered saline (PBS) and snap-frozen in liquid nitrogen before being stored at −80 ​°C for protein isolation.

Clinical study design and sample collection

The clinical samples from enrolled subjects used in this study were from another larger retrospective nest case–control study to investigate the biomarkers of HO in early stage. This study was approved by the ethics committee of Shanghai Jiao Tong University Affiliated Sixth People's Hospital (code: 2020-KY-086(K)). All study procedures were followed in accordance with the ethical standards of the World Medical Association and the sponsoring institution. Patients were recruited in the study before they received surgical treatment for elbow fractures at Shanghai Jiao Tong University Affiliated Sixth People's hospital, including patients diagnosed as the fracture which types had a high risk of developing tHO, such as isolated radial head fractures, isolated olecranon fractures, coronoid fractures, distal humerus fractures, terrible triad injuries, fracture-dislocation of elbow, floating elbow injuries, and elbow injuries concomitant with burns, or head injury. Subjects currently being treated for cancers or metabolic disease involving the bone were excluded. Informed consents had been signed by all participants before they enrolled in this study, and permissions from patients to publish the images of physical examination with their privacy well observed had been obtained. Clinical parameters of subjects including age and gender, body weight, type of injuries; presence of compound fracture, time to surgery were recorded when they received surgery in the hospital. Tissue samples were taken from subjects during their open reduction and internal fixation (ORIF) operation at Shanghai Jiao Tong University Affiliated Sixth People's hospital. Briefly, specimens were collected through the surgical approach. After exposing the surgical site, tissue at the edge between devitalized and healthy-appearing tissue that would otherwise be discarded as surgical waste was harvested without causing additional injury than the standard ORIF procedure.

After the surgery, the patients were followed-up for 6 months to get their radiological data and physical examination results. Radiological parameters included the presence or absence of HO on radiographs, as well as the location and size of HO were documented. The results of function and ROM of subjects’ elbows during the period of follow-up were recorded, too. The primary outcome was the case–control status of each participant which was classified into three categories: clinically relevant HO, clinically irrelevant HO or HO negative (control) according to the digital radiography and Mayo Elbow Performance Index during the period of follow up. The clinically relevant HO refers to HO resulting in functional limitation, including Class 2 or Class 3 of Hastings and Graham classification as well as the MEPI <70 [7]. Clinically irrelevant HO was defined as Class 1 of Hastings and Graham classification and the MEPI >70, which causing no functional limitation. And HO negative (control) indicated that no sign of heterotopic ossification was found in radiographic results. Samples were assigned to clinically relevant HO, clinically irrelevant HO or HO negative (control) group based on primary outcome of each subject.

Protein extraction and fractionation for DDA library generation

Samples harvested during the development of tHO were homogenized and lysed in SDT buffer (4% SDS, 100 ​mM DTT, 150 ​mM Tris–HCl, pH 8.0). The lysates were sonicated and boiled for 15 ​min. After samples were centrifuged at 14,000×g for 40 ​min, the protein content in the supernatant was extracted and quantified using the BCA Protein Assay Kit (Bio-Rad, San Francisco, California, USA). Proteins were separated by SDS-PAGE, and the gels were stained with Coomassie Blue. Inter-rater agreement for all samples was examined. An equal aliquot from each sample in this experiment was pooled into one sample for DDA library generation and quality control. Protein digestion was performed according to the FASP procedure [21]. The iRT-Kits (Biognosys, Zurich, Switzerland) were used to rectify the relative retention time differences between runs with a volume proportion of 1:3 for iRT standard peptides versus sample peptides.

Data-dependent acquisition and data-independent acquisition mass spectrometry assay

Data-dependent acquisitions (DDA) were used to generate a spectral library to form a query database for data-independent acquisition (DIA) mass spectra in subsequent analyses. The DDA spectra were analyzed using MaxQuant analytical software and were filtered to achieve an FDR of 1% at the peptide and protein levels. Then, each sample peptide was analyzed using an LC-MS/MS system (Thermo Scientific Q Exactive HF X, Waltham, Massachusetts, USA) operating in DIA mode. Quality control samples (pooled samples from equal aliquots of each sample in the experiment) were injected in the DIA mode, which was used to monitor the MS performance [22].

Data analysis

For DDA library data, the FASTA sequence database was searched using Spectronaut Pulsar X TM_12.0.20491.4 (Biognosys). The database was downloaded from the Ensembl website (http://asia.ensembl.org/index.html). All reported data were based on 99% confidence for protein identification, as determined by a false discovery rate (FDR=N (decoy)∗2/(N (decoy)+ N (target))) ≤1%. The spectral library was constructed by importing the original raw files and DDA search results into the Spectronaut Pulsar X TM. The DIA data were analyzed with Spectronaut Pulsar X TM by searching the above constructed spectral library. All results were filtered by setting the Q value cutoff at 0.01, equivalent to FDR<1% [23]. The mass spectrometry proteomics data were deposited in the ProteomeXchange Consortium (http://proteomecentral.proteomexchange.org) via the iProX partner repository [24]. The dataset was authorized using a PXD number: PXD025678.

Bioinformatics analysis

The protein sequences of differentially expressed proteins (DEPs) were examined in batches retrieved from the UniProtKB database (Release 2016_10) in FASTA format. The retrieved sequences were locally searched against the SwissProt database (rat) using the NCBI BLAST ​+ ​client software (ncbi-blast-2.2.28+-win32. exe) to identify homologous sequences from which the functional annotation could be transferred to the studied sequences. SIMCA-P 14.1, Umetrics, Umea, Sweden) was used to perform principal component analysis to identify and visualize the relationships of DEPs from different injury tissue samples. Gene ontology (GO) mapping and annotation were performed using Blast2GO (Version 3.3.5) [25] and InterProScan [26], and GO annotation results were plotted using R scripts. Ingenuity Pathway Analysis (Qiagen, Redwood City, CA, USA) was used to conduct the pathway enrichment analysis. Hierarchical clustering analysis and heat map visualization were performed using the R package. Temporal expression patterns of the proteins were analyzed by fuzzy C-means clustering using the cmeans function in R package e1071 with Euclidean distance, and the results were visualized using the Python package matplotlib and Blues color palette.

Weighted gene Co-expression network analysis

To identify differentially co-expressed protein modules, the “WGCNA” package in R was used to construct a co-expression network for the proteins that were used for the proteomic analyses. Canonical pathway analysis using Ingenuity Pathway Analysis (IPA) was used for functional annotation of the identified modules. The Eigen proteins of each module were used to measure the association between modules and pathological traits. The correlation between module eigenproteins (MEs) and pathological traits was calculated to identify the modules associated with the severity and development of heterotopic ossification. Then, STRING was utilized to perform the protein–protein interaction (PPI) network to screen for key proteins, and the results were visualized and further analyzed using Cytoscape (3.8.2). The cytoHubba plugin based on Cytoscape was used to identify highly connected hub proteins within the PPI network. The chosen proteins were used for further analysis and validation.

LC-PRM/MS-based quantitative validation of proteomics results

In order to validate the results obtained by DIA-based proteomic analysis, fifteen differential expressed proteins (DEPs) selected from the clinically related modules were quantified using a sensitive and rapid parallel reaction monitoring (PRM)-based LC-MS/MS method [27]. The methods optimized for collision energy, charge state, and retention times for the most significantly regulated peptides were generated experimentally using unique peptides of high intensity and confidence for each target protein. The raw data were analyzed using Skyline (MacCoss Lab, Seattle, Washington, USA) [28] where signal intensities of individual peptide sequences for each of the significantly altered proteins were quantified relative to each sample, and values were normalized to a standard reference.

MicroCT scan

Rat hind limbs were harvested and imaged 10-weeks post-injury by high resolution Micro-CT scanner (Skyscan 1176, software ​= ​Version 1.1 (build 6), Bruker, Kontich, Belgium), with parameter set at 18 ​μm resolution and 70 ​kV voltage. CT images were reconstructed with CTvox software and ectopic bone volume formation was calculated with CTan software (Version 1.15.4.0+, Bruker) at threshold Hounsfield units (HU) of 360 to determine the gross volume of mineralized tissues [29].

Quantitative real-time PCR assay

To extract the total RNA from rat and human tissue, EZ-press RNA Purification Kit (B0004D-100; EZBioscience, Roseville, MN, USA) was applied according to the manufacturer's protocol. We then used the cDNA Reverse Transcription Kit (EZBioscience, Roseville, MN, USA) to accomplish the reverse transcription of 1 ​μg of total RNA. Quantitative analysis of target gene mRNAs was performed using SYBR Green I Master Mix (EZBioscience, Roseville, MN, USA) and a LightCycler® 480 Real-time PCR system (Roche, Basel, Switzerland). The gene primers provided by Sangon Biotech (Shanghai, China) are listed in Table s1 with β-Actin as a housekeeping gene.

Western blot analysis

Tissue was lysed in RIPA buffer lysis system (Epizyme, Shanghai) supplemented with a proteinase inhibitor cocktail (Epizyme, Shanghai) on ice. The supernatant was collected to assess for total protein using BCA Protein Assay. Equal amounts of protein (20 ​μg) were loaded onto sodium dodecyl sulfatepolyacrylamide electrophoresis gels. Separated proteins in gels were transferred to polyvinylidene fluoride membranes. Subsequently, the membranes were blocked in 5% nonfat milk or bovine serum albumin (BSA) for 1 ​h and subsequently probed overnight at the temperature of 4 ​°C with diluted primary antibodies against β-Actin (1:5000, T0022, Affinity biosciences, ChangZhou, China), UBA1(1:2000, ab180125, Abcam Inc, Cambridge, UK), EIF3E (1:1000, ab134958, Abcam Inc, Cambridge, UK), RPL27 (1:500, 14980-1-AP, Proteintech, Chicago, USA), RPL17 (1:5000, 67223-1-Ig, Proteintech, Chicago, USA), RPS28(1:500, 14796-1-AP, Proteintech, Chicago, USA). Membranes were washed in wash buffer (1 ​× ​TBS 0.05% Tween-20). Following incubation with HRP-conjugated secondary antibodies for 1 ​h at room temperature, enhanced chemiluminescence reagent (Epizyme, Shanghai) was applied to develop the signal, which was detected by a ChemiDoc CRS imaging system (Bio-Rad, USA) [30].

Histology assays

Tissue specimens from rat hind limbs and clinical samples were fixed in 4% paraformaldehyde for 48 ​h and washed with phosphate buffered saline, then dehydrated in a graded ethanol series, vitrified with dimethylbenzene, and inserted in paraffin. Paraffin sections (4 ​μm) were deparaffinized in xylene, hydrated with gradient ethanol, and stained with standard H&E, SOFG, or Masson staining procedures. To perform the immunohistochemical staining, tissue slides were deparaffinized and rehydrated followed by antigen retrieval, endogenous peroxidase blocking and serum sealing. Then, the slides were incubated with antibodies against UBA1(1:100, ab180125, Abcam Inc, Cambridge, UK), EIF3E (1:50, ab134958, Abcam Inc, Cambridge, UK), RPL27 (1:50, 14980-1-AP, Proteintech, Chicago, USA), RPL17 (1:300, 67223-1-Ig, Proteintech, Chicago, USA), RPS28(1:50, 14796-1-AP, Proteintech, Chicago, USA) at 4 ​°C overnight. The next day, all sections were taken out and washed with PBS several times. Then the biotinylated secondary antibody (Servicebio, Wuhan, China) was used for 1 ​h at room temperature followed by the reaction with diaminobenzidine (Servicebio, Wuhan, China) and hematoxylin to develop color. Histological scores were calculated from the results of staining using Image J 6.0 (Media Cybernetics Corporation, USA) software.

Statistical analysis

The proteomic data represented three biological replicates for the AT, PAT, and control groups for each time point examined (2, 7, 14, and 28 dpi). Data were analyzed using SPSS (version 22.0; SPSS, Inc., Chicago, IL, USA). Categorical variables were analyzed using the chi-squared test. The differentially expressed proteins with a fold-change ≥1.5, <0.67, and p values ​< ​0.05, were considered statistically significant.

Results

Characterization of the temporal proteome during tHO in rat model

The mechanism underlying the development and severity of tHO remains unclear. In this study, we firstly built a valid animal model of traumatic heterotopic ossification which can induce different volume of ectopic bone. 10 weeks after surgery, the micro CT scan found that Achilles tenotomy stably induced significantly more amount of ectopic bone than partial Achilles tenotomy [20] (Fig. 1b and c). The histological analysis also proved that more trabecular structure and marrow cavities were found in injured tissue from Achilles tenotomy group than partial Achilles tenotomy group (Fig. 1d–g).

Fig. 1

(a) The graphical abstract of this study. The proteomics data are representative of n ​= ​3 biological replicates for AT and PAT rat models and control group (b, c) micro-CT images of rats of AT and PAT group 10 weeks post-injury, with three-dimensional reconstruction and quantification showing that significantly more ectopic bone was generated in AT group than PAT group (n ​= ​6/group, ∗∗∗∗P ​< ​0.0001, ∗Calculated using Student's t-test) (d, e, f, g) HE stains, Masson stains and Safranin-O stains and histological score semiquantitative evaluation of injury tissues of AT and PAT group at 10 weeks post modeling. n ​= ​6, Scale bars: 250 ​μm, ∗∗∗∗P ​< ​0.0001, ∗Calculated using Student's t-test.

To investigate the mechanism of tHO progression, injury site tissues specimens were collected from the sham (control), AT, and PAT tHO rat models at 2, 7, 14, and 28 days post-injury (dpi), and were subjected to DIA mode proteome analysis. This proteome profile identified 3807 quantifiable proteins in the samples from the AT, PAT and control groups, with no missing data at all time points (2, 7, 14, and 28 dpi). A total of 3547 of these proteins were detected in over 50% of samples and were used in the following bioinformatics analysis. Then we performed 12 pairwise comparison analyses among AT, PAT, and control groups, which were considered significantly altered if they exceeded the thresholds set at fold-change ≥1.5, or <0.67, and p values ​< ​0.05. The upregulated and downregulated proteins of all samples, are shown in Fig. 2a, Table 1.

Fig. 2

Data-independent proteomics analysis map the global protein changes in AT and PAT rat models during tHO development. (a) Venn diagrams show overlap between DEPs in AT and PAT groups during the development of tHO. The numbers of DEPs include both up- and down-regulated proteins in injury site tissue compared to normal tissue. AT: Achilles tenotomy group, PAT: partial Achilles tenotomy group, C: control group. (b) Principle component analysis of the data-independent acquisition (DIA) proteomic analysis. The AT, PAT and control group samples are marked with blue circles, grey rectangles, and yellow triangles, respectively. (c) Protein expression changes in samples from AT and PAT rat models at 2, 7, 14, and 28 days post injury. Heat map of selected proteins (|log2FC|>3) analyzed by hierarchical clustering analysis. (d) Volcano plot of protein expression changes in AT and PAT groups at day 2 post injury. Cut-offs using p ​< ​0.05 and log2-transformed ratio of >0.5 were applied to determine significantly altered proteins. Red dots represent the upregulated proteins and blue dots represent the downregulated proteins. (e) Top 10 Gene ontology (GO)-terms of biological process identified in pairwise comparisons between AT and PAT rat model samples at day 2 post injury (f, g) Top 10 canonical pathways and network of identified pathways in pairwise comparisons between AT and PAT rat model samples at day 2 post injury using Ingenuity Pathway Analysis (IPA). ​“Ratio” ​refers to the ratio of the number of molecules that map into the pathway versus the total number of molecules that define the canonical pathway by the IPA knowledge base. Z-score indicated the activation (>0) or inhibition (<0) of pathways.

Table 1

The number of up- and downregulated proteins identified in 12 pairwise comparisons.

Comparisons	Upregulated	Downregulated	Sum
2DPAT_vs_C	167	300	467
7DPAT_vs_C	314	467	781
14DPAT_vs_C	294	434	728
28DPAT_vs_C	192	381	573
2DAT_vs_C	184	365	549
7DAT_vs_C	461	386	847
14DAT_vs_C	390	363	753
28DAT_vs_C	479	354	833
2DPAT_vs_2DAT	174	132	306
7DPAT_vs_7DAT	130	258	388
14DPAT_vs_14DAT	67	307	374
28DPAT_vs_28DAT	130	464	594

Principal component analysis and hierarchical clustering of quantified proteins showed differences in the proteomes of AT and PAT groups during the development of tHO (Fig. 2b–c, Fig. s1). In general, proteins of both the AT and PAT groups were significantly altered after tenotomy operation, compared with the control group. However, the expression pattern at 2 dpi in the PAT groups could also be separated from that of the AT and PAT groups at other time points, as demonstrated by their distance in the dendrogram and concordance in the heat map. From 7 to 28 days post injury, the expression profiles for the AT and PAT groups gradually became concordant. GO annotation and canonical pathway analysis of the DEPs from pairwise comparison between AT and PAT group samples at 2,7,14,28 dpi were performed, which indicated that cytoskeleton organization, energy metabolism, mitochondrial dysfunction, acute phase response signaling, actin cytoskeleton signaling, ERK/MAPK signaling, and FAK signaling were the main pathways that differed between AT and PAT groups at day 2 post injury (Fig. 2e–g, Fig. s2). Overall, our observations suggest that AT and PAT groups exhibit obvious differences during the early phases of heterotopic ossification, which may result in the difference in the size of ectopic bone between AT and PAT groups.

The protein expression profile clusters during ectopic ossification of injury site tissue from AT and PAT groups

The analysis of protein expression profiles from a sequential series of samples in the tHO formation could provide a better understanding of the process of development of heterotopic ossification and identify the signaling pathways and mechanisms that are differentially regulated. We characterized the temporal dynamics of the AT and PAT proteomes by examining the significantly altered proteins by fuzzy C-means clustering. The analysis revealed eight distinct patterns of protein expression, which were generally similar between AT and PAT inflicted rat models (Fig. 3a and b). The expression of these proteins was increased in clusters 5, 6, and 8, slightly raised or not significantly increased in cluster 7, or decreased in clusters 1, 2, and 3; slightly decreased or unapparently reduced in cluster 4. Although clusters, 5, 6, and 8 contained proteins with increased expression, the different groups were classified into the same clusters and different clusters corresponded to different trends and intensities of change. Clusters 5 and 8 encompassed proteins that were dramatically increased (2- to 4-fold change), whereas cluster 6 contained proteins with moderate increases (∼1.5-fold change), and cluster 7 had proteins that were slightly increased or not significantly increased (0–0.5-fold change). Meanwhile, a distinct trend and intensity of protein expression was observed between AT and PAT groups in the same cluster. In cluster 8, the increase in protein expression in the AT group (∼4-fold change) was higher than that in the PAT group (∼2.5-fold change). In cluster 5, protein expression level was off at approximately 2-fold change of increase after 2 dpi in the AT group, whereas it was not until 7 dpi that protein expression in the PAT group reached a plateau. A similar situation was found for the downregulated clusters, where clusters 1, 2, and AT cluster 3 included intensely decreased proteins (2- to 3.5-fold change), PAT cluster 3 included moderately decreased proteins (∼1 fold change), and cluster 4 slightly decreased or unapparently reduced proteins (0–0.6 fold change). In the same cluster, proteins of the AT group exhibited a higher fold change of decrease than those of the PAT group in cluster 1 (∼4.5-fold change for AT and ∼3-fold change for PAT) and cluster 2 (∼2.5-fold change for AT and ∼2-fold change for PAT), and proteins in the PAT group showed a longer decline than AT (steady after day 2 post injury) before level off at 7 dpi. Next, we conducted canonical pathway analysis of the proteins from each cluster (Fig. 3c). It is noteworthy that clusters with similar trends and magnitudes of expression were involved in specific pathways or functions across the AT and PAT groups. For instance, in cluster 2 which included proteins that were dramatically downregulated, the proteins were involved in “mitochondrial dysfunction,” “oxidative phosphorylation,” “sirtuin signaling pathway,” which are related to cell energy metabolism. In cluster 6, which included moderately upregulated proteins, proteins were enriched in “EIF2 signaling,” “mTOR signaling,” “actin cytoskeleton signaling”, and “leukocyte extravasation signaling,” which are associated with inflammatory responses as well as cell proliferation and differentiation. These observations indicate that heterotopic ossification is characterized by particular biological functions and signaling pathways, whereas the trend and extent of changes in several specific pathways may influence the severity of heterotopic ossification. Furthermore, we investigated the expression profile of proteins in several pathways related to aberrant differentiation and osteogenesis during ectopic ossification (Fig. 3d). The proteins related to focal adhesion kinase signaling, such as ITGB1, were enriched in cluster 6, whose expression profile was moderately increased. EIF4E and HK1, associated with HIF signaling, were enriched in clusters 5 and 6. The proteins involved in TGF-β signaling, such as RAP1A and RHOA, were also enriched in cluster 6.

Fig. 3

Comparison of temporal expression profiles between AT and PAT group during the development of tHO. (a,b) Temporal patterns of tHO development proteome from AT and PAT groups were determined by Fuzzy c-means clustering analysis. Membership scores indicate the degree of relationship between proteins and each cluster. The number of proteins belonging to each cluster is indicated in the table below. (c) Canonical pathway analysis ​of ​the proteins in each fuzzy c-means cluster using IPA. Bolded terms indicate pathways that were strongly enriched and/or strongly contrasted between AT and PAT group. (d) Table of proteins in canonical signaling pathways of tHO found to be differentially expressed between AT and PAT group.

Co-expression analysis and hub protein selection

To further explore the proteins associated with the size of tHO, co-expression network analysis was conducted using the WGCNA package in R. Hierarchical clustering together with Dynamic Tree Cut algorithm identified ten modules after merging the correlated modules with a threshold power of β ​= ​7 (Fig. 4a). In our study, group labels were used as pathological traits to reveal the difference in the severity of ectopic ossification between AT and PAT rat models. To characterize the phenotype-related modules, we calculated the correlations of pathological traits with each Eigenprotein. The magenta and turquoise modules were found to have positive correlations with the AT groups, while they had weak or even negative relationships with the PAT and control groups (Fig. 4b). The GO enrichment annotation and canonical pathway analysis revealed the pathways affected in the turquoise and magenta modules (Fig. 4c and d), including cell migration-related pathways such as actin cytoskeleton signaling, paxillin signaling, and RhoA signaling; signaling associated with energy metabolism and mitochondrial functions such as mitochondrial dysfunction, oxidative phosphorylation, sirtuin signaling pathway, Huntington's disease signaling; pathways such as remodeling of epithelial adherens junctions, FAK signaling, PPARα/RXRα activation, ERK/MAPK signaling, and EIF2 signaling, which are involved in mechanotransductive signaling and cell differentiation. To verified the data of DIA proteomics profile, fifteen DEPs from magenta and turquoise module were randomly chosen and screened based on peptide identification with a Skyline cut-off score >0.95, including SEC61B, PSMD6, RPL4, PDIA3, TLN1, TPM2, MYLPF, VPS35, ACTN2, HSPA5, ACTN3, MYL3, RPSL1, APRT, DYNC1H1. These proteins evaluated by PRM and results showed that they shared similar trends in expression to those observed in DIA, which proved the reliability of proteomic profiles of magenta and turquoise modules (Fig. 5a, Fig. s3).

Fig. 4

Weighted gene co-expression network analysis (WGCNA) of AT and PAT rat models' proteome. (a) Determination of the soft-threshold power (β ​= ​7) for WGCNA after comparing various soft-threshold powers with their fitted scale-free indices and analyzing of the mean connectivity adapted to various soft-threshold powers; Sample clustering diagram. Number ​+ ​D represents different development stages of tHO; Hierarchical clustering tree of the DEPs clusters from the AT and PAT groups based on a dissimilarity measure (1- TOM). Each colored row contains a group of highly connected proteins, and similar rows were merged into 10 colored modules. (b) Heatmap of the correlation between modules and pathological traits. The abscissa represents groups with different pathological traits, and the ordinate represents different modules. The corresponding correlation and p-value are shown in each cell (c, d) Top 10 Gene ontology (GO)-terms analysis, canonical pathway analysis and the network of identified pathways of magenta and turquoise modules using Ingenuity Pathway Analysis (IPA).

Fig. 5

PRM validation of proteomic profile, Protein–Protein Interaction (PPI) Network Construction and hub protein identification. (a) Representative results of verifying the expression of proteins selected from magenta and turquoise modules by parallel reaction monitoring (PRM). (b) Constructed protein–protein interaction (PPI) networks of the proteins in the magenta and turquoise modules by Cytoscape (3.8.2). (c) The hub proteins of magenta and turquoise modules identified by Cytoscape (3.8.2) with cytoHubba plugin. The ranked score of each protein were recorded in the tables.

In addition, we also constructed a protein–protein interactions (PPI) network of the expressed proteins in the magenta and turquoise modules. The PPI network was generated using the STRING database in Cytoscape version 3.8.2 (http://cytoscape.org/) (Figs. 5a), and 11 scoring methods including the newly developed algorithm Maximal Clique Centrality (MCC) were employed by use of cytoHubba plugin [31]. The top 10 proteins with high ranked score were identified (Fig. 5b). The top 10 proteins in magenta module were UBA1, EIF3EL1, RPL27, APRT, PRPF8, VPS35, RPL4, IARS, IQGAP1, AP2A2. In turquoise module, the top 10 proteins were RPS28, RPL13A, RPL17, RPLP0, LOC108348142, RPL26, RPS23, LOC102555453, RPS3, RPL18A. The expression abundance of these proteins in the proteome profile are listed as follow (Table 2).

Table 2

1 Selected hub Proteins in magenta module as Identified by Proteomics Analysis

2 Selected hub Proteins in turquoise

module as Identified by Proteomics Analysis.

accession no.	symbol	name	2 ​d at/c	2 ​d at/pat	7 ​d at/c	7 ​d at/pat	14 ​d ​at/c	14 ​d ​at/pat	28 ​d ​at/c	28 ​d ​at/pat
ENSRNOP00000033950	Uba1	ubiquitin-like modifier activating enzyme 1	0.996319	1.057786	1.073135	1.122171	1.239392	1.498498	1.207784	1.337503
ENSRNOP00000029790	EIF3El1	eukaryotic translation initiation factor 3, subunit E-like 1	1.316914	0.937133	1.587955	1.448881	1.851214	2.168231	1.710788	1.40617
ENSRNOP00000028060	Rpl27	ribosomal protein L27	1.145059	1.168121	1.594847	1.432375	1.603811	1.735296	1.366063	1.390397
ENSRNOP00000061449	Aprt	adenine phosphoribosyl transferase	0.906954	1.107429	1.192763	1.44788	1.441801	1.535757	1.34019	1.234416
ENSRNOP00000005016	Prpf8	pre-mRNA processing factor 8	1.382218	1.106188	1.66131	1.438628	1.739696	1.855086	1.470597	1.323707
ENSRNOP00000024020	Vps35	VPS35 retromer complex component	1.691717	1.403596	1.372667	1.280299	1.957108	1.942958	2.095432	1.27129
ENSRNOP00000013462	Rpl4	ribosomal protein L4	1.298475	1.071218	1.79609	1.522275	1.974425	1.891286	1.713565	1.503293
ENSRNOP00000019862	Iars	isoleucyl-tRNA synthetase	1.197335	0.960149	1.374458	1.23537	1.534706	1.626175	1.553057	1.281874
ENSRNOP00000018021	Iqgap1	IQ motif containing GTPase activating protein 1	1.262384	1.382507	1.463243	1.346282	1.563769	1.306482	1.952925	1.524503
ENSRNOP00000060992	Ap2a2	adaptor-related protein complex 2, alpha 2 subunit	1.294847	1.403621	1.26517	1.409355	1.731262	1.579578	1.904217	1.26655

Current result of nest case–control study and validation of selected proteins

The tissue samples used in this study came from subjects had been recruited between September 2020 and June 2021. In total, forty-five tissue samples were originally collected from forty-five recruited participants in the nest case control study. Considering the consistency and comparability between clinical specimens and tendon tissue sample of rat tenotomy model, thirty-six ligament tissue samples were used in this research after quality control. Mean age of forty-five patients was 44.4 years (range, 18–72) and contained 23 males and 22 females. 24 right and 21 left sides elbows were affected. Distal humerus fractures and isolated radial head fractures, which included 28 patients, were the most common fracture type. This figure is followed by olecranon fractures (7 patients) and coronoid fractures (5 patients). A terrible traid was diagnosed in 5 patients and fracture-dislocations or compound fracture in 9 patients (Table 3).

Table 3

Patient demographics.

Demographics	n ​= ​45
Mean age (range), years	18-72（44.4）
Gender	males ​= ​23,females ​= ​22
Injury site	right elbow ​= ​24,left elbow ​= ​21
Fracture type, n
Distal humeral fracture	18
Isolated radial head fracture	10
Isolated olecranon fracture	7
Coronoid fractures	5
Terrible triad injury	5
Combined fracture	7
Fracture-dislocation	2
Time to surgery, n
Within 24 ​h	0
2day–7day	44
>7day	1

HO developed in 16 (35.5%) and clinically relevant HO in 4 (8.8%) surgically treated elbows. All clinically relevant HO was found developed in patients diagnosed with distal humeral fracture, with Class 3 of the Hastings and Graham classification developed in 2 patients, and Class 2 of the Hastings and Graham classification occurred in another 2 patients. There are 12 patients who developed clinically irrelevant HO, including 3 patients with terrible traid, 4 with isolated radial head fracture, 2 patients with distal humerus fracture, and 3 patients with fracture-dislocations (Table 4). The common site of clinically irrelevant HO contains lateral and medial equally, and 4 patients with elbow ankylosis or limitation of ROM were blocked by the ectopic bone developed in the front of anterior joint capsule (Fig. 6a and b). Tissue samples collected from these participants were selected and pooled into clinically relevant HO, clinically irrelevant HO or HO negative (control) group based on the primary outcome from nest case–control study and used for qPCR, western blot and immunohistochemistry assay.

Table 4

Types of elbow fracture and incidence of traumatic heterotopic ossification.

Elbow fracture types	HO,16 (35.5%)	Clinically relevant HO,4 (8.8%) Class 3 and Class 2 of the Hastings and Graham classification	Clinically irrelevant HO,12 (26.7%) Class 1 of the Hastings and Graham classification
Distal humeral fracture	6 (33.3%)	4 (22.2%)	2 (11.1%)
Isolated radial head fracture	4 (40%)	0 (0)	4 (40%)
Terrible traid injury	3 (60%)	0 (0)	3 (60%)
Fracture-dislocations	3 (42.9%)	0 (0)	3 (42.9%)

Fig. 6

X-ray images of clinically relevant HO, clinically irrelevant HO and normal elbow of patients after being surgically treated.

Some biological functions, such as cell proliferation and migration, cytoskeleton organization, and inflammatory response, play an important role in heterotopic ossification. Therefore, several representative proteins, whose differential expressions among control groups, partial Achilles tenotomy groups and Achilles tenotomy groups and involved in the above biological functions, were chosen for validation by qPCR, western blot and immunohistochemistry assay.

The qPCR results showed that, compared to the tissue of HO negative patients (control), five screened proteins, UBA1, EIF3E, RPL27, RPS28, RPL17, were significantly upregulated in the tissue of clinically relevant HO patients and clinically irrelevant patients in transcriptional level. And the mRNA expression levels of five identified proteins in the tissue of clinically relevant HO patients were generally higher than that of clinically irrelevant patients (Fig. 7a). In the tissue of rat models, the mRNA expression levels of five identified proteins were significantly upregulated in AT groups and PAT groups compared to control group, and the transcriptional levels of these five proteins in AT group (AT2d, AT7d, AT14d, AT28d) were generally higher than that of PAT group (PAT2d, PAT7d, PAT14d, PAT28d) (Fig. 7b). Western blot analysis and immunohistochemistry staining were applied to further verified the expression of the selected proteins. The results of immunoblot showed that five proteins, including UBA1, EIF3E, RPL27, RPS28, RPL17, exhibited the higher expression in the samples of clinically relevant HO patients and clinically irrelevant HO patients versus control group; Among these proteins, UBA1, EIF3E, RPL27, RPS28 were upregulated in the tissue of clinically relevant HO patients compared to clinically irrelevant HO patients. Meanwhile, the higher expression of RPL17 was found in the tissue samples of clinically relevant HO patients and clinically irrelevant HO patients than control group, but there is no obvious difference between the samples of clinically relevant HO patients and clinically irrelevant HO patients (Fig. 7c, e). Compared to control group, the protein expression of five identified proteins showed significantly upregulation in the tissue of AT and PAT rat models, And the protein expression levels of UBA1, EIF3E, RPL27, RPS28, RPL17 were generally upregulated in AT groups than PAT groups while the difference was more significant in the day2, day7 post injury than day 14, day28 post injury (Fig. 7d,f). Furthermore, the IHC results indicated that the expression levels of five screened proteins were significantly upregulated in AT groups (AT2d, AT7d, AT14d, AT28d) and PAT groups (AT2d, AT7d, AT14d, AT28d) compared to the tissue of control group. 10 weeks after modeling, the expression of five identified proteins were still strong but confined to the marrow cavities of ectopic bone (Fig. 8a–e, Fig. 9f). The expression levels of UBA1, EIF3E, RPL27, RPS28, RPL17 were significantly higher in the tissue of clinically relevant HO patients compared with the clinically irrelevant HO patients and HO negative patients (Fig. 9a–e, g). Collectively, results from qPCR, Western blot and immunohistochemistry staining assays were generally consistent with proteomic analysis data of rat HO model. These results proved the significance and applicability of selected hub proteins in the aspect of predicting the severity of traumatic heterotopic ossification in clinical course.

Fig. 7

qPCR and immunoblot assay compared the differential expression of the UBA1, EIF3E, RPL27, RPS28, RPL17 in tissue samples of rat models and clinical specimens. The relative gene expression of UBA1, EIF3E, RPL27, RPS28, RPL17 was evaluated in (a) clinical specimens and (b) tissue samples of AT, PAT and control rat models using qRT-PCR (c, d) Western blots for UBA1, EIF3E, RPL27, RPS28, RPL17 on whole-tissue lysate from clinical specimens and control, AT, PAT rat models. ​β-Actin served as a loading control. The relative expression level of UBA1, EIF3E, RPL27, RPS28, RPL17 in (e) clinical specimens and (f) tissue samples of control, AT, PAT rat models quantified using ImageJ and normalized to b-actin is shown. All tests in triplicate, ∗ represent significant change in the expression level compared to control group, ∗P ​< ​0.05, ∗∗P ​< ​0.01.

Fig. 8

Representative IHC staining of (a–e) UBA1, EIF3E, RPL27, RPS28, RPL17 in tendon tissue of control group and injury tissue from AT and PAT groups 2days, 7days, 14days, 28days, 10weeks after modeling. Original magnification is 20x. Inserts are approximately 4x magnified images of the boxed area. Scale bars: 250 ​μm.

Fig. 9

Representative IHC staining of (a–e) UBA1, EIF3E, RPL27, RPS28, RPL17 in ligament tissue of clinically relevant HO patients, clinically irrelevant HO patients and HO negative patients. Semiquantitative analysis of immunohistochemical staining of UBA1, EIF3E, RPL27, RPS28, RPL17 expression in tissue samples of control, AT and PAT rat models (f) and clinical specimens (g). ImageJ software was applied to calculate the average optical density. Original magnification is 20x. Inserts are approximately 4x magnified images of the boxed area. Scale bars: 250 ​μm n ​= ​3/group. ∗ represent significant change in the expression level compared to control group, ∗P ​< ​0.05, ∗∗P ​< ​0.01.

Discussion

Traumatic heterotopic ossification (tHO) remains a common complication in patients with orthopedic trauma, burns, traumatic brain injury, or spinal cord injury. The occurrence rate and severity of tHO do not appear to match the severity of trauma. Furthermore，according to the Hastings classification, the patients developed with Class 2 and class 3 of tHO will form bony ankylosis or blocking the joint activity by the large bone spurs developed in the elbow joint. But the tiny bone spurs formed in the elbow of the patients with class 1 will not cause symptoms or limiting the joint activity. For patients who may only develop clinically irrelevant tHO, the risk of prophylactic administration might overweigh its benefit due to side effects, including gastrointestinal toxicity and delayed fracture union. Recent studies have revealed a series of signaling pathways and proteins that play key roles during the development of tHO, such as TGF-β and its downstream BMP/SMAD pathways, the mTOR/AKT signaling axis, and the hypoxia inducible factor-1alpha (HIF-1α) pathway. However, few studies focus on the mechanism affecting the different volume of tHO, which determining the classification of tHO in the clinical course. This study was just designed to reveal the proteomic profile of development of different size of ectopic bone in traumatic heterotopic ossification in an attempt to obtain new insight into the complex mechanisms underlying the specificity of traumatic heterotopic ossification progression.

Although no consensus has been reached on the gold-standard model of trauma-induced HO yet, the Achilles tenotomy is a generally recognized method to induce tHO in animal studies. There are mainly two reasons: both the HO following Achilles tenotomy in animal models and HO following orthopaedic trauma in patients were induced by acute trauma, and endochondral ossification is the common pathological feature of these two types of heterotopic ossification [[32], [33], [34]]. In previous experiments, we unexpectedly found that partial Achilles tenotomy (PAT) could be used as a valid method to stably induce less ectopic bone than AT in rat models, without any other pharmacological or transgenic approaches. Therefore, the prerequisite of this study was to establish an animal model of distinct volume of ectopic ossification. The equal number of rats were performed Achilles tenotomy or partial Achilles tenotomy. 10 weeks after surgery, results of microCT scan of the rats’ hind limbs demonstrate the Achilles tenotomy stably induce significantly more volume of ectopic bone than partial Achilles tenotomy. These results correspond to the results of histological analysis which indicated that more trabecular structure and marrow cavities were found in injured tissue from Achilles tenotomy group than partial Achilles tenotomy group. In this way, the desired animal model was successfully built for the ensuing studies. Although the compound burn injury or central nerve system injury can also lead to increased severity of heterotopic ossification, the confounding factor beyond the local environment may undermine the reliability of proteomic analysis of the local heterotopic ossification. Accordingly, this model seemed to represent a usable model for investigate the severity of heterotopic ossification.

To further acquire the dynamic proteome profile and explore the mechanism during the progression of tHO, we performed AT and PAT in rats models and tissue samples were collected from tendon of sham (control) group, and at injury site from the AT, and PAT groups on 2, 7, 14, and 28 day post injury.

A data-independent acquisition mass spectrometry analysis of was applied on above tissue samples. A total of 3547 proteins detected in over 50 percent of the samples in the proteome analysis were used for further differential expression analysis. Notably, there were more DEPs in the AT group than in the PAT group. Interestingly, although the expression pattern at 2 dpi in the PAT groups differed from that in the AT and PAT groups at other time points according to the dendrogram and heatmap, the number of DEPs between AT and PAT samples at 2 dpi was the least among pairwise comparisons between AT and PAT samples across all time points.

The GO enrichment analysis of DEPs between AT and PAT samples at 2 dpi determined that the biological processes mainly included cytoskeleton organization, acute inflammatory response, and cellular respiration, while the canonical pathway analysis revealed the pathways involved in oxidative phosphorylation, mitochondrial dysfunction, sirtuin signaling pathway, ERK/MAPK signaling, and other pathways such as FAK signaling and Rho family GTPases signaling. The increased abundance of RhoA and decreased abundance of FBN1 in proteomes indicated that the increase in TGFβ signaling contributed to the development of tHO [35]. During the early stage of HO, a variety of progenitor cells migrate to the inflammatory milieu and proliferate [36], which was proven by the increase in RhoGDI signaling and FAK signaling in our proteome. Rho GTPase, a modulator of mechanotransduction that exerts its effects via the TGFβ receptor, was also found post-injury, where it modulates the cell cytoskeleton via the Rho GTPase and ROCK pathways to enhance chondrogenic/osteogenic differentiation [37]. Previous studies have shown that the activation of FAK and YAP/TAZ signaling induced by mechanical stimuli via mechanotransduction changes mesenchymal stem cells (MSCs) fate and initiates HO [38]. Furthermore, the ECM interaction receptor pathways activated by increased fibronectin and collagen subtypes would contribute to chondrogenesis/osteogenesis [39]. In addition, the PI3K/AKT and ERK/MAPK signaling pathways were found to be enriched in the cluster with increasing protein expression. It was previously reported that significant activation of the PI3K/AKT signaling pathway increases osteogenesis by directly enhancing mitochondrial respiration, as well as through the endothelial–mesenchymal transition process [40,41], while the ERK/MAPK pathways play a significant role in regulating chondrocyte differentiation [42].

In cluster 1 and cluster 2, which include intensely or moderately decreased proteins, the canonical pathway analysis revealed that the related pathways of these clusters were oxidative phosphorylation, mitochondrial dysfunction, TCA cycle, and sirtuin signaling pathway. This result indicated that proteins in cluster 1 and cluster 2 were mainly associated with energy generation and mitochondrial metabolism. It is noteworthy that the expression pattern of proteins in the AT group exhibited a higher fold change of decrease than those in the PAT group in both clusters, and proteins from the PAT group showed a more moderate decline than those in the AT group (steady after 2 dpi) before level off at 7 dpi. Previous studies have shown that during the initial phase of HO, uncoupled aerobic respiration in transient brown adipocyte-like cells results in a high energy expenditure, as well as hypoxia in the aberrant wound repaired tissue [43,44]. Meanwhile, the functional proteins associated with metabolic pathways were inhibited, and the enzymes involved in oxidative phosphorylation and TCA cycle pathways were markedly downregulated in tHO. Increased ratios of ATP:ADP, acetyl-CoA:CoA, and NADH:NAD have been shown to inhibit citrate synthase, and residual citrate is used for bone matrix construction [45]. The levels of cellular reactive oxygen species, which are mainly produced by mitochondria, can regulate the differentiation of local MSCs, and decreased mitochondrial metabolism and ROS generation can enhance pathways such as Wnt signaling, MAPK signaling, and BMP signaling to induce MSC differentiation into osteocytes and chondrocytes [46,47].

While the traditional DEP analysis and expression pattern clustering provided plentiful information, we applied weighted co-expression network analysis (WGCNA) to identify the proteins that specifically correlated with the differences in bone mass between the AT and PAT groups. After applying group labels as pathological traits to represent the different bone masses, we identified a magenta and a turquoise module that were associated with a larger volume of ectopic bone (AT group), while having a weak or negative correlation with less bone mass (PAT group) during the aberrant post-injury repair process. This result indicates that key proteins within the above modules may act as potential markers of the severity of heterotopic ossification. We identified the proteins in the magenta and turquoise modules that were enriched in canonical pathways associated with cell migration and proliferation as well as several other categories, including pathways associated with EIF2 signaling, actin cytoskeleton signaling, integrin signaling, RhoGDI signaling, RhoA signaling, ILK signaling, and Huntington's disease signaling.

To identify the key proteins in these modules, a PPI network was constructed, and hub proteins were identified with Cytoscape software and cytoHubba plugin, respectively. Meanwhile, a cluster of early stage traumatic heterotopic ossification tissue samples from a nest case–control study contained clinically relevant tHO, clinically irrelevant tHO and control patients were used to validate the reliability and applicability of the selected proteins. Top scored proteins of magenta and turquoise modules were chosen for validation. The WB assay confirmed that five representative proteins, including UBA1, EIF3E, RPL17, RPL27, RPS28, were differentially expressed among clinically relevant tHO patients, clinically irrelevant tHO patients and control patients.

UBA1, the Ubiquitin-activating enzyme E1, belongs to the ubiquitin-activating E1 family of enzymes and serves as a key regulatory role in the ubiquitin-proteasome system (UPS) and selective autophagy. The ubiquitin-proteasome system is known to have effects in functions such as cell cycle progression, DNA damage repair, transcription, translation, vesicle transport and apoptosis [48]. UBA1 was found involved in the activation of NF-κB signaling through Lys63 (K63)-linked polyubiquitin chains and led to inflammation and cell proliferation in the restenosis and acute myelocytic leukemia [49,50]. NF-κB/MAPK signaling can activate ACVR1 and induce heterotopic ossification in FOP patients [51]. UBA1 mutation was also link with deficient bone development in some neural muscular diseases [52]. Furthermore, UBA1 was also found involve in VEGFA signaling and may promote the endothelial cell proliferation and angiogenesis, which is essential for the bone formation [53,54]. Considering these findings, we would suppose that the up-regulated expression of UBA1, as detected by this study, may play a vital role in the regulation of the abnormal ectopic bone formation.

EIF3E, which is the “e” subunit of the translation initiation factor eIF3 (Eukaryotic initiation factors), plays a key role in regulating total protein synthesis, controlling cell growth, size and proliferation, and involving in tumorigenesis [55,56]. Recent studies found that EIF3E is an important regulator factors of HIF signaling. It has been proven that the HIF pathway and its downstream signaling such as vascular endothelial growth factor (VEGF) and bone morphogenetic proteins (BMPs) were activated by hypoxic microenvironment induced by the HO lesion. Hif1α expression coincides with the expression of chondrogenesis factor such as Sox9 [(sex determining region Y)-box 9] and pharmacologic inhibition of Hif1α can significantly decreased extraskeletal bone formation [57]. Importantly, EIF3E can promote the vascular remodeling and development via Int6/eIF3e-HIFs pathway [58]. Based on this studies, it could probably be suggested that the higher expression of EIF3E, as revealed in this study, may promote the expression of HIF signaling and increase the progression of heterotopic ossification [59].

Three ribosomal protein, L17, L27, S28, were identified as potential biomarkers of traumatic heterotopic ossification, too. Ribosomal proteins (RPs) involved in many key biosynthesis progresses such as ribosome biogenesis, peptide bond formation, and protein synthesis rate. During the wound healing and bone development, high volume protein synthesis requiring ribosomal activity is essential for growth and bone matrix production. Recent studies proved that disruption of ribosome-related genes contribute to the pathogenesis ankylosing spondylitis [60]. Furthermore, previous studies proved that RPL17 was correlated to angiogenesis which is essential for the HO formation [61]. Thus the increased expression level of ribosomal proteins may indicate the initiation and formation of aberrant ectopic ossification. The previous bioinformatics analysis about the late stage heterotopic ossification and proteomic analysis of HO ​+ ​tissue samples also identified ribosomal proteins such as RPL17, RPS18 exhibit a higher expression and act as key proteins of the molecular mechanism of heterotopic ossification [62,63].

However, there are still some limitations in our present study. Firstly, the sample size of the nest case–control study is small, thus we need larger group and more specimens to further validate the results of this manuscript. While we evaluated tissue-level expression of five screened proteins that were associated with HO prognosis, the expression levels of these potential biomarkers in the serum and other body fluids needed to be explored. Secondly, five proteins and several signaling pathways correlated with the progression of HO have been uncovered in this study, but their detail functions in the HO formation and interactions with classical signaling pathways of HO still needed to be further investigated in vivo and in vitro experiments.

Conclusion

In conclusion, we hope that the new established animal model and proteomic database generated in this study could serve as a solid foundation for the comprehensive investigation and the mechanism elucidation of the progression of traumatic heterotopic ossification. And the identified 5 proteins (UBA1, EIF3E, RPL17, RPL27, RPS28) were associated with the severity of traumatic heterotopic ossification. It is hypothesized that the results of our study may provide novel insight for the prevention, diagnosis, and treatment of traumatic heterotopic ossification.

Authors' contributions

ZYW and SG planned and executed experiments, ZYW wrote the manuscript, SG analyzed results and composed figures, HWW and CZ critically reviewed the manuscript, ZY and JRY recruited participants and collected specimens. BZ oversaw the project, critically reviewed and edited the manuscript.

Funding

This work was supported by the National Natural Science Foundation of China (grant numbers 81772426, 81974330); the Program of Science and Technology Commission of Shanghai (grant number 15dz2353800); the Interdisciplinary Program of (grant number YG2019QNA23); and Affiliated Sixth People's Hospital's own funds (LYZY-0262).

Declaration of competing interest

All authors declare they have no competing interests.