Transcriptional profiling of intramembranous and endochondral ossification after fracture in mice

Abstract

Bone fracture repair represents an important clinical challenge with nearly 1 million non-union fractures occurring annually in the U.S. Gene expression differs between non-union and healthy repair, suggesting there is a pattern of gene expression that is indicative of optimal repair. Despite this, the gene expression profile of fracture repair remains incompletely understood. In this work, we used RNA-seq of two well-established murine fracture models to describe gene expression of intramembranous and endochondral bone formation. We used top differentially expressed genes, enriched gene ontology terms and pathways, callus cellular phenotyping, and histology to describe and contrast these bone formation processes across time. Intramembranous repair, as modeled by ulnar stress fracture, and endochondral repair, as modeled by femur full fracture, exhibited vastly different transcriptional profiles throughout repair. Stress fracture healing had enriched differentially expressed genes associated with bone repair and osteoblasts, highlighting the strong osteogenic repair process of this model. Interestingly, the PI3K-Akt signaling pathway was one of only a few pathways uniquely enriched in stress fracture repair. Full fracture repair involved a higher level of inflammatory and immune cell related genes than did stress fracture repair. Full fracture repair also differed from stress fracture in a robust downregulation of ion channel genes following injury, the role of which in fracture repair is unclear. This study offers a broad description of gene expression in intramembranous and endochondral ossification across several time points throughout repair and suggests several potentially intriguing genes, pathways, and cells whose role in fracture repair requires further study.

Keywords: Animal models; Bone; Fracture repair; Transcriptome.

Figure 1.. Experimental overview of fracture models, time points, and RNA-seq pipeline.

(A) 97 female C57BL/6J mice were injured with stress or full fracture, sacrificed at 4 hours, 1, 3, 5, 7, or 14 days post-injury, and processed for RNA-seq or histology. (B) RNA-seq analysis began with RNA extraction from pulverized bone tissue. RNA was sequenced and reads were mapped to mm 10 genome. Differential expression analysis was performed to create lists of differentially expressed genes (DEGs) for each time point. Comparison to previously generated qPCR data of stress and full fracture callus was used to validate RNA-seq data. Callus component analysis using transcriptional profiling of callus cells was performed at each time point. Pathway analysis and GO annotation were performed on the DEG lists of each time point.

Figure 2.. Progression of healing of stress fracture and full fracture over time.

H&E staining of paraffin sections of stress fracture or full fracture callus. Scale bars represent 250 μm for stress fracture images and 1 mm for full fracture images. Ct = Cortical Bone, Ma = Marrow, Mu = Muscle, ▲ = stress fracture or full fracture line, dotted line = callus.

Figure 3.. Validation of RNA-seq data.

Principle component analysis (PCA) of all samples in stress fracture (A) and full fracture (B). Samples within time points cluster together and time point clusters were grouped closely with adjacent time points. (C) Correlation of RNA-seq data with analogous qPCR data from previous published reports of stress fracture (blue triangle) and full fracture (grey circle) and concurrent qPCR of full fracture samples (red circle). Log2 fold change (injured vs. control) from RNA-seq was plotted versus log2FC from qPCR (+/− standard deviation). Data is plotted from multiple time points for each injury. A linear regression confirmed strong correlation between RNA-seq and qPCR data. Genes from published reports are detailed in Supplementary Table 1. Select genes from concurrent qPCR are labeled with arrows.

Figure 4.. DEG trends across time.

Total number of DEGs at each time point are shown for stress fracture (A) and full fracture (B). Following stress fracture (S.Fx), few DEGs occurred early, but DEGs reached robust levels at later time points and peaked at day 5. Following full fracture (Fx), DEGs were immediately expressed at high levels that persisted through our experimental window, and peaked at day 7. Comparisons of the overlap of DEGs at each time point for stress fracture (C) and full fracture (D) are shown in 5-way Venn diagrams. Total number of DEGs at each time point are shown in parentheses. The largest lists were those unique to one time point or shared between adjacent time points.

Figure 5.. Gene ontology analysis of DEGs.

(A) Pipeline for gene ontology (GO) term analysis. (B) Comparison of the 115 GO terms uniquely enriched in stress fracture (dark blue), enriched in both injuries (light blue), or uniquely enriched in full fracture repair (grey). (C) Graphical list of all 115 enriched GO terms across all injuries and time points. Colors on chart match location of term from panel B. Go terms uniquely enriched in early full fracture time points were broken out into the smaller sub-panel for figure compactness (Black arrows). Blank spaces indicate GO term was not statistically enriched at time point or injury. H4 – 4 hour, D1 – day 1, D3 – day 3, D5 – day 5, D7 – day 7, D14 – day 14. S.Fx – stress fracture, Fx- full fracture

Figure 6.. Pathway analysis of DEGs.

(A) Pipeline for pathway analysis. (B) Comparison of pathways enriched in stress and full fracture repair. (C) Graphical list of all 60 enriched pathways across injuries and time points. Colors on chart match location of pathway on panel B. Blank spaces indicate no significance of pathway at that time point/fracture condition. H4 – 4 hour, D1 – day 1, D3 – day 3, D5 – day 5, D7 – day 7, D14 – day 14. S.Fx – stress fracture, Fx- full fracture

Figure 7.. Callus component phenotyping heat maps.

Heat maps (log2 fold change) were generated for each time point and injury of curated genes lists associated with (A) Neutrophils, (B) Macrophages, (C) Monocytes, (D) T Cells, (E) B Cells, (F) Osteoblasts, (G) Osteoclasts, (H) Endothelial Cells, and (I) Chondrocytes. Gene lists were generated using literature sources in Supplementary Table 2.

Figure 8.. Immune cell infiltration in fracture callus.

Immunohistochemistry with antibodies for Gr-1, F4/80, and CD45 was used to stain for neutrophils, macrophages, and leukocytes, respectively. Cartoon depictions of the fracture callus show Region of interest (ROI) (Black box) in reference to (A) stress fracture (S.Fx) or (B) full fracture (Fx) location. All images are 20X magnifications and black scale bars are 100 μm. Ct = Cortical Bone, Ma = Marrow, Ps = Periosteum

Figure 9.. Ion channels are downregulated following full fracture.

The Log2 FC of Calcium, Potassium, and Sodium Voltage Gated Channels which were DEGs are displayed across all time points for stress fracture (S.FX) and full fracture (FX). These channels were disproportionally down-regulated throughout full fracture repair. Due to these down-regulated DEGs, Go terms and pathways such as “ion transport”, “cellular calcium homeostasis”, “potassium ion transmembrane transport’, “calcium signaling pathway”, and “cardiac muscle contraction” were enriched during full fracture repair. A full list of these genes is provided in Supplementary Table 4.