Transformation of Mature Osteoblasts into Bone Lining Cells and RNA Sequencing-Based Transcriptome Profiling of Mouse Bone during Mechanical Unloading

A Ram Hong; Kwangsoo Kim; Ji Yeon Lee; Jae-Yeon Yang; Jung Hee Kim; Chan Soo Shin; Sang Wan Kim

doi:10.3803/EnM.2020.35.2.456

Endocrinol Metab > Volume 35(2); 2020 > Article

Hong, Kim, Lee, Yang, Kim, Shin, and Kim: Transformation of Mature Osteoblasts into Bone Lining Cells and RNA Sequencing-Based Transcriptome Profiling of Mouse Bone during Mechanical Unloading

Original Article

Endocrinology and Metabolism 2020;35(2):456-469.

Published online: June 30, 2020

DOI: https://doi.org/10.3803/EnM.2020.35.2.456

Transformation of Mature Osteoblasts into Bone Lining Cells and RNA Sequencing-Based Transcriptome Profiling of Mouse Bone during Mechanical Unloading

A Ram Hong^1,^*

, Kwangsoo Kim^2,^*

, Ji Yeon Lee³, Jae-Yeon Yang³, Jung Hee Kim³, Chan Soo Shin³, Sang Wan Kim³

¹Department of Internal Medicine, Chonnam National University Medical School, Gwangju, Korea

²Seoul National University Hospital Biomedical Research Institute, Seoul, Korea

³Department of Internal Medicine, Seoul National University College of Medicine, Seoul, Korea

Corresponding author: Sang Wan Kim, Department of Internal Medicine, Seoul Metropolitan Government Seoul National University Boramae Medical Center, Seoul National University College of Medicine, 20 Boramae-ro 5-gil, Dongjak-gu, Seoul 07061, Korea,
Tel: +82-2-870-2223, Fax: +82-2-870-3863, E-mail: swkimmd@snu.ac.kr

*These authors contributed equally to this work.

Received October 28, 2019 Revised March 27, 2020 Accepted April 03, 2020

This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (https://creativecommons.org/licenses/by-nc/4.0/) which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article has been corrected. See Endocrinol Metab. 2021 Nov 18; 36(6): 1314.

Abstract

Background

We investigated RNA sequencing-based transcriptome profiling and the transformation of mature osteoblasts into bone lining cells (BLCs) through a lineage tracing study to better understand the effect of mechanical unloading on bone loss.

Methods

Dmp1-CreERt2(+):Rosa26R mice were injected with 1 mg of 4-hydroxy-tamoxifen three times a week starting at postnatal week 7, and subjected to a combination of botulinum toxin injection with left hindlimb tenotomy starting at postnatal week 8 to 10. The animals were euthanized at postnatal weeks 8, 9, 10, and 12. We quantified the number and thickness of X-gal(+) cells on the periosteum of the right and left femoral bones at each time point.

Results

Two weeks after unloading, a significant decrease in the number and a subtle change in the thickness of X-gal(+) cells were observed in the left hindlimbs compared with the right hindlimbs. At 4 weeks after unloading, the decrease in the thickness was accelerated in the left hindlimbs, although the number of labeled cells was comparable. RNA sequencing analysis showed downregulation of 315 genes in the left hindlimbs at 2 and 4 weeks after unloading. Of these, Xirp2, AMPD1, Mettl11b, NEXN, CYP2E1, Bche, Ppp1r3c, Tceal7, and Gadl1 were upregulated during osteoblastogenic/osteocytic and myogenic differentiation in vitro.

Conclusion

These findings demonstrate that mechanical unloading can accelerate the transformation of mature osteoblasts into BLCs in the early stages of bone loss in vivo. Furthermore, some of the genes involved in this process may have a pleiotropic effect on both bone and muscle.

Keywords: Bone; Osteoporosis; Genes; Gene expression profiling

INTRODUCTION

The skeleton enables vertebrate animals that live on land to engage in physical activities such as locomotion. During physical activity or exercise, mechanical forces are loaded onto bones through gravity and by the contraction of skeletal muscles [1-3]. In particular, skeletal muscles are attached to bones, and thereby cause those bones to move when they contract. As a consequence, loss of muscle mass and strength reduces the mechanical forces that are exerted on bones, consequently contributing to bone loss [4,5]. Studies in humans have revealed significant changes in bone turnover markers and loss of cancellous and cortical bone as a result of mechanical unloading, as occurs during prolonged bed rest or long-term space flight [6,7]; and rodent studies have demonstrated that bone formation and muscle mass were decreased after tail suspension, whereas bone resorption was increased [8-10]. A representative clinical example of bone loss by mechanical unloading is disuse osteoporosis, which is observed in disabled patients after stroke or brain or spinal cord injury. In addition to pathological conditions, the physiological aging process is also associated with decreases in muscle mass and strength [11,12] and consequent bone loss [13,14]. However, the cellular mechanisms underlying bone loss induced by mechanical unloading have not been fully elucidated.

Bone lining cells (BLCs) are inactive osteoblasts lining the bone surface that can be recruited to form bone as a source of active osteoblasts. BLCs can be reactivated by anabolic agents. We previously showed that intermittent parathyroid hormone (PTH) or anti-sclerostin antibody treatment can convert quiescent periosteal lining cells into active osteoblasts, and our osteoblast lineage tracing studies revealed that intermittent exposure to PTH can also delay transformation of mature osteoblasts into BLCs [15-17]. Mechanical loading also might activate BLCs. Mechanical loading experiments with rats revealed that the lining cell surface area decreased after a single mechanical loading session [18]. Therefore, the idea that mechanical unloading also affects BLCs could be proposed; however, no direct data exist regarding this possibility.

Bone-muscle interactions are not simply a mechanical coupling in nature, but instead involve complex molecular and biochemical cross-talk since both bone and muscle are endocrine organs [19-21]. Moreover, both osteoblasts and myocytes differentiate from common progenitors (i.e., mesenchymal stem cells), and the incidence of osteoporosis and sarcopenia increase concomitantly with aging. Therefore, identifying pleiotropic genes that regulate both bone and muscle mass may reveal potential therapeutic targets for the treatment of these diseases [22,23].

In the present study, we investigated the transformation of mature osteoblasts into BLCs in response to mechanical unloading on bone using an osteoblast lineage tracing study. Furthermore, we present a comprehensive database of genes regulated by mechanical unloading through RNA-sequencing (RNA-seq)-based transcriptome profiling.

METHODS

Mice

To trace cells of the osteoblast lineage, we used a temporally-controlled transgene expression system [17]. Dmp1-cyclic recombinase (Cre)ERt2 mice were crossed with Rosa26R mice, in which Escherichia coli β-galactosidase expression is induced by Cre-mediated recombination using a universally expressed promoter [24]. Male and female mice were used for the studies. To assess changes in mRNA expression in response to mechanical unloading, wild-type (C57BL/6) male mice were used under the same experimental protocol. These studies were approved by the Institutional Animal Care and Use Committee of Seoul National University (approval number: SNU-160513-3).

Tamoxifen administration

Cre-mediated recombination in Dmp1-CreERt2:Rosa26R mice was induced by 4-hydroxy-tamoxifen (4-OHTam, Takeda, Osaka, Japan) injection. The 4-OHTam (2.5 mg) was dissolved in 100 μL of dimethylformamide (Thermo Fisher Scientific, Waltham, MA, USA) and this stock solution was diluted to 2.5 mg/mL in corn oil (Sigma-Aldrich, St. Louis, MO, USA). The Dmp1-CreERt2:Rosa26R mice were intraperitoneally injected with 1 mg of 4-OHTam three times a week starting from postnatal week 7.

Mechanical unloading

Hindlimb unloading was achieved by botulinum toxin injection and/or Achilles tenotomy. A 50-μL volume of toxin (2 U/100 g) was injected into the quadriceps, triceps surae, and calf muscles of the left hindlimb when the mice were at postnatal weeks 8 and 10. The same dose of normal saline was injected into the right hindlimb, which served as a non-unloaded control. Achilles tenotomy was performed on the left hindlimb at postnatal week 9. The animals were euthanized at postnatal weeks 8, 9, 10, and 12 (2 days, 1, 2, and 4 weeks, respectively, after the last 4-OHTam injection in Dmp1-CreERt2:Rosa26R mice) (Fig. 1). Mice euthanized after postnatal week 8 (without intervention) served as controls for the longitudinal analysis.

In vivo measurement of bone mass

Bone mineral density (BMD; g/cm²) and bone mineral content (BMC; g) of the whole body and hindlimb were measured in vivo by dual-energy X-ray absorptiometry (DXA) using a PIXImus scanner (GE Lunar, Madison, WI, USA) at postnatal weeks 8, 9, 10, and 12 before euthanasia. Mice were anesthetized with a 3:1 mixture of Zoletil (Virbac Laboratories, Carros, France) and Rompun (Bayer, Leverkusen, Germany) and placed prone on the platform of a PIXImus densitometer for BMD and BMC measurements.

Ex vivo measurement of bone microarchitecture by micro-computed tomography

Whole tibiae were harvested from the mice after euthanasia, fixed in formaldehyde solution for 48 hours, and then placed in 70% ethanol and stored at 4°C until imaging. The proximal tibia from each mouse was scanned using high-resolution micro-computed tomography (μCT; SkyScan 1173, Bruker microCT, Kontich, Belgium) at 90 kV and 88 μA with an isotropic voxel size of 7.1 μm using a 1.0-mm aluminum filter. For the metaphyseal tibia, a 1.5-mm section (starting 500 μm below the growth plate) was analyzed. Scanned images were reconstructed using NRecon v.1.6 software (Bruker microCT) by correcting for beam hardening and ring artifacts. Data were analyzed using a CT analyzer v.1.6 (Bruker microCT). For trabecular bone regions, bone volume/total volume (BV/TV), trabecular thickness (Tb.Th), trabecular separation (Tb.Sp), and trabecular number (Tb.N) were assessed.

Histological analysis, 5-bromo-4-chloro-3-indolyl-β-d-galactopyranoside staining, and cell thickness measurement

Lineage tracing was achieved by 5-bromo-4-chloro-3-indolyl-β-d-galactopyranoside (X-gal) staining for β-galactosidase activity and quantitative analysis of X-gal(+) cells [15,16,25]. Briefly, both femurs were dissected from each mouse and all soft tissue was removed. Each sample was rinsed twice with phosphate-buffered saline (PBS), fixed in 10% formalin for 30 minutes at 4°C, washed three times with PBS, and then stained overnight at 37.8°C in X-gal solution containing 1 mg/mL X-gal (Takeda), 5 mM potassium ferricyanide, 5 mM potassium ferrocyanide, 2 mM MgCl₂, 0.01% sodium deoxycholate, and 0.02% nonidet P-40. Samples were decalcified with buffered ethylenediaminetetraacetic acid (EDTA) for 2 weeks and then embedded and processed for paraffin sectioning. Sections were counterstained with eosin. X-gal(+) cells in eight to 16 microscopic fields at ×400 magnification from four to six comparable sections were counted. Periosteal X-gal(+) cells from each femur were analyzed with a Leica Application Suite camera (DM 2500, Leica Microsystems, Buffalo Grove, IL, USA) and associated software (LAS v.3.8).

Detection of serum levels of bone turnover markers

Blood was collected by cardiac puncture before mice were euthanized. Aliquots of serum were stored at −80°C until use. Serum levels of N-terminal propeptide of type I procollagen (P1NP) were measured by enzyme-linked immunosorbent assay (Immunodiagnostics Systems, Boldon, UK).

Statistical analysis for the lineage tracing study

All numerical data are presented as mean±standard error of the mean. We used the Student t test to compare the effects of unloading between or within groups. All statistical analyses were performed using SPSS for Windows version 21 (IBM Corp., Armonk, NY, USA). A P<0.05 was considered to indicate statistical significance.

mRNA library preparation and sequencing

For transcriptome profiling, we extracted total RNA from femoral diphyseal bones of wild-type (C57BL/6) mice at postnatal weeks 8 (n=3), 10 (n=2), and 12 (n=3). The femur was dissected immediately after sacrifice and cleaned of adherent tissues. After removing the distal and proximal epiphyses, the bone marrow was flushed out with 1 mL of RNAlater solution (Ambion, Austin, TX, USA), and RNA was extracted from the femur by homogenization in 1 mL of TRIzol reagent (Invitrogen, Carlsbad, CA, USA) using a tissue homogenizer. RNA then was extracted from the homogenate and purity and quantity were assessed with a spectrophotometer (Nanodrop 8000, Thermo Fisher Scientific). The quality of the total RNA was assessed with a Bioanalyzer (Model 2100, Agilent Technologies, Santa Clara, CA, USA). Each RNA sample with a RNA integrity number ≥7 was considered as being of sufficient quality for sequencing library preparation.

The libraries were prepared for 100-bp paired-end sequencing using a TruSeq RNA Library Preparation kit (Illumina, San Diego, CA, USA). Briefly, mRNA was purified and fragmented from 2 μg of total RNA using oligo (dT) magnetic beads. The fragmented mRNAs were synthesized as single-stranded cDNAs via random hexamer priming, which were then applied as a template for second-strand synthesis to prepare double-stranded cDNA. After sequential end repair, A-tailing, and adapter ligation, the cDNA libraries were amplified by polymerase chain reaction (PCR). Library quality was evaluated with the BioAnalyzer 2100 and quantification was performed with the KAPA library quantification kit (Kapa Biosystems, Wilmington, MA, USA) according to the manufacturer’s protocol. Following cluster amplification of denatured templates, paired-end sequencing (2×100 bp) was performed on a HiSeq2500 instrument (Illumina).

Raw read processing

Raw FASTQ files were trimmed using Trimmomatic with the following parameters: -threads 1 ILLUMINACLIP, LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36 [26]. These trimmed reads were then aligned to the University of California at Santa Cruz Mus musculus mm10 reference genome using STAR [27]. After alignment, featureCounts was used to count the reads [28], and the count estimates were normalized by upper-quartile normalization.

Selection of differentially expressed genes

Paired RNA-seq data were generated for the unloaded and control limb of each mouse. Three pairs were sequenced for postnatal weeks 8, 10, and 12. However due to the poor quality of one mouse at week 10, a pair was removed. Log₂ fold change (log₂FC [=log₂(unloaded limb expression+1)/(control limb expression+ 1)]) were calculated for all genes of each mouse. Differentially expressed gene (DEG) analysis was performed using eight log₂ FC values obtained from each gene. For convenience, the log₂ FC of each gene of mouse i in week j was denoted as LFCi,j.

Due to the presence of a different number of mice, genes with min_i_∈{1,2}(|LFC_i_,₁₀|)>1.5 and genes with min_i_∈{1,2,3}(|LFC_i_,₁₂|)>1.2 and average_i_∈{1,2,3}(|LFC_i_,₁₂|)>1.5 were identified as DEGs for week 10 and 12, respectively. In the DEG lists for weeks 10 and 12, DEGs at week 8 determined by the same criteria as week 12 were removed, because it was difficult to consider them to have been influenced by mechanical unloading. As week 10 included only two mice, and the other weeks were independent, DEGs were determined by utilizing the FC values without any other statistical selection methods, such as the (paired) t test or DESeq.

Gene set analysis and network analysis

The DEGs were subjected to enrichment analyses of gene ontology (GO) terms and Kegg pathways using Database for Annotation, Visualization and Integrated Discovery (DAVID) [29]. Significantly enriched gene sets were identified based on the Benjamini-Hochberg adjusted P value (q<0.05) of the Fisher exact test.

The protein-protein interaction (PPI) network of DEGs was imported from Cytoscape stringApp v1.5.0 (https://cytoscape.org/) [30]. The EAGLE algorithm included in ClusterViz v1.0.3 [31] was applied to identify functional modules. The subnetwork with the highest modularity value was expanded by adding its first neighbors, using CluePedia v1.5.5 [32]. Those genes in the expanded subnetwork were used for gene set analysis in ClueGO v2.5.5 [33].

Quantitative real-time PCR

Differential mRNA expression of some genes identified by RNA-seq was confirmed in MC3T3-E1, C2C12, and Ocy454 cell lines by quantitative real-time (qRT)-PCR on a LightCycler 96 system (Roche Diagnostics, Manheim, Germany). Data were normalized to the levels of glyceraldehyde 3-phosphate dehydrogenase and β-actin. Relative expression levels were calculated with the 2^(−ΔΔCT) method. The PCR primer sequences are shown in Supplemental Table S1.

RESULTS

Hindlimb unloading induced bone loss in vivo and ex vivo

We analyzed the bone phenotype of mice with hindlimb unloading and controls by measuring body weight-adjusted total body BMD and hindlimb BMC by DXA at different time points and found that the total body BMD of mice euthanized at 4 weeks after unloading was lower than that of controls (postnatal week 8) (Fig. 2A). Furthermore, hindlimb BMC was significantly reduced in the left hindlimbs compared with the right hindlimbs (non-unloaded control) in mice euthanized at postnatal weeks 9, 10, and 12 (Fig. 2B). To investigate the changes in bone microarchitecture induced by hindlimb unloading, we performed μCT analysis of trabecular bone regions in the proximal tibia. Representative three-dimensional μCT images from mice euthanized at postnatal week 12 (4 weeks after unloading) are shown in Fig. 2C. There were significant changes in the trabecular microarchitecture, in addition to bone mass. At 1 week after unloading, there were significant decreases in the proximal tibial BV/TV, Tb.Th, and Tb.N in the left hindlimbs (P=0.001, P=0.014, and P<0.001, respectively) compared with the right hindlimbs (Fig. 2D-G). Similar trends were observed at 2 and 4 weeks after unloading. At the latter time point, Tb.Sp was increased in the left hindlimbs compared with the right hindlimbs (P= 0.007), which was not the case at 1 and 2 weeks after unloading.

Unloading induced early decreases in the number and thickness of X-gal(+) cells on the periosteal surface of the femur

The 10-kb Dmp1-CreERt2 transgene is active in osteocytes and in a subset of osteoblasts on the surface of calvariae and long bones [15,16,34]. In this system, 4-OHTam administration induced CreERt2-mediated gene recombination. Since 4-OHTam activity is detected for only 24 hours post-injection [35], this inducible system allows genetic labeling of CreERt2-expressing cells at a specific time.

To evaluate the effect of mechanical unloading on the transformation of mature osteoblasts into quiescent BLCs, mice were euthanized at postnatal weeks 8 (before hindlimb unloading), 9, 10, and 12 (Fig. 1). In the mice of postnatal week 8, euthanasia was performed 2 days after the last 4-OHTam injection. Many labeled osteoblasts and osteocytes were observed in the femur by X-gal staining (postnatal week 8) (Fig. 3A). At 4 weeks after the last 4-OHTam injection, cells were visible on the periosteal surface of the femur, but were fewer in number (postnatal week 12) (Fig. 3A). The number and thickness of LacZ+ osteoblastic descendants on the periosteal surface of the femur were quantified 2 days (postnatal week 8), 1 week (week 9), 2 weeks (week 10), and 4 weeks (week 12) after the last 4-OHTam injection and compared between the left and right hindlimbs (non-unloaded control). At 1 week after unloading, there were no significant differences in the number (right: 15.5±2.2/mm; left: 9.6±1.4/mm; P=0.790) or thickness (right: 2.7±0.2 μm; left: 2.4±0.3 μm; P=0.273) of X-gal(+) cells between the left and right hindlimbs (postnatal week 9) (Fig. 3). However, at 2 weeks after unloading, a significant decrease in the number (right: 10.9±0.8/mm; left: 7.3±2.4/mm; P=0.038) and a subtle change in the thickness of X-gal(+) cells (right: 3.5±0.2 μm; left: 2.9±0.2 μm; P=0.071) were observed in the left hindlimb relative to the right hindlimb (postnatal week 10) (Fig. 3). At 4 weeks after unloading, the decrease in thickness was greater on the left side than on the right side (right: 2.7±0.2 μm; left: 2.0±0.3 μm; P=0.030), although there were comparable numbers of cells on both sides (right: 9.3±1.9/mm; left: 8.1±1.1/mm; P= 0.760) (postnatal week 12) (Fig. 3).

Unilateral hindlimb unloading decreased serum P1NP levels

To investigate whether unilateral hindlimb unloading induced changes in bone turnover, we measured the serum levels of P1NP in mice at each time point (Fig. 4). Serum P1NP levels decreased in mice euthanized at postnatal weeks 9 and 12, but increased in mice euthanized at postnatal week 10 compared with those euthanized at postnatal week 8 (postnatal week 8: 8.11±0.34 ng/mL; week 9: 2.94±0.35 ng/mL; week 10: 10.37± 0.88 ng/mL; week 12: 3.78±0.42 ng/mL; week 8 vs. 9, P<0.001; week 8 vs. 10, P=0.067; week 8 vs. 12, P<0.001; week 9 vs. 10, P<0.001; week 9 vs. 12, P=0.177; week 10 vs. 12, P<0.001).

DEG analysis

An analysis of the datasets from postnatal weeks 10 and 12 revealed 3211 DEGs (Fig. 5, Supplemental Table S2). At postnatal week 10, 99 genes were upregulated, while 1,364 were downregulated, and at postnatal week 12, 1,001 were upregulated and 1,091 were downregulated (Fig. 5A, B). We integrated the two datasets and divided the DEGs into up-up, up-down, down-down, and down-up (to illustrate, if gene A was always upregulated in both datasets, it was classified as “up-up”; if gene B was upregulated at postnatal week 10 and downregulated at postnatal week 12, it was classified as “up-down”). As shown in the Venn diagram in Fig. 5C, there were five up-up, seven up-down, 17 down-up, and 315 down-down genes (Supplemental Tables S3-6).

Gene set analysis of DEGs using DAVID and PPI network analysis

We performed a gene set analysis of the 315 down-down genes using DAVID (Fig. 6). When we separately analyzed the gene sets in terms of biological process, cellular component, molecular function, and Kegg pathway, several genes—including those involved in muscle contraction, calcium signaling, and glucose metabolism—were found to be significantly affected by mechanical unloading (Supplemental Table S7).

We constructed a final PPI network with 216 nodes and 552 edges (average node degree, 5.111; average local clustering coefficient, 0.223). A total of seven subnetworks were discovered through the EAGLE algorithm, and subnetwork 1 showed the highest modularity (modularity, 6.795). Subnetwork 1 had 56 nodes and 299 edges, with four major down-down DEGs including AMPD1 (adenosine monophosphate deaminase 1), Bche (butyrylcholinesterase), NEXN (nexilin, F-actin binding protein), and Xirp2 (encoding Xin-actin binding repeat containing 2) (Fig. 7). CluePedia-based network expansion of 56 genes in subnetwork 1 found additional 47 neighboring genes, and a total of 103 genes were used for the gene set analysis. The gene set analysis identified potential interactions among the DEGs, such as glucose metabolism and calcium signaling, which were similar to the results of the DAVID analysis (Supplemental Table S8).

mRNA expression levels of down-down genes were upregulated during osteoblastogenic/osteocytic and myogenic differentiation

Of the 315 down-down genes, we selected nine that had high fold change values and were considered to have biological functions related to bone and muscle (i.e., Xirp2, AMPD1, Mettl11b [methyl transferase-like protein 11b], NEXN, CYP2E1, Bche, Ppp1r3c [protein phosphatase 1 regulatory subunit 3c], Tceal7 [transcription elongation factor A like 7], and Gadl1 [glutamate decarboxylase-like protein 1]). To verify the RNA-seq results, we confirmed the mRNA expression levels of these genes by qRT-PCR in MC3T3-E1 and C2C12 cells. Xirp2, AMPD1, Mettl11b, NEXN, CYP2E1, Bche, Ppp1r3c, Tceal7, and Gadl1 mRNA expression was upregulated during both osteoblastogenic and myogenic differentiation (Fig. 8A, B). Since osteocytes comprise the majority of cells extracted from femur specimens, we further confirmed mRNA expression levels of these genes during osteocytic differentiation using Ocy454 osteocytic cells that expressed SOST [36]. We observed significant increases in the levels of Xirp2, AMPD1, Mettl11b, Ppp1r3c, Tceal7, and Gadl1 during osteocytic differentiation (Fig. 8C).

DISCUSSION

The present findings provide direct evidence for the effect of unloading on the fate of mature osteoblasts differentiating into BLCs. Unloading accelerated this process in a direction opposite to that of intermittent PTH or anti-sclerostin antibody treatment. Furthermore, the RNA-seq data revealed that mechanical unloading led to significant changes in the transcriptome profiles following periods of unloading compared with the non-unloaded control. We detected 315 differentially expressed protein-coding genes that were continuously downregulated during unloading. Notably, the mRNA expression of Xirp2, AMPD1, Mettl11b, NEXN, CYP2E1, Bche, Ppp1r3c, Tceal7, and Gadl1 increased in both differentiated osteoblast/osteocytic and myoblast cell lines.

We used two independent approaches—lineage tracing and RNA-seq-based transcriptome profiling—to investigate the cellular and molecular mechanisms by which mechanical unloading induces bone loss. We achieved unilateral hindlimb unloading in mice using a combination of a botulinum toxin injection and Achilles tenotomy to induce bone loss in vivo. Mature osteoblasts can either remain on the bone surface as quiescent BLCs or be embedded in the bone matrix as osteocytes. Previous lineage-tracing studies showed that the majority of LacZ+ cells (i.e., mature active osteoblasts) were plump and cuboidal at 2 days after the last 4-OHTam injection on the periosteal surface [15-17]. However, they became flat and consequently hardly visible on the bone surface over time because they were transformed into quiescent BLCs. Using a transgenic lineage-tracing model, we observed a significant decrease in the number of labeled cells on the bone surface at 2 weeks, and a decrease in their thickness at 4 weeks after hindlimb unloading compared with the non-unloaded control. These results imply that mechanical unloading accelerates the transformation of active osteoblasts into quiescent BLCs in vivo. We also observed significant decreases in serum P1NP levels following unloading, indicating that bone formation was reduced. This provides the first evidence that mechanical unloading also affects BLCs—specifically, in the opposite manner to mechanical loading. It is difficult to explain the increased P1NP level at week 10. In addition to the botulinum toxin injections at week 8 and 10, we performed Achilles tenotomy on the left hindlimbs to induce unloading at week 9. This procedure may have led to greater bone loss than occurred in response to the botulinum toxin injections, and thereby may have stimulated a compensatory mechanism in the contralateral limbs to activate osteoblasts related to the increased P1NP level. In our study, we did not evaluate unloading-induced changes in the number of osteocytes (either LacZ+ osteocytes or total osteocytes), proliferation, or apoptosis of osteoblast lineage cells. In previous studies, intermittent PTH or anti-sclerostin antibody treatment did not affect the number of labeled osteocytes in long bones [15,17].

To clarify the molecular mechanisms underlying unloading-induced bone loss, we analyzed the transcriptome profile obtained using the same unloading protocol as that in the lineage-tracing study. Since bone is a tissue with a low cell/matrix ratio, contamination by tissues with a high cell/matrix ratio (e.g., muscle or bone marrow) can lead to variations in gene expression patterns, the extent of such contamination depends on the individual mouse. To minimize contamination effects, we determined the time sequence of DEGs in paired bones from the same mouse, in which the left hindlimb was subjected to unloading while the right hindlimb served as a control. When we performed principal component analysis (Supplemental Fig. S1) of femur gene expression patterns, there was a high correlation between the unloaded and control femurs at each time point, indicating low heterogeneity in our samples. Since we ultimately analyzed relative changes in DEGs between unloaded and control limbs over time, and not absolute gene expression levels at each time point, the results of our RNA-seq analysis are reliable.

In the present study, genes that were differentially expressed between unloaded limbs and control limbs were analyzed at 2 and 4 weeks after unloading. Changes in gene expression patterns occur before changes in bone phenotypes; therefore, our experiments could not capture early responses after unloading, and a combination of a delayed response and the induction of compensatory mechanisms may have influenced our data. Moreover, hindlimb unloading was achieved not by a single procedure, but by repeating a procedure at 1-week intervals, which could complicate the interpretation of our findings. Given these limitations, we focused on gene sets that showed similar expression patterns (e.g., up-up or down-down) at 2 and 4 weeks after unloading. Indeed, most genes exhibited the down-down pattern during unloading.

To identify pleiotropic genes that potentially have biological functions in both bone and muscle, we reviewed the gene expression atlases of the Mouse Genome Database (www.informatics.jax.org) and the European Bioinformatics Institute (www.ebi.ac.uk). Among genes that are highly expressed in bone and muscle, we excluded some that are abundantly expressed in other tissues. We also reviewed previously published RNA-seq data by Ayturk et al. [37] from the tibia, muscle, and bone marrow of 10-week-old mice; although they did not perfectly match our own data, they were a useful reference for narrowing down the number of candidate genes that we searched.

We ultimately selected the following nine genes for in vitro validation by qRT-PCR and confirmed that they were upregulated during osteoblastogenic/osteocytic and myogenic differentiation: Xirp2, AMPD1, Mettl11b (methyl transferase-like protein 11b), NEXN, CYP2E1, Bche, Ppp1r3c, Tceal7, and Gadl1. Xirp2 is expressed in cardiac and skeletal muscle, and a hypomorphic Xirp2 allele is associated with cardiac hypertrophy [38]. A mutation in NEXN has been linked to hypertrophic cardiomyopathy [39]. CYP2E1 induces oxidative stress, resulting in impaired insulin action [40]. Similarly, Ppp1r3c modulates glycogen metabolism, especially glycogen synthesis [41]. AMPD1 participates in the regulation of glucose metabolism by modulating AMP-activated protein kinase activation in muscle, and AMPD1 deficiency contributes to the improvement of insulin resistance [42]. Tceal7 encodes an ovarian tumor suppressor that represses myoblast proliferation and enhances the differentiation of myoblasts into myotubes [43]. Bche, which can cleave acetylcholine, was detected in MC3T3-E1 cells and may modulate cell-matrix interactions in bone by inducing the degradation of acetylcholine in early stages of osteoblast differentiation [44]. The physiological functions of Gadl1 and Mettl11b have not yet been determined.

It remains unclear how mechanical unloading affects the gene expression patterns of active mature osteoblasts, resulting in rapid differentiation of these cells into quiescent BLCs. The best way to compare transcriptome profiling between BLCs and active mature osteoblasts is to perform single cell sequencing. However, at present, it is technically very difficult to isolate BLCs and active osteoblasts in vivo. Therefore, we cannot be certain that the gene sets affected by unloading are similar to those in mature osteoblasts and osteocytes. Further studies are needed to address this issue.

In conclusion, our results suggest possible cellular mechanisms by which mechanical unloading induces bone loss and provide gene sets for dissecting the regulatory networks that govern this process. The identified genes might exert pleiotropic effects on bone and muscle.

Supplementary Information

enm-2020-35-2-456-suppl1.pdf

enm-2020-35-2-456-suppl2.pdf

enm-2020-35-2-456-suppl3.xls

enm-2020-35-2-456-suppl4.pdf

enm-2020-35-2-456-suppl5.pdf

enm-2020-35-2-456-suppl6.xls

enm-2020-35-2-456-suppl7.pdf

enm-2020-35-2-456-suppl8.xls

enm-2020-35-2-456-suppl9.xlsx

ACKNOWLEDGMENTS

This work was supported by grants from the National Research Foundation of Korea (2012R1A1A2042211), and from Promising-Pioneering Researcher Program through Seoul National University (2016).

NOTES

CONFLICTS OF INTEREST

No potential conflict of interest relevant to this article was reported.

AUTHOR CONTRIBUTIONS

Conception or design: A.R.H., K.K., S.W.K. Acquisition, analysis, or interpretation of data: A.R.H., K.K., J.Y.L., J.Y.Y., S.W.K. Drafting the work or revising: A.R.H., K.K., S.W.K. Final approval of the manuscript: A.R.H., K.K., J.Y.L., J.Y.Y., J.H.K., C.S.S., S.W.K.

Fig. 1

Experimental design. For a lineage tracing study, Dmp1-CreERt2:Rosa26R mice were injected with 1 mg of 4-hydroxy-tamoxifen (4-OHTam) three times on postnatal week 7. To determine the fate of mature osteoblasts in vivo, animals were euthanized at postnatal weeks 8, 9, 10, and 12 (2 days, 1, 2, and 4 weeks after the last 4-OHTam treatment). Mechanical unloading was achieved by botulinum toxin and/or achilles tenotomy in the left hindlimb. The right hindlimb served as the control. RNA sequencing-based transcriptome profiling was performed with wild-type (C57BL/6) mice using the same experimental protocol.

Fig. 2

(A, B) Changes in bone mass and microarchitecture during mechanical unloading. In vivo measurements of bone mass using dual-energy X-ray absorptiometry. (A) Body weight-adjusted total body bone mineral density (BMD) and (B) hindlimb bone mineral content (BMC). (C-G) Ex vivo measurements of microarchitecture using micro-computed tomography (μCT). (C) The data are representative three-dimensional μCT images of the trabecular bone regions of proximal tibiae from mice euthanized on postnatal week 12 (4 weeks after unloading, left: control hindlimb; right: unloading hindlimb). (D) Bone volume/total volume (BV/TV) (%), (E) trabecular thickness (Tb.Th), (F) trabecular number (Tb.N), (G) trabecular separation (Tb.Sp). Data are expressed as mean±SEM (n=8 to 10 per group). ^aP<0.05.

Fig. 3

Effects of mechanical unloading on the transformation of mature osteoblasts into bone lining cells on the periosteal surface of the femur. 5-Bromo-4-chloro-3-indolyl-β-d-galactopyranoside (X-gal) staining was performed in both hindlimbs on 2 days (postnatal 8 weeks), 1 week (9 weeks), 2 weeks (10 weeks), or 4 weeks (12 weeks) after the last 4-hydroxy-tamoxifen injection, respectively. Unloading was performed on the left hindlimbs, while the right hindlimbs served as controls. (A) Data are representative of the experimental analyses performed on sections. Quantitative analysis of the number (B) and thickness (C) of X-gal(+) cells on the periosteum. Data are expressed as mean±SEM (n=4 to 6 per group). ^aP<0.05 between controls; ^bP<0.05 between unloaded hindlimbs; ^cP<0.05 between unloaded hindlimbs and controls at each time point.

Fig. 4

Effects of mechanical unloading on serum levels of N-terminal propeptide of type I procollagen (P1NP). Data are expressed as mean±SEM (n=15 per group). ^aP<0.05.

Fig. 5

Scatter plot of log₂ fold change (FC) and Venn diagram for differentially expressed gene (DEG) analysis. RNA sequencing-based analysis of differentially expressed mRNA after mechanical unloading at postnatal weeks 10 and 12 mice. DEGs were defined using the following criteria: |log₂(FC)|>log₂(1.2) and |avg(log₂(FC))|>log₂(1.5) at postnatal weeks 8 and 12; and each log₂(FC)>log₂(1.5) at postnatal week 10. Scatter plot of DEG analysis at postnatal weeks 10 (A) and 12 (B). Red dots: upregulated (Up) DEGs, blue dots: downregulated (Dn) DEGs, and gray dots: non-DEGs. (C) Venn diagram for DEG analysis. Each number indicates the number of upregulated and downregulated genes at postnatal weeks 10 and 12 (n=3/group in postnatal weeks 8 and 12; n=2/group in postnatal week 10).

Fig. 6

Gene set analysis for 315 genes downregulated at both postnatal weeks 10 and 12 using the Database for Annotation, Visualization and Integrated Discovery (DAVID) web-based tool. P values were calculated using the Fisher exact test. Significantly enriched gene sets were determined by a Benjamini-Hochberg adjusted P value (q value) below 0.05.

Fig. 7

Protein-protein interaction network analysis of differentially expressed genes (DEGs) in subnetwork 1. The STRING database was used to analyze the protein-protein interaction network based on the proteins corresponding to selected DEGs.

Fig. 8

mRNA expression levels of nine genes downregulated at both postnatal weeks 10 and 12 on (A) osteoblastogenic, (B) myogenic, and (C) osteocytic differentiation using quantitative real-time polymerase chain reaction. ^aP<0.05 between un-differentiated and differentiated cell lines.

REFERENCES

1. Lanyon LE, Hampson WG, Goodship AE, Shah JS. Bone deformation recorded in vivo from strain gauges attached to the human tibial shaft. Acta Orthop Scand 1975;46:256-68.
[CROSSREF] [PUBMED]

2. Usui T, Maki K, Toki Y, Shibasaki Y, Takanobu H, Takanishi A, et al. Measurement of mechanical strain on mandibular surface with mastication robot: influence of muscle loading direction and magnitude. Orthod Craniofac Res 2003;6(Suppl 1):163-7. 179-82.
[CROSSREF] [PUBMED]

3. Hong AR, Kim SW. Effects of resistance exercise on bone health. Endocrinol Metab (Seoul) 2018;33:435-44.
[CROSSREF] [PUBMED] [PMC]

4. Globus RK, Bikle DD, Morey-Holton E. The temporal response of bone to unloading. Endocrinology 1986;118:733-42.
[CROSSREF] [PDF]

5. Vandamme K, Holy X, Bensidhoum M, Deschepper M, Logeart-Avramoglou D, Naert I, et al. Impaired osteoblastogenesis potential of progenitor cells in skeletal unloading is associated with alterations in angiogenic and energy metabolism profile. Biomed Mater Eng 2012;22:219-26.
[CROSSREF] [PUBMED]

6. Inoue M, Tanaka H, Moriwake T, Oka M, Sekiguchi C, Seino Y. Altered biochemical markers of bone turnover in humans during 120 days of bed rest. Bone 2000;26:281-6.
[CROSSREF] [PUBMED]

7. Vico L, Collet P, Guignandon A, Lafage-Proust MH, Thomas T, Rehaillia M, et al. Effects of long-term microgravity exposure on cancellous and cortical weight-bearing bones of cosmonauts. Lancet 2000;355:1607-11.
[CROSSREF]

8. LeBlanc A, Marsh C, Evans H, Johnson P, Schneider V, Jhingran S. Bone and muscle atrophy with suspension of the rat. J Appl Physiol (1985) 1985;58:1669-75.
[CROSSREF] [PUBMED]

9. Fitts RH, Metzger JM, Riley DA, Unsworth BR. Models of disuse: a comparison of hindlimb suspension and immobilization. J Appl Physiol (1985) 1986;60:1946-53.
[CROSSREF] [PUBMED]

10. Judex S, Garman R, Squire M, Busa B, Donahue LR, Rubin C. Genetically linked site-specificity of disuse osteoporosis. J Bone Miner Res 2004;19:607-13.
[CROSSREF] [PUBMED]

11. Doherty TJ. The influence of aging and sex on skeletal muscle mass and strength. Curr Opin Clin Nutr Metab Care 2001;4:503-8.
[CROSSREF] [PUBMED]

12. Palmer IJ, Runnels ED, Bemben MG, Bemben DA. Muscle-bone interactions across age in men. J Sports Sci Med 2006;5:43-51.
[PUBMED] [PMC]

13. Raisz LG, Seeman E. Causes of age-related bone loss and bone fragility: an alternative view. J Bone Miner Res 2001;16:1948-52.
[CROSSREF] [PUBMED]

14. Raisz LG, Rodan GA. Pathogenesis of osteoporosis. Endocrinol Metab Clin North Am 2003;32:15-24.
[CROSSREF] [PUBMED]

15. Kim SW, Pajevic PD, Selig M, Barry KJ, Yang JY, Shin CS, et al. Intermittent parathyroid hormone administration converts quiescent lining cells to active osteoblasts. J Bone Miner Res 2012;27:2075-84.
[CROSSREF] [PUBMED] [PMC]

16. Jang MG, Lee JY, Yang JY, Park H, Kim JH, Kim JE, et al. Intermittent PTH treatment can delay the transformation of mature osteoblasts into lining cells on the periosteal surfaces. J Bone Miner Metab 2016;34:532-9.
[CROSSREF] [PUBMED] [PDF]

17. Kim SW, Lu Y, Williams EA, Lai F, Lee JY, Enishi T, et al. Sclerostin antibody administration converts bone lining cells into active osteoblasts. J Bone Miner Res 2017;32:892-901.
[CROSSREF] [PUBMED] [PMC]

18. Boppart MD, Kimmel DB, Yee JA, Cullen DM. Time course of osteoblast appearance after in vivo mechanical loading. Bone 1998;23:409-15.
[CROSSREF] [PUBMED]

19. Dallas SL, Prideaux M, Bonewald LF. The osteocyte: an endocrine cell... and more. Endocr Rev 2013;34:658-90.
[CROSSREF] [PDF]

20. Pedersen BK. Muscle as a secretory organ. Compr Physiol 2013;3:1337-62.
[CROSSREF] [PUBMED]

21. Brotto M, Bonewald L. Bone and muscle: interactions beyond mechanical. Bone 2015;80:109-14.
[CROSSREF] [PUBMED] [PMC]

22. Kramer I, Baertschi S, Halleux C, Keller H, Kneissel M. Mef2c deletion in osteocytes results in increased bone mass. J Bone Miner Res 2012;27:360-73.
[CROSSREF] [PUBMED]

23. Huang J, Hsu YH, Mo C, Abreu E, Kiel DP, Bonewald LF, et al. METTL21C is a potential pleiotropic gene for osteoporosis and sarcopenia acting through the modulation of the NF-κB signaling pathway. J Bone Miner Res 2014;29:1531-40.
[CROSSREF] [PUBMED] [PMC]

24. Soriano P. Generalized lacZ expression with the ROSA26 Cre reporter strain. Nat Genet 1999;21:70-1.
[CROSSREF] [PUBMED] [PDF]

25. Chung UI, Lanske B, Lee K, Li E, Kronenberg H. The parathyroid hormone/parathyroid hormone-related peptide receptor coordinates endochondral bone development by directly controlling chondrocyte differentiation. Proc Natl Acad Sci U S A 1998;95:13030-5.
[CROSSREF] [PUBMED] [PMC]

26. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 2014;30:2114-20.
[CROSSREF] [PUBMED] [PMC] [PDF]

27. Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics 2013;29:15-21.
[CROSSREF] [PUBMED] [PDF]

28. Liao Y, Smyth GK, Shi W. FeatureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics 2014;30:923-30.
[CROSSREF] [PUBMED] [PDF]

29. Dennis G Jr, Sherman BT, Hosack DA, Yang J, Gao W, Lane HC, et al. DAVID: database for annotation, visualization, and integrated discovery. Genome Biol 2003;4:P3.
[CROSSREF] [PUBMED]

30. Doncheva NT, Morris JH, Gorodkin J, Jensen LJ. Cytoscape StringApp: network analysis and visualization of proteomics data. J Proteome Res 2019;18:623-32.
[CROSSREF] [PUBMED]

31. Wang J, Zhong J, Chen G, Li M, Wu FX, Pan Y. ClusterViz: a Cytoscape app for cluster analysis of biological network. IEEE/ACM Trans Comput Biol Bioinform 2015;12:815-22.
[CROSSREF] [PUBMED]

32. Bindea G, Galon J, Mlecnik B. CluePedia Cytoscape plugin: pathway insights using integrated experimental and in silico data. Bioinformatics 2013;29:661-3.
[CROSSREF] [PUBMED] [PMC] [PDF]

33. Bindea G, Mlecnik B, Hackl H, Charoentong P, Tosolini M, Kirilovsky A, et al. ClueGO: a Cytoscape plug-in to decipher functionally grouped gene ontology and pathway annotation networks. Bioinformatics 2009;25:1091-3.
[CROSSREF] [PUBMED] [PMC] [PDF]

34. Powell WF Jr, Barry KJ, Tulum I, Kobayashi T, Harris SE, Bringhurst FR, et al. Targeted ablation of the PTH/PTHrP receptor in osteocytes impairs bone structure and homeostatic calcemic responses. J Endocrinol 2011;209:21-32.
[CROSSREF] [PUBMED] [PMC]

35. Maes C, Kobayashi T, Selig MK, Torrekens S, Roth SI, Mackem S, et al. Osteoblast precursors, but not mature osteoblasts, move into developing and fractured bones along with invading blood vessels. Dev Cell 2010;19:329-44.
[CROSSREF] [PUBMED] [PMC]

36. Wein MN, Spatz J, Nishimori S, Doench J, Root D, Babij P, et al. HDAC5 controls MEF2C-driven sclerostin expression in osteocytes. J Bone Miner Res 2015;30:400-11.
[CROSSREF] [PUBMED] [PMC]

37. Ayturk UM, Jacobsen CM, Christodoulou DC, Gorham J, Seidman JG, Seidman CE, et al. An RNA-seq protocol to identify mRNA expression changes in mouse diaphyseal bone: applications in mice with bone property altering Lrp5 mutations. J Bone Miner Res 2013;28:2081-93.
[CROSSREF] [PUBMED] [PMC]

38. McCalmon SA, Desjardins DM, Ahmad S, Davidoff KS, Snyder CM, Sato K, et al. Modulation of angiotensin II-mediated cardiac remodeling by the MEF2A target gene Xirp2. Circ Res 2010;106:952-60.
[CROSSREF]

39. Wang H, Li Z, Wang J, Sun K, Cui Q, Song L, et al. Mutations in NEXN, a Z-disc gene, are associated with hypertrophic cardiomyopathy. Am J Hum Genet 2010;87:687-93.
[CROSSREF] [PUBMED] [PMC]

40. Armoni M, Harel C, Ramdas M, Karnieli E. CYP2E1 impairs GLUT4 gene expression and function: NRF2 as a possible mediator. Horm Metab Res 2014;46:477-83.
[CROSSREF] [PUBMED] [PDF]

41. Lee SK, Moon JW, Lee YW, Lee JO, Kim SJ, Kim N, et al. The effect of high glucose levels on the hypermethylation of protein phosphatase 1 regulatory subunit 3C (PPP1R3C) gene in colorectal cancer. J Genet 2015;94:75-85.
[CROSSREF] [PUBMED] [PDF]

42. Cheng J, Morisaki H, Toyama K, Sugimoto N, Shintani T, Tandelilin A, et al. AMPD1: a novel therapeutic target for reversing insulin resistance. BMC Endocr Disord 2014;14:96.
[CROSSREF] [PUBMED] [PMC] [PDF]

43. Shi X, Garry DJ. Myogenic regulatory factors transactivate the Tceal7 gene and modulate muscle differentiation. Biochem J 2010;428:213-21.
[CROSSREF] [PUBMED] [PDF]

44. En-Nosse M, Hartmann S, Trinkaus K, Alt V, Stigler B, Heiss C, et al. Expression of non-neuronal cholinergic system in osteoblast-like cells and its involvement in osteogenesis. Cell Tissue Res 2009;338:203-15.
[CROSSREF] [PUBMED] [PDF]