A proteomics approach to study mouse long bones: examining baseline differences and mechanical loading-induced bone formation in young-adult and old mice

Protein and gene expression in bone correlate moderately

Paired tibias from young-adult and old mice at baseline were analyzed by proteomics and RNA-seq (Figure 1A, 1H). Femurs from these mice exhibited the expected age-related changes in bone morphology (Figure 1B–1G). After filtering, 1903 proteins were detected by proteomics, and 16273 genes were detected by RNA-seq (Supplementary Figure 1). At both the protein and RNA levels, multidimensional scaling (MDS) showed that young-adult samples separated from old samples. Of the detected proteins, 93% (1773/1903) were also detectable by RNA-seq (Figure 1I). Using the average peptide spectral matches (PSMs) and average counts per million (CPMs), the abundance of the 1773 targets detected by both proteomics and RNA-seq were correlated (Figure 1J), as reported [31]. The correlation was moderately positive (Spearman ρ = 0.40, p < 0.001), consistent with results in other tissues [31, 32]. Differential expression analysis between young-adult and old bone was performed separately for proteomics and RNA-seq. 183 proteins and 2290 genes met the p-value cutoff to be age-related differentially expressed proteins (DEPs) and differentially expressed genes (DEGs), respectively (Figure 1K).

Figure 1. RNA-seq and proteomics were used to characterize cortical bone from young-adult and old mice at baseline. (A) Untreated 5-month-old (young-adult) and 22-month-old (old) female C57BL/6N mice were sacrificed. Paired right and left tibial mid-diaphyses were isolated, removed of marrow, and snap frozen. From the right tibias, proteins were extracted using 4% SDS. Proteins from 5 tibias per age were analyzed by proteomics using a tandem mass tag (TMT)-11. From the left tibias, RNA was isolated using TRIzol. RNA from 7 tibias per age was sequenced. (B) MicroCT of the distal right femurs from these mice confirmed the expected age-related differences in the cortical bone. (C, D) The distal cortical bone area and cortical thickness were lower with age. (E) MicroCT also confirmed age-related changes in the trabecular bone of the distal femur. (F, G) The bone volume per total volume (BV/TV) and trabecular number were lower with age. (H) Proteomics and RNA-seq raw data were analyzed, and differential expression analysis was performed separately. For both methods, a Benjamini-Hochberg-adjusted p-value cutoff of 0.05 was used to identify differentially expressed genes (DEGs) and differentially expressed proteins (DEPs). Downstream analyses included correlations, overlaps, weighted gene co-expression network analysis (WGCNA), gene ontology (GO) analysis, pathway analysis, and COMPBIO analysis. (I) 93% (1773/1904) of proteomics hits (PSM≥3) were detectable by RNA-seq (non-zero CPM for all samples). (J) The abundance of the 1773 targets detected by both proteomics and RNA-seq (after PSM and CPM filtering) were correlated (Spearman). (K) Comparing young-adult and old bone at baseline, 183 proteomics targets and 2290 RNA-seq targets met the p-value cutoff to be DEPs and DEGs, respectively. Abbreviations: SDS: Sodium dodecyl sulfate; BV/TV: Bone Volume/Total Volume; DEG: Differentially Expressed Gene; DEP: Differentially Expressed Protein; CPMs: Counts per million; PSMs: Peptide spectral matches.

Many of the most significantly up- and downregulated proteins and genes between young-adult and old bone have been associated with bone phenotypes in GWAS

We compiled the top 15 most significantly differentially expressed (by FDR) proteins and genes (both up- and downregulated) between young-adult and old bone at baseline (Table 1). Several of these targets have been associated with phenotypes related to bone [33–49] or aging [50, 51] in genome-wide association studies (GWAS). The only top-15 target identified by both platforms was MMP13, which was more highly expressed in old bone at both the protein and transcript level. MMP13 is critical for osteocyte perilacunar remodeling [52] and maintains cartilage homeostasis [53].

The most significantly upregulated protein in old bone at baseline was CTSS, a cysteine protease known as cathepsin S that regulates extracellular matrix (ECM) remodeling and antigen presentation. CTSS interacts with osteocalcin (Bglap) and has been shown to control osteoblast differentiation and bone turnover [54]. PYCARD, which promotes apoptosis and inflammation, and OSTF1, which is known as osteoclast stimulating factor 1 and directly induces osteoclast differentiation and bone resorption [55], were almost two-fold higher in old bone. Additional proteins that were higher in old bone include IPO5, an important nuclear transport receptor; CD36, which is important for both osteoblast [56] and osteoclast [57] function; DYSF, which is known to act as a calcium sensor in muscle to facilitate membrane repair [58, 59]; and BTF3, which plays a role in c-Jun transcriptional activity.

The most significantly downregulated protein in old bone at baseline was SOST, commonly called sclerostin. SOST is a well-known Wnt antagonist in bone that is highly expressed by mature osteocytes and is the target of the most recently approved osteoporosis drug romozosumab [60]. BASP1 is also lower in old bone, and although little is known about BASP1 in the context of bone, it was detected in a proteomic profiling of osteoblast differentiation [61]. COL5A1, a target of TGF-beta in bone [62], and COL11A1, which regulates bone microarchitecture during development [63], are both lower in old bone. Four other proteins that were lower in old bone included: CD44, which inhibits inflammatory bone loss [64]; MBL1, a mannose binding lectin that may suppress osteoclastogenesis [65]; TNN, the murine version of tenascin W [66] that is induced in bone cells by loading and Wnt signaling [67]; and S100A1, a calcium binder that is important in mineralization [68].

The most significantly upregulated genes in old bone at baseline play a role in osteoclastic resorption, including Tubb6 [69], Pstpip1 [70], and Dock5 [71, 72]. Tyrob affects the myeloid lineage (precursors of osteoclasts) and has been linked to a human disease involving bone cysts [73, 74]. Cebpb is known to be important in osteoblast differentiation [75]. H2-k1 and H2-d1 are orthologs of the human HLA-A gene, and while HLA-A is expressed by osteoblast-lineage cells [76], the higher expression may reflect more inflammation and immune cell infiltration in older bone [77]. H6pd, of the pentose phosphate shunt, and Pik3ap1, a player in B-cell development, were also upregulated in old bone. Last, Acp5, better known as tartrate-resistant acid phosphatase (TRAP), was also higher in old bone. TRAP is expressed by bone-resorbing osteoclasts and is also expressed by osteocytes during perilacunar remodeling [78].

The most significantly downregulated gene in old bone at baseline was Hdac9, an inhibitor of osteoclastogenesis [79]. A number of other downregulated genes are related to TGF-beta or bone morphogenetic protein (BMP) signaling. Ltbp1, known as latent TGF-beta binding protein 1, directly modulates TGF-beta activity and is substrate of matrix metalloproteinases (MMPs) [80]. The lower expression of Ltbp1 may relate to lower levels of the identified BMP regulators Smad9 [81], Chrdl1 [82], and Adamts17 [83]. Nduf4l2, which has been suggested to be important for the metabolic transition of osteoblasts into osteocytes, was lower in old bone [84]. Similar to findings at the protein level, Sost was also a downregulated DEG (FC: -3.9, FDR: 9.5E-05). Other genes more lowly expressed in old bone included Magi2 and Pip4k2a, which may play a role in PI3K-AKT signaling, and Tcf7l2, which is a Wnt effector transcription factor [85].

Integrated analysis of GWAS BMD hits with baseline proteomics and RNA-seq from young-adult and old bone identified eight targets including Asrgl1 and Timp2

Between proteomics and RNA-seq, 71 shared targets were differentially expressed between bones from young-adult and old mice at baseline. To further narrow these candidates to those most relevant to human disease, we intersected them with the hits identified in a recent GWAS that identified genetic determinants of BMD [39], resulting in a list of eight targets (Figure 2A). These included the well-known bone factors Sost, Col1a1, Col1a2, and Mepe, as well as Ncam1 [86] and Itgb5 [87, 88] which have also been studied in the context of bone. We also identified two targets Asrgl1 (Asparaginase and isoaspartyl peptidase 1) and Timp2 (Tissue inhibitor of metalloproteinases 2) that have been less described in the context of bone.

Figure 2. Integrated analysis of GWAS BMD hits with baseline proteomics and RNA-seq identified eight targets including Asrgl1 and Timp2. (A) Intersection of baseline proteomics DEPs, RNA-seq DEGs, and a GWAS database [39] revealed 8 targets. Of these, Sost, Asrgl1 and Timp2 were further evaluated histologically. (B) SOST expression was validated using IHC in tibias of young-adult and old mice (transverse sections cut at the proximal diaphysis). SOST expression was restricted to osteocytes. (C) Fewer SOST-positive and (D) total osteocytes (per 600 mm x 600 mm ROI) were present in old bone. (E) The percentage of SOST-positive osteocytes and (F) the areal density of SOST-positive osteocytes were not lower in old bone. (G) Asrgl1 expression was characterized using RNA ISH. (H) Most of the periosteal surface was covered by cells expressing Asrgl1, and this percentage of cells decreased with age. (I) Periosteal Asrgl1 expression was lower in old bone. (J) Slightly less than half of osteocytes (OCY) within cortical bone expressed Asrgl1, and this percentage of OCY was not lower in old bone. (K) The areal density of Asrgl1-positive osteocytes trended to be lower in old bone but did not reach significance. (L) Timp2 expression was characterized using RNA ISH. (M) Most of the periosteal surface was covered by cells expressing Timp2, and this percentage of cells trended to decrease with age. (N) Periosteal Timp2 expression was two times lower in old bone. (O) Most osteocytes within cortical bone expressed Timp2, and this percentage of OCY decreased with age. (P) The areal density of Timp2-positive osteocytes decreased with age. Abbreviations: OCY: Osteocyte; CB: Cortical bone; Ma: Marrow; Mu: Muscle; Pos: Positive; Ps: Periosteal; Surf: Surface). Data shown as mean +/−SD. p-values calculated by unpaired, 2-tailed t-test; n = 6–7 mice per age.

To validate our proteomics findings, we assessed SOST expression histologically at the protein level in a separate set of mice. As expected, SOST expression was restricted to osteocytes, particularly the more mature osteocytes further from the surface (Figure 2B). Fewer SOST-positive and total osteocytes per transverse section were present in old bone (Figure 2C, 2D), consistent with the measured lower SOST at the protein level. Contrary to our expectations, the percentage of SOST-positive osteocytes and the areal density of SOST-positive osteocytes was not lower in old bone (Figure 2E, 2F), suggesting that the decrease in measured SOST at the protein level is due to the lower absolute number of osteocytes.

To better characterize Asrgl1 and Timp2, we used RNA ISH to confirm mRNA expression in bone cells and characterize the expression pattern in young-adult and old bone. Periosteal Asrgl1-expressing cells cover less than a quarter (23%) of the surface in young-adult bone, and this percentage is lower in old bone (12%, p = 0.004, Figure 2G, 2H). Consistent with this, the total periosteal surface expression of Asrgl1 (based on DAB-positive area) decreased with age (young: 5440, old: 4735 µm²/mm, ^*p = 0.040, Figure 2I). Slightly less than half of osteocytes also express Asrgl1 in young-adult bone (43%), and this percentage of positive osteocytes was not significantly lower with age (38%, p = 0.22, Figure 2J). The areal density of Asrgl1-positive osteocytes trended to be lower but did not reach significance (340 vs. 264 cells/mm², p = 0.066, Figure 2K). While Asrgl1 expression is detected in fewer than one half of bone cells, Timp2 is expressed by most of the cells on the periosteal surface in young-adult and old bone (96% and 89%, respectively; Figure 2L, 2M). The total periosteal surface expression of Timp2 was decreased by half with age (young: 5115 vs. old: 2446 µm²/mm, p = 0.007, Figure 2N). Most osteocytes in young-adult bone also express Timp2 (88%), but this percentage decreases with age (76%, p = 0.025, Figure 2O), as does the areal density of Timp2-positive (858 vs. 649 cells/mm², p = 0.019, Figure 2P). In summary, we observed reduced expression of Asrgl1 and Timp2 in bones of old mice using three methods—proteomics, RNA-seq, and RNA ISH.

An age-related module indicating baseline differences in TGF-beta and Wnt signaling was identified through co-expression analysis of the RNA-seq data

A co-expression network analysis (WGCNA) of the RNA-seq data identified 44 unique co-expression modules (Figure 3A). Three modules had a significant association with age, but only Module 40, with 2452 genes, correctly separated all young-adult and old samples by hierarchical clustering (Figure 3B and Supplementary Figure 2). We then used COMPBIO to analyze the genes in Module 40, limiting the analysis to the 1055 genes that had an FDR <0.05 in the DEG analysis (young-adult vs. old). Within Module 40, we identified numerous themes that were enriched. The top three themes were TGF-beta signaling, Wnt signaling, and Pi3K/AKT and MAP/ERK, and notably these themes were interconnected in the network (Figure 3C). Other top themes included histone H3 acetylation, cell-cell adherens junctions, and endochondral ossification (Figure 3D). PANTHER pathway analysis of all 2452 genes reinforced differences in Wnt signaling and TGF-beta signaling (Figure 3E). It also identified other pathways including gonadotropin-releasing hormone receptor, Alzheimer disease-presenilin, integrin signaling pathway, and angiogenesis. The genes driving this enrichment were examined for all nine pathways, and nearly all of the genes were more lowly expressed in old bone, suggesting that these pathways are less active in aged bones. KEGG analysis of all 2452 also identified differences in TGF-beta signaling and Wnt signaling but further identified differences that included: axon guidance; hippo signaling; ECM-receptor interaction; parathyroid hormone synthesis, secretion, and action; and AGE/RAGE signaling pathway in diabetic complications (Figure 3F). The diversity of pathways identified in this age-related module underscores the complexity of aging processes occurring simultaneously in bone.

Figure 3. A single age-related module with perfect separation was identified using weighted gene co-expression network analysis (WGCNA) of the RNA-seq data. (A) Gene counts across all samples were clustered into 44 modules and correlated with age status. Three modules (black and orange) were significantly different between ages using a Bonferroni correction (p = 0.0011 = 0.05/44). (B) Only Module 40 (orange) displayed perfect separation using hierarchical clustering. (C, D) COMPBIO analysis of the subset of 1055 genes from Module 40 that had an FDR <0.05 between young-adult and old samples identified TGF-beta signaling and Wnt signaling as the top themes. (E, F) PANTHER Pathway and KEGG analyses of all 2452 genes in Module 40 also revealed Wnt signaling and TGF-beta signaling, among other pathways. Examination of the fold-changes of the genes from the PANTHER analysis showed that all 9 pathways are reduced with aging. Abbreviation: FE: fold-enrichment.

Baseline age-related differences in ECM/MMPs and TGF-beta signaling were identified in both the proteome and transcriptome

To investigate the age-related processes by both proteomics and transcriptomics, we used COMPBIO to analyze the top 500 (by p-value) proteins (Figure 4A) and genes (Figure 4B) that were differentially expressed at baseline between old and young-adult bone. At the protein level, the top themes were nuclear lamin/progeria, ECM/MMPs, and ER-associated degradation (ERAD) (Figure 4C). The most robustly interconnected group, which contained ECM/MMP, also included themes related to actin cytoskeleton, osteoblasts, and TGF-beta signaling. At the RNA level, the top themes were ECM/MMPs, dendrites, and TGF-beta-signaling (Figure 4D). Again, the ECM/MMP-containing group was the most robustly interconnected and additionally included actin cytoskeleton and cell-cell adherens junctions. Notably, only proteomics detected age-related differences related to nuclear lamin/progeria, ERAD, ubiquitination, and chaperone proteins.

Figure 4. Baseline age-related differences in ECM/MMPs and TGF-beta signaling were detected in both the proteome and transcriptome. COMPBIO analysis of (A) the proteome and (B) the transcriptome indicated baseline changes in ECM/MMPs. (C) At the protein level, the top 10 themes included Lamin/Progeria, ECM/MMPs, and ER-Associated Degradation. (D) At the RNA level, the top 10 themes included ECM/MMPs, Dendrites, and TGF-beta signaling. Examination of age-related differences in individual (E) proteins and (F) genes related to TGF-beta signaling and perilacunar remodeling (PLR) showed reduced TGF-beta levels but higher Mmp13. (G) RNA ISH for Tgfb2 showed predominant expression in osteocytes. (H) Periosteal Tgfb2 expression was minimal in young-adult bone and even lower in old bone. (I) Over half of osteocytes (OCY) within cortical bone expressed Tgfb2, and this percentage was lower in old bone. (J) The areal density of Tgfb2-positive osteocytes was about half in old versus young bone.

Collectively, the most prominent signals involved ECM/MMPs and TGF-beta signaling. Previously, osteocyte-intrinsic TGF-beta signaling has been shown to control perilacunar remodeling (PLR) [89], an important process for maintaining bone quality and fracture resistance [52]. Specifically, TGF-beta signaling is thought to control the expression of MMPs and other matrix-remodeling enzymes such as cathepsin K (CTSK) and tartrate resistant acid phosphatase (ACP5). Therefore, we examined baseline changes in PLR factors in our data. By proteomics, TGF-beta 1 (TGFB1) was the only TGF-beta type detectable, and it was two times lower in old bone (Figure 4E). MMP2 was significantly lower in old bone while MMP13 expression was significantly higher. Other proteins related to perilacunar remodeling such as CTSK and ACP5 trended higher in old bone but did not reach significance. At the transcript level, all three TGF-beta types were detected, and in young-adult bone, Tgfb1 and Tgfb2 were the most highly expressed at baseline (Supplementary Figure 3). In old bone, Tgfb2 expression was nearly four-fold lower (Figure 4F) while Tgfb1 expression was also significantly lower but by a much smaller magnitude. At the transcript level, the directional changes of the PLR genes were the same as observed by proteomics, with Mmp13, Ctsk, and Acp5 being significantly higher in old bone.

To better localize age-related TGF-beta signaling factors, we used RNA ISH to examine the expression of Tgfb2. We found that Tgfb2 is predominantly expressed by osteocytes in young-adult and aged bone (Figure 4G). On the periosteal surface, less than 10% of periosteal surface cells expressed Tgfb2; this percentage significantly decreased with age (9% vs. 3%, p = 0.005, Figure 4H). On the other hand, over 50% of osteocytes in young-adult bone expressed Tgfb2, and this percentage also decreased in old bone (58% vs. 39%, p = 0.0005, Figure 4I). The areal density of Tgfb2-expressing osteocytes was approximately half in old bone (495 vs. 272 cells/mm², p = 0.0005, Figure 4J). Therefore, both proteomic and transcriptomic approaches revealed baseline differences in TGF-beta signaling, particularly related to Tgfb2, which may be involved in PLR.

The proteome differed more with age than mechanical loading

We next sought to study the proteomes of bones from young-adult and old mice following a mechanical loading stimulus. Tibias from a separate set of young-adult and old mice were loaded in vivo for 1 or 5 days (Figure 5A). Time points were selected to sample distinct phases of the loading response: early mechanosensation (day 1) and active bone formation (day 5). Paired loaded and non-loaded tibias were analyzed by proteomics (Figure 5B, 5C). After filtering, 2300 and 2140 proteins were detected by proteomics at days 1 and 5, respectively (Supplementary Figure 1). Multidimensional scaling (MDS) showed that young-adult samples separated from old samples regardless of loading status at day 1 (Figure 5D). As expected, at day 5, when bone formation is actively occurring [11], samples separated marginally better based on loading status, but the strongest proteomic differences were still due to age status (Figure 5E).

Figure 5. Experimental design for proteomics loading experiment. (A) 20 female C57BL/6N mice (5 mice/day/age) were subjected to daily in vivo axial tibial compression for either 1 (Day 1) or 5 (Day 5) bouts and sacrificed 6 hours after their final bout of loading. (B) The loaded (right) and non-loaded (left) tibial mid-diaphyses were isolated, removed of marrow, and snap frozen. Proteins were extracted in 4% SDS. Proteins for all 40 samples were then analyzed by proteomics using a tandem mass tag (TMT)-11 design. (C) Proteomics raw data were analyzed, and differential expression analysis was performed using ProteoQ. An unadjusted p-value cutoff of 0.05 and a fold-change threshold of 1.1 were used to identify loading-regulated proteins for downstream analyses, which included gene ontology (GO), PANTHER pathways, COMPBIO, and correlation with previously published RNA-seq data [11]. (D) At Day 1, multidimensional scaling (MDS) showed that the strongest differences were between ages rather than with loading status. (E) A Day 5, MDS showed slightly better separation between loaded and non-loaded samples, but the strongest separation was still between ages.

The top 10 upregulated and downregulated proteins following loading were compiled for each age and day combination (Table 2). Using the same DEP criterion of an FDR <0.05 to compare loaded and non-loaded samples, only five proteins (SMPD3, FKBP7, SSR3, EIF4G1, and EIF2S3X) reached significance across the entire loading experiment, all upregulated in young-adult bones at day 5. SMPD3, a sphingomyelinase that produces ceramide [90], was the most significantly upregulated protein. Ceramides are bioactive lipids that play a role in triggering proliferation and differentiation [91]. SMPD3 has GWAS associations with both BMD and height and is known to be important in bone for development [92] and fracture healing [93]. FKBP7 is a calcium-binding molecular chaperone. In old bones, neither protein was significantly upregulated following loading at day 5. Given that Smdp3 and Fkbp7 were also previously identified as upregulated DEGs (7- and 3-fold up, respectively) in young-adult mice but downregulated DEGs in old mice (3- and 2-fold down, respectively) [11], they may play a role in the diminished loading-induced bone-formation with aging. SSR3, a signal sequence receptor important for protein translocation across the ER membrane, was also significantly upregulated at the protein level. It was also a loading DEG [11] and has been linked to BMD and height, although little is known about its direct role in bone cells. The last upregulated DEPs were EIF4G1 and EIF2S3X, which both play a role in translation, and may reflect increased matrix synthesis during the loading-induced bone formation response. Notably, in young-adult mice at day 5, we also detected a trend toward reduced MMP13, consistent with our recent findings using gene microarray and IHC [94].

Proteomes following loading showed distinct pathway, protein class, process enrichment

As input for pathway analysis, we defined loading-responsive proteins using the following criteria: a linear fold-change of at least 1.1 (up or down) and an unadjusted p-value < 0.05. At day 1, tibias from young-adult mice had 150 upregulated and 109 downregulated loading-responsive proteins, whereas tibias from old mice had 38 upregulated and 38 downregulated loading-responsive proteins (Figure 6A). At day 5, tibias from young-adult mice had 88 upregulated and 137 downregulated loading-responsive proteins, whereas tibias from old mice had 59 upregulated and 120 downregulated loading-responsive proteins. We next examined the pathways that changed at the protein level following mechanical loading using combined upregulated and downregulated protein lists with PANTHER pathways (Figure 6B). We extended this pathway analysis with a more expansive and directional (upregulated and downregulated proteins separately) analysis: Reactome pathways (Figure 7A), gene ontology (GO) terms (Supplementary Figure 4), protein classes (Supplementary Figure 5), and COMPBIO themes (Supplementary Figure 6).

Figure 6. Tibias from old mice had fewer protein-level changes following loading compared to young-adult mice. (A) Volcano plots for each day and age combination following the loading experiment. The total number of loading regulated proteins is shown above with the upregulated number in red and the downregulated number in blue. (B) Loading-related proteins were input into PANTHER to identify enriched PANTHER pathways. (C) Loading-related PANTHER pathways were temporally arranged.

Figure 7. Temporal map of PANTHER Reactome Pathways for proteomics loading experiment. (A) Loading-responsive protein lists for each age and day combination (up and down lists separate) were input to PANTHER to identify enriched pathways. (B) Reactome pathways were arranged temporarily for upregulated or downregulated protein enrichment at each day and age combination. Both upregulated and downregulated lists were enriched for pathways.

In young-adult bone at day 1, we observed enrichment of integrin signaling, nicotinic acetylcholine receptor signaling, Parkinson disease, and TCA cycle (Figure 6C). Consistent with the TCA cycle signal, Reactome pathway analysis indicated that upregulated proteins were enriched for pathways related to complex I and the formation of ATP by chemiosmotic coupling (Figure 7B). Protein class analysis showed enrichment of ATP synthase and dehydrogenase classes while GO processes related to mitochondrial electron transport. Actin-binding cytoskeletal proteins and ECM structural proteins showed protein class and GO enrichment, but interestingly, downregulated proteins were enriched for collagen-related processes. As expected at this early time point, we observed a prominent transcriptional signal.

In young-adult bone at day 5, no proteins were significantly enriched for any PANTHER pathway, but Reactome analysis showed that upregulated proteins were enriched for collagen biosynthesis. COMPBIO also identified an osteoblast theme within its top ten. The upregulated proteins were enriched for processes related to translation, consistent with the GO, protein class, and COMPBIO enrichments. Downregulated proteins at day 5 showed enrichment for amino acid regulation of mTORC1 and regulation of insulin-like growth factor.

In old bone at day 1, no proteins were significantly enriched for any PANTHER pathway, but Reactome analysis indicated an immune response. Protein class analysis indicated enrichment of calcium-binding proteins and ribosomal proteins while top COMPBIO themes related to ubiquitin degradation, actin filament dynamics, and RHO GTPase activity.

In old bone at day 5, we observed enrichment of integrin signaling and glycolysis. Reactome pathway analysis indicated that downregulated proteins were enriched for glycolysis and gluconeogenesis, consistent with the GO process enrichment and the protein class enrichment of isomerases and dehydrogenases. Downregulated proteins also showed enrichment for amino acid regulation of mTORC1, as in the young-adult samples. Upregulated proteins were enriched for collagen-related pathways, consistent with the GO and protein class enrichment related to the ECM. COMPBIO analysis reinforced these findings and also indicated changes in nuclear lamins/envelope. In contrast to the GO analysis in young-adult samples at this time point, we observed enrichment of processes related to transcription rather than translation, potentially suggesting a delay in the bone formation cascade with age.

Mice

Female C57BL/6N mice were obtained at 5 months (young-adult) and 22 months (old) from the aged rodent colony at the National Institute on Aging, which is managed by Charles River Laboratories (Figure 1A). Female mice were selected because: (1) osteoporosis is more prevalent in females [1]; (2) male mice often fight, which can impact the skeleton and confound the effects of loading [122]; and (3) to facilitate comparison with our previous RNA-seq study, which used female mice [26–30, 11]. Across all experiments, 58 mice (34 young-adult and 24 old) were used in this study. Mice were housed in groups of up to five animals of the same age and kept on a 12-hour light/dark cycle under standard conditions with ad libitum access to water and chow (Purina 5053 and 5058). All animal work was approved by and in compliance with the Washington University IACUC. All included mice were healthy throughout the experiment. For the baseline experiment, we performed paired proteomics and RNA-seq on tibias of young-adult (5-month) and old (22-month) C57BL/6N female mice not subjected to any interventions (Figure 1A). For the loading experiment, we used proteomics to compare the loading responses between tibias of young-adult and old mice following 1 or 5 days of axial loading (Figure 5A). MicroCT of the right femurs from the mice used in the baseline comparison confirmed the expected age-related changes in bone, specifically decreases in distal cortical bone area, distal cortical thickness, trabecular bone volume per total volume (BV/TV), and trabecular number (Figure 1B–1G and Supplementary Figure 7).

In vivo mechanical loading

In the loading experiment, mice were loaded for either 1 or 5 bouts of daily loading. Mice were anesthetized (3% isofluorane) and subjected to loading each morning for the specified number of days. With the mouse prone, the right leg (tibia) was placed vertically in the loading fixture with the knee positioned superiorly in a semi-spherical cup (10 mm diameter) attached to the system actuator, and the foot held in a static fixture inferiorly (20° of dorsiflexion). A preload (−0.5 N) was applied, and tibias were subjected to axial compression for 60 cycles/day (4Hz haversine waveform) using the Electropulse 1000 materials testing system (Instron). This loading protocol is anabolic for both cortical and trabecular bone in young-adult mice [123]. We used a strain-matched study design. Based on prior strain gauging analyses [124] and consistent with our RNA-seq study, age-specific peak forces of -8N and -7N were selected for the 5- and 22-month-old mice, respectively, to engender average peak compressive periosteal strains of −2200 µε at the cortical mid-shaft [28]. Corresponding tensile strains on the anterio-medial surface were approximately 1200 µε [124]. After each loading bout, buprenorphine (0.1 mg/kg subcutaneously) was given to mitigate pain from loading [125], and mice were returned to their cages to resume unrestricted activity. The left tibias served as non-loaded, contralateral controls. Six hours after their final loading bout, mice were euthanized by CO₂ asphyxiation.

Tibial isolation

For proteomics and RNA-seq analyses, right and left tibias were stripped of muscle, cut at the distal tibiofibular junction and 2 mm distal to the tibial plateau, placed in ice-cold PBS, centrifuged to remove the bone marrow [125, 126], flushed with PBS, and snap frozen in liquid nitrogen. Samples were stored at -80°C until protein or RNA extraction. For the histology assays, the left tibias were cut at the ankle and above the tibial plateau and trimmed of muscle with care not to disturb the periosteum; samples were immediately fixed in 10% neutral-buffered formalin (NBF).

Protein extraction by homogenization and sonication

Preliminary experiments were performed to compare the protein yield and peptide detection of several tissue homogenization methods (Supplementary Figure 8). Frozen tibial samples were homogenized with the cryoPREP (Covaris CP02) pulverizer in small impactor bags (tissueTUBES, Covaris # 520071), as previously described in other tissues [127]. With the bag partly submerged in liquid nitrogen, the bone was transferred into the bag, and the cap was lightly screwed on (to allow air to escape during the impact). The CryoPrep was activated on impact level 5. After impact, 100 µL of SDS (4 % SDS, 100 mM Tris-HCl, pH 8.0) at room temperature (RT) were added to the bag. A pipette was used to spread and agitate the buffer within the bag to cover all bone chips. Eventually, the bone lysate became ‘stringy’ and was pushed into a corner of the bag. Using a pipette (with cut tip), the lysate was transferred to a sonication vial (Covaris #520130). After incubating at RT for 5 min, the samples were sonicated in the S220X focused-ultrasonicator (Covaris) for 4 min with the following settings: PIP: 500 Watt, cycles per burst: 500, duty factor: 10%. Sonicated samples were transferred to 1.7 mL Eppendorf tubes and spun at 14,000 × g for 15min at 8°C to pellet debris. The supernatant was collected, transferred to a new 1.7 mL tube, and 2 µL were reserved for protein concentration determination, according to BCA kit instructions (Pierce #23227); the remaining lysate was reduced with 100 mM DTT followed by heating at 95°C for 10 min. Samples were stored at −80°C.

Peptide preparation

The samples were digested using a modification of the filter-aided sample preparation (FASP) method, as previously described [128]. The reduced samples were mixed with 200 µL of 100 mM Tris-HCL buffer, pH 8.5 containing 8M urea (UA buffer) and transferred onto the top chamber of a 30,000 MWCO cutoff filtration unit (Millipore #MRCF0R030) and spun in a microcentrifuge at 14,000 × g for 10 min. An additional 200 µL of UA buffer was added to the top chamber of the filter unit and the filter was centrifuged at 14,000 × g for 15–20 min. The flow through was discarded, and the proteins were alkylated by adding 100 µL of 50 mM Iodoacetamide (Pierce #A39271) in UA buffer to the top chamber of the filtration unit and gyrating at 550 rpm in the dark at RT for 30min using a thermomixer (Eppendorf #05-400-205). The filter was spun at 14,000 × g for 15 min, and the flow through was discarded. Unreacted iodoacetamide was washed through the filter with two sequential additions of 200 µL of UA buffer and centrifugation at 14,000 × g for 15–20 min after each addition. The urea buffer was exchanged into digestion buffer (DB) consisting of 50 mM ammonium bicarbonate buffer at pH 8. Two sequential additions of 200 uL of DB with centrifugation after each addition to the top chamber were performed. The top filter units were transferred to a new collection tube, and 1 micro unit of LysC (Wako Chemicals #129-02541) was added. Samples were digested at 37°C for 2hr. After LysC digestion, 1 µg of sequencing-grade trypsin (Promega #V5113) was added, and samples were digested overnight at 37°C. The filters were spun at 14,000 × g for 15min to collect the peptides in the flow through. The filter was washed with 50 µL of 100 mM ammonium bicarbonate buffer, and the wash was collected with the peptides. In preparation for desalting, peptides were washed 3 times in 1mL ethyl acetate followed by acidification to 1% (vol/vol) trifluoroacetic acid (TFA) final concentration. The peptides were desalted using porous graphite carbon two micro-tips (Glygen BIOMEK NT3CAR) on a Beckman robot (Biomek NX), as previously described [129]. The peptides were eluted with 60% (vol/vol) acetonitrile in 0.1% TFA (vol/vol) and dried in a SpeedVac (Thermo Scientific Savant DNA 120) after adding TFA to 5% (vol/vol). The peptides were dissolved in 20 µL of 1% (vol/vol) acetonitrile in water. An aliquot (10%) was removed for quantification using the Pierce Quantitative Fluorometric Peptide Assay kit (Thermo Scientific, Cat. No. 23290). The remainder of the enriched samples and reference pool sample were transferred into a 0.5 mL Eppendorf tube, dried in the Speed-Vac and dissolved in 12 µL of HEPES buffer (100 mM, pH 8.0) (Sigma, H3537).

The peptides were labeled with tandem mass tag (Thermo Scientific TMT10 or TMT11-131C) reagents according to the manufacturer protocol. For the baseline comparison between old and young-adult bones, the samples were analyzed in a single plex using TMT-10 (Supplementary Figure 9). For the loading comparison at each time point (Figure 5A), young-adult and old samples (paired tibias from n = 5 mice per age) were analyzed on two age-specific runs separately. An 11th bridge sample containing pooled young-adult and old samples per time point was created to facilitate quantitative comparisons between ages. The labeled samples that were assigned to each ten- or eleven-plex were pooled, dried, and dissolved in 120 µL of 1% FA. The TMT-labeled samples were desalted as described above for the unlabeled peptides. The eluates were transferred to autosampler vials (Sun-Sri #200046), dried, and stored at −80ºC.

Proteomics nano-Liquid Chromatography–Mass Spectrometry (nano-LC-MS)

The samples in formic acid (1%) were loaded (2.5 µL) onto a 75 µXm i.d. × 50cm Acclaim PepMap 100 C18 RSLC column (Thermo Fisher Scientific) on an EASY nanoLC (Thermo Fisher Scientific) at a constant pressure of 700bar at 100% A (0.1%FA). Prior to sample loading, the column was equilibrated to 100% A for a total of 11 µL at 700bar. Peptide chromatography was initiated with mobile phase A (1% FA) containing 5% B (100%ACN, 1%FA) for 1 min, then increased to 15% B over 108 min, to 25% B over 87 min, to 35% B over 40 min, to 70% B in 6 min, to 95% B over 2 min and held at 95% B for 18 min, with a flow rate of 300 nL/min. The data were acquired in data-dependent acquisition (DDA) mode. The full-scan mass spectra were acquired with the Orbitrap mass analyzer with a scan range of m/z = 375 to 1500 and a mass resolving power set to 70,000. Twelve data-dependent high-energy collisional dissociations were performed with a mass resolving power set to 35,000, a fixed lower value of m/z 100, an isolation width of 1.2Da, and a normalized collision energy setting of 32. The maximum injection time was 60ms for parent-ion analysis and 120ms for product-ion analysis. The target ions that were selected for MS/MS were dynamically excluded for 20sec. The automatic gain control (AGC) was set at a target value of 3e6 ions for full MS scans and 1e5 ions for MS2. Peptide ions with charge states of 1 or ≥7 were excluded for HCD acquisition.

Proteomics analysis

PSM files were imported into ProteoQ (https://github.com/qzhang503/proteoQ) for normalization, quantitation, and analysis. Data were further explored using principal component analysis and multidimensional scaling. At baseline, differential expression analysis was performed between young-adult and old samples. For the loading comparison, differential expression analysis was performed separately within an age for each time-point between loaded and non-loaded samples. Linear modeling was performed using the contrast fit approach in limma [130] to assess the statistical significance in protein abundance differences between indicated groups of contrasts. To facilitate the high-confidence assignment of proteins, we used a strict criterion of peptide spectral matches (PSM) ≥3 to call each protein.

RNA isolation and library prep

Frozen bones were homogenized using a mikro-dismembrator (Braun Biotech International), and total RNA was extracted using TRIzol (Ambion) with the RNeasy Kit (Qiagen # 74004). Total RNA integrity was determined using the 4200 Tapestation, and concentrations were measured with the Qubit fluorometer. The median RNA integrity number (RIN) was 5.8 with a range of 5.1–6.8 (Supplementary Figure 1). No samples were excluded. Library preparation was performed with 500 ng of total RNA, and ribosomal RNA was blocked using FastSelect reagents (Qiagen) during cDNA synthesis. RNA was fragmented in reverse transcriptase buffer with FastSelect reagent, and mRNA was reverse transcribed to yield cDNA using SuperScript III RT enzyme (Life Technologies) and random hexamers.

RNA sequencing and aligning

A total of 14 samples were subjected to RNA-seq (Figure 1B) by the Washington University Genome Technology Access Center. Samples were indexed, pooled, and sequenced on an Illumina NovaSeq 6000. Basecalls and demultiplexing were performed with bcl2fastq2 (Illumina), and the reads were aligned to Ensembl release 76 using STAR (2.5.1a) [131]. Gene counts were derived from uniquely aligned, unambiguous reads using Subread:featureCount (1.4.6-p5) [132].

RNA-seq analysis

To adjust for differences in library size between samples, normalization factors from gene counts were calculated using EdgeR [133]. Ribosomal genes and genes not expressed >1 count-per-million (CPM) in at least six samples were excluded from initial analyses. The size factors and matrix of counts were then imported into limma [130], and weighted likelihoods based on the observed mean-variance relationship were calculated using voomWithQualityWeights [134]. For the downstream analysis, genes that were not expressed (CPM = 0) in any samples were filtered out. Statistical model fitting and robustness were then assessed with multidimensional scaling (Supplementary Figure 1). Next, differential expression analysis between young-adult and old samples was performed using limma’s moderated t-tests, and the results were filtered for genes with Benjamini-Hochberg adjusted p-values < 0.05.

Co-expression network analysis for module construction

For the RNA-seq data, weighted gene co-expression network analysis (WGCNA) was used to generate co-expression networks for genes that differed between young-adult and old samples [135]. The Limma voomWithQualityWeights-moderated log2 counts-per-million generated for the previous differential expression analysis were used as inputs. All genes were correlated across each other using Pearson correlations and clustered by expression similarity into unsigned modules using a power threshold empirically determined from the data. To identify modules associated with age, the eigengene of each module was determined and correlated with age status (young-adult vs. old).

Gene ontology, pathway, protein class, and COMPBIO analyses

For the gene ontology (GO), pathway, and protein class analyses, the differentially expressed genes (DEG), differentially expressed proteins (DEP), or loading-regulated protein lists were input into PANTHER (version 14.1) [136] to identify enrichment. GO Slim was used for the GO analysis, and both PANTHER and Reactome pathways were analyzed in the pathway analysis. GO terms, pathways, and protein classes were arranged temporally by age and then alphabetically. For the COMPBIO (COmprehensive Multi-omics Platform for Biological InterpretatiOn, https://becker.wustl.edu/resources/software/compbio/) analysis, the described gene lists were input with the default parameters to determine the most prominent themes. The COMPBIO platform uses an ontology-free approach to generate a comprehensive and contextual map of the core biological concepts and themes associated with input entities. Specifically, it assembles the maps using contextual language processing algorithms to scan all PubMed abstracts to identify enriched concepts associated with the input entities. The platform utilizes conditional probability analysis to compute the statistical enrichment of biological concepts (processes/pathways) over those that occur by random sampling. Related concepts built from the input list of differentially expressed genes are clustered into themes (e.g., biological pathways/processes, cell types and structures) and further interconnected into groups.

Bone fixation for histology and RNAscope in situ hybridization

After 32hrs of fixation in NBF at RT, tibias were washed in PBS and decalcified in 14% ethylenediaminetetraacetic acid (EDTA) for 18 days at 4°C. EDTA was changed daily for the first 3 days and then every other day thereafter. Full decalcification was confirmed with X-ray. The tibias were processed for transverse (cut at the proximal diaphysis) paraffin sectioning (5 µm) by the Musculoskeletal Research Center Histology Core.

RNA in situ hybridization (ISH) using RNAscope

The manufacturer’s protocol for the RNAscope^® 2.5 HD Detection Reagents (ACD #322310) was followed except where noted. Paraffin sections were baked overnight at 60°C and deparaffinized. Hydrogen peroxide (ACD #322335) was applied for 10 min to quench endogenous peroxides, and pepsin (Sigma #R2283) was applied for pre-treatment at 37°C for 45 min. The remaining protocol was followed for the target probe hybridization (Supplementary Table 1), amplification, DAB reaction, and counterstaining. A negative and positive control probe were included with each batch, and for each round, the control results were as expected (Supplementary Figure 10). Slides were cover slipped with VectaMount (Vector #H-5000-60).

Immunohistochemistry (IHC)

For MMP13, sections were deparaffinized and rehydrated. Antigen retrieval was accomplished using UNI-TRIEVE (Innovex Biosciences #NB325) at 60°C (water bath) for 30min. For the remaining staining, the Universal Animal IHC Kit (Innovex #329ANK) was used with slight modifications. Endogenous peroxidases were blocked with PEROX-BLOCK for 15 min at RT, and background staining was blocked using Fc-Block for 45 min followed by Background Buster for 45 min. Primary anti-MMP13 antibody (1:50 dilution in PBS, Abcam #ab39012) or normal rabbit IgG (Cell Signaling Technologies #2729) was incubated overnight at 4°C in a humidity chamber. The Multivalent Secondary Antibody was incubated for 10 min, and next, the Peroxidase (HRP) was incubated for 10 min. Finally, the DAB substrate was applied for 5 min, and enhanced with DAB enhance for 3 min. Sections were counterstained with hematoxylin (Vector #H-3404) and cover slipped using Advantage Mounting Medium (Innovex NB300). For SOST, the VectaStain ABC Goat IgG Kit (Vector Labs #PK-40005) and ImmPACT DAB (Vector #SK-4105) kits were used. Sections were deparaffinized and rehydrated, and no antigen retrieval was performed. Endogenous peroxidases were blocked with BLOX-ALL for 10 min, and background staining was blocked using 3% BSA with 1 drop of Normal Rabbit serum for 1hr. Primary anti-SOST antibody was prepared in the same blocking buffer (1:200) and incubated overnight at 4°C in a humidity chamber. After washing, the rabbit anti-goat secondary antibody was incubated for 30 min, followed by the ABC Elite reagent for 30 min, and finally the DAB reaction for 4min. Slides were counterstained with hematoxylin (Vector Labs #H-3404), dehydrated, cleared with Xylenes, and cover slipped with VectaMount (Vector #H-5000-60).

Histological analysis and cell counting

Slides were imaged on a NanoZoomer 2.0-HT system (Hamamatsu) in brightfield at 40X. Images were blinded and analyzed with either BIOQUANT (OSTEO II) or ImageJ (Fiji) to quantify positive cells staining on the periosteal surface and intracortically (osteocytes). For the surface percent positive analysis, BIOQUANT was used to analyze the entire bone perimeter. If any part of the cell was positive or displayed puncta, that cell length was traced as positive. For the total surface expression analysis, ImageJ was used to manually isolate the periosteum, and the default color deconvolution (H-DAB) was used to isolate the periosteal DAB signal. Finally, after converting the DAB channel image to binary, the pixels were counted, converted to an area, and normalized per surface length. For the osteocyte analysis, in BIOQUANT, a 600 µm × 600 µm region of interest (ROI) was defined on the tibial apex, and positive and negative osteocytes (OCY) within the ROI were counted. The bone area within the ROI was also measured, and using these measures, we calculated both the percent positive (posOCY/totOCY) and areal density (posOCY/ROI Area) measures. On the surface, only the periosteum was analyzed because the marrow often separated from the endosteal surface, preventing reliable analysis of this surface.

Statistical analysis

For both proteomics and RNA-seq in the baseline experiment, we defined DEGs and DEPs as targets that had a Benjamani-Hochberg adjusted p-value < 0.05. For the loading experiments, we compared loaded versus non-loaded samples within each age, and very few proteins reached the FDR <0.05. We performed a sensitivity analysis with a variety of parameters, and regardless of the criteria used, the overall trends in the number of loading-responsive proteins were consistent for each age and day combination (Supplementary Figure 11). Therefore, to facilitate downstream analysis on the most differentially expressed protein targets with loading, we relaxed to the following criteria: a linear fold-change cutoff of >1.1 (up or down) and an unadjusted p-value < 0.05. For the pathway analyses, a Fisher’s exact test was used with an FDR of 0.05. For the co-expression network analysis, the significance threshold was set using a Bonferroni correction (p = 0.05/44 = 0.0011). For the RNA ISH and other histological analyses, after observing consistent staining with the probes, we analyzed one slide per target (n = 7 per age). If the periosteum lifted from the slide, that slide was excluded from the periosteal surface analysis, resulting in n = 6 per age for some targets. An unpaired t-test was used to compare measures between young-adult and old bones (GraphPad Prism 9.0).