Abstract
We aimed to investigate differentially expressed genes (DEGs) in different stages after femoral fracture based on rat models, providing the basis for the treatment of sport-related fractures. Gene expression data GSE3298 was downloaded from Gene Expression Omnibus (GEO), including 16 chips. All femoral fracture samples were classified into earlier fracture stage and later fracture stage. Total 87 DEGs simultaneously occurred in two stages, of which 4 genes showed opposite expression tendency. Out of the 4 genes, Rest and Cst8 were hub nodes in protein-protein interaction (PPI) network. The GO (Gene Ontology) function enrichment analysis verified that nutrition supply related genes were enriched in the earlier stage and neuron growth related genes were enriched in the later stage. Calcium signaling pathway was the most significant pathway in earlier stage; in later stage, DEGs were enriched into 2 neurodevelopment-related pathways. Analysis of Pearson's correlation coefficient showed that a total of 3,300 genes were significantly associated with fracture time, none of which was overlapped with identified DEGs. This study suggested that Rest and Cst8 might act as potential indicators for fracture healing. Calcium signaling pathway and neurodevelopment-related pathways might be deeply involved in bone healing after femoral fracture.
1. Introduction
As the 2008 Beijing Olympics were successfully held in Beijing, sports developed rapidly in China. More and more inhabitants, professional or amateur, take part in daily physical activities. However, improper movement may cause injury. The intense sports (like pugilism, football, and basketball) and hazardous sports (like motorcycle race, drift motion, and bungee jumping) are all high-risk sports. Collisions with the ground, objects, and other players are common, and unexpected dynamic force on limbs and joints can cause injury [1]. In human, the femur fracture is one of the most common injuries resulted from improper movement [2]. The femur is the only bone in the thigh with the formation of long, slender, and cylindrical bone and is capable of walking, running, or jumping [3]. The femoral fracture is involved in the femoral head, femoral neck, or the shaft of the femur, accounting for 1-2% of all fractures in children and adolescents [4, 5].
For a long time, femoral fractures have been treated by using traction and/or casting [6]. More recently, surgery has gained popularity [7]. However, femur fracture is still difficult to manage because of the multifocal fractures of the femur [8, 9]. Although numerous surgical operations have been described to manage this injury, evidence for which to choose is lacking and individual approach is strongly emphasized during the treatment of these injuries [9, 10]. It is necessary for us to study the differences of gene expression at different stages after femoral fracture, in the purpose of finding the indicator of fracture healing.
Rats grow rapidly to attain their adult size. At 4 weeks, femur growth is near its maximum rate. At the age of 10 weeks, the linear growth of femur has slowed due to mitosis and hypertrophy in the chondrocytes of the physis [11, 12]. Based on rat model, changes in mRNA gene expression of femoral heading have been reported [12, 13]. Briefly, in 4-week-old female Sprague-Dawley (SD) rats, at 0 (intact), 0.1, 0.4, 1, 2, 3, 4, and 6 weeks after fracture, mRNA gene expression in the femoral heading after unilateral midshaft femoral fracture was identified, including 8,002 genes, about half increasing and half decreasing. These upregulated genes were related to cartilage, blood vessels, osteoprotegerin, osteomodulin, and most ribosomal proteins. Meanwhile, downregulated genes were related to bone, growth promoting cytokines, G proteins, GTPase-mediated signal transduction factors, cytokine receptors, mitosis, integrin-linked kinase, and the cytoskeleton. The relevant microarray data were deposited in GEO (Gene Expression Omnibus) database (ID: GSE3298) [12, 13].
In this present study, based on the microarray data of GSE3298, 2 weeks after femoral fracture was chosen as a split point, and thus the earlier stage and later stage were grouped. We aimed to identify DEGs at different stages of femoral fracture healing by bioinformatics methods, in order to provide the basis for the treatment of sport-related fractures.
2. Materials and Methods
2.1. Microarray Data
The mRNA expression profiling data was obtained from the research of Meyer et al., which were displayed in GEO (http://www.ncbi.nlm.nih.gov/geo/) database (ID: GSE3298) [12]. Briefly, female SD rats, aged 4 weeks at surgery, were subjected to a unilateral, simple, transverse, and middiaphyseal femoral fracture and stabilized with an intramedullary rod. At 0 (intact), 0.1, 0.4, 1, 2, 3, 4, and 6 weeks after fracture, the femoral head with the proximal physis was collected from fractured and intact femora. The RNA was extracted, processed to biotin labeled cRNA, and hybridized to Affymetrix Rat 230 2.0 GeneChip microarrays. The full microarray data has been deposited in the NCBI GEO as series GSE3298.
2.2. Data Preprocess
The microarray data in CEL files were downloaded from GEO database, including 16 chips, converted into fluorescence intensity values and standardized via the robust multiarray average (RMA) method [14]. For genes corresponding to multiple probe sets that had a plurality of expression values, the expression values of those probe sets were summed.
2.3. Differentially Expressed Gene Analysis
Considering the different healing level in different periods after fracture, 2 weeks was set as the split point. Chips data were divided into 2 groups: earlier stage (0.1, 0.4, 1 and 2 weeks after fracture) and later stage (2, 3, 4, and 6 weeks after fracture). The LIMMA package in R language was used to identify DEGs between earlier stage and later stage [15]. The P value <0.05 and the were used as the cut-off criterion.
2.4. Construction of Interaction Network
For genes differentially expressed in two stages, the STRING (Search Tool for the Retrieval of Interacting Genes) [16] database was used to analyze their interaction network. For genes with consistent expression in two stages, BisoGenet [17] software was performed to map these genes to STRING database or BOND database for interaction network analysis. The value <0.05 was chosen as cut-off criterion.
2.5. Pathway Enrichment Analysis
For function analysis of DEGs, the DEGs of two stages were, respectively, inputted into DAVID (Database for Annotation, Visualization, and Integrated Discovery) [18, 19] for KEGG (Kyoto Encyclopedia of Genes and Genomes) [20] pathway and GO (Gene Ontology) [21] enrichment analysis. The count number larger than 5 and value less than 0.01 were chosen as cut-off criterion.
2.6. Correlation Analysis
A Pearson correlation coefficient was calculated between expression level of every expressed gene after fracture and fracture time via cor.test in R language [22]. The was chosen as cut-off criterion. Then, DAVID tool was used to identify function classification associated with these significant genes.
3. Results
3.1. Differentially Expressed Genes
After standardization, there were 31,042 probes corresponding to 30,641 genes. In earlier stage, total 1,004 DEGs had been identified, including 301 upregulated genes and 703 downregulated genes. In later stage, total 986 DEGs were obtained, including 446 upregulated genes and 540 downregulated genes. The most significant DEGs from two stages were displayed in Tables 1 and 2.
Among DEGs from two stages, 87 DEGs occurred in both earlier stage and later stage, including 26 upregulated genes, 57 downregulated genes, and 4 differentially regulated genes. Briefly, one DEG (GenBankAcc: BF402112) was upregulated in earlier stage and downregulated in later stage, and 3 DEGs (GenBankAcc: AF037203 (Rest), AI071395, and NM019258 (Cst8)) were downregulated in earlier stage and upregulated in later stage (Figure 1).
3.2. Interaction Network of DEGs
The obtained 87 DEGs were mapped to STRING in order to construct the interaction network. Among the 26 upregulated DEGs, only MOBKL3 was the hub node that interacted with other genes (Figure 2(a)). Meanwhile, among the 57 downregulated DEGs, 5 DEGs showed interaction with other genes in rats (Figure 2(b)), such as syt and stx families.
(a)
(b)
In addition, the protein-protein interaction (PPI) network of Rest was constructed via STRING tool and displayed in Figure 3, suggesting that Rest protein might interact with 11 proteins in rats. The PPI network of Cst8 was built as well, in which Cst8 was the hub protein connected with 10 proteins (Figure 4).
3.3. Enrichment Analysis of DEGs
For function analysis, all DEGs were inputted into DAVID for GO function and KEGG pathway enrichment analysis. was set as significant difference.
GO function enrichment analysis of DEGs in earlier stage showed 170 significant GO terms, which were divided into 20 clusters, including material transportation in cells, regulation of biological process, structure development, neurodevelopment, and the blood pressure regulation. The most significant GO term was GO: 0051179 (localization), of which the fold enrichment was 1.576. The top 10 GO terms were shown in Table 3 (upper). Similarly, total 111 significant GO terms were obtained from GO analysis of DEGs in later stage and were divided into 13 clusters, including neurons and synapses development, ion transport, regulation of gene expression, and hormone secretion. The most significant GO term was GO: 0045202 (synapse), of which the fold enrichment was 2.53. The top 10 GO terms of later stage were shown in Table 3 (lower).
Additionally, KEGG pathway enrichment analysis of DEGs in earlier stage showed 5 significant pathways (Table 4, upper). Calcium signaling pathway was the most significant one (fold enrichment: 2.69). Meanwhile, DEGs in later stage were enriched into 2 significant pathways, mainly focused on neurodevelopment (Table 4, lower).
3.4. Correlation Analysis
Among the expressed genes after fracture, a Pearson correlation coefficient was calculated between gene expression level and fracture time via cor.test in R language. With the strict cut-off of , total 3,300 genes significantly associated with fracture time were collected, including negative correlation (1,714 genes) and positive correlation (1,586 genes) (Table 5). None of the 3,300 significant correlation genes was overlapped with DEGs identified using LIMMA package. The function annotation of these significant correlation genes showed relationship with illness, cancer, and immune system, indicating that surgical approach did not cause serious damage to health of animals. Besides, the correlation analysis of DEGs from two stages did not show significant correlation with fracture time.
4. Discussion
In daily life, sport-related fractures are common in adolescents, particularly in males [23]. Femoral fracture, a common sport injury, has great impact on human physical exercise ability and improper treatment can lead to nerve injury, infection, pain, or dyskinesia [5]. For professional athletes, femoral fracture is very popular and the outcome of treatment affects their athletic career [24]. It is necessary for us to identify DEGs after femoral fracture and to explore the key gene of the bone healing, which will provide theoretical basis for future treatment of these sport-related fractures.
In this study, the chip data were divided into earlier stage and later stage based on 7 time points after femoral fracture. In earlier and later stages, 1,004 and 986 DEGs were identified by comparing with control group, respectively. For example, among DEGs in early stage, Oprm1 was opioid receptor [25], the reduced expression of which in dorsal root ganglion neurons was found to be associated with bone cancer pain in mouse models [26]. Tmem200a was a transmembrane protein which might inhibit overgrowth of myelocyte [27]. Meanwhile, among DEGs in later stage, Bcl2l 1 encoded Bcl-2-like 1 protein, a critical regulator of programmed cell death, belongs to Bcl-2 protein family [28]. Consistently, it is reported that Bcl-2 plays an important role in regulating the apoptosis of osteoclast and osteocyte [29]. Furthermore, Wt1 (Wilms tumor 1) might act as a novel oncogene facilitating development of the lethal metastatic phenotype in prostate cancer [30].
Among DEGs between two stages, there was no significant difference in the number of DEGs and total 87 DEGs were shared by two stages, indicating different expression profiles between two stages. There were 4 DEGs oppositely regulated in earlier and later stages, which might act as indicators for femur healing. Among the 4 DEGs, Rest, similar to Tmem200a, might inhibit overgrowth of myelocyte combined with myc gene [31]. Rest gene is a transcriptional repressor of diverse neuronal genes, the downregulation of which contributed to the proper development of neurons [32]. Similarly, in the current study, Rest was downregulated in earlier stage but upregulated in later stage. Moreover, Rest is involved in the differentiation from pluripotent cell to neural stem cell and from stem cell to mature neurons [33]. Cyst8 belongs to cystatin family of proteins [34]. Many members of the cystatin superfamily such as gelatin could protect matrix metalloproteinases without affecting their biological activities, which are critical for tissue modeling [35]. Total 57 DEGs were downregulated in both two stages, of which interaction network showed that 5 genes were interacted with other reported genes in rats, such as syt and stx families. Syt1 was a key factor controlling neurotransmitters release via binding to calcium ion [36]. Consistently, this study showed that calcium signaling pathway was also enriched in early stage, suggesting the critical role of calcium signaling in bone healing after femoral fracture. Besides, Syt1 might control neural signal transmission combined with SNAP-25 [37, 38] and STX1A [39]. STX1A was involved in vesicle fusion process which is critical for calcium-dependent neurotransmitters release. Importantly, it has been reported that increase of Syt-1 might play a role in impairment of learning and memories attributed to aging in mouse model [40].
Pearson’s correlation coefficient analysis between gene expression and fracture time indicated that significant correlation genes between gene expression and fracture time were not overlapped with identified DEGs, which demonstrated that rats underwent surgical operation without other infections and injuries.
GO analysis of DEGs from two stages was enriched into different GO terms. Briefly, in earlier stage, abundant DEGs were related to material transportation and synthesis in cells, and a few genes were enriched in synapse growth, while, in later stage, in contrast, the majority of DEGs were related to synapse growth and a small number of genes were related to transporter activity. These discrepancies suggested that fracture healing involved distinct functions in earlier and later stages. Besides, system development was enriched in both earlier and later stages, revealing its importance in the whole process of fracture healing. KEGG pathway analysis showed that neuroactive ligand-receptor interaction was needed in two stages, indicating its important role in fracture healing. Meanwhile, in earlier stage, DEGs were significantly enriched into calcium signaling pathway and neuroactive ligand-receptor interaction pathway. For later stage, neuroactive ligand-receptor interaction becomes the most important pathway, and axon guidance pathway was also enriched. The two pathways were closely associated with neurodevelopment. These findings indicated difference of physical growth between two stages. In earlier stage, more nutrients for vegetative growth were needed to repair fracture; in later stage, nervous systems were repaired to restore movement ability, which were consistent with general understanding.
5. Conclusions
In conclusion, based on rat model, identification of DEGs after femoral fractures was useful for investigation of the proper treatment and providing foundation for exercise capacity recovery. However, further genetic studies were still needed to confirm our observation.
Conflict of Interests
The author declares that there is no conflict of interests regarding the publication of this paper.