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% ([@b1-mmr-15-06-3969]). The prevalence of osteoporotic fractures, hip fractures in particular, is increasing in many regions of the world ([@b2-mmr-15-06-3969]). 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 ([@b3-mmr-15-06-3969]). 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 ([@b4-mmr-15-06-3969]). 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 ([@b5-mmr-15-06-3969]). *Lrp5* deficiency results in decreased osteoblast proliferation and function, which induces a low bone mass phenotype ([@b6-mmr-15-06-3969]). 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 ([@b7-mmr-15-06-3969]). 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 ([@b8-mmr-15-06-3969]). 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 ([@b9-mmr-15-06-3969]). In addition, this previous study identified a set of differentially expressed genes (DEGs) in the two mouse models using microarray analysis ([@b9-mmr-15-06-3969]). 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* ([@b9-mmr-15-06-3969]) 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 ([@b9-mmr-15-06-3969]) 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) ([@b10-mmr-15-06-3969]) of Bioconductor (version 3.4; <http://www.bioconductor.org/>). Subsequently, the org.Hs.eg.db (version 3.4.0) ([@b11-mmr-15-06-3969]) and illuminaHumanv3.db (version 1.26.0) ([@b12-mmr-15-06-3969]) 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>) ([@b13-mmr-15-06-3969]), 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 ([@b14-mmr-15-06-3969]). Only the genes with FDR values \<0.05 and log~2~ 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>) ([@b15-mmr-15-06-3969]) 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 ([@b16-mmr-15-06-3969]). 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 ([@b17-mmr-15-06-3969]). 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 ([@b18-mmr-15-06-3969]). 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 ([@b19-mmr-15-06-3969]), betweenness centrality ([@b20-mmr-15-06-3969]) and subgraph centrality methods ([@b21-mmr-15-06-3969]). The scores obtained from the degree, betweenness and subgraph methods were calculated using the CytoNCA plug-in (version 2.1.6) ([@b22-mmr-15-06-3969]) 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) ([@b23-mmr-15-06-3969]).

### Coexpression associations of hub genes with DEGs

The Pearsons correlation coefficient (PCC) method ([@b24-mmr-15-06-3969]) 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](#f1-mmr-15-06-3969){ref-type="fig"}).

### 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](#f2-mmr-15-06-3969){ref-type="fig"}). 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](#f2-mmr-15-06-3969){ref-type="fig"}).

In the *Lrp5*^−/−^ vs. WT group, the upregulated genes were implicated in muscle contraction (GO: 0006936) and muscle system process (GO: 0003012) ([Fig. 2B](#f2-mmr-15-06-3969){ref-type="fig"}). The downregulated genes were markedly associated with the regulation of transcription (GO: 0006355) and RNA metabolic processes (GO: 0051252) ([Fig. 2B](#f2-mmr-15-06-3969){ref-type="fig"}). 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](#f3-mmr-15-06-3969){ref-type="fig"}). 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](#tI-mmr-15-06-3969){ref-type="table"}). Following the removal of duplicates in [Table I](#tI-mmr-15-06-3969){ref-type="table"}, 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](#f4-mmr-15-06-3969){ref-type="fig"}).

The 66 hub genes were significantly associated to the five signaling pathways ([Table II](#tII-mmr-15-06-3969){ref-type="table"}). 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](#tII-mmr-15-06-3969){ref-type="table"}).

### 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](#f5-mmr-15-06-3969){ref-type="fig"}). 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](#f5-mmr-15-06-3969){ref-type="fig"}).

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 ([@b25-mmr-15-06-3969]). 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 ([@b26-mmr-15-06-3969]). 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 ([@b27-mmr-15-06-3969],[@b28-mmr-15-06-3969]). *Sdc2* functions as an integral membrane protein and mediates cell-to-matrix interactions via its ECM protein receptor ([@b29-mmr-15-06-3969]). *Sdc2* is a crucial determinant of chemosensitivity in osteoblasts, and it stimulates the mitogenic activity of granulocyte-macrophage colony-stimulating factor ([@b30-mmr-15-06-3969]). *Fkbp10* is a part of the FKBP-type peptidyl-prolyl cis/trans isomerase family and interacts with collagens ([@b31-mmr-15-06-3969]). A homozygous splicing mutation in *Fkbp10* leads to osteogenesis imperfecta with a mineralization defect via a reduction in bone collagen content ([@b32-mmr-15-06-3969],[@b33-mmr-15-06-3969]). 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 ([@b34-mmr-15-06-3969]). *Ifit1* is known to be an important innate immune bottleneck ([@b35-mmr-15-06-3969]). During the response of osteoblasts to immune cytokine interferon-β, the expression of *Ifit1* is induced ([@b36-mmr-15-06-3969]). *Ifit2* and *Oasl2* are involved in innate immunity ([@b37-mmr-15-06-3969],[@b38-mmr-15-06-3969]). 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 ([@b9-mmr-15-06-3969]). 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* ([@b9-mmr-15-06-3969]), 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.

DEGs

:   differentially expressed genes

PPIs

:   protein-protein interactions

BMP-2

:   bone morphogenetic protein-2

PTH

:   parathyroid hormone

*LRP5*

:   lipoprotein receptor-related protein 5

*KREMEN2*

:   kringle containing transmembrane protein 2

DKK1

:   dickkopf homolog 1

GEO

:   Gene Expression Omnibus database

RMA

:   robust microarray analysis

FDR

:   false discovery rate

GO

:   Gene Ontology

KEGG

:   Kyoto Encyclopedia of Genes and Genomes

DAVID

:   Database for Annotation, Visualization and Integrated Discovery

STRING

:   Search Tool for the Retrieval of Interacting Genes

PCC

:   pearson correlation coefficient

CAMs

:   cell adhesion molecules

*Thbs2*

:   thrombospondin 2

*Sdc2*

:   syndecan 2

*Fkbp10*

:   FK506 Binding Protein 10

![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.](MMR-15-06-3969-g00){#f1-mmr-15-06-3969}

![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.](MMR-15-06-3969-g01){#f2-mmr-15-06-3969}

![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.](MMR-15-06-3969-g02){#f3-mmr-15-06-3969}

![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.](MMR-15-06-3969-g03){#f4-mmr-15-06-3969}

![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).](MMR-15-06-3969-g04){#f5-mmr-15-06-3969}

###### 

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

  A, Subgraph      
  ---------------- -----------
  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

###### 

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

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