Skip to main content
Advertisement
Browse Subject Areas
?

Click through the PLOS taxonomy to find articles in your field.

For more information about PLOS Subject Areas, click here.

  • Loading metrics

Effects of Age and Estrogen on Skeletal Gene Expression in Humans as Assessed by RNA Sequencing

  • Joshua N. Farr ,

    Contributed equally to this work with: Joshua N. Farr, Matthew M. Roforth

    Affiliations Division of Endocrinology, Department of Medicine, Mayo Clinic College of Medicine, Rochester, MN, 55905, United States of America, Robert and Arlene Kogod Center on Aging, Mayo Clinic College of Medicine, Rochester, MN, 55905, United States of America

  • Matthew M. Roforth ,

    Contributed equally to this work with: Joshua N. Farr, Matthew M. Roforth

    Affiliations Division of Endocrinology, Department of Medicine, Mayo Clinic College of Medicine, Rochester, MN, 55905, United States of America, Robert and Arlene Kogod Center on Aging, Mayo Clinic College of Medicine, Rochester, MN, 55905, United States of America

  • Koji Fujita,

    Affiliations Division of Endocrinology, Department of Medicine, Mayo Clinic College of Medicine, Rochester, MN, 55905, United States of America, Robert and Arlene Kogod Center on Aging, Mayo Clinic College of Medicine, Rochester, MN, 55905, United States of America

  • Kristy M. Nicks,

    Affiliations Division of Endocrinology, Department of Medicine, Mayo Clinic College of Medicine, Rochester, MN, 55905, United States of America, Robert and Arlene Kogod Center on Aging, Mayo Clinic College of Medicine, Rochester, MN, 55905, United States of America

  • Julie M. Cunningham,

    Affiliation Department of Experimental Pathology and Laboratory Medicine, Mayo Clinic College of Medicine, Rochester, MN, 55905, United States of America

  • Elizabeth J. Atkinson,

    Affiliation Division of Biomedical Statistics and Informatics, Department of Health Sciences Research, Mayo Clinic College of Medicine, Rochester, MN, 55905, United States of America

  • Terry M. Therneau,

    Affiliation Division of Biomedical Statistics and Informatics, Department of Health Sciences Research, Mayo Clinic College of Medicine, Rochester, MN, 55905, United States of America

  • Louise K. McCready,

    Affiliation Division of Endocrinology, Department of Medicine, Mayo Clinic College of Medicine, Rochester, MN, 55905, United States of America

  • James M. Peterson,

    Affiliation Division of Endocrinology, Department of Medicine, Mayo Clinic College of Medicine, Rochester, MN, 55905, United States of America

  • Matthew T. Drake,

    Affiliations Division of Endocrinology, Department of Medicine, Mayo Clinic College of Medicine, Rochester, MN, 55905, United States of America, Robert and Arlene Kogod Center on Aging, Mayo Clinic College of Medicine, Rochester, MN, 55905, United States of America

  • David G. Monroe,

    Affiliations Division of Endocrinology, Department of Medicine, Mayo Clinic College of Medicine, Rochester, MN, 55905, United States of America, Robert and Arlene Kogod Center on Aging, Mayo Clinic College of Medicine, Rochester, MN, 55905, United States of America

  • Sundeep Khosla

    khosla.sundeep@mayo.edu

    Affiliations Division of Endocrinology, Department of Medicine, Mayo Clinic College of Medicine, Rochester, MN, 55905, United States of America, Robert and Arlene Kogod Center on Aging, Mayo Clinic College of Medicine, Rochester, MN, 55905, United States of America

Abstract

Precise delineation of the specific genes and pathways altered with aging and estrogen (E) therapy may lead to new skeletal biomarkers and the development of novel bone therapeutics. Previous human bone studies, however, have been limited by only examining pre-specified genes and pathways. High-throughput RNA sequencing (RNAseq), on the other hand, offers an unbiased approach to examine the entire transcriptome. Here we present an RNAseq analysis of human bone samples, obtained from iliac crest needle biopsies, to yield the first in vivo interrogation of all genes and pathways that may be altered in bone with aging and E therapy in humans. 58 healthy women were studied, including 19 young women (mean age ± SD, 30.3 ± 5.4 years), 19 old women (73.1 ± 6.6 years), and 20 old women treated with 3 weeks of E therapy (70.5 ± 5.2 years). Using generally accepted criteria (false discovery rate [q] < 0.10), aging altered a total of 678 genes and 12 pathways, including a subset known to regulate bone metabolism (e.g., Notch). Interestingly, the LEF1 transcription factor, which is a classical downstream target of the Wnt/β-catenin signaling pathway, was significantly downregulated in the bones from the old versus young women; consistent with this, LEF1 binding sites were significantly enriched in the promoter regions of the differentially expressed genes in the old versus young women, suggesting that aging was associated with alterations in Wnt signaling in bone. Further, of the 21 unique genes altered in bone by E therapy, the expression of INHBB (encoding for the inhibin, beta B polypeptide), which decreased with aging (by 0.6-fold), was restored to young adult levels in response to E therapy. In conclusion, our data demonstrate that aging alters a substantial portion of the skeletal transcriptome, whereas E therapy appears to have significant, albeit less wide-ranging effects. These data provide a valuable resource for the potential identification of novel biomarkers associated with age-related bone loss and also highlight potential pathways that could be targeted to treat osteoporosis.

Trial Registration

ClinicalTrials.gov NCT02349113

Introduction

Aging is the single largest risk factor for bone loss in both sexes [1]. While virtually all current therapies target osteoclast-mediated bone resorption, age-related bone loss results, in large part, from a defect in the number and/or function of osteoblasts—the cells within basic multicellular units (BMUs) responsible for forming new bone. Thus, reflecting the age-related defect in bone formation, histologically measured mean wall thickness, a measure of the work done by osteoblasts in BMUs, declines with age in both sexes [2]. However, while serum bone formation markers steadily decline with age in men [3,4], they generally increase in older women [4]. This is because the marked estrogen (E) deficiency in postmenopausal women leads to an overall increase in bone turnover, resulting in more BMUs, even though there is a relative reduction in bone formation at the cellular level [5]. Ultimately, this imbalance between bone resorption and formation leads to net bone loss. Therefore, impaired bone formation is a hallmark of age-related bone loss in both sexes.

Despite this understanding, directly identifying the underlying mechanisms for impaired bone formation with aging and E deficiency in humans is a significant gap in knowledge. Accordingly, the results of studies aimed at these issues may lead to novel approaches to prevent or reverse age-related bone loss. In addition, such studies may lead to the identification of new skeletal biomarkers to permit better targeting of therapies to individual patients. However, the specific genes and pathways in human bone that are regulated by aging or E remain unclear. These genes and pathways must be defined more precisely in order to develop novel therapeutic approaches to combat age-related bone loss.

To address this issue, we have developed and validated an approach to obtain and analyze small needle bone biopsies (1–2 mm diameter) from the posterior iliac crest of humans [6,7]. Using this approach, we have to date obtained bone samples from 60 women, including 20 young women as well as 40 old women (20 per group) receiving either no therapy or 3 weeks of short-term E therapy. Previously, we coupled this methodology to customized, in-house quantitative polymerase chain reaction (QPCR) analyses of nearly 300 genes related to bone metabolism in this cohort of women [6,7]. A limitation of these studies, however, was that we only examined pre-specified pathways and genes using QPCR. High-throughput RNA sequencing (RNAseq), on the other hand, offers an unbiased approach to examine the entire transcriptome. Here we present a high-throughput RNAseq analysis of our previously characterized human bone samples [6,7] to yield the first in vivo interrogation of all potential genes and pathways in bone that may be altered with aging and in response to E therapy in women.

Materials and Methods

Study subjects

This study was approved by the Mayo Clinic Institutional Review Board (IRB), and full informed written consent was obtained from all subjects. As described previously [6,7], we recruited a total of 60 healthy women, including 20 young women (mean age ± SD, 30.0 ± 5.4 years, range 22 to 40 years) as well as 40 old women (n = 20 per group) receiving either no therapy (72.9 ± 6.5 years, 65 to 88 years) or 3 weeks of E therapy (70.5 ± 5.2 years, 65 to 88 years). All old women were ≥65 years of age and postmenopausal (i.e., absence of menses >1 year; serum follicle stimulating hormone [FSH] >20 IU/L), whereas all young women were premenopausal. Old women receiving E were treated with 100 μg/d of transdermal 17β-estradiol (Alora, Watson Pharma Inc., Corona, CA); patches were changed every 3–4 days for the 3 week study duration. All potential subjects were rigorously screened for coexisting disease and excluded if they had low body stores of vitamin D (serum total 25-hydroxyvitamin D [25(OH)D] of <20 ng/ml), any clinical abnormalities in other serum biochemistries (i.e., calcium, phosphorus, alkaline phosphatase, aspartate transaminase, creatinine, parathyroid hormone [PTH], thyroid stimulating hormone [TSH], and FSH), or any disorders associated with altered skeletal structure or function. This included the presence of chronic renal impairment, chronic liver disease, severe neuropathic disease, unstable cardiovascular disease, malignancy, chronic gastrointestinal disease, hypo- or hyperparathyroidism, acromegaly, Cushing’s syndrome, hypopituitarism, severe chronic obstructive pulmonary disease, alcoholism, or diabetes. In addition, subjects with a history of pathological fractures (e.g., due to Paget’s disease, myeloma, metastatic malignancy) were excluded. Further, subjects were excluded if undergoing treatment with any medication that can affect bone metabolism such as corticosteroids (>3 months at any time or >10 days within the previous year), anticonvulsant therapy (within the previous year), pharmacological doses of thyroid hormone (causing TSH to decline below normal), adrenal or anabolic steroids, aromatase inhibitors, calcitonin, calcium supplementation >1200 mg/d (within the preceding 3 months), bisphosphonates, denosumab, E or selective E receptor modulator (within the past year), sodium fluoride, or teriparatide. Clinical details in the medical records were reviewed to determine if subjects met study criteria. These details were confirmed during an interview with each subject.

The young and untreated old women were studied under a common IRB protocol. The women treated with E were studied under a separate IRB protocol; since this group received an intervention, it is now considered a clinical trial. As such, this component of the study was registered at ClinicalTrials.gov (registration number NCT02349113). The study was not registered before the enrollment of patients started because when this study enrolled patients (2011–2012), the requirements for registration of all studies with an intervention, even for non-therapeutic (i.e., mechanistic) studies, were less stringent. The authors do confirm that all ongoing and related trials for this drug/intervention are registered. Subjects were recruited and studied between February 1, 2011 and September 13, 2012. A flowchart for enrollment and allocation of the subjects in the E arm of the study is provided in Fig 1.

thumbnail
Fig 1. Flow diagram for enrollment and allocation of the subjects in the E arm of the study.

https://doi.org/10.1371/journal.pone.0138347.g001

Study protocol

All studies were performed and samples collected in the outpatient Clinical Research Unit at the Mayo Clinic (Rochester, Minnesota, USA) with study subjects maintaining their usual diet and calcium intake. After an overnight fast, morning blood was drawn from potential candidates and serum was analyzed for determination of serum biochemistries (see Study subjects). Blood samples were stored at -80°C. Screening procedures also included a medical history and physical examination. Anthropometric data were collected on all subjects wearing light-weight clothing and no shoes. Weight was obtained using an electronic scale (Model 5002, Tronic, Inc., White Plains, NY) and height was measured using a customized stadiometer (Mayo Section of Engineering). Body mass index (BMI, kg/m2) was calculated as the ratio of weight to height2. Small needle biopsies of bone were obtained from all subjects as detailed below.

Obtaining and processing needle biopsies of bone

Using an approach slightly modified from that used clinically by hematologists to obtain bone marrow aspirates and biopsies, we obtained small needle bone biopsies from the posterior iliac crest of all subjects using an 8G needle under local anesthesia (1% lidocaine) and monitored intravenous (IV) sedation (1–3 mg of IV midazolam and 50–100 μg of fentanyl), as described previously [6,7]. All 60 biopsies were performed without any complications. Each biopsy is 1–2 mm in width and ~1–2 cm in length, and contains both cortical and trabecular bone (see [8] for a μCT analysis of a similar bone sample). Of note, we found that the standard practice of snap freezing and storing the bone samples for subsequent homogenization results in significant RNA degradation, with higher variability in RNA integrity number (RIN) values. Thus, we immediately homogenized (Tissue Tearor, Cole-Parmer) the bone samples in lysis buffer (QIAzol; QIAGEN) and subsequently stored them at -80°C until the time of further processing. This procedure permitted isolation of intact RNA (RIN values consistently >9) from 58 of the 60 (19 young, 19 old, and 20 E therapy) bone biopsy samples; these 58 samples were included in the RNAseq analysis of this report.

RNA isolation and whole transcriptome RNA sequencing (RNAseq)

As described previously [6,7,8], total RNA was isolated using the RNeasy Micro Kit (Qiagen, Valencia, CA) and treated with the Turbo DNA-freeTM Kit (Life Technologies, Grand Island, NY) to remove potential contaminating DNA that may lead to false-positive amplification. RNA quality and purity was confirmed with a Nanodrop spectrophotometer (Thermo Scientific). The resulting intact RNA was used for either RNAseq analysis or QPCR analyses (see below). RNAseq was performed using RNA from the bone samples, using methods as extensively described previously from our laboratory [7]. The only exception to this was that because ample total RNA was obtained from the bone biopsies, there was no need to amplify the mRNA species, as done previously [8].

Briefly, first strand cDNA was generated from ~100 ng of total RNA using DNA/RNA chimeric primers and reverse transcriptase, creating a cDNA/RNA hybrid, followed by second strand cDNA synthesis containing a DNA/RNA duplex. The resulting double-stranded cDNA products were modified by random priming and extension to create double-stranded products suitable for generating sequencing libraries. The double-stranded products then underwent blunt-end repair. Adapter molecules were ligated to the 5’ and 3’ ends of each fragment to facilitate PCR amplification of the fragments to produce the final library. Unique indices were created for each sample and incorporated at the adaptor ligation step for loading multiple samples per flow cell. Three distinct libraries were loaded per flow cell and sequenced on an Illumina HiSeq 2000 using TruSeq SBS sequencing software (version 3) and SCS data collection software (version 1.4.8). Base calling was performed using Illumina RTA (version 1.12.4.2). The RNAseq data associated with this manuscript are accessible through GEO Series accession number GSE72815.

QPCR gene expression analyses

QPCR primers were designed using the Primer Express® program (version 3.0.1, Applied Biosystems) and are available on request. QPCR reactions were run using the ABI Prism 7900HT Real time System (Applied Biosystems) with SYBR Green (Qiagen) as the detection method. Normalization for variations in input RNA was performed based on a panel of 10 housekeeping genes (18S, G6PDH, GAPDH, GUSB, L13A, RPII, TBP, TUBA1B, Β2M, ACTB) from which the 3 most stable housekeeping genes were selected using the geNorm algorithm (http://medgen.ugent.be/~jvdesomp/genorm/) [9,10], while the PCR Miner algorithm [11] was used to correct for variations in amplification efficiencies.

Serum biochemical assays

Total 25(OH)D (within assay CV of 2.4% and between assay CV of 6.8%) was measured using liquid chromatography-tandem mass spectrometry (API 5000; Applied Biosystems-MDS Sciex, Foster City, CA). Creatinine was measured by isotope dilution mass spectroscopy; the estimated glomerular filtration rate (eGFR) was calculated using the Modification of Diet in Renal Disease (MDRD) equation [12]. Bone formation was assessed by serum amino-terminal propeptide of type I collagen (PINP; μg/L) as measured by double-antibody radioimmunoassay (within assay CV of 7.5%, between assay CV of 6.5%; DiaSorin, Stillwater, MN). Bone resorption was evaluated by serum cross-linked C-telopeptide of type I collagen (CTx; ng/mL) using an enzyme-linked immunosorbent assay (ELISA, within assay CV of 4.6%, between assay CV of 8.0%; Roche Diagnostics, Indianapolis, IN). Serum sclerostin levels were measured using two different techniques: 1) by ELISA (pg/mL, within- and between-assay CVs of <6%; Biomedica, Alpco, Salem, NH); and 2) using an electrochemiluminescense assay (pg/mL, within assay CV of 6%, between assay CV of 10%; Meso-Scale Discoveries [MSD] 96-Well MULTI-ARRAY Human Sclerostin Assay, Gaithersburg, MD), as previously described [13]. The MSD assay has a detection limit of 1 pg/mL (with a range of 1 to 10,000 pg/mL), and rigorous epitope mapping has demonstrated that this assay detects the full intact sclerostin molecule [14].

Statistical analyses

Clinical characteristics and serum biochemistries were compared between groups using t-tests or ANOVA as appropriate. The QPCR analysis of the bone biopsies obtained from this cohort has been described in detail [6,7]. As previously described [8], the RNAseq data were analyzed using several steps: alignment, quality control, obtaining genomic features per sample and summarizing the data across samples for subsequent comparisons between samples (see below). Pair-end reads were aligned by TopHat 2.0.6 [15] against the hg19 genome build using the bowtie1 aligner option [16]. Initial assessment of quality control was made using RSeQC software [17] to estimate the distance between paired end reads, evaluate the sequencing depth for alternate splicing analysis, determine the rate of duplicate reads, evaluate coverage of reads across genes, etc. Gene counts were generated using HTseq software (http://www-huber.embl.de/users/anders/HTSeq/doc/overview.html) and gene annotation files were obtained from Illumina (http://cufflinks.cbcb.umd.edu/igenomes.html). Genes with median gene counts of <10 in both young and old groups were considered “non-expressed”, as previously described [8]. Conditional quantile normalization using the cqn Bioconductor package [18] was applied to the RNAseq data to reduce variability introduced by GC content, gene size, and total gene counts per sample. Significance in differential expression was determined using the edgeR Bioconductor package [19] assuming a negative binomial error structure accounting for the cqn derived offset. We used a correction for the false discovery rate (FDR, q-value) in order to account for multiple comparisons, specifying a q-value of < 0.10, as is generally accepted [16,17].

Pathway enrichment analysis was conducted using the Ingenuity Pathway Analysis software (IPA, Ingenuity Systems, Inc., Redwood City, CA, USA; http://www.ingenuity.com) to identify significant canonical pathways in which the Differentially Expressed Genes (DEGs) in the tested samples were enriched. The IPA program applies Fisher’s exact test to calculate a p-value that represents the probability of the DEGs in the pathway being found together due to random chance. The Benjamini-Hochberg FDR (q < 0.10) correction was also applied in order to account for multiple comparisons in the IPA. In addition, IPA was used to determine networks associated with aging in the tested samples based on DEGs. For network generation, a dataset with gene identifiers and corresponding fold change values was uploaded into IPA. Default settings were used to identify differentially regulated molecules. These molecules were overlaid onto a global molecular network (contained in the Ingenuity Knowledge Base) and algorithmically connected. The functional analysis in IPA determined the biological functions and diseases associated with each network.

The transcription factor binding site analysis was conducted using the C3 module (motif gene sets) from the Molecular Signatures Database (MSigDB; http://www.broadinstitute.org/gsea/msigdb/collections.jsp#C3) [20]. This analysis identifies common cis-regulatory DNA motifs that are statistically enriched in the query gene list. The transcription factor binding sites are defined in the TRANSFAC (version 7.4; http://www.gene-regulation.com) database. These motifs are catalogued and represent known or likely DNA regulatory elements contained within the DNA sequence surrounding the transcriptional start site of a particular gene of interest [21]. Upstream regulator analysis, a component of the IPA software package, was performed as described in the IPA software documentation.

Results

Descriptive characteristics

Clinical characteristics of the young, old, and E groups are shown in Table 1. In this cohort of healthy women, the premenopausal women were, as expected, much younger (mean age ± SD, 30.3 ± 5.4 years) than both groups (old and E therapy) of postmenopausal women. By contrast, the old (age 73.1 ± 6.6 years) and E (age 70.5 ± 5.2 years) groups did not differ in age. Whereas the young women were taller than both groups of old women, the old and E groups were similar in height. In addition, none of the groups differed in weight or BMI. All women had sufficient levels of total 25(OH)D and creatinine, levels of which were in the normal range and did not differ among the groups. By contrast, eGFR was significantly higher in young as compared to both groups of old women. At baseline, the young women tended to have lower serum levels of CTx (a marker of bone resorption), but similar serum levels of PINP (a marker of bone formation), as compared to both groups of old women. Further, as compared to the old non-treated women, PINP and CTx levels were not different at baseline in the old E-treated women. As previously described [7], short-term E therapy markedly increased serum PINP (by 55.4%, p < 0.001), reflecting the effects of E in enhancing bone formation; this increase is followed by a subsequent decrease in PINP levels [7] due to the coupling of bone formation and bone resorption, which is uniformly reduced by E therapy. The latter is reflected by the decrease in serum CTx (by 28.6%, p < 0.001) levels relative to baseline after 3 weeks of E treatment (Table 1). By contrast, neither parameter changed in the non-treated old women over the same time course.

thumbnail
Table 1. Clinical characteristics and biochemical markers in the young, old, and estrogen (E)-treated women.

https://doi.org/10.1371/journal.pone.0138347.t001

High-throughput RNAseq of bone biopsies

We next extended our previous QPCR analyses of the bone biopsies [6,7] to include identification of genes and pathways altered with aging and in response to E therapy using high-throughput RNAseq to compare the skeletal transcriptome of the young versus old and old versus E-treated women. Using fairly stringent criteria (median count ≥ 10, FDR [q] < 0.10), 678 genes were altered with aging (in old relative to young subjects), including 446 up- and 232 down-regulated genes, respectively (S1 Table). Additional adjustment for height and eGFR, two variables that differed between the young and old subjects, resulted in essentially no change in the findings (data not shown). By contrast, using the same criteria, only 21 genes were altered by E therapy (in E-treated relative to non-treated old subjects), including 18 up- and 3 down-regulated genes, respectively (S2 Table).

QPCR confirmation of RNAseq results

Among the genes in S1 and S2 Tables, a number have been previously shown to be altered in human bone biopsies with aging [6] or E therapy [7] using QPCR, implying that our RNAseq analysis was accurate and reliable. In order to validate this further, we examined the association of gene expression levels between RNAseq and QPCR in a randomly selected subset of genes significantly altered with aging (q < 0.10 by RNAseq, n = 46). As evident in S1 Fig, a very strong correlation (r = 0.95, p < 0.001) was observed for gene expression fold changes between RNAseq and QPCR.

Ingenuity Pathway Analysis (IPA)

In order to interrogate the specific canonical cellular pathways altered with skeletal aging, we performed IPA on the young versus old RNAseq dataset. It should be noted, however, that IPA was not run on the old versus E treatment dataset due to the comparatively fewer genes altered in response to E therapy. The IPA analysis identified 73 canonical pathways that were altered with aging at the p < 0.05 level (S3 Table). After applying a q < 0.10, a total of 12 pathways remained significant (Table 2). Although a number of pathways that were altered with aging at the p < 0.05 level did not meet FDR criteria (q < 0.10), several of these pathways have known functions in bone metabolism and therefore warrant consideration as potential targets to treat osteoporosis; Table 3 lists 10 cellular pathways that fell in this category.

thumbnail
Table 2. Pathways altered (q < 0.10) in the young versus old RNAseq dataset based on pathway analysis using the Ingenuity Pathway Analysis software (see Statistical analyses).

The associated p-value, false discovery rate (q), and ratio for each pathway, as well as each gene fold change (in parenthesis) are provided.

https://doi.org/10.1371/journal.pone.0138347.t002

thumbnail
Table 3. Selected pathways of interest with known function in bone in the young versus old RNAseq dataset based on pathway analysis using the Ingenuity Pathway Analysis software (see Statistical analyses).

The associated p-value, false discovery rate (q), and ratio for each pathway, as well as each gene fold change (in parenthesis) are provided.

https://doi.org/10.1371/journal.pone.0138347.t003

Particularly noteworthy, and consistent with our previous study using QPCR [6], was the significant alteration (q = 0.021) of the Notch signaling pathway with aging in the bone biopsies (Table 2). Indeed, expression levels of key regulatory genes from the Notch pathway including NOTCH3 and NOTCH4 (both single-pass transmembrane cell surface receptors), as well as JAG2 and DLL4 (two Notch receptor ligands that cause proteolytic cleavage of Notch receptors upon binding) and two critical direct downstream Notch pathway targets (HES1 and HEY1) were significantly up-regulated in bone biopsies of old versus young subjects (Fig 2A). It is also noteworthy that the Wnt/β-catenin signaling pathway was altered with aging, although as noted above, this pathway fell into the category of bone-related pathways with p < 0.05, but q > 0.10 (Table 3). In this pathway, transcriptional levels of the genes for SOX4, SFRP5, and LEF1 were decreased, whereas levels for SOX17, SOX7, FZD4, CDH3, SFRP1, SOX18, SOX13, and CCND1 were all increased in bone biopsies of old versus young women (Fig 2B). By contrast, however, there were no changes in any of the lipoprotein receptor-related proteins previously linked to bone (e.g., LRP4, LRP5, and LRP6) or the associated DKK family of Wnt antagonists with aging. Finally, as evident in Table 3, additional pathways with known roles in bone metabolism, including Inhibition of Matrix Metalloproteases, eNOS signaling, Nitric Oxide Signaling in the Cardiovascular System, PDGF Signaling, Gαi Signaling, and Oncostatin M Signaling were all altered in bone biopsies of old versus young subjects (at least at the p < 0.05 level, but only the first 2 at the q < 0.10 level).

thumbnail
Fig 2. (A) Notch pathway and (B) Wnt/β-catenin pathway genes are altered with aging in human bone biopsies.

Genes in the Notch pathway and Wnt/β-catenin pathway are significantly altered in old relative to young women revealed by RNAseq based on pathway analysis using the Ingenuity Pathway Analysis software (see Statistical analyses). Values are presented as median fold changes (95% CIs) for old relative to young subjects. p < 0.01; p < 0.001.

https://doi.org/10.1371/journal.pone.0138347.g002

Network analyses

We next utilized IPA to identify networks based on the 678 network-eligible molecules that were altered (median count ≥ 10, q < 0.10) in bone with aging S4 Table lists all identified networks, as well as their associated scores and network molecules). The top 10 highest scored networks are listed in Fig 3A. Particularly noteworthy, genes identified among the highest scored networks (Fig 3B) included LEF1, which plays an important role in Wnt/β-catenin signaling [22,23], as well as CDKN1A (p21), a key mediator of cellular senescence [24]. Other molecules implicated in bone metabolism identified among the highest scored networks were Notch signaling components (NOTCH3, NOTCH4, HEY1, JAG2, and DLL4) (Fig 3C) and SEMA3A; the latter regulates bone metabolism via sensory innervations [25].

thumbnail
Fig 3. Networks derived from age-related differentially expressed genes (DEGs) in bone.

(A) Top 10 scoring networks and associated network functions derived from the 678 DEGs altered (median count ≥ 10, q < 0.10) in the young versus old dataset determined using the Ingenuity Pathway Analysis software (see Statistical Analysis). (B) Network of cell morphology, hematological system development and function, and protein synthesis (score = 33); key molecules implicated in bone metabolism include LEF1 and CDKN1A (p21). (C) Network of cardiovascular system development and function, cellular development, and organismal development (score = 31); key molecules implicated in bone metabolism include Notch signaling components (NOTCH3, NOTCH4, HEY1, JAG2, and DLL4) and SEMA3A. Green indicates downregulated genes, and red indicates upregulated genes.

https://doi.org/10.1371/journal.pone.0138347.g003

Identification of enriched cis-regulatory DNA motifs

In addition to the IPA-derived canonical pathways associated with the DEGs from the young versus old RNAseq dataset, it is also of interest to identify common cis-regulatory DNA elements that are enriched with aging, which then can provide clues regarding shared molecular mechanisms involved in the regulation of these DEGs. Therefore, in order to identify putative cis-regulatory elements that are enriched in the promoter regions of the 678 DEGs from the young versus old RNAseq dataset, we queried the C3 module (motif gene sets) from the Molecular Signatures Database [MSigDB; [20]]. This analysis identified statistically enriched cis-regulatory DNA sequences contained within 2-kilobases (kb) surrounding the transcription start site for each gene (e.g. from -2kb to +2kb). Table 4 summarizes the top 20 cis-regulatory DNA motifs based on p-value (a complete list of the top 100 motifs can be found in S5 Table). Particularly noteworthy was the identification of two related DNA motifs for the LEF1 transcription factor, a classical downstream target of the Wnt/β-catenin signaling pathway [22,23]. This is consistent with our IPA and QPCR analyses (Table 3 and Fig 2B, respectively) and suggests that aging is associated with alterations in Wnt signaling in bone. The specific DEGs from the young versus old RNAseq dataset containing each of the LEF1 cis-regulatory motifs are shown in S6 Table.

thumbnail
Table 4. Top 20 transcription factor binding sites enriched in the DEGs from the young versus old RNAseq dataset.

The ratio (the percentage of genes in the overlap between the 678 DEGs and the MSigDB gene set), P-value and false discovery rate (FDR) are provided.

https://doi.org/10.1371/journal.pone.0138347.t004

Upstream Regulator Analytic (URA) tool

Another interesting feature of the IPA software package is the Upstream Regulator Analytic (URA) function. This tool enables identification of putative upstream regulators based on the DEGs identified in the canonical pathway analysis, and predicts whether these regulators are activated or inhibited based on the directionality of the DEGs. The URA tool reports an Activation z-score, where z-scores ≥2 are considered activated and z-scores ≤-2 are considered inhibited. Table 5 summarizes the top 10 activated and inhibited upstream regulators in the young versus old RNAseq dataset. The top activated upstream regulator was the growth factor TGFβ1. Interestingly, the top repressed upstream regulator was SMAD7, which is an inhibitory SMAD involved in ression of TGFβ signaling [26]. A complete list of upstream regulators with z-scores ≥2 or ≤-2 is provided in S7 Table.

thumbnail
Table 5. Top 10 activated and repressed upstream regulators in the young versus old dataset.

https://doi.org/10.1371/journal.pone.0138347.t005

Bone and serum sclerostin levels

Since we have previously reported age [6] and E [7] effects on sclerostin (SOST) mRNA levels by QPCR in these women, we next examined SOST gene expression levels in the RNAseq dataset, as well as serum levels of sclerostin in the young, old, and E-treated subjects using two different sclerostin assays. Consistent with our previous QPCR data [6], [7], there was no effect of age on SOST mRNA levels by RNAseq (fold change in the old versus young = 0.99, p = 0.981), whereas E treatment did significantly reduce SOST mRNA levels (fold change in the E-treated versus untreated old women = 0.59, p = 0.026), although this change did not meet FDR criteria (q = 0.46). Interestingly, although serum sclerostin levels were significantly higher in the old versus young subjects using the Biomedica assay, they were no different in the old versus young subjects using the MSD assay (Fig 4A). By contrast, consistent with both the QPCR and RNAseq data, serum sclerostin levels using both assays were reduced in old women by E therapy (Fig 4B). Thus, the MSD assay more accurately reflected changes in bone SOST mRNA levels with aging and E treatment, showing concordant changes in circulating sclerostin levels under both conditions.

thumbnail
Fig 4. Effects of age and estrogen (E) on serum sclerostin levels.

Serum sclerostin levels in the young versus old subjects (A) and in the old versus E-treated subjects (B) by the Biomedica and Meso Scale Discovery (MSD) assays. Data are mean ± SEM; note the difference in scales for the two sclerostin assays. *p < 0.05; p < 0.01.

https://doi.org/10.1371/journal.pone.0138347.g004

Genes altered by E therapy

Finally, we focused on the 21 genes in bone tissue altered in old women in response to E therapy as revealed by RNAseq (S2 Table). Particularly noteworthy were the marked increases in expression of LGR5 (by 3.9-fold) and PPARGC1A (by 3.7-fold), while the expression of C19orf80 (coding for Betatrophin) decreased (by 0.5-fold) in response to E therapy. Furthermore, the expression of INHBB (inhibin, beta B), which significantly decreased with aging (Fig 5A; q < 0.001), was restored to young adult levels (Fig 5B; q < 0.001) in response to E therapy. Correspondingly, expression of INHBA (inhibin, beta A) followed the same pattern with aging (Fig 5A; p < 0.01; q = 0.118) and in response to E therapy (Fig 5B; p < 0.001; q = 0.179), although these gene expression changes did not meet FDR criteria (q < 0.10).

thumbnail
Fig 5. Effects of age and estrogen (E) on gene expression of INHBA and INHBB.

Bone INHBB and INHBA gene expression levels by RNAseq in the old relative to young women (A) and in the E-treated relative to untreated old women (B). Data are shown as median fold changes (95% CIs). p < 0.01; p < 0.001.

https://doi.org/10.1371/journal.pone.0138347.g005

Discussion

It is important to more precisely define the specific genes and pathways in human bone that are altered with aging and E deficiency for several reasons. First, direct identification of the underlying mechanisms for impaired bone formation with aging in humans is a significant gap in knowledge. Second, such studies may lead to the development of new skeletal biomarkers, which could result in the identification of circulation-based biomarkers to define patients most at risk for bone loss. Third, because virtually all current skeletal therapies target osteoclast-mediated bone resorption, the results of such studies may lead to novel approaches to stimulate bone formation in patients and prevent or reverse age-related bone loss. However, what has been lacking is an unbiased assessment of the potential genes and pathways regulated by aging and E in human bone.

Based on these considerations, the genes identified in this report using RNAseq that were altered in human bone with aging and in response to short-term E therapy may represent novel skeletal biomarkers. Furthermore, our results demonstrate that aging significantly altered a total of 12 canonical pathways in bone tissue, including a subset known to regulate bone metabolism (e.g., Notch). Interestingly, a number of the identified pathways and genes were previously unrecognized as changing in bone with aging. Perhaps unexpectedly, many of these pathways and genes have critical functions in other tissues and disease processes, although their role(s) in the pathogenesis of skeletal fragility requires additional investigation.

Consistent with our previous report using QPCR [6], our RNAseq analysis revealed that genes associated with Notch signaling were significantly altered in bone with aging. Based on transgenic mouse models, some genetic mutations in individual Notch surface receptors or certain Notch receptor ligands can be embryonically lethal [27,28,29], whereas others have no apparent skeletal effects [30,31,32,33,34]. More recently, Notch signaling has been shown to impair osteoblast maturation and differentiation through both HES and HEY family members, which inhibit Runx2 transcriptional activity through direct physical interactions, as demonstrated in various in vitro culture systems (e.g., MC3T3 cells [35], COS7 cells [36], and both CHO and ST2 cells [37]). Moreover, mutant mice overexpressing HEY1 develop an osteopenic phenotype [38]. Thus, given the crucial functions of Notch signaling in osteoblast maturation and differentiation, our finding that the Notch signaling pathway was altered in human bone with aging and specifically, that levels of both HES1 and HEY1 were significantly up-regulated in bone of old relative to young subjects, indicate that increased Notch signaling may contribute to impaired bone formation with aging. Moreover, these findings suggest that bone formation in vivo may be stimulated by targeted ression of the Notch signaling pathway [37].

Although originally discovered in rare human mutations affecting bone [39,40,41,42,43,44], mouse genetics has confirmed the importance of canonical Wnt signaling in the regulation of bone homeostasis [45]. Overall, activation of this pathway leads to increased bone formation, whereas inhibition leads to reduced bone formation. Thus, components of the Wnt pathway are currently targets of therapeutic interventions to restore bone mass and strength in aging humans. Perhaps surprisingly, we did not find that any of the genes encoding for the 19 known mammalian Wnt proteins were altered in human bone with aging or in response to E therapy. Similarly, none of the lipoprotein receptor-related proteins previously linked to bone (e.g., LRP4, LRP5, and LRP6) or the associated DKK family of Wnt antagonists were regulated in aging bone. Interestingly, however, another antagonist of Wnt signaling, SFRP1, was significantly up-regulated in bone of old relative young women. Capable of forming a complex with Frizzled receptors to thereby inhibit Wnt signaling [46], SFRP1 up-regulation may play a causal role in ressing bone formation with aging. Importantly, LEF1 mRNA levels were significantly reduced with aging, and LEF1 was found among the highest scored IPA networks (Fig 3B). In addition, two related DNA motifs for the LEF1 transcription factor were identified among the top cis-regulatory DNA elements enriched with aging. Since LEF1 is a classical Wnt target gene [22,23,47], these findings are consistent with reduced Wnt action in bone with aging, perhaps related to increased SFRP1 expression.

Mounting data in both animals and humans suggest that inhibitors of sclerostin (encoded by SOST), an osteocyte-secreted soluble Wnt/β-catenin signaling pathway antagonist, can stimulate bone formation. Using the Biomedica assay, previous studies in humans have demonstrated that serum levels of sclerostin increase with advancing age [6,48]. However, a concern with the Biomedica sclerostin assay is that it uses an anti-sclerostin antibody that may detect differential forms or proteolytic fragments of sclerostin [49]. By contrast, the MSD assay uses an antibody which, through rigorous epitope mapping, appears to detect only the full intact sclerostin molecule [14]. Thus, the Biomedica assay may overestimate true circulating sclerostin levels, which may at least in part explain why we found that sclerostin levels were, on average, 5- to 7-fold higher across groups of women when measured with the Biomedica assay as compared to the MSD assay. Importantly, these findings are consistent with a previous report [50]. Further, in contrast to the Biomedica assay results, we found no significant effect of aging on circulating sclerostin levels according to the MSD assay, which essentially mirrored the lack of changes in human bone SOST mRNA levels in the young versus old women (by both QPCR and RNAseq). In addition, according to both sclerostin assays, we found that E therapy in postmenopausal women reduced circulating sclerostin levels, which is not only consistent with previous studies [7,51], but also with the E-induced decrease in human bone SOST mRNA levels (by both QPCR and RNAseq). Collectively, these findings suggest that changes in serum sclerostin levels with aging and in response to E therapy, at least according to the MSD assay, reflect changes that occur in bone at the tissue level.

In addition to changes with aging, another goal of the present study was to utilize RNAseq analysis of human bone samples to interrogate potential gene targets that may be altered in vivo in response to E therapy. This analysis yielded a total of 21 genes (3 down- and 18 up-regulated, respectively) that were altered in E-treated women versus non-treated old women of a similar age. Interestingly, the expression C19orf80 (coding for Betatrophin), which has roles in β-cell proliferation [52] and autophagy [53], decreased in response to E therapy. Further, it is noteworthy that the single most up-regulated gene (by 3.9-fold) in bone with E therapy, the orphan receptor LGR5, has been implicated in Wnt-β-catenin signaling [54], while the second most up-regulated gene (by 3.7-fold), PPARGC1A, is involved in mitochondrial biogenesis [55], antioxidant defense [56], and also integrates the mammalian clock with energy metabolism [57]. Thus, these E-regulated genes have critical functions in other non-skeletal tissues and disease processes, although elucidating their role(s) in the pathogenesis of bone fragility will require further work.

Our RNAseq analysis of human bone samples also revealed that INHBB expression decreased with aging but was restored to young adult levels in response to E therapy. It should be noted that the INHBB-encoded protein product can either homodimerize to form the activin B protein complex or heterodimerize with the INHBA-encoded protein product to form the activin AB protein complex [58]. While activins are known to be present at high levels in bone and likely play a major role in the regulation of bone metabolism [59], existing data are conflicting in that the effects of activins on bone formation and resorption are different among species and in vitro/in vivo models [60]. These findings may be explained by differences in cellular context, whereby the net effect of altered expression of mRNAs that encode for activins in bone, and its change in response to E, depends on the relative expression of activin pathway signaling components (e.g., receptors, follistatin antagonist, intracellular SMADs). Nevertheless, targeted inhibition of activin signaling in vivo using soluble activin type II receptors has been shown to enhance bone formation and inhibit bone resorption in both animals [61,62,63] and humans [64]. Thus, the potential biological significance of the decrease in INHBB mRNA levels with age and the restoration of these levels with E therapy remain to be more clearly defined.

One potential concern regarding RNAseq is validity of the technique. Beyond the fact that a number of the genes in bone that were significantly altered by RNAseq with aging or in response to E therapy were are also altered when assessed by QPCR [6,7], our direct comparison between RNAseq and QCPR showed a very strong correlation between the techniques, which provides additional ort for the high accuracy and reliability of RNAseq. Hence, we believe that RNAseq is a valid approach for transcriptional profiling. Given the advantages of RNAseq in that it provides an unbiased approach to profiling the entire transcriptome, application of this technique enabled us to identify several sets of genes and pathways altered in human bone biopsies with aging and in response to E therapy that were previously unrecognized. We do, however, acknowledge a limitation of our approach in that the changes we observed may not represent causality. Thus, some of the altered signaling pathways in bone with aging may not be the cause of bone loss, or alternatively could reflect compensatory changes that occur with aging. Clearly, further studies are needed to define the potential causal roles of the observed altered pathways in mediating bone loss.

A further limitation of our study is that the bone biopsies contain heterogeneous populations of cells, including bone marrow elements. As such, we do not know in which precise population(s) the observed changes in gene expression patterns are occurring, although we are in the process of untangling this complexity by using highly purified populations of osteoblasts and osteocytes from human bone biopsy samples [8]. This may also explain why we observed relatively few changes in the biopsies following E therapy, despite the major effects of E on bone metabolism. It is possible that changes occurring in bone cells (osteoblasts, osteocytes, or osteoclasts) following E therapy may have been dampened due to the heterogeneous nature of the biopsy samples, and future studies using more enriched cell populations may show a larger number of regulated genes in response to E. Alternatively, in future studies, RNAseq could be performed on marrow alone followed by subtracting from the bone plus marrow results, as previously described [65]. We should note, however, that the dramatic effects of E on bone (increase in PINP and decrease in CTX) may not necessarily require changes in a large number of genes. For example, we found that E reduced both bone sclerostin mRNA levels (by RNAseq and QPCR) and serum sclerostin levels (using 2 different commercial assays). This change in sclerostin production could well account for most of the effects of E on bone, given the marked effects of the sclerostin antibody in bone turnover and bone mass in humans [66], including very recent evidence that, in addition to inhibiting bone formation, sclerostin may also lead to an increase in bone resorption due to stimulation of RANKL production by osteocytes [67]. However, the identification of the further downstream gene expression changes following sclerostin inhibition by E will likely require more refined techniques, as noted above.

In conclusion, our findings demonstrate that aging alters a subset of the skeletal transcriptome, whereas E therapy appears to have significant, albeit less wide-ranging effects. To date, our results provide the most comprehensive assessment of the effects of aging and E therapy on the human skeleton’s gene expression profile, and serve as a valuable resource for the identification of novel biomarkers associated with age-related bone loss to better target therapies to specific patients. In addition, our study validates, in humans, several pathways associated with age-related bone loss in various animal models. Thus, the coordinated alteration of these pathways with aging emphasizes their potential importance in the pathogenesis of osteoporosis in humans. However, these signaling pathways are integrated with other networks and this overlap leads to enormous complexity which remains to be understood. Notwithstanding, interventions that modulate these pathways may lead to increased bone formation in humans and improved skeletal health.

Supporting Information

S1 Fig. QPCR confirms RNAseq data in human bone biopsies.

Spearman’s correlation (r = 0.95) showing strong association of gene expression levels between RNAseq and QPCR in a subset of genes significantly (p < 0.05, q < 0.10) altered with aging (by RNAseq, n = 46).

https://doi.org/10.1371/journal.pone.0138347.s001

(TIF)

S1 Protocol. Mayo Clinic Institutional Review Board approved protocol for young versus old study.

https://doi.org/10.1371/journal.pone.0138347.s002

(PDF)

S2 Protocol. Mayo Clinic Institutional Review Board approved protocol for estrogen study.

https://doi.org/10.1371/journal.pone.0138347.s003

(PDF)

S1 Table. Genes regulated by aging in human bone.

Complete gene list of the genes altered based on RNAseq analysis (median count ≥ 10, p < 0.05, false discovery rate [q] < 0.10) in human bone samples of young women relative to old women.

https://doi.org/10.1371/journal.pone.0138347.s004

(XLSX)

S2 Table. Genes regulated by short-term estrogen (STE) therapy in human bone.

Complete gene list of the genes altered based on RNAseq analysis (median count ≥ 10, p < 0.05, false discovery rate [q] < 0.10) in bone samples of old women treated with STE therapy relative to old non-treated women.

https://doi.org/10.1371/journal.pone.0138347.s005

(XLSX)

S3 Table. Canonical pathways regulated by aging in human bone.

Complete list of 73 significant (p < 0.05) pathways altered in human bone samples with aging by Ingenuity Pathway Analysis. Ratio refers to the ratio of genes altered in bone with aging divided by the total number of genes in that pathway.

https://doi.org/10.1371/journal.pone.0138347.s006

(XLSX)

S4 Table. Networks regulated by aging in human bone.

Complete list of 25 networks altered in human bone samples with aging by Ingenuity Pathway Analysis, as well as corresponding diseases/functions, scores, and molecules.

https://doi.org/10.1371/journal.pone.0138347.s007

(XLSX)

S5 Table. Transcription factor binding sites regulated by aging in human bone.

List of top 100 transcription factor binding sites of genes differentially expressed in the young versus. old dataset.

https://doi.org/10.1371/journal.pone.0138347.s008

(XLSX)

S6 Table. Genes associated with LEF1 transcription factor binding sites.

Lists of genes associated with two separate LEF1 transcription factor binding sites: 1) TTGTTT (n = 93); and 2) CTTTGA (n = 60), as well as the 20 genes common to both lists.

https://doi.org/10.1371/journal.pone.0138347.s009

(XLSX)

S7 Table. Complete list of upstream regulators.

List of all upstream regulators of genes altered with aging (median count ≥10, q < 0.10) in the young versus old dataset.

https://doi.org/10.1371/journal.pone.0138347.s010

(XLSX)

Acknowledgments

We thank the women for their participation in this study. This work was orted in part by an investigator initiated grant from Merck and by P01 AG004875 from the National Institute of Aging (NIA). It was also made possible by the Mayo Clinical Center for Clinical and Translational Science (CCaTS), Grant Number UL1 TR000135 from the National Center for Advancing Translational Science (NCATS), a component of the National Institutes of Health (NIH). Its contents are solely the responsibility of the authors and do not necessarily represent the official view of NIH. The funders had no role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript. There are no other relevant declarations. Dr. Farr is orted by a career development award from the Robert and Arlene Kogod Center on Aging (Mayo Clinic). Dr. Khosla had full access to all the data in the study and takes responsibility for the integrity of the data and the accuracy of the data analysis.

Author Contributions

Conceived and designed the experiments: SK MMR JNF DGM JMC. Performed the experiments: MMR KF LKM MTD. Analyzed the data: EJA JMP JNF SK. Wrote the paper: JNF MMR KF KMN EJA TMT LKM JMP MTD DGM SK.

References

  1. 1. Khosla S (2013) Pathogenesis of age-related bone loss in humans. J Gerontol A Biol Sci Med Sci 68: 1226–1235. pmid:22923429
  2. 2. Lips P, Courpron P, Meunier PJ (1978) Mean wall thickness of trabecular bone packets in the human iliac crest: changes with age. Calcif Tissue Res 26: 13–17. pmid:737547
  3. 3. Fatayerji D, Eastell R (1999) Age-related changes in bone turnover in men. J Bone Miner Res 14: 1203–1210. pmid:10404022
  4. 4. Khosla S, Melton LJ 3rd, Atkinson EJ, O'Fallon WM, Klee GG, Riggs BL (1998) Relationship of serum sex steroid levels and bone turnover markers with bone mineral density in men and women: A key role for bioavailable estrogen. J Clin Endocrinol Metab 83: 2266–2274. pmid:9661593
  5. 5. Riggs BL, Khosla S, Melton LJ 3rd (2002) Sex steroids and the construction and conservation of the adult skeleton. Endocr Rev 23: 279–302. pmid:12050121
  6. 6. Roforth MM, Fujita K, McGregor UI, Kirmani S, McCready LK, Peterson JM, et al. (2014) Effects of age on bone mRNA levels of sclerostin and other genes relevant to bone metabolism in humans. Bone 59: 1–6. pmid:24184314
  7. 7. Fujita K, Roforth MM, Demaray S, McGregor U, Kirmani S, McCready LK, et al. (2014) Effects of estrogen on bone mRNA levels of sclerostin and other genes relevant to bone metabolism in postmenopausal women. J Clin Endocrinol Metab 99: E81–88. pmid:24170101
  8. 8. Fujita K, Roforth MM, Atkinson EJ, Peterson JM, Drake MT, McCready LK, et al. (2014) Isolation and characterization of human osteoblasts from needle biopsies without in vitro culture. Osteoporos Int 25: 887–895. pmid:24114401
  9. 9. Radonic A, Thulke S, Mackay IM, Landt O, Siegert W, Nitsche A (2004) Guideline to reference gene selection for quantitative real-time PCR. Biochem Biophys Res Commun 313: 856–862. pmid:14706621
  10. 10. Vandesompele J, De Preter K, Pattyn F, Poppe B, Van Roy N, De Paepe A, et al. (2002) Accurate normalization of real-time quantitative RT-PCR data by geometeric averaging of multiple internal control genes. Genome Biol 3: research0034.0031-0030-0034.0011. pmid:12184808
  11. 11. Zhao S, Fernald RD (2005) Comprehensive algorithm for quantitative real-time polymerase chain reaction. J Comput Biol 12: 1047–1064. pmid:16241897
  12. 12. Levey AS, Coresh J, Greene T, Stevens LA, Zhang YL, Hendriksen S, et al. (2006) Using standardized serum creatinine values in the modification of diet in renal disease study equation for estimating glomerular filtration rate. Ann Intern Med 145: 247–254. pmid:16908915
  13. 13. van Lierop AH, Witteveen JE, Hamdy NA, Papapoulos SE (2010) Patients with primary hyperparathyroidism have lower circulating sclerostin levels than euparathyroid controls. Eur J Endocrinol 163: 833–837. pmid:20817762
  14. 14. van Lierop AH, Hamdy NAT, Hamersma H, van Bezooijen RL, Power J, Loveridge N, et al. (2011) Patients with sclerosteosis and disease carriers: human models of the effect of sclerostin on bone turnover. J Bone Miner Res 26: 2804–2811. pmid:21786318
  15. 15. Trapnell C, Pachter L, Salzberg SL (2009) TopHat: discovering splice junctions with RNA-Seq. Bioinformatics 25: 1105–1111. pmid:19289445
  16. 16. Langmead B, Trapnell C, Pop M, Salzberg SL (2009) Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol 10: R25. pmid:19261174
  17. 17. Wang L, Wang S, Li W (2012) RSeQC: quality control of RNA-seq experiments. Bioinformatics 28: 2184–2185. pmid:22743226
  18. 18. Hansen KD, Irizarry RA (2012) Removing technical variability in RNA-seq data using conditional quantile normalization. Biostatistics 13: 204–216. pmid:22285995
  19. 19. Robinson MD, McCarthy DJ, Smyth GK (2010) edgeR: a bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26: 139–140. pmid:19910308
  20. 20. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. (2005) Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci USA 102: 15545–15550. pmid:16199517
  21. 21. Xie X, Lu J, Kulbokas EJ, Golub TR, Mootha V, Lindblad-Toh K, et al. (2005) Systematic discovery of regulatory motifs in human promoters and 3' UTRs by comparison of several mammals. Nature 434: 338–345. pmid:15735639
  22. 22. Giese K, Amsterdam A, Grosschedl R (1991) DNA-binding properties of the HMG domain of the lymphoid-specific transcriptional regulator LEF-1. Gene Dev 5: 2567–2578. pmid:1752444
  23. 23. Love JJ, Li X, Case DA, Giese K, Grosschedl R, Wright PE (1995) Structural basis for DNA bending by the architectural transcription factor LEF-1. Nature 376: 791–795. pmid:7651541
  24. 24. Baker DJ, Weaver RL, van Deursen JM (2013) p21 both attenuates and drives senescence and aging in BubR1 progeriod mice. Cell Rep 3: 1164–1174. pmid:23602569
  25. 25. Fukuda T, Takeda S, Xu R, Ochi H, Sunamura S, Sato T, et al. (2013) Sema3A regulates bone-mass accrual through sensory innervations. Nature 497: 490–493. pmid:23644455
  26. 26. Hayashi H, Abdollah S, Qiu Y, Cai J, Xu YY, Grinnell BW, et al. (1997) The MAD-related protein Smad7 associates with the TGFbeta receptor and functions as an antagonist of TGFbeta signaling. Cell 89: 1165–1173. pmid:9215638
  27. 27. Hamada Y, Kadokawa Y, Okabe M, Ikawa M, Coleman JR, Tsujimoto Y (1999) Mutation in ankyrin repeats of the mouse Notch2 gene induces early embryonic lethality. Development 126: 3415–3424. pmid:10393120
  28. 28. Gale NW, Dominguez MG, Noguera I, Pan L, Hughes V, Valenzuela DM, et al. (2004) Haploinsufficiency of delta-like 4 ligand results in embryonic lethality due to major defects in arterial and vascular development. Proc Natl Acad Sci U S A 101: 15949–15954. pmid:15520367
  29. 29. Xue Y, Gao X, Lindsell CE, Norton CR, Chang B, Hicks C, et al. (1999) Embryonic lethality and vascular defects in mice lacking the Notch ligand Jagged1. Hum Mol Genet 8: 723–730. pmid:10196361
  30. 30. Krebs LT, Xue Y, Norton CR, Shutter JR, Maguire M, Sundberg JP, et al. (2000) Notch signaling is essential for vascular morphogenesis in mice. Genes Dev 14: 1343–1352. pmid:10837027
  31. 31. Swiatek PJ, Lindsell CE, del Amo FF, Weinmaster G, Gridley T (1994) Notch1 is essential for postimplantation development in mice. Genes Dev 8: 707–719. pmid:7926761
  32. 32. Domenga V, Fardoux P, Lacombe P, Monet M, Maciazek J, Krebs LT, et al. (2004) Notch3 is required for arterial identity and maturation of vascular smooth muscle cells. Genes Dev 18: 2730–2735. pmid:15545631
  33. 33. Hrabe de Angelis M, McIntyre J 2nd, Gossler A (1997) Maintenance of somite borders in mice requires the Delta homologue DII1. Nature 386: 717–721. pmid:9109488
  34. 34. Jiang R, Lan Y, Chapman HD, Shawber C, Norton CR, Serreze DV, et al. (1998) Defects in limb, craniofacial, and thymic development in Jagged2 mutant mice. Genes Dev 12: 1046–1057. pmid:9531541
  35. 35. Zamurovic N, Cappellen D, Rohner D, Susa M (2004) Coordinated activation of notch, Wnt, and transforming growth factor-beta signaling pathways in bone morphogenic protein 2-induced osteogenesis. Notch target gene Hey1 inhibits mineralization and Runx2 transcriptional activity. J Biol Chem 279: 37704–37715. pmid:15178686
  36. 36. Garg V, Muth AN, Ransom JF, Schluterman MK, Barnes R, King IN, et al. (2005) Mutations in NOTCH1 cause aortic valve disease. Nature 437: 270–274. pmid:16025100
  37. 37. Hilton MJ, Tu X, Wu X, Bai S, Zhao H, Kobayashi T, et al. (2008) Notch signaling maintains bone marrow mesenchymal progenitors by suppressing osteoblast differentiation. Nat Med 14: 306–314. pmid:18297083
  38. 38. Salie R, Kneissel M, Vukevic M, Zamurovic N, Kramer I, Evans G, et al. (2010) Ubiquitous overexpression of Hey1 transcription factor leads to osteopenia and chondrocyte hypertrophy in bone. Bone 46: 680–694. pmid:19857617
  39. 39. Gong Y, Slee RB, Fukai N, Rawadi G, Roman-Roman S, Reginato AM, et al. (2001) LDL receptor-related protein 5 (LRP5) affects bone accrual and eye development. Cell 107: 513–523. pmid:11719191
  40. 40. Little RD, Carulli JP, Del Mastro RG, Dupuis J, Osborne M, Folz C, et al. (2002) A mutation in the LDL receptor-related protein 5 gene results in the autosomal dominant high-bone-mass trait. Am J Hum Genet 70: 11–19. pmid:11741193
  41. 41. Boyden LM, Mao J, Belsky J, Mitzner L, Farhi A, Mitnick MA, et al. (2002) High bone density due to a mutation in LDL-receptor-related protein 5. N Engl J Med 346: 1513–1521. pmid:12015390
  42. 42. Brunkow ME, Gardner JC, Van Ness J, Paeper BW, Kovacevich BR, Proll S, et al. (2001) Bone dysplasia sclerosteosis results from loss of the SOST gene product, a novel cystine knot-containing protein. Am J Hum Genet 68: 577–589. pmid:11179006
  43. 43. Balemans W, Patel N, Ebeling M, Van Hul E, Wuyts W, Lacza C, et al. (2002) Identification of a 52 kb deletion downstream of the SOST gene in patients with van Buchem disease. J Med Genet 39: 91–97. pmid:11836356
  44. 44. Loots GG, Kneissel M, Keller H, Baptist M, Chang J, Collette NM, et al. (2005) Genomic deletion of a long-range bone enhancer misregulates sclerostin in Van Buchem disease. Genome Res 15: 928–935. pmid:15965026
  45. 45. Baron R, Kneissel M (2013) WNT signaling in bone homeostasis and disease: from human mutations to treatments. Nat Med 19: 179–192. pmid:23389618
  46. 46. Chim CS, Pang R, Fung TK, Choi CL, Liang R (2007) Epigenetic dysregulation of Wnt signaling pathway in multiple myeloma. Leukemia 21: 2527–2536. pmid:17882284
  47. 47. Khosla S, Westendorf JJ, Oursler MJ (2008) Building bone to reverse osteoporosis and repair fractures. J Clin Invest 118: 421–428. pmid:18246192
  48. 48. Modder UI, Hoey KA, Amin S, McCready LK, Achenbach SJ, Riggs BL, et al. (2011) Relation of age, gender, and bone mass to circulating sclerostin levels in women and men. J Bone Miner Res 26: 373–379. pmid:20721932
  49. 49. Clarke BL, Drake MT (2013) Clinical utility of serum sclerostin measurements. BoneKEy Reports 2: 361. pmid:24578825
  50. 50. Durosier C, van Lierop A, Ferrari S, Chevalley T, Papapoulos S, Rizzoli R (2013) Association of circulating sclerostin with bone mineral mass, microstructure, and turnover biochemical markers in healthy elderly men and women. J Clin Endocrinol Metab 98: 3873–3883. pmid:23864703
  51. 51. Modder UIL, Clowes JA, Hoey K, Peterson JM, McCready L, Oursler MJ, et al. (2011) Regulation of circulating sclerostin levels by sex steroids in women and men. J Bone Miner Res 26: 27–34. pmid:20499362
  52. 52. Yi P, Park JS, Melton DA (2013) Betatrophin: a hormone that controls pancreatic beta cell proliferation. Cell 153: 747–758. pmid:23623304
  53. 53. Tseng YH, Ke PY, Liao CJ, Wu SM, Chi HC, Tsai CY, et al. (2014) Chromosome 19 open reading frame 80 is upregulated by thyroid hormone and modulates autophagy and lipid metabolism. Autophagy 10: 20–31. pmid:24262987
  54. 54. Carmon KS, Gong X, Lin Q, Thomas A, Liu Q (2011) R-spondins function as ligands of the orphan receptors LGR4 and LGR5 to regulate Wnt/beta-catenin signaling. Proc Natl Acad Sci U S A 108: 11452–11457. pmid:21693646
  55. 55. Puigserver P, Wu Z, Park CW, Graves R, Wright M, Spiegelman BM (1998) A cold-inducible coactivator of nuclear receptors linked to adaptive thermogenesis. Cell 92: 829–839. pmid:9529258
  56. 56. St-Pierre J, Drori S, Uldry M, Silvaggi JM, Rhee J, Jager S, et al. (2006) Suppression of reactive oxygen species and neurodegeneration by the PGC-1 transcriptional coactivators. Cell 127: 397–408. pmid:17055439
  57. 57. Liu C, Li S, Liu T, Borjigin J, Lin JD (2007) Transcriptional coactivator PGC-1alpha integrates the mammalian clock and energy metabolism. Nature 447: 477–481. pmid:17476214
  58. 58. Makanji Y, Zhu J, Mishra R, Holmquist C, Wong WP, Schwartz NB, et al. (2014) Inhibin at 90: from discovery to clinical application, a historical review. Endoc Rev 35: 747–794.
  59. 59. Lotinun S, Pearsall RS, Horne WC, Baron R (2012) Activin receptor signaling: a potential therapeutic target for osteoporosis. Curr Mol Pharmacol 5: 195–204. pmid:21787285
  60. 60. Nicks KM, Perrien DS, Akel NS, Suva LJ, Gaddy D (2009) Regulation of osteoblastogenesis and osteoclastogenesis by the other reproductive hormones, Activin and Inhibin. Mol Cell Endocrinol 310: 11–20. pmid:19615428
  61. 61. Pearsall RS, Canalis E, Cornwall-Brady M, Underwood KW, Haigis B, Ucran J, et al. (2008) A soluble activin type IIA receptor induces bone formation and improves skeletal integrity. Proc Natl Acad Sci USA 105: 7082–7087. pmid:18460605
  62. 62. Lotinun S, Pearsall RS, Davies MV, Marvell TH, Monnell TE, Ucran J, et al. (2010) A soluble activin receptor Type IIA fusion protein (ACE-011) increases bone mass via a dual anabolic-antiresorptive effect in Cynomolgus monkeys. Bone 46: 1082–1088. pmid:20080223
  63. 63. Fajardo RJ, Manoharan RK, Pearsall RS, Davies MV, Marvell T, Monnell TE, et al. (2010) Treatment with a soluble receptor for activin improves bone mass and structure in the axial and appendicular skeleton of female cynomolgus macaques (Macaca fascicularis). Bone 46: 64–71. pmid:19781677
  64. 64. Ruckle J, Jacobs M, Kramer W, Pearsall AE, Kumar R, Underwood KW, et al. (2009) Single-dose, randomized, double-blind, placebo-controlled study of ACE-011 (ActRIIA-IgG1) in postmenopausal women. J Bone Miner Res 24: 744–752. pmid:19049340
  65. 65. Ayturk UM, Jacobsen CM, Christodoulou DC, Gorham J, Seidman JG, Seidman CE, et al. (2013) 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 28: 2081–2093. pmid:23553928
  66. 66. McClung MR, Grauer A, Boonen S, Bolognese MA, Brown JP, Diez-Perez A, et al. (2014) Romosozumab in postmenopausal women with low bone mineral density. N Engl J Med 370: 412–420. pmid:24382002
  67. 67. Tu X, Delgado-Calle J, Condon KW, Maycas M, Zhang H, Carlesso N, et al. (2015) Osteocytes mediate the anabolic actions of canonical Wnt/β-catenin signaling in bone. Proc Natl Acad Sci USA 112: E478–E486. pmid:25605937