Cross-Species RNA-Seq Study Comparing Transcriptomes of Enriched Osteocyte Populations in the Tibia and Skull

Local site-specific differences between bones in different regions of the skeleton account for their different properties and functions. To identify mechanisms behind these differences, we have performed a cross-species study comparing RNA transcriptomes of cranial and tibial osteocytes, from bones with very different primary functions and physiological responses, collected from the same individual mouse, rat, and rhesus macaque. Bioinformatic analysis was performed to identify 32 genes changed in the same direction between sites and shared across all three species. Several well-established key genes in bone growth and remodeling were upregulated in the tibias of all three species (BMP7, DKK1, FGF1, FRZB, SOST). Many of them associate or crosstalk with the Wnt signaling pathway. These results suggest Wnt signaling-related candidates for different control of regulatory mechanisms in bone homeostasis in the skull and tibia and indicate a different balance between genetically determined structure and feedback mechanisms to strains induced by mechanical loading at the different sites.

Introduction

Patterning during embryogenesis accounts for the development of all the specialized tissues and organs in multicellular organisms (1). Skeletal patterning accounts for the shape and structure of different bones and their behavior during life. However, mechanical forces provide additional potent influences on most parts of the skeleton, both before (2, 3) and after birth (4). This process of skeletal adaptation tunes bone strength to prevailing functional demands, adding material above a genetically determined baseline structure that forms and is retained even in the absence of loading stimuli.

The physical demands on bones, and therefore their optimal structures vary with their site in the skeleton. These form/function relationships vary widely throughout the skeleton (5), but we can consider examples at the opposite ends of the spectrum. Long bones need to be rigid levers for efficient locomotion, while the skull provides protection for the vulnerable brain. Structure of long bones in the limbs is regulated by the loads they experience during physiological activity (6), while the skull must retain mass and strength to resist damaging events despite the absence of any significant or regular loading. We have shown that in a human subject, the peak strains experienced in the skull even under extreme physiological circumstances do not reach 1/10th of the magnitude or rate those in the tibia of the same person (7). Such low strains would be perceived as disuse in a long bone, yet skull bone integrity is retained even in wasting conditions that deplete material in other parts of the skeleton. As all the bones in the skeleton experience the same circulating milieu of hormones and nutrients, it seems therefore that site-specific differences in regulation of bone formation and resorption must exist to account for regional differences. It might be considered that any differences between the skull and long bones are related to the embryological origins of the bone at the two sites, rather than their functional differences. This is unlikely because the clavicle is an example of a bone that is formed partly through endochondral ossification (like the long bones) and partly through intramembranous ossification (like the skull). The clavicles are bones which are known to adapt to function and lose bone mass through disuse, like long bones and unlike the skull (8), suggesting that intramembranous ossification is not invariably associated with a lack of sensitivity to disuse. Furthermore, the lateral end of the clavicle, which is the part formed by intramembranous ossification is affected by a bone wasting condition that has no effect on bone mass in the skull (9).

In order to identify the biological basis for site specific regulation of bone mass and structure we have used RNA sequencing to compare the transcriptomes of samples of skull and tibial bones from three species: mouse, rat, and rhesus macaque. Others have tried to address this question in the past, using bone samples (that included bone marrow) from young rats (10, 11) and mice (12), identifying several 100 genes that appeared regulated differently between sites. The cross-species design allowed us to focus our results on what we believe are broad biological mechanisms, and resulted in only a relatively small number of differences between sites that were shared across species and are therefore likely to reflect mechanisms in humans. The method we designed was to focus attention on enriched osteocytes populations, because osteocytes are widely accepted to be the mechanosensors in bone (13), responding rapidly to mechanical loading (14), and so influencing bone formation and resorption. More importantly, evidence has suggested that osteocytes, constituting more than 90% of bone cells, use the Wnt signaling pathway involving sclerostin and Dkk1, and crosstalk with the prostaglandin pathway in response to mechanical loading (15–17). Our current study has provided further evidence to support this notion but from a site-specific differences perspective.

Materials and Methods

Experimental Animals

Four-month old female C57BL/6 mice (Charles River, Harlow, UK), 6-month old female Wistar rats (Charles River, Harlow, UK), and 4–6 year old female rhesus macaques (Macaca mulatta) were used in this study. These ages represent bones which have undergone their periods of maximum growth rate for each species, and so transcriptome analysis should not be confounded by rapid growth of the skeleton (18). Long bones (tibias) and calvariae were dissected free of soft tissue after animals were euthanased. All procedures complied with the United Kingdom Animals (Scientific Procedures) Act 1986 (PPL 40/3499).

Bone Tissue Preparation

To obtain highly enriched populations of osteocytes, tibias of mouse and rat were cut to expose the metaphysis and centrifuged briefly to remove bone marrow (1,500 g, 30 s). Calvarial samples were obtained from the rodents with scissors, removing the sutures and avoiding any marrow spaces (diploe). Samples were obtained from the Macaque bones with a hacksaw irrigated with normal saline. All samples were dissected free of periosteum and accessible surfaces were scraped with a scalpel. They were then immersed briefly in a 1 mg/ml solution of collagenase (Sigma Aldrich) for 3–5 min at 37°C to remove any adherent surface cells. Samples were then washed in saline before being snap frozen by complete submersion in liquid nitrogen and then stored at −80°C. Bone samples were pulverized using a Mikro-Dismembrator-S (Braun Biotech International GmbH Melsungen, Germany), in which the bone is shaken in a robust PTFE mill chamber with an 8 mm tungsten carbide ball (both cooled in liquid nitrogen before use, and again after adding the bone pieces, before agitation). The mill with tissue sample was placed within the Mikro-Dismembrator-S and set to shake at 2,500 rpm for 45 s. The weight of the fine bone powder was recorded and it was stored for RNA extraction.

RNA Extraction From Bone

One milliliter of Trizol Reagent (Ambion) was added per 125 mg of pulverized tissue (typical tissue amounts 350–500 mg); samples were incubated in Trizol reagent for 10 min at room temperature. Samples were centrifuged at 500 g for 5 min at 4°C, the supernatant was removed carefully to ensure no debris was transferred. 0.3 ml of chloroform was added to each 1 ml of Trizol reagent, the sample was thoroughly mixed and allowed to incubate at room temperature for 5 min. The samples were centrifuged at 12,000 g for 20 min at 4°C. The colorless upper phase was collected and transferred to a separate tube, to which an equal volume of 70% ethanol was added and incubated at room temperature for 10 min.

The samples were then applied to spin cartridges from TRIzol Plus RNA purification kit (Ambion) following the manufactures instructions using the optional On-Column PureLink DNase (Invitrogen) treatment step. Samples were quantified using spectrophotometry (NanoDrop Thermo) and quality measured using a Bioanalyzer (Agilent). Only RNA with a RNA integrity number (RIN) above 8 was stored in −80°C and used for RNA-Seq. Samples used for analysis were one long bone and skull pair pooled from two rats, two long bone, and skull pairs each pooled from 3 mice, and 3 long bone, and skull pairs each from individual macaques.

Library Preparation and Sequencing

cDNA library preparation and sequencing were performed by Eurofins Genomic (Ebersberg, Germany). From the total RNA sample, poly(A)+ RNA was enriched and randomized primer was used for first strand cDNA synthesis. Sequencing was performed in the 1 × 125 bp run mode, 6 samples/channel on Illumina HiSeq 2500 and generated 30 million reads/sample on average. The %Q30 of each samples were all above 86% and detailed quality control metrics were shown in Supplementary Table 1. Complete RNA sequencing data was submitted to the Gene Expression Omnibus under accession record GSE151971.

RNA-Seq Data Processing and Functional Gene Annotation

RNA-seq data analysis pre-processing and functional gene annotation was performed by Bioinformatics Consultants (Stockholm, Sweden). Reference genome sequences and gene annotations for Rattus norvegicus (rn6) and Mus musculus (mm10) were downloaded from UCSC. For Macaca mulatta we used the MacaM genome assembly (https://www.unmc.edu/rhesusgenechip/index.htm) and its corresponding gene annotations (19). Exons of genes with multiple isoforms were merged according to their genomic coordinates using a custom Perl script. HISAT2 (v. 2.0.5) with default settings was used for mapping the sequencing reads to each respective genome (20). Sequencing reads were then counted on gene models with the HTSeq-count program (the –s flag set to “no”) (21). Differential gene expression was analyzed with the R package edgeR v. 3.16.5 for each species/individual separately (22). Any gene with zero read counts in at least one sample was removed, in order to lower the chance of quantifying expression from genes with incomplete annotation in one or more genomes. The trimmed mean of M-values (TMM) normalization was applied to raw read counts prior to testing for differential expression (23). Finally, differential expression was defined as genes with a false discovery rate <20% and an absolute fold change higher than 1.1 (24, 25). In the heat-map plots, raw read counts were TMM normalized, and a variance stabilizing transformation was applied (26). Heat-maps were rendered with the heatmap.2 function in the gplots R package. Three-way ortholog pairs across the three genomes were determined through their gene symbols with the same direction of gene expression. For functional gene annotation analysis, gene Ontology enrichment analysis was conducted on all identified three-way ortholog pairs with GOrilla (Gene Ontology enRIchment anaLysis and visuaLizAtion tool) (27), PANTHER (protein annotation through evolutionary relationship) classification system (28), and DAVID (Database for Annotation, Visualization and Integrated Discovery) (29). Gene Set Enrichment Analysis (GSEA) (30) was also conducted with GSEA v.2.2.3, ran with MSigDB v. 5.2 (gene set definitions) (31).

Results

Pipeline Design

Orthologs across the three genomes were initially determined using OrthoDB (32). A principal component analysis (PCA) was conducted on the 12 samples using normalized expression of 1:1:1 orthologs. The result of this analysis showed that samples generally clustered according to species and not by sample site (Figure 1A), which indicates that the phylogenetic divergence is stronger than site differences. Similar results can also be seen through pairwise comparisons in a heat-map where samples from the same species but different sites have more similar expression profile than samples of the same sites from a different species (Figure 1B). This could be a consequence of analysis based on the limited 14,308 three-way orthologs that can be mapped 1:1:1 between Rat-Mouse-Macaque, due to incomplete genome annotations at the time OrthoDB was compiled. Therefore, to achieve enough power to detect transcriptome differences between different sites of skeleton across species, the differential expression (DE) analysis was carried out for each species/individual separately and orthologs were then determined using a simpler approach where the ortholog relationship is found through their gene symbols and same direction of gene expression. To this end, we propose a pipeline to determine transcriptome differences across species, which comprises of sequencing reads mapped using HISAT2, DE analysis with the R package edgeR v. 3.16.5, and gene symbols based three-way ortholog pairs determination, followed by Gene Ontology (GO) enrichment analysis conducted with GOrilla, Panther classification system, and DAVID, and Gene Set Enrichment Analysis with GSEA v.2.2.3 (Figure 1C).

FIGURE 1

Figure 1. Pipeline design to discover cross species transcriptome signature of osteocytes in tibia and calvaria using RNA-Seq. (A) A principal component analysis (PCA) suggested that samples cluster according to species but not to tissue sites (data are based on ortholog triplets determined with OrthoDB). (B) Pairwise comparisons in a heat-map suggested samples from the same species but different sites have much more similar expression profile than samples of the same sites from a different species. (C) The pipeline to determine cross-species transcriptome signature: Raw RNA-seq data was aligned to relevant genomes using HISAT2. Differential expression (DE) was analyzed with the R package edgeR v. 3.16.5. Three-way ortholog pairs were determined based on gene symbols and the same direction of expression. Gene Ontology (GO) enrichment analysis was conducted with GOrilla, Panther classification system, and DAVID, in addition to Gene Set Enrichment Analysis with GSEA. After raw RNA-Seq data was mapped to relevant genomes, DE expression was performed for each species separately.

Differential Expression (DE) Analysis

DE expression was analyzed with the R package edgeR v. 3.16.5 for each species separately and differentially expressed genes were defined as genes with a false discovery rate <20% and an absolute fold change higher than 1.1 in mouse and macaque. This is to privilege sensitivity over specificity and highlight broad trends firstly (24). Therefore, 1,187 and 302 differentially expressed genes between tibias and calvariae were found in macaque and mouse, respectively (Figures 2A,B). By comparing gene symbols, there were 64 genes differentially expressed in both macaque and mouse with the same direction of change (Figure 2C). A conventional DE analysis with edgeR was not applied to the pair of rat samples as count dispersion could not be estimated. Instead, log2([cpm calvaria]/[cpm tibia]) >2 was selected as a threshold, which gave 946 DE genes in rat. By comparing gene symbols, there were 32 genes across all three species sharing same direction of gene expression (Figure 2D), with specifically interested genes annotated in Figures 2A,B. The list of genes with all quantitative data were ranked in a descending order in Table 1, based on fold changes of gene expression between calvariae and tibias.

FIGURE 2

Figure 2. Differential expression (DE) and orthologs analysis. After raw RNA-Seq data was mapped to relevant genomes, DE expression was performed for each species separately. For macaque and mouse, genes with a false discovery rate <20% and an absolute fold change higher than 1.1 between samples from tibias and calvariae were defined as having significantly differential expression. (A,B) The distribution of transcriptomes were shown in the volcano plots built based on Log2 FC and FDR, with significant altered genes marked in red and specifically interested genes annotated. (C) There were 1,187 and 302 significantly DE genes in macaque and mouse, respectively. There are 64 genes sharing the same direction of expression in both macaque and mouse. (D) For rat, log2([cpm calvaria]/[cpm tibia]) > 2 was regarded as significantly DE, which produced 946 DE genes in rat. By comparing gene symbols, there were 32 genes common to all three species that shared same direction of gene expression.

TABLE 1

Table 1. Genes sharing same direction of expression across three species (fold changes: calvaria over tibia).

Gene Ontology (GO) Enrichment Analysis

To identify enriched GO terms in this list of shared DE genes, GO enrichment analysis using the web-based GOrilla application was performed first based on the 64 DE genes shared between macaque and mouse. DE genes in rat were then taken into account for quantitative analysis. Results suggested that a group of genes (HOXA7, HOXA11, HOXC4, HOXC8, HOXC9, HOXC10, ZIC1) are enriched in pattern specification process at the level of biological process (Figures 3A,B). These enriched genes are mainly binding sequence-specific DNA as their molecular function (Figure 3C) and located in nuclear part of the cellular component (Figure 3D).

FIGURE 3

Figure 3. Gene Ontology (GO) enrichment analysis. GOrilla application was firstly performed on the 64 DE genes shared between macaque and mouse. DE genes in rat were then taken into account when quantitative analysis was performed. (A) At the level of biological process, a group of genes (HOXA7, HOXA11, HOXC4, HOXC8, HOXC9, HOXC10, ZIC1) were enriched in pattern specification process. (B) These genes were significantly up-regulated in tibias in all three species, apart from ZIC1. (C) These enriched DE genes mainly bind sequence-specific DNA as their molecular function and (D) located in nuclear part of the cellular component. P-value color scale for (A,C,D).

Gene function analysis with the PANTHER classification system suggested these DE genes mainly function in the biological processes of organismal development and morphogenesis (Table 2), not surprisingly in skeletal system development and morphogenesis (Figures 4A,B). This result is also consistent with the GSEA analysis which suggests in the tibia that there are multiple gene sets involved in skeletal development and morphogenesis. PANTHER protein class analysis confirms the results by GOrilla analysis that many of these genes belong to the family of homeobox transcription factor and are DNA binding proteins.

TABLE 2

Table 2. Gene function in biological process analysis with the PANTHER classification.

FIGURE 4

Figure 4. Gene function analysis. Gene function analysis with the PANTHER classification system and DAVID. PANTHER analysis suggested a list of DE genes that were up-regulated in osteocytes from tibias across species in biological processes of skeletal system development (A) and morphogenesis (B). (C) Functional gene annotations using DAVID suggested that 10 cross-species DE genes (FOSB, BMP7, DKK1, FGF1, FRZB, HOXA11, HOXA7, MEPE, SOST, and TNFSF10) were implicated genes in bone mineral density related diseases, which ranks the top of the list. (D) Among these DE genes, only FOSB was down-regulated in tibias in all three species.

Using the DAVID bioinformatics resources, functional gene annotations suggest that many of these DE genes are disease implicated genes. Bone mineral density ranks the top of the disease list and includes 10 shared DE genes (FOSB, BMP7, DKK1, FGF1, FRZB, HOXA11, HOXA7, MEPE, SOST, and TNFSF10) (Figure 4C), many of which associate with the Wnt signaling pathway. Among these DE genes, all but FOSB were up-regulated in tibias in all three species (Figure 4D).

Discussion

The results of our experiments show a clearly defined set of genes that are regulated differently in the calvaria and tibia across mouse, rat, and macaque. Given the different physiological functions and responses to stimuli of bone cells at these two sites, there is a compelling case that these 32 candidates account for the bulk of this site specificity, although the mechanisms by which they do this are not completely clear. The results of the comparisons of the macaque paired samples are likely to be closer to human physiology than rodents because both are primates and therefore more related. However, the analysis of differences across species allows us to identify candidates that have a high likelihood of being broadly based regulators of mammalian site-specific differences rather than those only in any one of the separate species. Specifically, our analysis shows that samples cluster according to species more than to bone sample site (Figure 1A). We can interpret this to mean that the impact of phylogenetic divergence on these samples is stronger than any site-specific differences. This difference could be due to fundamental regulatory differences in the skeleton of different species, e.g., different functional units within compact bone. However, because of species differences, we believe the genes we have identified are representative of site specific rather than other differences. Using a single species for analysis could therefore affect the power to detect important general site-specific differences in expression, as any calvaria-tibia difference would have to be larger than the gene expression divergence between species. That inference is perhaps unsurprising, but it has important implications on likely relevance of our data when compared with studies performed in only one species (10, 11).

The bioinformatic analysis of the data provided robust results within and across species. Even where a conventional differential expression (DE) analysis could not be applied to the pair of rat samples (because count dispersion could not be estimated), we were able to use a different method to assess significant DE, identifying 946 DE rat genes. The segregation of our data into patterning genes (HOXA7, HOXA11, HOXC4, HOXC8, HOXC9, HOXC10, ZIC1) and 8 others known to be related to the regulation of bone mass and architecture (FOSB, BMP7, DKK1, FGF1, FRZB, MEPE, SOST, and TNFSF10).

We show that some of these genes (i.e., HOXC8, HOXC9, HOXC10) are almost exclusively expressed in tibias, at levels >64-fold higher than in calvariae. This is consistent with other studies that suggested HOX genes are pivotal in patterning the axial and appendicular skeleton, particularly for embryonic long bone development (33–39). This is consistent with their participation in cartilage differentiation in the process of endochondral ossification which is the way that long bones grow, but is not involved in the intramembranous ossification that leads to the formation of the flat bones of the skull (39–45). The involvement of a large group of HOX transcription factors indicates that the site-specific difference in regulating response to mechanical loading by osteocytes is partly a genetically determined baseline structure since embryonic stage.

Using the DAVID bioinformatics resources to identify DE genes involved in diseases, we found bone mineral density diseases ranking the top of the list and include implicated genes: BMP7, DKK1, FGF1, FRZB, MEPE, SOST, FOSB, and TNFSF10 (or TRAIL). Most of these genes code for secreted proteins and are up-regulated in osteocytes from tibias except FOSB. This list includes not only anabolic factors that would be expected to increase bone density (BMP7, FGF1, FRZB, and FOSB) but also catabolic factors (DKK1, SOST, MEPE, TNFSF10) with roles in bone remodeling, skeletogenesis and patterning (46–52). More interestingly, many of these genes, including FRZB (53, 54), BMP (55), FGF (56), DKK1 (57, 58), SOST (59), and MEPE (60–62), have long been shown to regulate skeletal development in association or crosstalk with the canonical Wnt signaling (55, 63). The identification of these bone mineral density and Wnt signaling-related, soluble proteins suggests osteocytes at different sites response to mechanical loading differently via potential paracrine and/or endocrine agents. Upon response to strain, osteocytes regulate bone remodeling with precisely tuned balance between anabolic and catabolic effects, targeting both osteoblasts and osteoclasts. Intriguingly, Sost^−/− mice were shown to have resistance to mechanical unloading-induced bone loss and exhibited high bone mass (64), which underscores our results that Sost is less expressed in the skull (with its insensitivity to disuse) but is more expressed in the tibia which requires significant mechanical loading to maintain mass. Moreover, this is also reflected by pathological conditions caused by disrupting transcription of the SOST gene and subsequently Wnt signaling, such as sclerosteosis and the van Buchem syndrome (VBCH). Both these two closely related diseases feature generalized sclerosis, including the enlargement of the skull which endures lower physical mechanical strains (65–67). All these indicate that the canonical Wnt signaling pathway dictates skeletal site differences through regulating bone's adaptive response to mechanical loading (16, 68). This could be achieved by osteocytes through a Sost-dependent feedback mechanism of maintaining quiescent bone-lining cells which are the main source of osteoblasts in adulthood (69, 70). However, other possible pathways could not be overlooked, for example, zinc finger protein of cerebellum (Zic) family member transcription factor Zic1, one of only two down-regulated gene in tibias in our study. Although not annotated by DAVID as a bone mineral density disease implicated gene, ZIC1 has been suggested a site-specific expression pattern by osteocytes and to act as a link between mechanosensing and Wnt signaling (71). The lower expression of ZIC1 in tibias suggests a potential negative feedback pathway in osteocyte in response to mechanical stimuli.

We recognize that the lack of additional technical validation and low number of biological replicates impacts on the strength of the conclusions that can be drawn from the study and the potential targets we have identified need to be confirmed and explored using in vitro models and transgenic mice in future. We can conclude that a relatively small number of genes are differentially expressed between two skeletal sites in multiple species: results that are consistent with some but not all differences found by others exploring site-specific differences in the skeleton and our data cast new light onto possible mechanisms for those differences as well as potential future clinical benefits. The approach we have used appears to have advantages over microarray studies in a single species, allowing greater refinement of data to identify key regulators of the site-specific characteristics of the skeleton. This has been underscored by a very recent study using scRNA-Seq, suggesting that osteoblasts isolated from calvaria and cortical long bone in mice have similar transcriptomes (12), although the authors suggest that changes after isolation of the cells may have contributed to the lack of differences identified. To impact upon bone pathology, therapies to make long bones behave more like the flat bones of the skull in respect of their susceptibility to loss in response to aging, disuse, and hormonal or nutritional changes could provide powerful methods to maintain bone strength throughout life and provide a strong biological explanation for the promising development of sclerostin antibodies for the treatment of osteoporosis (72, 73).

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Ethics Statement

The animal study was reviewed and approved by the local Research Ethics Committee of the University of Sheffield (Sheffield, UK). All the procedures were performed under a British Home Office project license (PPL 40/3499) and in compliance with the UK Animals (Scientific Procedures) Act 1986.

Author Contributions

TS and GR: study design. NW and CN: study conduct and data collection. NW, NL, and TS: data analysis and interpretation. NW and TS: drafting manuscript. NW, CN, GR, and TS: approving final version of manuscript.

Funding

This research was supported by a Frontier Engineering Award (EP/K03877X/1) from the Engineering and Physical Sciences Research Council (EPSRC). SB's research was funded by the Medical Research Council (MRC) (MR/P012922/1 and MR/P023967/1).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We would like to thank Professor Stuart Baker of the Institute of Neuroscience Newcastle University for provision of the macaque bone samples.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fendo.2020.581002/full#supplementary-material

References

1. Wolpert L. Positional information and pattern formation. Curr Top Dev Biol. (2016) 117:597–608. doi: 10.1016/bs.ctdb.2015.11.008

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Gomez C, David V, Peet NM, Vico L, Chenu C, Malaval L, et al. Absence of mechanical loading in utero influences bone mass and architecture but not innervation in Myod-Myf5-deficient mice. J Anat. (2007) 210:259–71. doi: 10.1111/j.1469-7580.2007.00698.x

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Sharir A, Stern T, Rot C, Shahar R, Zelzer E. Muscle force regulates bone shaping for optimal load-bearing capacity during embryogenesis. Development. (2011) 138:3247–59. doi: 10.1242/dev.063768

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Meakin LB, Price JS, Lanyon LE. The contribution of experimental in vivo models to understanding the mechanisms of adaptation to mechanical loading in bone. Front Endocrinol. (2014) 5:154. doi: 10.3389/fendo.2014.00154

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Currey JD. Mechanical properties of bone tissues with greatly differing functions. J Biomech. (1979) 12:313–9. doi: 10.1016/0021-9290(79)90073-3

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Rubin CT, Lanyon LE. Kappa delta award paper. osteoregulatory nature of mechanical stimuli: function as a determinant for adaptive remodeling in bone. J Orthop Res. (1987) 5:300–10. doi: 10.1002/jor.1100050217

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Hillam RA, Goodship AE, Skerry TM. Peak strain magnitudes and rates in the tibia exceed greatly those in the skull: an in vivo study in a human subject. J Biomech. (2015) 48:3292–8. doi: 10.1016/j.jbiomech.2015.06.021

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Teodoro Ezequiel Guerra M, Isabel Pozzi M, Busin G, Crestana Zanetti L, Antonio Lazzarotto Terra Lopes J, Orso V. Densitometric study of the clavicle: bone mineral density explains the laterality of the fractures. Rev Bras Ortop. (2014) 49:468–72. doi: 10.1016/j.rboe.2014.07.002

PubMed Abstract | CrossRef Full Text | Google Scholar

9. DeFroda SF, Nacca C, Waryasz GR, Owens BD. Diagnosis and management of distal clavicle osteolysis. Orthopedics. (2017) 40:119–24. doi: 10.3928/01477447-20161128-03

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Kingsmill VJ, McKay IJ, Ryan P, Ogden MR, Rawlinson SC. Gene expression profiles of mandible reveal features of both calvarial and ulnar bones in the adult rat. J Dent. (2013) 41:258–64. doi: 10.1016/j.jdent.2012.11.010

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Rawlinson SC, McKay IJ, Ghuman M, Wellmann C, Ryan P, Prajaneh S, et al. Adult rat bones maintain distinct regionalized expression of markers associated with their development. PLoS ONE. (2009) 4:e8358. doi: 10.1371/journal.pone.0008358

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Ayturk UM, Scollan JP, Ayturk DG, Suh ES, Vesprey A, Jacobsen CM, et al. Single cell Rna sequencing of calvarial and long bone endocortical cells. J Bone Miner Res. (2020). doi: 10.1002/jbmr.4052. [Epub ahead of print].

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Aarden EM, Burger EH, Nijweide PJ. Function of osteocytes in bone. J Cell Biochem. (1994) 55:287–99. doi: 10.1002/jcb.240550304

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Skerry TM, Bitensky L, Chayen J, Lanyon LE. Early strain-related changes in enzyme activity in osteocytes following bone loading in vivo. J Bone Miner Res. (1989) 4:783–8. doi: 10.1002/jbmr.5650040519

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Lara-Castillo N, Kim-Weroha NA, Kamel MA, Javaheri B, Ellies DL, Krumlauf RE, et al. In vivo mechanical loading rapidly activates beta-catenin signaling in osteocytes through a prostaglandin mediated mechanism. Bone. (2015) 76:58–66. doi: 10.1016/j.bone.2015.03.019

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Galea GL, Lanyon LE, Price JS. Sclerostin's role in bone's adaptive response to mechanical loading. Bone. (2017) 96:38–44. doi: 10.1016/j.bone.2016.10.008

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Bonewald LF, Johnson ML. Osteocytes, mechanosensing and Wnt signaling. Bone. (2008) 42:606–15. doi: 10.1016/j.bone.2007.12.224

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Kilborn SH, Trudel G, Uhthoff H. Review of growth plate closure compared with age at sexual maturity and lifespan in laboratory animals. Contemp Top Lab Anim Sci. (2002) 41:21–6.

PubMed Abstract | Google Scholar

19. Zimin AV, Cornish AS, Maudhoo MD, Gibbs RM, Zhang X, Pandey S, et al. A new rhesus macaque assembly and annotation for next-generation sequencing analyses. Biol Direct. (2014) 9:20. doi: 10.1186/1745-6150-9-20

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. (2015) 12:357–60. doi: 10.1038/nmeth.3317

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Anders S, Pyl PT, Huber W. HTSeq–a python framework to work with high-throughput sequencing data. Bioinformatics. (2015) 31:166–9. doi: 10.1093/bioinformatics/btu638

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. (2010) 26:139–40. doi: 10.1093/bioinformatics/btp616

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Robinson MD, Oshlack A. A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biol. (2010) 11:R25. doi: 10.1186/gb-2010-11-3-r25

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Komljenovic A, Li H, Sorrentino V, Kutalik Z, Auwerx J, Robinson-Rechavi M. Cross-species functional modules link proteostasis to human normal aging. PLoS Comput Biol. (2019) 15:e1007162. doi: 10.1371/journal.pcbi.1007162

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Wu C, Chen X, Shu J, Lee CT. Whole-genome expression analyses of type 2 diabetes in human skin reveal altered immune function and burden of infection. Oncotarget. (2017) 8:34601–9. doi: 10.18632/oncotarget.16118

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Law CW, Chen Y, Shi W, Smyth GK. voom: precision weights unlock linear model analysis tools for RNA-seq read counts. Genome Biol. (2014) 15:R29. doi: 10.1186/gb-2014-15-2-r29

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Eden E, Navon R, Steinfeld I, Lipson D, Yakhini Z. GOrilla: a tool for discovery and visualization of enriched GO terms in ranked gene lists. BMC Bioinform. (2009) 10:48. doi: 10.1186/1471-2105-10-48

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Mi H, Huang X, Muruganujan A, Tang H, Mills C, Kang D, et al. PANTHER version 11: expanded annotation data from gene ontology and reactome pathways, and data analysis tool enhancements. Nucleic Acids Res. (2017) 45:D183–D9. doi: 10.1093/nar/gkw1138

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Huang da W, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protocol. (2009) 4:44–57. doi: 10.1038/nprot.2008.211

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci USA. (2005) 102:15545–50. doi: 10.1073/pnas.0506580102

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Liberzon A, Birger C, Thorvaldsdottir H, Ghandi M, Mesirov JP, Tamayo P. The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Syst. (2015) 1:417–25. doi: 10.1016/j.cels.2015.12.004

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Kriventseva EV, Tegenfeldt F, Petty TJ, Waterhouse RM, Simao FA, Pozdnyakov IA, et al. OrthoDB v8: update of the hierarchical catalog of orthologs and the underlying free software. Nucleic Acids Res. (2015) 43:D250–6. doi: 10.1093/nar/gku1220

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Gonzalez-Martin MC, Mallo M, Ros MA. Long bone development requires a threshold of Hox function. Dev Biol. (2014) 392:454–65. doi: 10.1016/j.ydbio.2014.06.004

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Boulet AM, Capecchi MR. Targeted disruption of hoxc-4 causes esophageal defects and vertebral transformations. Dev Biol. (1996) 177:232–49. doi: 10.1006/dbio.1996.0159

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Davis AP, Capecchi MR. Axial homeosis and appendicular skeleton defects in mice with a targeted disruption of hoxd-11. Development. (1994) 120:2187–98.

PubMed Abstract | Google Scholar

36. Carpenter EM, Goddard JM, Davis AP, Nguyen TP, Capecchi MR. Targeted disruption of Hoxd-10 affects mouse hindlimb development. Development. (1997) 124:4505–14.

PubMed Abstract | Google Scholar

37. Wahba GM, Hostikka SL, Carpenter EM. The paralogous Hox genes Hoxa10 and Hoxd10 interact to pattern the mouse hindlimb peripheral nervous system and skeleton. Dev Biol. (2001) 231:87–102. doi: 10.1006/dbio.2000.0130

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Wellik DM, Capecchi MR. Hox10 and Hox11 genes are required to globally pattern the mammalian skeleton. Science. (2003) 301:363–7. doi: 10.1126/science.1085672

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Hostikka SL, Gong J, Carpenter EM. Axial and appendicular skeletal transformations, ligament alterations, and motor neuron loss in Hoxc10 mutants. Int J Biol Sci. (2009) 5:397–410. doi: 10.7150/ijbs.5.397

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Yueh YG, Gardner DP, Kappen C. Evidence for regulation of cartilage differentiation by the homeobox gene Hoxc-8. Proc Natl Acad Sci USA. (1998) 95:9956–61. doi: 10.1073/pnas.95.17.9956

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Ye L, Fan Z, Yu B, Chang J, Al Hezaimi K, Zhou X, et al. Histone demethylases KDM4B and KDM6B promotes osteogenic differentiation of human MSCs. Cell Stem Cell. (2012) 11:50–61. doi: 10.1016/j.stem.2012.04.009

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Morgan BA, Izpisua-Belmonte JC, Duboule D, Tabin CJ. Targeted misexpression of Hox-4.6 in the avian limb bud causes apparent homeotic transformations. Nature. (1992) 358:236–9. doi: 10.1038/358236a0

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Mao L, Ding J, Zha Y, Yang L, McCarthy BA, King W, et al. HOXC9 links cell-cycle exit and neuronal differentiation and is a prognostic marker in neuroblastoma. Cancer Res. (2011) 71:4314–24. doi: 10.1158/0008-5472.CAN-11-0051

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Taher L, Collette NM, Murugesh D, Maxwell E, Ovcharenko I, Loots GG. Global gene expression analysis of murine limb development. PLoS ONE. (2011) 6:e28358. doi: 10.1371/journal.pone.0028358

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Gao RT, Zhan LP, Meng C, Zhang N, Chang SM, Yao R, et al. Homeobox B7 promotes the osteogenic differentiation potential of mesenchymal stem cells by activating RUNX2 and transcript of BSP. Int J Clin Exp Med. (2015) 8:10459–70.

PubMed Abstract | Google Scholar

46. Klein-Nulend J, Louwerse RT, Heyligers IC, Wuisman PI, Semeins CM, Goei SW, et al. Osteogenic protein (OP-1, BMP-7) stimulates cartilage differentiation of human and goat perichondrium tissue in vitro. J Biomed Mater Res. (1998) 40:614–20. doi: 10.1002/(SICI)1097-4636(19980615)40:4&<;614::AID-JBM13&>;3.0.CO;2-F

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Beederman M, Lamplot JD, Nan G, Wang J, Liu X, Yin L, et al. BMP signaling in mesenchymal stem cell differentiation and bone formation. J Biomed Sci Eng. (2013) 6:32–52. doi: 10.4236/jbise.2013.68A1004

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Bandyopadhyay A, Tsuji K, Cox K, Harfe BD, Rosen V, Tabin CJ. Genetic analysis of the roles of BMP2, BMP4, and BMP7 in limb patterning and skeletogenesis. PLoS Genet. (2006) 2:e216. doi: 10.1371/journal.pgen.0020216

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Takei Y, Minamizaki T, Yoshiko Y. Functional diversity of fibroblast growth factors in bone formation. Int J Endocrinol. (2015) 2015:729352. doi: 10.1155/2015/729352

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Su N, Jin M, Chen L. Role of FGF/FGFR signaling in skeletal development and homeostasis: learning from mouse models. Bone Res. (2014) 2:14003. doi: 10.1038/boneres.2014.3

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Vinarsky V, Krivanek J, Rankel L, Nahacka Z, Barta T, Jaros J, et al. Human embryonic and induced pluripotent stem cells express TRAIL receptors and can be sensitized to TRAIL-induced apoptosis. Stem Cells Dev. (2013) 22:2964–74. doi: 10.1089/scd.2013.0057

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Inoue D, Kido S, Matsumoto T. Transcriptional induction of FosB/DeltaFosB gene by mechanical stress in osteoblasts. J Biol Chem. (2004) 279:49795–803. doi: 10.1074/jbc.M404096200

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Enomoto-Iwamoto M, Kitagaki J, Koyama E, Tamamura Y, Wu C, Kanatani N, et al. The Wnt antagonist Frzb-1 regulates chondrocyte maturation and long bone development during limb skeletogenesis. Dev Biol. (2002) 251:142–56. doi: 10.1006/dbio.2002.0802

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Thysen S, Cailotto F, Lories R. Osteogenesis induced by frizzled-related protein (FRZB) is linked to the netrin-like domain. laboratory investigation. J Tech Methods Pathol. (2016) 96:570–80. doi: 10.1038/labinvest.2016.38

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Wang Y, Li YP, Paulson C, Shao JZ, Zhang X, Wu M, et al. Wnt and the Wnt signaling pathway in bone development and disease. Front Biosci. (2014) 19:379–407. doi: 10.2741/4214

CrossRef Full Text | Google Scholar

56. Maruyama T, Mirando AJ, Deng CX, Hsu W. The balance of WNT and FGF signaling influences mesenchymal stem cell fate during skeletal development. Sci Signal. (2010) 3:ra40. doi: 10.1126/scisignal.2000727

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Li J, Sarosi I, Cattley RC, Pretorius J, Asuncion F, Grisanti M, et al. Dkk1-mediated inhibition of Wnt signaling in bone results in osteopenia. Bone. (2006) 39:754–66. doi: 10.1016/j.bone.2006.03.017

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Lieven O, Knobloch J, Ruther U. The regulation of Dkk1 expression during embryonic development. Dev Biol. (2010) 340:256–68. doi: 10.1016/j.ydbio.2010.01.037

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Li X, Zhang Y, Kang H, Liu W, Liu P, Zhang J, et al. Sclerostin binds to LRP5/6 and antagonizes canonical Wnt signaling. J Biol Chem. (2005) 280:19883–7. doi: 10.1074/jbc.M413274200

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Nampei A, Hashimoto J, Hayashida K, Tsuboi H, Shi K, Tsuji I, et al. Matrix extracellular phosphoglycoprotein (MEPE) is highly expressed in osteocytes in human bone. J Bone Mineral Metab. (2004) 22:176–84. doi: 10.1007/s00774-003-0468-9

PubMed Abstract | CrossRef Full Text | Google Scholar

61. David V, Martin A, Hedge AM, Rowe PS. Matrix extracellular phosphoglycoprotein (MEPE) is a new bone renal hormone and vascularization modulator. Endocrinology. (2009) 150:4012–23. doi: 10.1210/en.2009-0216

PubMed Abstract | CrossRef Full Text | Google Scholar

62. Staines KA, Poulet B, Farquharson C, Pitsillides AA. Sclerostin/MEPE axis in OA: lessons from long bone development. Bone Abstracts. (2013) 1:PP27. doi: 10.1530/boneabs.1.PP27

CrossRef Full Text | Google Scholar

63. Liu F, Kohlmeier S, Wang CY. Wnt signaling and skeletal development. Cell Signal. (2008) 20:999–1009. doi: 10.1016/j.cellsig.2007.11.011

PubMed Abstract | CrossRef Full Text | Google Scholar

64. Lin C, Jiang X, Dai Z, Guo X, Weng T, Wang J, et al. Sclerostin mediates bone response to mechanical unloading through antagonizing Wnt/beta-catenin signaling. J Bone Miner Res. (2009) 24:1651–61. doi: 10.1359/jbmr.090411

PubMed Abstract | CrossRef Full Text | Google Scholar

65. Van Hul W, Balemans W, Van Hul E, Dikkers FG, Obee H, Stokroos RJ, et al. Van Buchem disease (hyperostosis corticalis generalisata) maps to chromosome 17q12-q21. Am J Hum Genet. (1998) 62:391–9. doi: 10.1086/301721

PubMed Abstract | CrossRef Full Text | Google Scholar

66. Maupin KA, Droscha CJ, Williams BO. A comprehensive overview of skeletal phenotypes associated with alterations in Wnt/beta-catenin signaling in humans and mice. Bone Res. (2013) 1:27–71. doi: 10.4248/BR201301004

PubMed Abstract | CrossRef Full Text | Google Scholar

67. Boudin E, Fijalkowski I, Hendrickx G, Van Hul W. Genetic control of bone mass. Mol Cell Endocrinol. (2016) 432:3–13. doi: 10.1016/j.mce.2015.12.021

PubMed Abstract | CrossRef Full Text | Google Scholar

68. Moustafa A, Sugiyama T, Prasad J, Zaman G, Gross TS, Lanyon LE, et al. Mechanical loading-related changes in osteocyte sclerostin expression in mice are more closely associated with the subsequent osteogenic response than the peak strains engendered. Osteoporos Int. (2012) 23:1225–34. doi: 10.1007/s00198-011-1656-4

PubMed Abstract | CrossRef Full Text | Google Scholar

69. Matic I, Matthews BG, Wang X, Dyment NA, Worthley DL, Rowe DW, et al. Quiescent bone lining cells are a major source of osteoblasts during adulthood. Stem Cells. (2016) 34:2930–42. doi: 10.1002/stem.2474

PubMed Abstract | CrossRef Full Text | Google Scholar

70. Regard JB, Zhong Z, Williams BO, Yang Y. Wnt signaling in bone development and disease: making stronger bone with Wnts. Cold Spring Harb Perspect Biol. (2012) 4:a007997. doi: 10.1101/cshperspect.a007997

PubMed Abstract | CrossRef Full Text | Google Scholar

71. Kalogeropoulos M, Varanasi SS, Olstad OK, Sanderson P, Gautvik VT, Reppe S, et al. Zic1 transcription factor in bone: neural developmental protein regulates mechanotransduction in osteocytes. FASEB J. (2010) 24:2893–903. doi: 10.1096/fj.09-148908

PubMed Abstract | CrossRef Full Text | Google Scholar

72. Bhattacharyya S, Pal S, Chattopadhyay N. Targeted inhibition of sclerostin for post-menopausal osteoporosistherapy: A critical assessment of the mechanism of action. Eur J Pharmacol. (2018) 826:39–47. doi: 10.1016/j.ejphar.2018.02.028

PubMed Abstract | CrossRef Full Text | Google Scholar

73. McClung MR. Sclerostin antibodies in osteoporosis: latest evidence and therapeutic potential. Ther Adv Musculoskelet Dis. (2017) 9:263–70. doi: 10.1177/1759720X17726744

PubMed Abstract | CrossRef Full Text | Google Scholar