Gene markers of fracture healing in early stage and the regulatory mechanism during the process using microarray analysis.

BACKGROUND: The aim of this study was to explore crucial markers and uncover the regulatory mechanisms of fracture healing in the early stage.
METHODS: Gene expression profile of GSE45156 was downloaded, in which 3 fractured samples and 3 unfractured samples were used in our present study. Based on the threshold value, differentially expressed genes (DEGs) were selected between two kinds of samples using limma package in R. Enrichment analysis of these DEGs was performed by DAVID software. Furthermore, protein-protein interaction (PPI) network was established integrating information in STRING database, and visualized by Cytoscape software.
RESULTS: We identified a set of 960 DEGs including 509 up-regulated and 451 downregulated genes. Biological processes involving RNA splicing and cell cycle were significantly enriched for the up-regulated genes such as Snrpd2, Eftud2, Plk1 and Bub1b, whereas skeletal system development and bone development processes were predominant for down-regulated genes like Ubc. In the constructed PPI network, all the five genes were the predominant nodes, of which Snrpd2 was linked to Eftud2, while Bub1b was to interact with Plk1.
CONCLUSION: Five candidate genes crucial for indicating the process of fracture in early stage were identified. Eftud2, Snrpd2, Bub1b and Plk1 might function through the involvement of cell-cycle-related BP, while Ubc might influence the protein degradation during bone development. However, more experimental validations are needed to confirm these results.

Introduction

Fracture is a serious problem worldwide and nearly 10% of all the fractures will develop into nonunion. The risk factors might be the inadequate functional stability, excessive soft tissue damage and infection. Fracture healing is commonly considered as a complex process that involves the coordinated action of multiple proteins and genes to repair the structural integrity without the formation of scar. Although standard clinical management of fracture healing exists, patients with extensive trauma would take a long duration to union, which might cause substantial morbidity and poor outcomes.

Initial inflammatory response, generation of soft and hard callus, bone union and bone remodeling are the primary events during the repair. Numerous cell types are involved during the fracture healing including hematopoietic cells, which are essential for hematoma formation and the inflammatory response, as well as endothelial cells.5, 6 Mesenchymal cells are another crucial cell type whose proliferation and differentiation contribute to the bone generation in the process of repair. The pluripotent mesenchymal progenitors are able to differentiate to osteoblasts or chondrocytes. In severe fractures with substantial soft tissue injury, mesenchymal progenitors can also facilitate to fracture healing. Due to the significant role of mesenchymal cells in fracture healing, multiple studies have been investigated to explore markers of this kind of cell. Among the established markers, alpha smooth muscle actin (αSMA) was identified as a marker of the mesenchymal progenitor cell facilitating to osteochondral events in fracture healing by a recent study. Moreover, the study labeled the αSMA-expressing cell populations which were involved in fracture repair, and further examined their expression profiles in the early stage of fracture healing. Though samples of 2 days and 6 days after fracture were both contained, their emphasis was obviously on 6 days, and consequently they found the increased genes relating to matrix production and remodeling as well as the decreased genes involving in Notch signaling process. However, crucial genes in early stage (2 days) after fracture as well as the potential regulatory interactions between them were scarce.

Therefore, we re-analyzed the microarray data of GSE45156 deposited by Matthews et al, and just utilized the samples of 2 days after fracture, to investigate the gene alterations in the initial of the fracture and the corresponding repair process. Bioinformatics methods such as differential analysis, functional enrichment analysis and protein–protein interaction (PPI) network were performed in our study, to explore the crucial gene markers and uncover the underlying regulatory mechanisms during the early stage of fracture healing.

Methods

The study did not involve any human or animals, so the ethical approval was not required.

Microarray data resource

Gene expression profile of GSE45156, which was deposited by Matthews et al in the public Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo) database, was utilized in our study. The expression profiles of six samples consisting of 3 Fractured Day 2 αSMA-labeled cells (fractured samples) and 3 Unfractured Day 2 αSMA-labeled cells (unfractured samples, control samples) were used in this study. The platform for gene expression profile analysis was GPL6885 Illumina MouseRef-8 v2.0 expression beadchip (Illumina, San Diego, California, USA).

Preprocessing of the expression profile data

Two annotation files (GSE45156_family.soft and GSE45156_series_matrix.txt) were downloaded from the platform, and the normalized gene expression data were extracted. If more than one probe corresponded to one single gene, the averaged probe value was considered as the final expression value of the specific gene.

Differential analysis

Differentially expressed genes (DEGs) between fractured sample and control sample were selected using limma (Linear Models for Microarray Analysis, http://www.bioconductor.org/packages/release/bioc/html/limma.html) package in R. The thresholds for DEG selection were false discovery rate (FDR) < 0.05 and |log2 (fold change)| > 0.5.

Function and pathway enrichment analysis

The Gene Ontology (GO, http://www.geneontology.org/) project provides ontologies to attributes of gene function in three domains, including molecular function, biological process and cellular component. Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.genome.jp/kegg/pathway.html) is a service of establishing manually curated pathway maps that contain current knowledge on biological networks. The DAVID (Database for Annotation, Visualization and Integration Discovery, http://david.abcc.Ncifcrf.gov/) online tool was recruited to perform the GO and KEGG pathway enrichment analysis, to ascertain the functions and pathways that might be disturbed by the identified DEGs, with the selection criterion of FDR < 0.1.

Construction of PPI network

Proteins often work together to exert a specific function during a biological process. To further explore the potential interactions among the DEGs from protein level, the STRTING (Search Tool for the Retrieval of Interacting Genes, http://string-db.org/) database was searched to screen the interplayed proteins of the protein production of the DEGs. The interactions with the confidence score >0.7 that have been undergone experimental validation, text mining, co-expressed analysis and recorded in relevant databases were extracted to construct the PPI network, visualizing by Cytoscape software (http://cytoscape.org/). In the PPI network, a protein (the protein product of a DEG) serves as a node and the degree of a node reflects the number of the interplayed proteins with the specific node. Hub nodes were deemed as the nodes with high degrees.

Results

DEGs between fractured and control samples and their altered functions and pathways

Based on the aforementioned selection criteria, a subset of 960 DEGs (509 were up-regulated and 451 were down-regulated) were identified. As presented in Table 1, the up-regulated DEGs such as Snrpd2 and Eftud2 were significantly enriched in RNA splicing (GO: 0008380) and mRNA processing (GO: 0006397) BP, while Plk1 and Bub1b were predominantly associated with cell cycle relating process (GO: 0007049, GO: 0022403 and GO: 0022402). Notably, Snrpd2 and Eftud2 were also significantly enriched in the pathway of Spliceosome (mmu03040) (Table 2).

Table 1

Significant enriched biological functions of the differentially expressed genes (DEGs).

Term	Count	FDR	Genes
Up-regulated DEGs
GO: 0006412 ∼ translation	33	5.34E-09	Tufm, Rpl18, Rps27a, Mrps16, Ptrh1 et al.
GO: 0008380 ∼ RNA splicing	25	8.35E-08	Raly, Snrpd3, Snrpd2, Eftud2, Sfrs1 et al.
GO: 0006397 ∼ mRNA processing	26	3.98E-06	Raly, Snrpd3, Snrpd2, Sf3b5, Eftud2 et al.
GO: 0006396 ∼ RNA processing	34	4.46E-06	Raly, Elac2, Pus1, Snrpd3, Snrpd2 et al.
GO: 0016071 ∼ mRNA metabolic process	26	7.02E-05	Raly, Snrpd3, Snrpd2, Eftud2, Sfrs9 et al.
GO: 0006259 ∼ DNA metabolic process	31	9.13E-05	Rad23b, Clspn, Mre11a, Dbf4, Hus1 et al.
GO: 0007049 ∼ cell cycle	38	1.67E-04	E2f1, Plk1, Bub1b, E2f3, Tsg101 et al.
GO: 0006260 ∼ DNA replication	17	8.63E-04	Ssrp1, Ssbp1, Dbf4, Brca2, Mcm10 et al.
GO: 0022403 ∼ cell cycle phase	24	4.54E-03	Plk1, Bub1b, Ran, Dbf4, Mre11a et al.
GO: 0022402 ∼ cell cycle process	26	9.47E-03	Plk1, Bub1b, Tsg101, Mre11a, Dbf4 et al.
GO: 0043039 ∼ tRNA aminoacylation	9	1.6E-02	Iars, Nars, Rars, Farsb, Khsrp, Lars, Gars, Vars, Kars
Down-regulated DEGs
GO: 0001568 ∼ blood vessel development	27	8.16E-09	Mef2c, Emcn, Edn1, Tiparp, Foxo1 et al.
GO: 0001944 ∼ vasculature development	27	1.43E-08	Mef2c, Emcn, Edn1, Tiparp, Foxo1 et al.
GO: 0048514 ∼ blood vessel morphogenesis	22	1.09E-06	Mef2c, Klf5, Bmp4, Emcn, Flt1 et al.
GO: 0001525 ∼ angiogenesis	17	2.07E-05	Klf5, Bmp4, Emcn, Flt1, Egfl7 et al.
GO: 0001501 ∼ skeletal system development	21	2.69E-03	Ubc, Mef2c, Bmp4, Tiparp, Edn1 et al.
GO: 0051270 ∼ regulation of cell motion	12	1.92E-02	Selp, Pdgfb, Egfl7, Sema3f, Pecam1 et al.
GO: 0007167 ∼ enzyme linked receptor protein signaling pathway	19	2.12E-02	Bmp4, Flt1, Pdgfb, Ltbp4, Tiparp et al.
GO: 0042127 ∼ regulation of cell proliferation	28	2.54E-02	Pdgfb, Edn1, Clu, Pparg, Foxo1 et al.
GO: 0030334 ∼ regulation of cell migration	11	2.95E-02	Selp, Pdgfb, Egfl7, Pecam1, Gab1, Robo4, Thbs1, Sele, Pik3r1, Vcl, Cxcl10
GO: 0007389 ∼ pattern specification process	19	3.63E-02	Bmp4, Irx3, Mfng, Flt1, Edn1
GO: 0060348 ∼ bone development	12	4.87E-02	Bmp4, Mef2c, Hdac5, Pthlh, Hoxb4, Ubc, Foxc1, Gnas, Npr2, Igf2, Igfbp3, Bmp6

GO: Gene Ontology; Count: Numbers of the genes enriched in a specific GO function term; FRD: false discovery rate.

Table 2

Dysfunctional pathways of the up-regulated differentially expressed genes.

Term	Count	FDR	Genes
mmu03040: spliceosome	18	1.10E-04	Snrpd3, Eftud2, Snrpd2, Sfrs1, Sf3b5,Sf3a1, Sf3a3, Ddx46, Sfrs9, Usp39 et al.
mmu03030: DNA replication	8	6.59E-02	Rfc5, Rpa1, Rfc3, Ssbp1, Rfc2, Pold1, Mcm4, Fen1

Count: Numbers of the genes enriched in a specific pathway term; FRD: false discovery rate.

On the other hand, the down-regulated DEGs such as Ubc were significantly enriched in skeletal system development (GO: 0001501) and bone development (GO: 0060348) BP (Table 1). Unfortunately, based on the pre-defined criteria, no significant pathways were enriched for the down-regulated DEGs.

PPI network of the identified DEGs

A PPI network of the DEGs was established, in which Mrto4 (degree = 43), Snrpd2 (degree = 35), Eftud2 (degree = 38), Bub1b (degree = 36), Plk1 (degree = 40), Polr1a (degree = 31), Ubc (degree = 79) and Rps27a (degree = 42) were the hub nodes (Fig. 1). Notably, Snrpd2 was linked to Eftud2, whereas Bub1b was interacted with Plk1.

Fig. 1

Protein–protein interaction network of the differentially expressed genes. White round denotes protein products of the up-regulated genes and the gray represents protein products of the down-regulated genes. Lines between two proteins indicate their predicted interactions.

Discussion

Large amount of cellular activities are involved in the complex process of fracture healing. Through re-analyzing the expression profile data, GSE45156, we identified a set of 960 DEGs including 509 up-regulated genes and 451 down-regulated genes, which might be crucial for fracture healing. BPs such as RNA splicing and cell cycle were significantly enriched for the up-regulated DEGs such as Snrpd2, Eftud2, Plk1 and Bub1b, whereas skeletal system development and bone development processes were the predominant ones for the down-regulated DEGs like Ubc. Notably, all the five genes were the remarkable nodes in the PPI network, in which Snrpd2 was predicted to link to Eftud2, while Bub1b was to interact with Plk1.

RNA splicing, which is known as a crucial regulatory mechanism for gene expression, is a process mediated by large protein-RNA complexes. Different RNA splicing is proposed as the causative factor for clinical behavior of a tumor. The Eftud2 encoded protein is a component of the splicesome involved in the process of RNA splicing. Knockdown of the protein results in the increase of the apoptotic cells and mutation of such protein is tightly correlated with genetic diseases such as mandibulofacial dysostosis. Snrpd2 is a member of the small nuclear ribonucleoprotein (SNR) core protein family that accounts for the pre-mRNA splicing and SNR biogenesis. Numerous SNR proteins such as Snrpd2 were related to the inositide-dependent phospholipase C β1 (PI-PLCβ1), which mediates cell division and cell proliferation in mouse erythroid leukemia cells. All these corroborate that alternative splicing exerts significant roles in the regulation of cellular differentiation and apoptosis. Additionally, multiple variants of the splice factors including Eftud2 and Snrpd2 have been discovered in the inflammatory breast cancer cell lines, suggesting their potential roles in the inflammation. In our present study, up-regulated Eftud2 and Snrpd2 were both enriched in RNA splicing BP, consistent with the functional enrichment analysis in fatal familial insomnia in brain tissues, further supporting the involvement of the two genes in this process. At present, regulatory relationships between these two genes have not been reported, however they were interacted in the PPI network of our study, providing a predicted regulation between them. Given that the process of fracture healing is highly associated with cell-cycle-related BP and inflammatory response,25, 26 it might be speculated that the involvement of Eftud2 and Snrpd2 in RNA splicing might be related to cell cycle and inflammation in the early stage of fracture healing.

Regulation of cell cycle is often involved during the process of fracture healing. TP508, a synthetic 23-amino acid peptide, is ascertained to promote the fracture healing via the expression modulation of proteins that are mainly related to cell cycle and cellular proliferation functions. In eukaryotic cells, PLK1 (Polo-Like Kinase 1) acts as a crucial mediator of cell division for it could trigger G2/M transition.28, 29 Expression of cell-cycle-related genes such as PLK1 is frequently measured to determine whether alteration of other genes would affect the proliferation of bone marrow.30, 31 The kinase encoded by Bub1b involves in spindle checkpoint function and plays significant roles in the control of appropriate chromosome segregation. Expression of Bub1b shows a cell-cycle-related pattern, and it is undetectable in G1 but abundant in G2/M in synchronized cells. Based on the pathway analysis, DEGs relating to osteoblasts such as Bub1b were pronouncedly enriched in the cell cycle pathway, in accordance with our findings. It is reported that in response to Ochratoxin A, one potent renal carcinogen in rodents, gene expressions of Bub1b and Plk1 were both altered. Notably, proteins of Bub1b and Plk1 were indirectly linked in the Bub1b subnetworks of glioblastoma tumors. Though mechanisms of tumor and fracture healing were distinct, this finding may also be implicit for the regulatory interactions between genes of Bub1b and Plk1 because expressions of extensive genes correlating to cell proliferation are changed during fracture healing. Combining the results that Bub1b and Plk1 were both remarkably enriched in cell cycle BP with that they were interplayed in the PPI network, it might be inferred these two genes might exert important regulatory roles in the process of fracture healing in the early stage via modulating the cell-cycle-related process.

Ubc is an ubiquitin gene of ubiquitin C. Dysregulation of Ubc expression is discovered in early fracture tissue in domestic sheep by EST sequencing. Up-regulation of several proteins including versican and chondroitin sulfate proteoglycan have been proposed as the candidate markers for osteoblasts differentiation. Considering that Ubc is tightly associated with the protein degradation, the down-regulated Ubc might contribute to the fracture healing for the decrease of protein degradation in the bone development process.

Despite the comprehensive analysis, several limitations should be discussed in this study. First, all the analyses were based on the microarray data downloaded from the public GEO database. The sample size was relatively small and only 6 samples were recruited in our study. Moreover, though the pooled periosteum or callus tissue from several sex-matched mice would dispel the heterogeneity across individuals, a portion of vital DEGs might also be undetectable due to the pooled samples. Last but not the least, all the predictive results in our present study need to be further confirmed through experimental validations, which will be considered in our future studies. Nevertheless, the study is also of great value in identification of crucial genes and revelation of mechanisms for fracture.

Conclusions

In conclusion, five candidate genes crucial for indicating the process of fracture in early stage were identified including Eftud2, Snrpd2, Bub1b, Plk1 and Ubc. Eftud2, Snrpd2, Bub1b and Plk1 might function through the disturbance of cell-cycle-related BP, while Ubc might influence the protein degradation during bone development. However, more experimental validations are warranted.

Competing interests

The authors declare that they have no competing interests.