Muscle transcriptomic circuits linked to periarticular physiology in end-stage osteoarthritis

Abstract

The ability of individuals with end-stage osteoarthritis (OA) to functionally recover from total joint arthroplasty is highly inconsistent. The molecular mechanisms driving this heterogeneity have yet to be elucidated. Furthermore, OA disproportionately impacts females, suggesting a need for identifying female-specific therapeutic targets. We profiled the skeletal muscle transcriptome in females with end-stage OA (n = 20) undergoing total knee or hip arthroplasty using RNA-Seq. Single-gene differential expression (DE) analyses tested for DE genes between skeletal muscle overlaying the surgical (SX) joint and muscle from the contralateral (CTRL) leg. Network analyses were performed using Pathway-Level Information ExtractoR (PLIER) to summarize genes into latent variables (LVs), i.e., gene circuits, and link them to biological pathways. LV differences in SX versus CTRL muscle and across sources of muscle tissue (vastus medialis, vastus lateralis, or tensor fascia latae) were determined with ANOVA. Linear models tested for associations between LVs and muscle phenotype on the SX side (inflammation, function, and integrity). DE analysis revealed 360 DE genes (|Log₂ fold-difference| ≥ 1, FDR ≤ 0.05) between the SX and CTRL limbs, many associated with inflammation and lipid metabolism. PLIER analyses revealed circuits associated with protein degradation and fibro-adipogenic cell gene expression. Muscle inflammation and function were linked to an LV associated with endothelial cell gene expression highlighting a potential regulatory role of endothelial cells within skeletal muscle. These findings may provide insight into potential therapeutic targets to improve OA rehabilitation before and/or following total joint replacement.

Keywords: osteoarthritis; skeletal muscle; transcriptomics.

Figure 1.

Methodological workflow utilized for the complementary single-gene (Limma-Voom) and network-based (PLIER) analysis pipelines. Transcripts were filtered for low expression, protein-coding biotype, then by study design. Batch effects and outliers were assessed: two surgical (SX) samples were outliers, and these samples along with corresponding contralateral (CTRL) samples were removed from further analyses. Limma-voom was used to detect single-gene differences in muscle gene expression between the SX and CTRL limbs. For PLIER, counts were converted to counts per million then z-scored. An input matrix linking genes to biological pathways from published prior knowledge was constructed resulting in 4,515 available pathways. Singular value decomposition (SVD) was performed in PLIER and generated 14 latent variables (LVs). SX vs. CTRL and TKA vs. THA differences were assessed with analysis of variance (ANOVA). For SX only, PLIER scores for each LV were assessed using linear models to determine associations between LVs and individual phenotypes of interest, including gene expression as measured by qPCR, histologically assessed fibrosis, and muscle function. PLIER, Pathway-Level Information ExtractoR; THA, total hip arthroplasty; TKA, total knee arthroplasty.

Figure 2.

Volcano plot of differentially expressed genes in skeletal muscle between SX and CTRL limbs of participants with OA. In total, 360 genes were significantly different between the SX and CTRL limbs with 340 higher and 20 lower on the SX limb, respectively. Significance was determined at a threshold of FDR ≤ 0.05 and |Log₂FD| ≥ 1. Labeled genes represent the five highest log₂ fold-difference (Log₂FD) in the positive and negative directions. These genes are further explored in Table 2. CTRL, contralateral; FDR, false discovery rate; OA, osteoarthritis; SX, surgical.

Figure 3.

Annotated latent variables (LVs) and the top annotations associated with them. Prior knowledge inputs provided to PLIER consisted of canonical pathways (KEGG, REACTOME, MIPS), gene ontology (GO), fiber type-specific expression profiles, and vastus lateralis-derived mononuclear cell expression profiles (, –46). Where available, the top pathway association was assigned as the annotation for each LV. FAP cells, fibro-adipogenic cells; FBN1⁺, fibrillin-1 positive; KEGG, Kyoto Encyclopedia of Genes and Genomes; MIPS, Munich Information Center for Protein Sequences; PLIER, Pathway-Level Information Extractor; PCV, postcapillary venule.

Figure 4.

A: difference between the latent variable (LV) score on the SX and CTRL limb for LV5. Individual points are indicated for each participant and circles denote total knee arthroplasty (TKA) participants and triangles denote total hip arthroplasty participants (THA). Generally, LV5 was higher in muscle on the SX limb compared with the CTRL limb (P-adjusted = 0.006). B: heatmap illustrating the normalized expression of the top 50 genes in LV5 demonstrating the difference in expression across limbs (SX and CTRL) and surgery location (TKA or THA). Significance was set as an adjusted P value ≤ 0.05, Trends were determined as an adjusted P valued ≤ 0.10. CTRL, contralateral; SX, surgical.

Figure 5.

A: difference between the latent variable (LV) score on the SX and CTRL limb for LV9 annotating to fibro-adipogenic cell expression (FAP). Individual points are indicated for each participant, and circles denote total knee arthroplasty (TKA) participants and triangles denote total hip arthroplasty participants (THA). Generally, LV9 was higher in muscle on the SX limb compared with the CTRL limb (P-adjusted = 0.056, Trend). B: heatmap illustrating the normalized expression of the top 50 genes in LV9 demonstrating the difference in expression across limbs (SX and CTRL) and surgery location (TKA or THA). Significance was determined as an adjusted P value ≤ 0.05. Trends were determined as an adjusted P valued ≤ 0.10. CTRL, contralateral; SX, surgical.

Figure 6.

Associations of LV14 (annotated to endothelial cells) to markers of inflammation (A, qPCR-determined Fn14 expression, P value = 0.01, n = 18), knee extension strength (B, MVIC, P value = 0.08, n = 16), and knee extension power (C, power, P value = 0.06, n = 16). Sample sizes for strength and power are reduced because some participants were unable to complete the measurement. MVIC, maximum voluntary isometric contraction.