This article is mentioned in:

Introduction

Osteoporotic fracture is a common event in the elderly, resulting in substantial mortality, and the mortality rate of hip fracture for 6 months is ~10–20% (1). The prevalence of osteoporotic fractures, hip fractures in particular, is increasing in many regions of the world (2). Current therapies focus on the prevention and treatment of osteoporotic fractures; however, this may easily lead to complications, thus it remains a worldwide public health concern. Therefore, a greater understanding of the underlying molecular mechanisms of fracture healing in the osteoporotic bone is required, as well as identifying candidate biomarkers for osteoporotic fracture therapies.

Over the past few years, a number of remarkable achievements have been made in the genetic study of fracture healing in osteoporosis. One such study demonstrated that transgenesis of bone morphogenetic protein-2 promotes fracture healing in osteoporosis by inducing increased callus density and a larger cross-sectional callus area (3). During remodeling of fractured bone, parathyroid hormone (PTH) promotes the formation of osteoclasts to restore the mechanical strength and structure of bones, and polymorphisms in genes encoding PTH influence the genetic regulation of bone mineral density (4). Low density lipoprotein receptor-related protein 5 (LRP5) serves a significant functional role in skeletal homeostasis, and mutations in LRP5 induce a variety of bone density-associated diseases (5). Lrp5 deficiency results in decreased osteoblast proliferation and function, which induces a low bone mass phenotype (6). Kringle containing transmembrane protein 2 (KREMEN2), also known as KRM2, is a high-affinity transmembrane receptor of dickkopf homolog 1, and is thought to be a regulator of bone remodeling (7). It has been demonstrated that Krm2−/− mice develop a high bone mass phenotype and overexpression of Krm2 in type I collagen (Col1a1)-Krm2 transgenic mice induces severe osteoporosis with decreased levels of osteoblasts and elevated osteoclast differentiation (8). Using a model of fracture healing in Col1a1-Krm2 transgenic mice and Lrp5−/− mice, a previous study revealed that fracture healing is greatly damaged in Col1a1-Krm2 transgenic mice and Lrp5−/− mice; however, the Col1a1-Krm2 mice were more severely impaired than Lrp5−/− mice (9). In addition, this previous study identified a set of differentially expressed genes (DEGs) in the two mouse models using microarray analysis (9). However, DEG interactions and functions require further investigation in order to provide a more comprehensive understanding of the effect of osteoporosis on fracture healing.

In order to investigate the interactions and functions of DEGs in Col1a1-Krm2 transgenic mice and Lrp5−/− mice further, the microarray data obtained by Liedert et al (9) were analyzed in the present study. Following identification of DEGs, enrichment analysis was performed. In addition, protein-protein interactions (PPIs) of DEGs and hub genes in the PPI network were analyzed. Furthermore, coexpression associations between hub genes and additional DEGs were examined. These results may contribute to a greater understanding of the effect of osteoporosis on fracture healing, and provide novel information that facilitates the development of future clinical therapies for osteoporotic fractures.

Materials and methods

Affymetrix microarray data

The raw gene expression profile dataset GSE51686 (9) was obtained from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/). The data was generated by the (Mouse430_2) Affymetrix Mouse Genome 430 2.0 Array platform (GEO accession, GPL1261; Affymetrix, Inc., Santa Clara, CA, USA). This dataset contained 9 fracture callus tissue samples obtained from the femora at 10 days following osteotomy in the wild-type (WT) mice (n=3), Col1a1-Krm2 transgenic mice with severe osteoporosis (n=3), and Lrp5−/− mice with low bone mass (n=3), respectively. All mice were female and 26 weeks of age.

The CEL and probe annotation files for this dataset were downloaded. The raw expression data were preprocessed by background correction, quantile normalization and probe summarization using the robust microarray analysis algorithm in the affy package (version 3.3.2) (10) of Bioconductor (version 3.4; http://www.bioconductor.org/). Subsequently, the org.Hs.eg.db (version 3.4.0) (11) and illuminaHumanv3.db (version 1.26.0) (12) packages of Bioconductor were used to translate probe identifications (IDs) to gene symbols. If one gene symbol was matched by multiple probe IDs, the mean expression value was selected as the expression level of this gene.

Identification of DEGs

DEGs in Col1a1-Krm2 mice and Lrp5−/− mice compared with the WT controls were identified using the linear models for microarray data (LIMMA) package (version 3.30.3; http://www.bioconductor.org/packages/release/bioc/html/limma.html) (13), which is a commonly used tool for the identification of DEGs. The P-value for each gene was calculated using the unpaired t-test in LIMMA, which was then adjusted for the false discovery rate (FDR) using the Benjamini-Hochberg method (14). Only the genes with FDR values <0.05 and log2 fold change values ≥0.5 were selected as DEGs.

The Venny online tool (version 2.0; http://bioinfogp.cnb.csic.es/tools/venny/index.html) (15) was utilized to construct Venn diagrams for the upregulated and downregulated genes identified between the Col1a1-Krm2 vs. WT and Lrp5−/− vs. WT groups.

Enrichment analysis of DEGs

Functional Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses of upregulated and downregulated genes were performed using the Database for Annotation Visualization and Integrated Discovery (version 6.8; http://david.abcc.ncifcrf.gov/) database (16). The P-value was calculated using the modified Fishers exact test, and P<0.05 was considered to indicate a statistically significant difference. A gene count in each term ≥2 was set as the cut-off criteria. Additional parameters were set to the default values.

Construction of PPI networks

PPIs of DEGs were obtained from the Search Tool for the Retrieval of Interacting Genes database (version 10.0; http://string-db.org/), which integrates a variety of known and predicted protein associations (17). The combined score for each PPI was calculated, and a score of >0.4 was set as the cut-off criterion. Additional parameters were set to the default values. The PPI network was visualized using the Cytoscape software (version 3.4.0; http://cytoscape.org/), which is an open access software for visualizing biomolecular networks (18). In the network, ‘node’ represents a gene or protein, and ‘line’ represents an interaction between the two nodes. The degree of each node (number of interactions with other proteins) is equal to the number of nodes that interacted with this node.

Analysis of hub genes in the PPI network

Hub genes refer to the relatively key genes in the network. Hub genes were identified using three centricity methods in the PPI network, including the degree centrality (19), betweenness centrality (20) and subgraph centrality methods (21). The scores obtained from the degree, betweenness and subgraph methods were calculated using the CytoNCA plug-in (version 2.1.6) (22) in Cytoscape. High scores for the degree, betweenness and subgraph methods indicated that the nodes were more significant in the network. Hierarchical clustering of hub genes with higher scores was performed using the pvclust R package (version 1.3–2) (23).

Coexpression associations of hub genes with DEGs

The Pearsons correlation coefficient (PCC) method (24) was used to identify the coexpression associations of hub genes with other DEGs. Only coexpression associations with PCC values of >0.9 were selected for analysis. A PCC value of >0 indicated that the two genes were positively correlated, and a PCC value of <0 indicated that the two genes were negatively correlated.

Results

Statistical analysis

Based on the cut-off criteria, a total of 841 DEGs (335 upregulated and 506 downregulated) and 50 DEGs (16 upregulated and 34 downregulated) were identified in the Col1a1-Krm2 and Lrp5−/− mice when compared with WT mice, respectively. When compared with WT mice, 12 of these genes were upregulated and 25 were downregulated in the Col1a1-Krm2 and Lrp5−/− mice (Fig. 1).

Figure 1. Venn diagrams of differentially expressed genes between Col1a1-Krm2 vs. WT and Lrp5−/− vs. WT groups. (A) Venn diagram of upregulated genes. (B) Venn diagram of downregulated genes. Col1a1, type I collagen; Krm2, kringle containing transmembrane protein 2; Lrp5, low density lipoprotein receptor-related protein 5; Col1a1-Krm2, Col1a1-Krm2 transgenic mice; Lrp5−/−, mice deficient in Lrp5; WT, wild-type mice.

DEG function

To further reveal gene function in the two groups, GO and KEGG pathway enrichment analyses were performed. In the Col1a1-Krm2 vs. WT group, the upregulated genes were primarily associated with hemopoiesis (GO: 0030097), hemopoietic or lymphoid organ (GO: 0048534) and immune system development (GO: 0002520), as well as pathways associated with primary immunodeficiency (mmu05340) and nitrogen metabolism (mmu00910) (Fig. 2A). The downregulated genes were significantly associated with cell adhesion (GO: 0007155) and regulation of the smoothened signaling pathway, as well as the hedgehog signaling pathway (mmu04340) and cell adhesion molecules (mmu04514) (Fig. 2A).

Figure 2. GO and KEGG pathway terms enriched by differentially expressed genes between (A) Col1a1-Krm2 transgenic mice and wild-type mice and (B) Lrp5−/− mice and wild-type mice. ‘Up’ represents upregulated genes and ‘Down’ represents downregulated genes. Red bars indicate the P-values, and blue bars indicate the gene count. GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; Col1a1, type I collagen; Krm2, kringle containing transmembrane protein 2; Lrp5, low density lipoprotein receptor-related protein 5; BP, biological process; CC, cellular component; MF, molecular function.

In the Lrp5−/− vs. WT group, the upregulated genes were implicated in muscle contraction (GO: 0006936) and muscle system process (GO: 0003012) (Fig. 2B). The downregulated genes were markedly associated with the regulation of transcription (GO: 0006355) and RNA metabolic processes (GO: 0051252) (Fig. 2B). No significant pathways were enriched by the upregulated genes.

Analysis of PPI network

In order to determine interactions between DEGs, a PPI network was constructed. The network was composed of 551 nodes and 1,608 PPIs (Fig. 3). Based on the centricity methods, the top 40 nodes with the highest scores in the PPI network were selected as hub genes for further analysis, including 2–5-oligoadneylate synthase-like protein 2 (Oasl2), thrombospondin 2 (Thbs2), syndecan 2 (Sdc2), FK506 binding protein 10 (Fkbp10), interferon induced protein with tetratricopeptide repeats (Ifit) 1 and Ifit2 (Table I). Following the removal of duplicates in Table I, a total of 66 genes remained, which were clustered into two groups and used to distinguish the WT, Col1a1-Krm2 and Lrp5−/− samples in a heat map (Fig. 4).

Figure 3. Protein-protein interaction network of differentially expressed genes as determined using the Search Tool for the Retrieval of Interacting Genes database (http://string-db.org/). Each node represents a protein, and each line represents the interaction between the two proteins.

Figure 4. Heat map of hub genes in the protein-protein interaction network in WT, Lrp5−/− and Col1a1-Krm2 mice. Each row represents a single gene and each column represents a sample. The gradual color alteration from orange to blue represents the process from upregulation to downregulation of genes. WT, wild-type; Col1a1, type I collagen; Krm2, kringle containing transmembrane protein 2; Lrp5, low density lipoprotein receptor-related protein 5.

Table I. Top 40 nodes with a high score in the protein-protein interaction network.

Table I.

Top 40 nodes with a high score in the protein-protein interaction network.

A, Subgraph
Node	Score
Mmp2	132118.94
Oasl2	115994.88
Ifit1	95969.49
Trim30a	92368.98
Mmp9	89980.91
Usp18	87137.86
Ifit2	83711.57
Rsad2	78389.78
Parp9	65845.10
Psmb9	62329.73
Ifi44	61890.17
Ptprc	60581.35
Bst2	52788.61
Lgals3bp	52307.49
Parp12	50636.56
Col5a1	42741.73
Ifi204	32780.95
Thbs2	32724.00
Tlr3	32063.93
Lgals9	30456.22
Pcolce	27954.06
Mmp14	27788.79
Ifi203	27352.81
Pcna	26264.86
Nmi	25586.03
Frk	25365.84
Loxl1	24301.15
Lox	23572.33
Isg20	21952.91
Rhoc	20744.50
Mdk	20643.03
Fstl1	20386.14
Ddx41	19990.48
Sdc1	18627.42
Serpinh1	17851.99
Gpx8	17771.15
Rcn3	17400.74
Cd68	17032.79
Sdc2	14213.72
Fkbp10	14067.95
B, Degree
Mmp2	47.00
Ptprc	41.00
Mmp9	40.00
Oasl2	34.00
Pcna	31.00
Frk	31.00
Col5a1	27.00
Rhoc	27.00
Ddx41	25.00
Prkar2b	24.00
Spna1	23.00
Mapk13	23.00
Ifit1	22.00
Usp18	22.00
Psmb9	22.00
Gpx8	22.00
Acacb	22.00
Slc4a1	22.00
Trim30a	21.00
Pcolce	21.00
Mmp14	20.00
Obscn	20.00
Ifit2	19.00
Rsad2	19.00
Thbs2	19.00
Myh4	19.00
Alas2	19.00
Slc2a4	18.00
Actn3	18.00
Klf1	17.00
Parp9	16.00
Tlr3	16.00
Loxl1	16.00
Lox	16.00
Mdk	16.00
Fgfr2	16.00
Smc2	16.00
Ifi44	15.00
Cd68	15.00
Fkbp10	15.00
C, Betweenness
Ptprc	30306.94
Mmp2	30018.82
Pcna	26498.73
Oasl2	24770.62
Frk	20782.62
Rhoc	17367.36
Acacb	16584.40
Mmp9	16291.33
Spna1	13323.41
Mapk13	13268.12
Prkar2b	11615.84
Gpx8	11494.19
Slc2a4	11248.73
Ddx41	10606.64
Obscn	10462.83
Pfas	9557.37
Psmb9	9042.65
Col5a1	8806.07
Rps6ka1	8625.42
Pcolce	7343.19
Mmp14	6612.19
Actn3	6597.27
Dlg3	6507.87
Msn	6370.79
Myh4	6325.04
Alas2	6084.76
Fgfr2	5320.19
Glrx5	5243.41
Mdk	5180.35
Slc4a1	5082.71
Atp8a1	5063.36
Copb2	5019.16
Ppargc1a	4905.46
Cd68	4728.94
H2-Aa	4708.25
Cxcr4	4565.29
Ncf4	4559.99
Rps3	4274.92
Hk2	4171.70
Thbs2	3997.35

The 66 hub genes were significantly associated to the five signaling pathways (Table II). Matrix metalloproteinase (Mmp) 2 and Mmp9 were associated with the leukocyte transendothelial migration pathway, whereas Thbs2 and Sdc2 were associated with the extracellular matrix (ECM)-receptor interaction pathway. The protein tyrosine phosphatase receptor type C and Sdc2 were implicated in the cell adhesion molecule pathway (Table II).

Table II. Pathways enriched by the hub genes in the protein-protein interaction network.

Table II.

Pathways enriched by the hub genes in the protein-protein interaction network.

KEGG entry term: pathway	P-value	Gene count	Genes
mmu04670: Leukocyte transendothelial migration	1.82×10−5	7	Mapk13, Cxcr4, Mmp9, Ncf4, Msn, Actn3, Mmp2
mmu04910: Insulin signaling pathway	0.00420	5	Prkar2B, Slc2A4, Hk2, Acacb, Ppargc1A
mmu04512: ECM-receptor interaction	0.00740	4	Sdc1, Thbs2, Col5A1, Sdc2
mmu04514: Cell adhesion molecules	0.03803	4	Ptprc, Sdc1, H2-Aa, Sdc2
mmu04920: Adipocytokine signaling pathway	0.04178	3	Slc2A4, Acacb, Ppargc1A

[i] KEGG, Kyoto Encyclopedia of Genes and Genomes; ECM, extracellular matrix.

Analysis of the coexpression network

In order to investigate the coexpression associations between the selected hub genes and additional DEGs, a coexpression network was constructed. A total of 21 hub genes were determined to coexpress with additional DEGs (Fig. 5). A set of hub genes were observed to coexpress with each other, including Thbs2, Sdc2 and Fkbp10, as well as Oasl2, Ifit1 and Ifit2 (Fig. 5).

Figure 5. Coexpression network of hub genes and additional differentially expressed genes. Red-colored nodes represent hub genes in the protein-protein interaction network, and blue-colored genes represent differentially expressed genes that are not hub genes. Each node represents a protein, and each line represents the coexpression association between the two genes. Solid lines indicate the coexpression (pearson correlation coefficient) >0.9 (positive correlation), and the dotted lines indicate the coexpression (pearson correlation coefficient) <0.9 (negative correlation).

Discussion

In the present study, a set of 841 DEGs (335 upregulated and 506 downregulated) and 50 DEGs (16 upregulated and 34 downregulated) were identified in the Col1a1-Krm2 vs. WT and Lrp5−/− vs. WT groups, respectively. A number of DEGs demonstrated a high score in the PPI network, and were coexpressed in the coexpression network. These genes included Thbs2, Sdc2 and Fkbp10, as well as Oasl2, Ifit1 and Ifit2. Thbs2 and Sdc2 were associated with the ECM-receptor interaction pathway.

Thbs2 is a part of the thrombospondin family and mediates cell-to-cell and cell-to-matrix interactions (25). A previous review reported that disrupted Thbs2 expression increases cortical bone density, accelerates fracture healing, induces resistance to ovariectomy-induced bone loss and alters the pattern of load-induced bone formation (26). In Thbs2-null mice, marrow-derived osteoprogenitor cells are increased, and endosteal bone formation is promoted, indicating that Thbs2 modulates the proliferation of osteoprogenitor cells and bone remodeling (27,28). Sdc2 functions as an integral membrane protein and mediates cell-to-matrix interactions via its ECM protein receptor (29). Sdc2 is a crucial determinant of chemosensitivity in osteoblasts, and it stimulates the mitogenic activity of granulocyte-macrophage colony-stimulating factor (30). Fkbp10 is a part of the FKBP-type peptidyl-prolyl cis/trans isomerase family and interacts with collagens (31). A homozygous splicing mutation in Fkbp10 leads to osteogenesis imperfecta with a mineralization defect via a reduction in bone collagen content (32,33). There is no direct evidence to implicate Sdc2 and Fkbp10 in osteoporotic fracture healing, however, they are thought to coexpress with Thbs2. Therefore, Sdc2 and Fkbp10, as well as Thbs2 may serve key roles during the fracture healing process in osteoporosis, via their coexpression associations with each other.

In the present study, Oasl2, Ifit1 and Ifit2 demonstrated high scores in the PPI network and coexpressed with each other. Ifit1 and Ifit2 were interferon-induced proteins containing tetratricopeptide repeats (34). Ifit1 is known to be an important innate immune bottleneck (35). During the response of osteoblasts to immune cytokine interferon-β, the expression of Ifit1 is induced (36). Ifit2 and Oasl2 are involved in innate immunity (37,38). Only a limited number of studies have investigated the association between the Ifit1, Ifit2 and Oasl2 genes and fracture repair; however they present potential novel candidates for osteoporotic fracture repair therapies.

In the present study, the number of identified DEGs in the Col1a1-Krm2 vs. WT group was markedly higher than that observed in the Lrp5−/− vs. WT group, which was consistent with previous findings (9). According to the DEGs enrichment analysis, the DEGs in the Col1a1-Krm2 vs. WT group were primarily associated with immunity and cell adhesion. By contrast, the DEGs in the Lrp5−/− vs. WT group were significantly associated with muscle system processes (GO: 0003012) and the regulation of transcription (GO: 0006355). These results suggest that during the fracture repair process in osteoporosis, the DEGs induced by Krm2 overexpression or Lrp5 deficiency, and their functions, may be distinctly different.

Compared with the findings presented by Liedert et al (9), the present study identified the interactions and coexpression patterns among a set of genes, which was not determined previously. However, these predictions require validation in further studies. In a future study, the DEGs and their interactions will be determined in patients.

In conclusion, a series of DEGs, including Thbs2, Sdc2 and Fkbp10, as well as Oasl2, Ifit1 and Ifit2, demonstrated a significant role in the PPI network and were observed to form co-expression patterns. The results suggest that these genes may serve crucial roles during the fracture repair process in osteoporosis. Sdc2, Fkbp10, Oasl2, Ifit1 and Ifit2 were demonstrated to be novel genes associated with osteoporotic fracture healing.

Glossary

Abbreviations

Abbreviations:

References