Abstract
Osteoarthritis affects millions worldwide, yet effective treatments
remain elusive due to poorly understood molecular mechanisms. While
genome-wide association studies (GWAS) have identified hundreds of
osteoarthritis-associated loci, identifying the genes impacted at each
locus remains challenging. We investigate alternative splicing using
RNA-sequencing data from 101 human chondrocyte samples treated with
phosphate-buffered saline or fibronectin fragment, an osteoarthritis
trigger. We identified 590 differentially spliced genes between
conditions, with FN-f inducing splicing events similar to those in
primary osteoarthritis tissue. CRISPR/Cas9 mimicking of an SNRNP70
splicing event observed in osteoarthritis induced an
osteoarthritis-like expression pattern. Integration with genotyping
data revealed 7188 splicing quantitative trait loci (sQTL) affecting
3056 genes, including 738 and 343 condition-specific sQTLs for resting
and fibronectin fragment, respectively. Colocalization with
osteoarthritis GWAS identified 6 putative risk genes. Our study
highlights the significant impact of alternative splicing in
osteoarthritis and provides potential therapeutic targets.
Subject terms: RNA splicing, Cartilage
__________________________________________________________________
The authors identify thousands of genetic variants affecting RNA
splicing in primary human chondrocytes and link several of them to
osteoarthritis risk using genome editing and computational approaches.
Introduction
Osteoarthritis (OA) is a chronic, debilitating joint disease affecting
over 500 million people worldwide^[45]1. Despite its prevalence,
effective treatments to prevent disease progression remain elusive, due
in part to our limited understanding of the underlying molecular
mechanisms^[46]2. Genetic factors play a substantial role in OA risk,
with an estimated heritability of over 50%^[47]3. Genome-wide
association studies (GWAS) have identified hundreds of loci associated
with OA risk^[48]4,[49]5, yet understanding their functional impact
remains challenging, as most associated single nucleotide polymorphisms
(SNPs) reside in non-coding genomic regions^[50]6,[51]7.
One powerful method to determine the mechanisms through which
disease-associated variants act is to identify quantitative trait loci
(QTL) for a feature of interest (e.g., gene expression). Once mapped,
those QTLs can be colocalized with GWAS risk variants to identify
putative causal mechanisms. We and others have colocalized expression
QTL data sets with OA GWAS to identify putative causal
genes^[52]4,[53]8,[54]9. Despite these successes, many OA GWAS loci
remain unexplained. One possible reason is that the disease-associated
variants might influence traits other than gene expression. This may be
due to alternative splicing, as existing literature has described
substantial contributions of both expression QTLs and splicing QTLs to
disease risk; however, these disease-associated eQTLs and sQTLs often
influence disease risk via different genes^[55]10–[56]12. While OA
tissue is characterized by distinct changes in splicing
patterns^[57]13; sQTLs have yet to be mapped in cartilage or other
OA-relevant cell types.
To address this critical need, we mapped splicing QTLs in primary human
chondrocytes, the cell type responsible for maintaining the
extracellular matrix of articular cartilage. Splicing QTLs were mapped
in chondrocytes both in the resting state and in response to
fibronectin fragment (FN-f), a known OA trigger. FN-fs are present in
the cartilage and synovial fluid of OA patients^[58]14 and in cell
culture studies have been shown to initiate a range of catabolic
signaling pathways characteristic of chondrocytes isolated from OA
tissue, mimicking the OA-like state^[59]15–[60]17. This approach allows
us to capture the transition from a normal to an OA phenotype and
identify variants that exert effects in both homeostatic and
pathological settings.
In this work, we analyzed RNA-seq data from 101 human chondrocyte
samples treated with either PBS (resting) or FN-f, as well as 16 OA
donor samples. We identified hundreds of genes that were differentially
spliced between PBS-treated and either FN-f treated or OA chondrocytes.
Differential splicing events induced by FN-f mimicked those seen in OA
tissue, underscoring the accuracy of our FN-f model. We used CRISPR to
mimic a splicing event observed in both FN-f treated normal
chondrocytes and in OA chondrocytes and found that it induced OA-like
expression changes. We identified 7188 splicing quantitative trait loci
(sQTL) affecting 3056 genes, over 1000 of which exhibited
condition-specific effects. Colocalization of these sQTLs with OA GWAS
identified 6 putative OA risk genes, including putative candidate
PBRM1. This approach demonstrates the potential of sQTL analysis to
uncover genetic factors contributing to OA susceptibility. Our findings
both improve the mechanistic understanding of the disease and identify
genes that can be further studied and potentially targeted for
therapeutic intervention.
Results
FN-f induces alternative splicing in human chondrocytes
To determine how FN-f stimulation influences RNA splicing, we
reanalyzed our existing RNA-seq data sets obtained from primary human
chondrocytes isolated from 101 non-OA donors (donor information in
Supplementary Data [61]1) treated with either PBS or FN-f, a known OA
trigger. Initial quality control revealed minor batch effects
associated with the RNA extraction kit and FN-f batch (Supplementary
Fig. [62]1a). We corrected these technical confounders using the limma
package^[63]18, resulting in improved sample clustering by biological
condition (Supplementary Fig. [64]1b). Differential analysis using
LeafCutter identified 974 differential splice junctions corresponding
to 590 genes (adjusted p < 0.05, |∆PSI| > 0.15, Supplementary
Data [65]2). Hierarchical clustering of differentially spliced intron
junctions revealed distinct patterns between PBS resting and FN-f
treated samples (Fig. [66]1a). Over ninety percent of these
differential splicing events were previously annotated in GENCODE
release 45 (GRCh38.p14) as reported by LeafCutter. (Fig. [67]1b).
Furthermore, both our annotated and unannotated splice junctions
exhibited strong enrichment for canonical splicing motifs
(Supplementary Fig. [68]1c). Together, these results support the
validity of our splicing analysis.
Fig. 1. FN-f induces alternative splicing events in human chondrocytes.
[69]Fig. 1
[70]Open in a new tab
a Heatmap of differential splicing events between PBS and FN-f treated
chondrocytes (n = 101 biologically independent non-OA donors), with
each donor providing one chondrocyte culture split for PBS and FN-f
treatment. Rows show PSI values for 621 intron junctions from
LeafCutter clusters (|∆PSI| > 0.15, adjusted p < 0.05), derived from
590 genes. Significance determined using two-sided
Dirichlet-multinomial regression with FDR correction by LeafCutter.
Ancestry: AMR (Ad-mixed American), AFR (African), EUR (European), SAS
(South Asian). b Pie chart illustrates the distribution of differential
alternative splicing intron junction (n = 974) types. c Protein-protein
interaction network of 590 differentially spliced genes using
high-confidence interactions from STRING database v12.0 (interaction
score > 0.900). Network shows genes with protein-protein interactions
above the confidence threshold. Colors highlight pathway groups, bold
circles indicate genes also differentially spliced in OA tissue, text
indicates genes differentially spliced in both FN-f and OA conditions.
d Venn diagram showing overlap between FN-f-induced splicing changes
(n = 101 donors) and OA-associated splicing changes (n = 16 OA donors
versus n = 101 non-OA controls). Statistical significance determined
using two-sided Dirichlet-multinomial regression with FDR correction by
LeafCutter default method. e Box plots showing ∆PSI in FN-f response
for introns differentially spliced in OA. Each dot represents an
individual splice junction (n = 237 for increased PSI in OA, n = 221
for decreased PSI in OA). The analysis compares OA vs PBS splicing
changes (from 16 OA and 101 non-OA biological donors) with FN-f vs PBS
changes (from the same 101 non-OA donors). Box plots show median
(center line), 25th and 75th percentiles (box bounds), minimum and
maximum values within 1.5× interquartile range (whiskers), with
outliers as individual points. Statistical test: one-sided Wilcoxon
signed-rank test against theoretical median of zero (p = 2.23 × 10^−8
for increased PSI group; p = 1.23 × 10^−6 for decreased PSI group). No
adjustments were made for multiple comparisons.
We performed pathway and Gene Ontology (GO) enrichment analysis on the
differentially spliced genes to understand the potential functional
implications of these splicing changes. For FN-f-stimulated
chondrocytes as compared to PBS controls, this analysis revealed
significant enrichment in several pathways and processes, some of which
had possible connections to OA pathogenesis (Supplementary
Fig. [71]2a). Key enriched pathways included signal transduction, FGFR1
signaling in disease, deleted in colorectal cancer-mediated attractive
signaling, and osteoclast differentiation. Enriched GO terms included
anatomical structure morphogenesis, actin cytoskeleton organization,
and cell motility.
To further interrogate the relationship between differentially spliced
genes, we performed protein-protein interaction network analysis and
detected 64 nodes with high-confidence interactions (Interaction
score > 0.900). This network highlighted subnetworks corresponding to
critical pathways such as GPCR signaling, signal transduction,
osteoclast differentiation, hedgehog signaling, and steroid hormone
biosynthesis (Fig. [72]1c). Taken together, these results demonstrate
that FN-f induces large-scale changes in alternative splicing that
converge on specific pathways and processes, many of which have
possible relevance to OA.
FN-f induces OA-like RNA splicing in human chondrocytes
To determine if our FN-f-treated chondrocytes exhibited similar
splicing patterns to those found in OA tissue, we performed RNA-seq on
chondrocytes isolated from cartilage obtained from 16 donors who
underwent knee replacement surgery. Comparison between our PBS-treated
non-OA chondrocytes and OA chondrocytes revealed 466 significantly
differentially spliced intron junctions corresponding to 321 genes
(adjusted p < 0.05, |∆PSI| > 0.15, Supplementary Fig. [73]2b, c and
Supplementary Data [74]2). Genes that were differentially spliced in OA
compared to non-OA chondrocytes (treated with PBS) were enriched in
multiple KEGG pathways, including MCM complex, Rho GTPase cycle, and
fatty acid metabolism (Supplementary Fig. [75]2d and Supplementary
Data [76]3) and GO terms, including positive regulation of GTPase
activity, cell morphogenesis, and regulation of filopodium assembly.
Seventy-five intron junctions (corresponding to 75 distinct genes) were
significantly altered in both FN-f-stimulated and OA chondrocytes
(Fig. [77]1d, right panel). 75% of the splicing changes shared between
OA and FN-f treated samples showed a concordant directional effect, and
the effect sizes were significantly correlated (R^2 = 0.27, p < 0.001,
linear model regression, Supplementary Fig. [78]2e). Several of these
OA-associated differential spliced genes were also present in the
protein-protein interaction networks generated from FN-f induced
splicing events (Fig. [79]1c in bold), suggesting potential functional
connections. This exploratory network analysis highlighted genes with
established roles in cartilage biology or processes relevant to OA
pathogenesis. For instance, SNRNP70 is a splicing factor that has been
implicated in rheumatoid arthritis^[80]19; MALAT1 is a lncRNA involved
in inflammation and chondrocyte proliferation^[81]20, and NCOR2 is a
nuclear receptor co-repressor that is hypermethylated and downregulated
in OA tissue^[82]21.
Next, we sought to determine the extent to which splicing events in
FN-f-treated chondrocytes mirror those seen in OA. Differential
splicing events induced by FN-f were significantly correlated with
those seen in OA tissue (linear regression, p < 0.001, Supplementary
Fig. [83]2f). To investigate this further, we divided the 466 introns
that were differentially spliced between OA and non-OA (PBS-treated)
chondrocytes into those that showed increased or decreased percent
spliced in (PSI) in OA samples. On average, introns that exhibited
increased PSI in OA tissue also exhibited increased PSI in response to
FN-f (Fig. [84]1e). Conversely, introns that exhibited decreased PSI in
OA tissue also exhibited decreased PSI in response to FN-f
(Fig. [85]1e). Our FN-f model mimics the early events in OA onset, and
we wouldn’t expect it to fully recapitulate the changes seen in
end-stage OA tissue; however, the agreement between splicing changes in
OA and FN-f-treated chondrocytes, in terms of both broad directional
changes and specific high-confidence alterations, suggests that the
FN-f-treated cells provide a valuable model to interrogate OA-related
alternative splicing events. The identified set of consistently altered
splicing events provides a focused group of targets for future studies
aimed at understanding and potentially intervening in OA progression.
Splicing events that differed between FN-f-treated and OA samples were
enriched for GO terms, including ion transport processes, specific
protein modifications, and metabolic pathways (Supplementary
Data [86]3), which provide insight into which aspects of OA are and are
not fully modeled via FN-f treatment.
Forced splicing of SNRNP70 induces OA-like gene expression
We identified SNRNP70, a critical component of the U1 snRNP complex, as
differentially spliced in both FN-f-stimulated and OA chondrocytes
compared to PBS (resting). RNA-seq analysis revealed significant
changes in SNRNP70 splicing, particularly between exons 7 and 8
(Fig. [87]2a). The alternative exon 8 (chr19:49,102,114-49,103,587)
showed reduced inclusion in FN-f and OA conditions, with ∆PSI
differences of 22% (FN-f/PBS) and 23% (OA/PBS). This indicates
alternative exon 8 skipping in 67% and 68% of transcripts in FN-f and
OA conditions, respectively, compared to 45% in PBS (resting).
Exon-level expression analysis confirmed significant downregulation
(p < 0.01) of the alternative exon 8 in both conditions (Fig. [88]2b).
This differential splicing event was confirmed via qPCR (Supplementary
Fig. [89]3). Western blotting for SNRNP70 revealed a single band
representing the expected protein size, with reduced protein levels
under FN-f conditions matching the transcriptomic data (Supplementary
Fig. [90]4).
Fig. 2. Deletion of SNRNP70 alternative exon 8 mimics OA-related splicing
patterns.
[91]Fig. 2
[92]Open in a new tab
a RNA-seq signal track for PBS-treated, FN-f-treated, and OA
chondrocytes at the SNRNP70 gene locus. Below the signal tracks are
transcript annotations representing the MANE Select transcript
(ENST00000598441.6) and a transcript containing the alternative exon 8
(ENST00000401730.5). b Barplot represents log[2] fold change for
expressed exons with FN-f and OA compared to PBS. c Schematic of
CRISPR/Cas9 genome editing strategy for SNRNP70 alternative exon 8
deletion. d Barplot comparing log[2] fold change for expressed exons
after SNRNP70 alt ex8 deletion. e Barplot of KEGG and Reactome pathways
and GO terms enriched in differentially expressed genes between
alternative exon 8 (alt ex 8) heterozygous knockdown (Het-KD) and
wildtype (WT). Enrichment was calculated using cumulative
hypergeometric distribution of overlap genes to the total background
gene set, with statistical significance indicated by -log[10] of the
p-value. One-sided hypergeometric tests were used to identify
significant enrichment, without adjustment for multiple comparisons. A
total of 135 differentially expressed genes (padj <0.01,
|Log[2]FoldChange| > 2) were analyzed, with 71 upregulated and 64
downregulated genes in Het-KD compared to WT. f Box plots showing gene
expression changes (log[2]FoldChange) in FN-f vs PBS for genes that
were differentially expressed in the SNRNP70 alternative exon 8
heterozygous knockdown (Het-KD) vs wildtype (WT) comparison. The
analysis examines whether genes upregulated in Het-KD/WT also tend to
be upregulated in FN-f/PBS, and whether genes downregulated in
Het-KD/WT tend to be downregulated in FN-f/PBS. Each dot represents an
individual gene (n = 60 for genes upregulated in Het-KD vs WT, n = 57
for genes downregulated in Het-KD vs WT). Box plots show median (center
line), 25th and 75th percentiles (box bounds), 5th and 95th percentiles
(whiskers), and outliers as individual points. Statistical analysis:
one-sided Wilcoxon signed rank test with continuity correction
(upregulated genes: p-value = 0.008283; downregulated genes:
p-value = 0.01869). No adjustments were made for multiple comparisons.
To assess the functional role of SNRNP70 alternative exon 8, we
employed CRISPR-Cas9 to delete this region in primary human
chondrocytes (Fig. [93]2c). Our deletion strategy targeted guide RNAs
to both sides of a 1735 bp region encompassing the alternative exon 8
(Supplementary Fig. [94]5a). PCR screening of 46 single-cell-derived
colonies revealed 26% with heterozygous and 0% with homozygous deletion
of the alternative exon 8 (Supplementary Fig. [95]5b, c). Sanger
sequencing of Het-KD colonies confirmed the deletion, with an
additional cytosine deletion beyond the expected cut site
(Supplementary Fig. [96]5d).
To determine the phenotypic impact of heterozygous exon 8 deletion, we
performed RNA-seq on both edited and unedited colony-expanded cells.
The differential analysis confirmed decreased SNRNP70 alternative exon
8 inclusion (Fig. [97]2c) and identified 135 differentially expressed
genes (DESeq2, FDR-adjusted p-value < 0.01). Genes that were
upregulated in Het-KD cells were enriched for multiple GO terms and
KEGG pathways, including IL-17 signaling, chemokine receptor binding,
and rheumatoid arthritis (Fig. [98]2e and Supplementary Data [99]4).
Genes that were downregulated in Het-KD cells were enriched for
pathways including L1 and Ankyrin interactions and O-glycan
biosynthesis (Fig. [100]2e and Supplementary Data [101]4). Differential
splicing analysis revealed 259 differential splicing events (233
differential splicing genes) in Het-KD cells, although these changes
did not clearly correlate with those observed in FN-f treated cells
(Supplementary Fig. [102]5e).
Interestingly, the directionality of gene expression changes in Het-KD
cells significantly correlated with those in FN-f-stimulated
chondrocytes (p < 0.05 for both up- and downregulated genes,
Fig. [103]2f). Genes that were upregulated both in response to FN-f and
in response to alternative exon 8 deletion included CXCL8, CXCL2, IL11,
MMP12, and MMP13, which are associated with inflammation and
extracellular matrix degradation in OA^[104]22–[105]24. Downregulated
genes included PEG3, FGL2, SEMA3E, CSRNP3, and HLF, which have been
implicated in chondrocyte homeostasis^[106]25–[107]29. These results
suggest that the exclusion of SNRNP70 alternative exon 8 in response to
cartilage matrix damage modeled by FN-f treatment may be responsible
for some of the other expression changes associated with the OA
phenotype.
Genetic variants impact splicing in chondrocytes
We performed splicing QTL (sQTL) analysis on both PBS (resting) and
FN-f stimulated chondrocyte samples to determine the impact of genetic
differences on chondrocyte mRNA splicing. We tested the association of
each RNA splicing event with genetic variants within ±100 kb of the
start and end points of the splice intron junctions, using PSI values
adjusted for confounding factors (Supplementary Fig. [108]6a). We
identified a total of 5873 unique sQTLs corresponding to 2575 sGenes
(genes with at least one significant sQTL) in PBS resting, and 4606
sQTLs corresponding to 2132 sGenes in FN-f stimulated chondrocytes
(QTLtools^[109]30, FDR < 0.05, Fig. [110]3a and Supplementary
Data [111]5). Over half (3400 in PBS and 2631 in FN-f) of sQTLs were
not previously detected in any GTEx tissue (Supplementary
Fig. [112]6e). We attempted to perform sQTL analysis on our 16 OA donor
samples, but due to the very small sample size, we identified only a
single sQTL.
Fig. 3. Condition-specific sQTL discovery in resting and activated
chondrocytes.
[113]Fig. 3
[114]Open in a new tab
a Venn diagram showing the overlap between sGenes identified in PBS and
FN-f treated conditions. Linear mixed-effect modeling allowed the
detection of 677 condition-specific sGenes (477 PBS-specific and 200
FN-f-specific). b The distribution of distances between sSNPs and the
affected splice intron junctions (kb). c Percentage of sSNPs located
within the same gene (within a gene) or outside the gene. d–f Examples
of d shared (RHBDF2), e PBS-specific (SLC26A4), and f FN-f-specific
(MAPK8) sQTLs. Each panel contains locus zoom plots showing the
-log[10](p-value) of association between SNPs and splicing across
genomic positions for both PBS and FN-f conditions. Gene structures are
displayed beneath each locus zoom plot. Each panel also includes box
plots displaying the relationship between genotype and inverse normal
transformed PSI values for the indicated splice junction in both PBS
and FN-f conditions (n = 101 independent donor samples per condition).
Boxes indicate the interquartile range (25th–75th percentile) with a
line at the median (50th percentile). Whiskers extend to the minimum
and maximum values within 1.5× the interquartile range, with outliers
shown as individual points. The figure shows the lead SNP ID, effect
size (Beta), and p-value for each condition. Shared sQTLs show
significant associations in both PBS and FN-f conditions. PBS-specific
sQTLs show significant association only in PBS but not in FN-f
condition. FN-f-specific sQTLs show significant association only in
FN-f but not in PBS condition. Nominal p-values were calculated using
QTLtools, which performs two-sided tests for association by converting
correlation coefficients to p-values using a t-distribution of freedom
(n-2, where n = 101 samples).
Conditional analysis revealed 182 secondary signals, 13 tertiary
signals, and 2 quaternary signals for PBS-treated chondrocytes and 121
secondary signals and 1 tertiary signal for FN-f-treated chondrocytes
(Supplementary Fig. [115]6b). In total, 116 PBS sGenes and 82 FN-f
sGenes had two or more independent signals. For example, the CAST gene
exhibited multiple independent signals associated with the splice
intron junction chr5:96768431–96776404 in PBS resting, while the MICA
gene exhibited multiple independent signals with the splice intron
junction chr6:31410797–31411209 in FN-f stimulated chondrocytes
(Supplementary Fig. [116]7). To further characterize the genetic
architecture of splicing regulation, we systematically assessed allelic
heterogeneity across genes with multiple sQTLs. By examining linkage
disequilibrium (LD) patterns (R^2 < 0.7) between sQTLs in the same
gene, we identified 760 genes in PBS tissue (29.5% of 2575 sGenes) and
576 genes in FN-f tissue (27.0% of 2132 sGenes) in which two or more
independent variants influenced splicing in the same gene.
Additionally, we investigated the relationship between splicing
regulation and gene expression by analyzing the overlap between our
sQTLs and published eQTLs^[117]9 (Supplementary Data [118]5). We found
that the vast majority of our sQTLs (5557 in PBS and 4302 in FN-f)
showed no overlap (R^2 < 0.7) with eQTLs affecting the same gene,
suggesting that the genetic variants influencing splicing at these loci
did not also influence expression levels. While a smaller subset of
sQTLs (513 in PBS and 426 in FN-f) did overlap with eQTLs. Retesting
for sQTLs while conditioning on the genotype of the lead eQTL SNP
suggests that over 99% of these sQTLs that were in high LD with an eQTL
(512 out of 513 in PBS and 424 out of 426 in FN-f) also induced changes
in gene expression. In summary, approximately 10% of our detected sQTLs
also influenced gene expression.
We analyzed the genomic distribution of sQTLs (sSNPs) relative to their
associated splice intron junctions. Lead sSNPs were largely close to
splice junctions, with 43% and 45% of lead SNPs within 5 kb of the
affected splice junctions for PBS and FN-f, respectively
(Fig. [119]3b). Further analysis revealed that the majority of sSNPs
were located within the same gene as their associated splice intron
junction (Fig. [120]3c). Taken together, these results both support the
validity of our sQTLs and suggest that sSNPs largely act at the
affected splice junctions rather than via distal mechanisms.
To determine if any of the effects of variants were dependent on the
condition, we used a linear mixed-effect model and identified 1635
sGenes that exhibited significant genetic and condition interaction
effects (lme4^[121]31, ANOVA, p < 0.05). To compile a high-confidence,
large-effect size list of condition-specific sQTLs, we further refined
this list to include only those that were identified as an sQTL in only
one condition and exhibited an absolute beta difference between
conditions of greater than or equal to 0.2. In total, we identified 677
high-confidence condition-specific sGenes (Fig. [122]3a). 400 lead SNPs
had a stronger genetic effect on splicing in PBS-treated cells, and 277
lead SNPs had a stronger genetic effect on splicing in FN-f-treated
cells. Examples of shared, PBS-specific, and FN-f-specific sQTLs are
depicted in Fig. [123]3d–f. Several of these condition-specific sGenes
have known roles in cartilage biology and/or OA. For example,
Fig. [124]3e shows a PBS-specific condition sQTL for SLC26A4 involving
the splice junction chr7:107694480-107694621. SLC26A4, an anion
transporter, may influence chondrocyte metabolism in homeostatic
conditions, potentially affecting cartilage maintenance^[125]32.
Fig. [126]3f demonstrates an FN-f-specific sQTL for MAPK8, associated
with the splice junction chr10:48363284-48401612. MAPK8 (JNK1) is
involved in inflammatory signaling, and JNK inhibition is currently
being investigated as a treatment for OA^[127]33,[128]34. However,
knockout of Jnk1 or Jnk2 increased the senescence burden and the extent
of age-related osteoarthritis in mice^[129]35. The complex role of JNK
signaling in OA may result from high activity in a rare sub-population
of chondrocytes that serves an “inflammation-amplifying” role in OA
cartilage^[130]36.
RNA-binding proteins reveal putative sQTL mechanisms
RNA-binding proteins (RBPs) play crucial roles in post-transcriptional
regulation, influencing various cellular processes, including
inflammation and metabolic changes related to OA^[131]37,[132]38.
Recent studies have highlighted the dysregulation of RBPs in
chondrocytes, synovial fibroblasts, and osteoblasts, indicating their
significant impact on OA progression^[133]39,[134]40. To gain insight
into the putative mechanisms through which our sQTLs were influencing
splicing in chondrocytes, we leveraged publicly available enhanced
crosslinking and immunoprecipitation (eCLIP) data from the ENCODE
database^[135]41. This dataset encompasses binding sites for 140 RBPs
across three cell types: HepG2 (n = 85), K562 (n = 118), and SM-9MVZL
(n = 2). While these cell lines are not OA-specific, similar data is
not currently available in chondrocytes, and presumably, many of these
binding sites are shared across cell types.
Using the QTLtools fenrich module, we assessed the statistical
significance of RBP enrichment through permutation tests. For PBS
(resting) lead sQTLs, we identified 12 significantly enriched RBPs
(odds ratio > 1, empirical p < 0.05), including SF3B4, EFTUD2, and
PRPF8. The FN-f condition yielded 7 enriched RBPs, including UP3,
EFTUD2, and SF3B4 (Fig. [136]4a and Supplementary Data [137]6). We
found 255 sQTL lead variants that were located within the binding site
of at least one RNA-binding protein (RBP) and showed significant
correlation between the splicing pattern and RBP gene expression
(Pearson’s R^2 > 0.15, p < 0.01; Supplementary Data [138]6).
Fig. 4. RNA-Binding protein associations with sQTLs.
[139]Fig. 4
[140]Open in a new tab
a Volcano plot depicting the effect size (x-axis) and -log[10] p-values
(y-axis) of RBP sites enrichment in sQTLs identified in PBS- or
FN-f-treated chondrocytes. P-values and odds ratios were calculated
using QTLtools Fenrich module with 1000 permutations (“Methods”). The
enrichment analysis uses a two-sided empirical permutation test to
determine whether RBP binding sites are over- or under-represented in
regions containing sQTLs compared to control regions, without
adjustment for multiple comparisons. Data points are shown for
p-value < 0.05 and odds ratios > 1. b Analysis of AATF RBP binding and
expression correlation with the rs11871958 sQTL in the SNHG29 gene.
Locus zoom plots show -log[10](p-value) of SNP-splicing associations
across the genomic region in PBS and FN-f conditions (PBS:
p = 1.48 × 10^−21, β = −0.9; FN-f: p = 1.06 × 10^−24, β = −0.99, minor
allele frequency = 0.1832). These association statistics were
calculated using QTLtools. The middle track displays AATF eCLIP signal
(RPM) with the sQTL SNP location highlighted. Box plots show normalized
PSI values by rs11871958 genotype (n = 202 samples from 101
biologically independent donors, each providing PBS and FN-f samples).
Box plots display median (center line), 25th–75th percentiles (box),
minimum-maximum values within 1.5× interquartile range (whiskers), and
outliers (points). The scatter plot shows correlation between AATF gene
expression (VST-normalized DESeq values) and normalized PSI values
across all samples (R^2 = 0.408, p = 2.22 × 10^−28). Pearson
correlation was used with a two-sided linear regression model to
calculate the p-value, without adjustment for multiple comparisons. c
Analysis of EFTUD2 RBP binding and expression correlation with the
rs6809764 sQTL in the PCYT1A gene, following the same format as panel
b. QTLtools statistics for this sQTL: PBS (p = 4.65 × 10^−8, β = 0.44);
FN-f (p = 2.29 × 10^−7, β = 0.522, minor allele frequency = 0.3911).
The correlation between EFTUD2 gene expression and normalized PSI
values (n = 202 samples) shows R^2 = 0.164, p = 7.84 × 10^−9,
calculated using the same statistical methods.
An example of a putative RBP-mediated sQTL is rs11871958, which is
associated with alternative splicing of SNHG29 and overlaps an AATF RBP
binding site (Fig. [141]4b). The minor C allele of rs11871958 is
associated with decreased PSI of the associated intron junction
(chr17:16439414-16439528). Interestingly, AATF gene expression is
strongly correlated with the PSI of the associated intron junction,
lending further support to the role of AATF in this splicing event
(R^2 = 0.458, p = 2.233 × 10^−28). SNHG29, a long non-coding RNA, has
been previously reported to regulate the miR-223-3p/CTNND1 axis and
influence the Wnt/β-catenin signaling^[142]42, which has been
previously implicated in OA pathogenesis^[143]43.
Another example is rs6809764, which was identified as an sQTL for
PCYT1A and overlaps an EFTUD2 binding site (Fig. [144]4c). The minor G
allele of rs6809764 is associated with increased PSI of the intron
junction (chr3:196270541-196287615). Gene expression of EFTUD2
correlated with decreased junction usage (R^2 = 0.154,
p = 7.394 × 10^−9), which supports EFTUD2 as a possible mediator of
this sQTL. PCYT1A encodes an enzyme crucial for phosphatidylcholine
biosynthesis^[145]44. It has been observed to localize to the nuclear
envelope in hypertrophic chondrocytes of bone growth plates, coinciding
with a significant increase in cell volume^[146]44. This localization
suggests a potential role for PCYT1A in chondrocyte morphology changes
during hypertrophy, a process relevant to cartilage development and OA
pathogenesis^[147]45.
These findings provide insights into the regulatory mechanisms of
chondrocyte splicing events. However, we acknowledge the limitations of
using non-chondrocyte cell lines for RBP binding data. Further
biochemical studies in relevant cell types will be crucial to fully
elucidate the mechanisms through which these genetic variants alter
splicing patterns in cartilage tissue.
Colocalization analysis reveals putative OA risk genes
To understand the potential role of splicing in OA risk, we performed
colocalization analysis between our sQTLs and 100 independent OA GWAS
loci spanning 11 OA-related phenotypes reported by Boer et al.^[148]4.
We selected sQTLs whose lead SNP was in moderate to high LD (r^2 > 0.5)
with a lead SNP from any of the 11 OA GWAS data SNPs and tested for
colocalization using the coloc package (Supplementary Data [149]7). In
total, we identified 6 colocalized signals corresponding to 6 unique
sGenes across 4 OA phenotype subtypes (Table [150]1 and Fig. [151]5).
Table 1.
Chondrocyte sQTLs that colocalize with OA GWAS loci
Gene Intron Junction rsID Position Minor Allele MAF PPH3 PPH4 (PPH4)/
(PPH3 + PPH4) GWAS phenotype rsID (GWAS) sQTL Condition GTEX Tissue
Count
RNF144B chr6:18399699–18401957 rs1886248 chr6:18399163:C:T C 0.4059
0.00267 0.997 0.997 FingerOA rs9396861 Both 14
WWP2 chr16:69925484–69929448 rs7192245 chr16:69921902:T:C T 0.2228
0.179 0.714 0.780 KneeOA rs34195470 PBS 0
COLGALT2 chr1:183978520–184037581 rs74767794 chr1:184036994:A:G G 0.297
0.0302 0.921 0.9683 THR rs1327123 PBS 0
GCAT chr22:37808163–37809949 rs13057823 chr22:37792277:C:A A 0.223
0.0598 0.924 0.9392 THR rs12160491 Both 47
HMGN1 chr21:39347292–39347379 rs2249666 chr21:39347315:C:T T 0.451
0.0994 0.786 0.8877 TJR rs9981884 PBS 31
PBRM1 chr3:52658315–52662133 rs7628578 chr3:52594305:A:T T 0.4653 0.270
0.717 0.7263 THR rs3774354 Fn-f 10
[152]Open in a new tab
The table shows significant colocalization results between sQTL and
GWAS signals (PP H4 > 0.7) with the most significant splice intron
junction. The number of GTEx tissues (out of 50) where the same sQTL is
detected is also indicated. PP posterior probability, FingerOA
osteoarthritis of the finger, KneeOA osteoarthritis of the knee, THR
total hip replacement.
Fig. 5. Six sQTLs colocalize with OA GWAS signals.
[153]Fig. 5
[154]Open in a new tab
a Total hip replacement (THR) GWAS signal colocalizes with an sQTL
found in both conditions for PBRM1 associated with splice intron
junction chr3:52658315–52662133. b THR GWAS signal colocalizes with an
sQTL found in both conditions for GCAT associated with splice intron
junction chr22:37808163–37809949. c Total joint replacement (TJR) GWAS
signal colocalizes with an sQTL found in both PBS-treated chondrocytes
for HMGN1 associated with splice intron junction
chr21:39347292–39347379. d Finger OA GWAS signal colocalizes with an
sQTL found in both conditions for RNF144B associated with splice intron
junction chr6:18399699–18401957. e Knee OA GWAS signal colocalizes with
an sQTL found in PBS-treated chondrocytes for WWP2 associated with
splice intron junction chr16:69925484–69929448. f THR GWAS signal
colocalizes with an sQTL found in PBS-treated for COLGALT2 associated
with splice intron junction chr1:183978520–184037581.
To the best of our knowledge, these splicing events have not previously
been colocalized with OA GWAS signals; however, three of these genes
have been implicated in OA via colocalization with expression QTLs.
GCAT and HMGN1 have been previously associated with knee and/or hip OA
via colocalization with an osteoclast-specific eQTL^[155]46. RNF144B
was previously linked to OA via colocalization between OA GWAS variants
and eQTLs from chondrocytes, tibial nerves, and testis^[156]4,[157]9.
Our results build upon these previous studies by providing more
mechanistic detail into how these variants act, presumably via their
impact on alternative splicing.
Three of our colocalized genes have not been previously linked to OA
via colocalized eQTLs but do have supporting evidence for a role in OA.
WWP2 and COLGALT2 have been implicated in OA by their proximity to
mQTLs that colocalize with OA GWAS variants and control methylation
sites near WWP2 and COGALT2^[158]47. Methylation at these loci was
confirmed to affect WWP2 and COGALT2 via epigenome
editing^[159]48,[160]49. To the best of our knowledge, one of our
colocalized sGenes (PBRM1) has not been previously linked to OA via
genetic methods but was found to be hypomethylated in OA
cartilage^[161]50,[162]51.
These findings highlight the potential role of alternative splicing in
mediating OA risk and demonstrate the value of integrating sQTL data
with GWAS results to uncover the potential mechanisms underlying
genetic associations.
Discussion
This study establishes a comprehensive sQTL map in primary
chondrocytes, providing insights into the relationship between genetic
variation, alternative splicing, and OA risk. By integrating sQTL
analysis with OA GWAS data and functional validation, we have
identified potential mechanisms underlying OA pathogenesis and
highlighted putative OA risk genes for further investigation and future
therapeutic development.
The splicing alteration observed in both FN-f stimulated and OA
chondrocytes highlights the role of alternative splicing in OA
progression. Our identification of 974 and 466 differentially spliced
intron junctions in FN-f-stimulated and OA chondrocytes, respectively,
suggests broad splicing dysregulation in OA. The overlap between
FN-f-induced and OA-associated splicing events supports the relevance
of our ex vivo model and suggests that some of these splicing changes
may contribute to OA pathogenesis. However, while the overlap between
FN-f treated cells and cells from OA tissue is significant, their
correlation is modest, which may be due to the small number of OA
samples available for this study. In addition, the OA tissue is
obtained from patients undergoing knee replacement, which reflects the
end-stage of the disease, while the FN-f stimulus of normal
chondrocytes may induce splicing events important at earlier stages of
OA. Therefore, while FN-f treatment is a controllable OA-relevant
stimulus that provides insights into the genetics of OA risk, it does
not fully recapitulate the splicing changes observed in end-stage OA.
Our functional interrogation of SNRNP70 alternative splicing provides
evidence for the impact of splicing changes on OA-related gene
expression. The heterozygous deletion of SNRNP70 alternative exon 8 was
a mimic of the reduced retention of alternative exon 8 that occurs
during OA and in response to FN-f. While deletion of Snrnp70 itself is
embryonic lethal in mice^[163]52, it was anticipated that chondrocytes
could tolerate full deletion of alternative exon 8. The lack of
colonies with homozygous alternative exon 8 deletions may represent a
growth disadvantage as compared to intact or heterozygous deletion. The
exact mechanism by which this splicing event alters chondrocyte biology
is not entirely clear, but our data and the literature highlight
several possibilities. Alternative exon 8 is located in the RNA
recognition motif domain, which is critical for the protein RNA binding
function, and its inclusion could alter the structure or function of
the resulting protein. Alternatively, this splicing event could
influence overall SNRNP70 expression levels, as transcripts that
include alternative exon 8 are predicted by GENCODE v45 to induce
nonsense-mediated decay (NMD)^[164]53. While the strongest impact of
our CRISPR experiments was a sharp decrease in alternative exon 8
usage, we also observed small but significant decreases in the usage of
exons 9 and 10. This pattern is consistent with potential effects on
transcript stability or NMD processes following exon 8 deletion, though
the precise mechanism linking the deletion of alternative exon 8 to
reduced usage of downstream exons requires further investigation.
Alternative exon 8 deletion led to alterations in genes and pathways
crucial for OA progression, including inflammatory signaling and matrix
degradation, and suggests that alternative spliced genes may influence
cellular phenotypes relevant to OA. The identification of other
differentially spliced genes in both FN-f stimulated and OA
chondrocytes, such as STS, ANKRD36C, and PIK3CD, provides potential
targets for further investigation in OA. Interestingly, while the
strongest impact of our CRISPR experiments was a sharp decrease in
alternative exon 8 usage, we also observed small but significant
decreases in exon 9 and exon 10 usage. The direct impact of alternative
8 removal vs its downstream effect on exons 9 and 10 will require
further investigation.
Our sQTL analysis identified genetic variants influencing splicing in
both resting and activated chondrocytes. The prevalence of local
regulatory effects suggests that cis-acting splicing regulation
contributes to how genetic variation influences OA risk. We identified
hundreds of sQTLs that are influenced by condition (i.e., PBS vs FN-f),
highlighting the need to conduct QTL studies in the correct biological
context. What drives these condition-specific differences is
challenging to determine without further experiments, but could be
explained by differences in expression, splicing, nuclear localization,
or post-translational modification of the splicing factors that mediate
these events. Our RBP binding analysis identified potential regulatory
proteins that may mediate these sQTL effects, with 137 of 140 analyzed
RBPs showing detectable expression in chondrocytes. We identified 12
and 7 significantly enriched RBPs in PBS and FN-f conditions,
respectively, including key spliceosome components such as SF3B4,
EFTUD2, and PRPF8. The examples of AATF-mediated regulation of SNHG29
splicing and EFTUD2 potential role in PCYT1A splicing provide insights
into putative mechanisms underlying our sQTLs. While these associations
require experimental validation, and our analysis was limited by
reliance on non-chondrocyte eCLIP data, these findings illuminate the
complex regulatory networks that may drive chondrocyte splicing changes
and provide promising candidate genes for future mechanistic studies.
The colocalization analysis between our sQTLs and OA GWAS signals
identified 6 colocalized signals. Three of these colocalized genes have
been previously linked to OA via colocalized eQTLs, but our results
build on this knowledge by suggesting that their pathogenic mechanism
may act via alterations to mRNA splicing. While WWP2 has been
functionally linked to OA through knockout studies, both WWP2 and
COLGALT2 are implicated in OA by their proximity to colocalized mQTLs
and have been shown to exhibit allelic expression imbalance in
cartilage transcriptome studies, despite the absence of direct eQTL
associations^[165]47,[166]54.
One of our colocalized genes, PBRM1, provides a new potential OA
mediator that has not previously been implicated by hypothesis-driven
studies or genomic analysis. The rs3774354 GWAS signal was initially
assigned to ITIH1 based on genomic proximity^[167]4. Earlier studies of
the rs6976/rs11177 signal at the same locus implicated different genes,
including GNL3 and SPCS1^[168]55, which showed allelic expression
imbalance correlating with that OA risk genotype, while regulatory
activity was demonstrated for SNPs near both ITIH1 and GLT8D1^[169]56.
However, our colocalized sQTL suggests that altered splicing of PBRM1,
rather than effects on ITIH1, may be the mechanism driving OA risk at
this locus. As a SWI/SNF chromatin-remodeling complex
component^[170]57, PBRM1 influences gene expression patterns. Its role
in mesenchymal stromal cell differentiation is particularly intriguing
in the context of OA. Previous research has shown that PBRM1 affects
osteogenic differentiation and modulates SOX9 expression, a key
regulator of chondrogenesis and expression of chondrocyte matrix
genes^[171]58,[172]59. These observations suggest that PBRM1 might be
involved in maintaining the delicate balance between different cell
fates in joint tissues, a process often disrupted in OA.
Interestingly, GCAT and HMGN1 have been associated with hipOA, knee
and/or hip OA and/or all OA in osteoclast-specific eQTL^[173]46. Our
findings extend this by providing evidence for their involvement in
chondrocyte splicing regulation. GCAT, encoding glycine
C-acetyltransferase, plays a crucial role in glycine metabolism in
mitochondria^[174]60. Previous studies have shown that epigenetic
downregulation of GCAT can lead to respiratory defects in aged
fibroblasts, and glycine supplementation has been found to restore
these defects^[175]60. While the direct link between GCAT and OA
pathogenesis remains to be fully elucidated, our identification of GCAT
through sQTL analysis suggests that alterations in mitochondrial
function and glycine metabolism may contribute to OA risk, possibly
through effects on cellular energy production or extracellular matrix
synthesis. HMGN1 is a known non-histone chromatin remodeler and has
been implicated in the regulation of chondrogenesis^[176]61. Furusawa
et al. demonstrated that HMGN1 directly regulated the expression of
SOX9. This suggests that alternative splicing may modulate HMGN1
function and regulate chondrocyte-specific gene expression.
RNF144B showed different patterns of condition specificity in our sQTL
analysis compared to previous eQTL studies. We observed strong RNF144B
sQTLs in both PBS- and FN-f-treated conditions, whereas the eQTL for
RNF144B was highly specific to the FN-f-treated condition^[177]9. This
discrepancy highlights how sQTL and eQTL analyses can provide
complementary information for mRNA production and processing. Given
RNF144B role in regulating DNA damage-induced apoptosis, its
differential splicing could affect chondrocyte survival and turnover in
OA.
WWP2 has been previously implicated in OA through methylation studies.
In particular, Roberts et al. identified a differentially methylated
region (DMR) associated with the OA risk-conferring allele of
rs34195470^[178]48. This DMR was found to act as a
methylation-sensitive transcriptional repressor, with increased DNA
methylation resulting in increased WWP2 expression, particularly of the
WWP2-FL and WWP2-N transcripts. Our colocalization of a WWP2 sQTL with
an OA GWAS signal provides a direct connection to OA and adds another
layer to its regulatory complexity, potentially indicating that OA risk
variants alter splicing via changes to DNA methylation. Further work is
necessary to explore this hypothesis.
COLGALT2 has also been previously associated with OA
risk^[179]49,[180]62. Previous research has identified that the OA risk
allele of rs11583641 corresponds to reduced DNA methylation at nearby
CpGs and increased COLGALT2 expression^[181]63. Our sQTL findings now
implicate alternative splicing as another potential mechanism by which
COLGALT2 variation may influence OA risk, further highlighting the
complex regulatory landscape of this gene in OA pathogenesis.
While our study represents a significant advance in understanding the
role of splicing in OA, it has several limitations that should be
addressed in future research. Our sample size of 101 donors, while
substantial for an sQTL study in chondrocytes, may not have captured
all splicing-related OA GWAS genes. Future studies with larger sample
sizes and expanded diversity will be crucial to fully elucidate the
splicing landscape in OA. Additionally, while chondrocytes are highly
relevant to OA, investigating splicing changes in other cell types
involved in joint biology will provide a more comprehensive
understanding of OA pathogenesis. Our FN-f stimulation model, while
valuable, cannot fully recapitulate the complex joint environment.
Advanced 3D culture systems, such as cartilage organoids or
joint-on-a-chip models, could provide more physiologically relevant
contexts for validating our findings. Furthermore, our study focused on
a single time point of FN-f stimulation and collection of end-stage OA
cartilage. Future studies considering other time-points or from
patients at different stages of OA progression may capture splicing
changes do not present here. While our functional validation of SNRNP70
provides valuable insights, further experimental work is needed to
understand the detailed functional consequences of the colocalized
sQTLs we identified. Potential approaches for these investigations
could include CRISPR-Cas9 mediated modulation of splice sites, minigene
assays to confirm splicing effects, and functional assays to assess the
impact on chondrocyte biology and OA-related processes.
Despite these limitations, our study advances understanding of the role
of alternative splicing in OA. We present an sQTL analysis in
chondrocytes, and we built upon this resource by performing
colocalization analysis and functional validation. The identification
of these risk genes and splicing events provides potential targets for
therapeutic intervention and biomarker development in OA. As we
continue to unravel the complex genetic and molecular underpinnings of
OA, integrative approaches that consider multiple layers of gene
regulation will be instrumental in developing more effective diagnostic
tools and targeted therapies for this debilitating disease.
Methods
Ethical statement
Non-OA tissues were obtained from deceased donors with no known history
of arthritis from the Gift of Hope Organ and Tissue Donor Network
through the Rush Medical Center, Chicago, IL. OA tissues were obtained
from patients undergoing knee joint replacement surgery at the
University of North Carolina Hospitals. The University of North
Carolina Institutional Review Board (IRB) reviewed the study protocol
and determined that research involving de-identified tissues from
deceased donors and surgical waste material does not constitute human
subjects research under federal regulations, and therefore did not
require IRB approval (Study Number 14-0189). Therefore, the requirement
for written informed consent for the use of tissues in this study was
waived for both OA and non-OA tissues. The families of deceased tissue
donors do consent that tissue can be used for research purposes. All
samples were de-identified, and no individually identifiable protected
health information was accessible to investigators. Essential
demographic information, including sex and age ranges can be found in
Supplementary Data [182]1.
Dataset and sample information
Primary human chondrocytes were isolated from talar cartilage tissue
obtained from 101 non-osteoarthritis (non-OA) donors (101 PBS resting
and 101 FN-f stimulated) and from 16 donors undergoing knee joint
replacement due to OA. The complete dataset is available in the
Database of Genotypes and Phenotypes (dbGaP, accession:
phs003581.v2.p1).
Chondrocytes were isolated through enzymatic digestion of cartilage
tissues and maintained in a confluent monolayer culture to retain the
chondrocyte phenotype. Cultured cells were treated with
phosphate-buffered saline (PBS) as a resting condition or 1 µM
recombinant human fibronectin fragment (FN7-10; also known as FN-f) for
18 h each, simulating OA-like conditions. OA samples underwent the same
procedures, except for FN-f treatment. DNA extraction was performed
using the QIAamp DNA mini kit (Qiagen, #51304), and genotyping
performed using the Infinium Global Diversity Array-8 v.10 Kit
(Illumina #20031669). Genotype data were processed with comprehensive
quality control using PLINK (v1.90b3.45)^[183]64, including removal of
variants with excessive missing data (> 10%), deviation from
Hardy–Weinberg equilibrium (p < 1 × 10^−6), and low minor allele
frequency (< 1%). Sample relatedness was assessed through
identity-by-descent analysis, retaining only unrelated individuals.
Population structure was evaluated using principal component analysis
with EIGENSTRAT (v8.0.0) in combination with 1000 Genomes reference
data^[184]65,[185]66. Following quality control, data were imputed
against the TOPMed reference panel with Engle2 (v2.4), yielding
approximately 9.7 million high-quality autosomal variants^[186]67. RNA
extraction used the RNeasy kit (Qiagen #74104) with on-column DNase
digestion. RNA integrity was assessed with the Agilent TapeStation
4150. The New York Genome Center performed RNA-seq library preparation
and sequencing.
No sex-based analysis was performed because the study design was not
explicitly set up to examine sex differences. The disproportionate sex
distribution between groups would have limited statistical power for
meaningful comparisons.
RNA-seq processing and quality control
RNA-seq data libraries averaged approximately 101 million paired-end
reads (2 × 100 bp) per sample. Low-quality reads and adapters were
trimmed using TrimGalore! (v0.6.7)^[187]68, followed by FastQC
(v0.11.9)^[188]69 quality control. Trimmed fastqs were aligned to the
GENCODE release 45 (GRCh38.p14) reference genome using the STAR aligner
(v2.7.7a)^[189]70. Gene expression levels were estimated using Salmon
(v1.9)^[190]71 with –seqBias, –gcBias, and –validateMappings flags,
assembly of GENCODE version 45 transcript sequences. The
tximeta^[191]72 R package facilitated gene-level scaled transcript
analysis. RNA signal tracks for PBS resting, FN-f stimulated
chondrocytes, and OA donors were created using deepTools
(v3.5.4)^[192]73 and merged by condition. We evaluated genotype
consistency between RNA-seq and genotyping array data using VerifyBamID
(v1.1.3)^[193]74. Two genotyping sample swaps were corrected. Samples
with FREEMIX and CHIPMIX scores above 0.2 after correction were
omitted. This process was also applied to OA donor samples. The final
dataset included 101 PBS (resting) and 101 FN-f treated chondrocyte
samples from non-OA donors, and 16 untreated OA donor samples.
Splicing quantification of introns usage
We quantified alternative splicing events using LeafCutter
(v0.2.9)^[194]75, which allows for the detection of both annotated and
novel splicing events. Exon-exon junctions were extracted from uniquely
mapped reads using RegTools (v1.0.0)^[195]76 with default settings. To
determine how many of our differential splice junctions contained a
canonical motif, we collated the motif annotations reported by the STAR
aligner. LeafCutter analysis was performed with a minimum threshold of
100 supporting reads per intron and a maximum intron length of 500,000
base pairs, using the GENCODE release 45 (GRCh38.p14) reference genome
and GTF as reference. LeafCutter defines intron clusters representing
alternative splicing choices by grouping overlapping introns.
Differential splicing events were annotated by LeafCutter based on
overlap with the GENCODE v45 transcriptome as annotated, including
cryptic 5′, cryptic 3′, cryptic unanchored, and novel annotated pairs.
We calculated the PSI for each intron based on its relative usage
within its cluster. Intron excision ratios underwent standardization
across individuals and quantile normalization across samples. Global
alternative splicing principal components were calculated using the
prepare_phenotype_table.py script from LeafCutter, using raw counts of
intron excision from the intron clustering step.
Differential splicing analysis
Differential splicing analysis was performed using LeafCutter_ds.R,
considering only introns used in at least 10% of samples and with at
least 10 samples per group having sufficient coverage. We compared
PBS-resting vs. FN-f-stimulated chondrocytes and PBS-resting vs.
osteoarthritic (OA) chondrocytes. PSI values were normalized using a
custom method that accounts for cluster-wise proportions and handles
missing data. The batch correction was determined by correlations
between technical confounders and the top 10 principal components of
global splicing in resting and FN-f stimulated chondrocytes. We
identified that the RNA extraction kit batch and fragment batch were
the primary sources of batch effects in our data. To address these
issues, along with donor-specific variations, we applied the
removeBatchEffect() function from the limma package^[196]18. For the
final differential splicing analysis, we employed a mixed-effects
model:
[MATH:
model.m<
mi>atrix~condi<
/mi>tion+covariates+(1∣DonorID<
mo>) :MATH]
1
Covariates included RNA shipping date, RNA extraction kit batch, and
fragment batch. Significantly spliced intron junctions were identified
using a threshold of adjusted p-value (padj) < 0.05 (Benjamini–Hochberg
procedure for multiple test correction) and absolute PSI difference
(|∆PSI|) > 0.15.
Characterization of differentially spliced genes
Significantly spliced intron junctions (padj <0.05 and |∆PSI| > 0.15)
were identified in comparisons of PBS-resting vs. FN-f stimulated
chondrocytes and PBS-resting vs. OA chondrocytes. We performed pathway
(KEGG, Reactome source) and GO analyses on the associated genes using
the findMotifs.pl script from the HOMER^[197]77 software suite (v4.11).
Enriched GO terms (p < 0.01) and KEGG and Reactome pathways (p < 0.01)
were identified, with GO terms further refined based on semantic
similarity using rrvgo (v1.14.2)^[198]78. Protein-protein interaction
networks were constructed using the STRING database (v12.0)^[199]79,
considering only high-confidence interactions (combined score > 0.9).
Comparative splicing analysis in OA and FN-f chondrocytes
We compared differentially spliced intron junction-associated genes
between PBS-resting versus FN-f stimulated and PBS-resting versus OA
chondrocytes. Both comparisons used the same thresholds (padj < 0.05
and |∆PSI| > 0.15). To assess the relationship between OA-associated
and FN-f-induced splicing changes, we examined the PSI in OA for genes
showing differential splicing in response to FN-f in chondrocytes. The
significance of these comparisons was evaluated using the Wilcoxon
test, with p < 0.01 considered statistically significant.
RNA isolation and RT-qPCR analysis of SNRNP70 splicing
For RT-qPCR validation experiments, a separate set of primary human
chondrocytes was isolated from non-OA donors following the same
protocols described in our dataset information section. Cells were
treated with PBS or FN-f for 18 h under identical conditions to our
RNA-seq experiments. Total RNA was extracted using the RNeasy kit
(Qiagen #74104) with gDNA eliminator columns, followed by the Turbo
DNA-free kit (cat AM1907). RNA was reverse transcribed using the
High-Capacity RNA-to-cDNA kit (cat #4387406) according to the
manufacturer’s instructions. qPCR was performed using Power-up SYBR
Green Master Mix (cat A25742) on a 96-well format Applied Biosystems
QuantStudio™ 6 Flex machine (Applied Biosystems, Foster City, CA, USA).
Custom SYBR green primers specific to the SNRNP70 alternative splicing
junctions were designed. These included primers targeting the junction
between exon 7 and alternative exon 8 (Forward:
5′-CGAGACATGCACTCCACCAC-3′, Reverse: 3′-CTGAGGTGTTGTAAGGGGAGC-5′), the
junction between exon 7 and regular exon 8 (Forward:
5′-CGAGACATGCACTCCGCTTA-3′, Reverse: 3′-CTCCAGCCCTTCACGGTTC-5′), and a
constitutive region spanning exons 3–5 (Forward:
5′-GCGCATGGAGAGGAAAAGAC-3′, Reverse: 3′-CCCCTGAGCATTGGGATCAT-5′).
Expression levels were normalized to the constitutive exon region
(exons 3–5). PCR cycling conditions were 95 °C for 2 min; 40 cycles of
95 °C for 15 s, 60 °C for 60 s. To confirm amplicon specificity and
expected product sizes, PCR products were visualized on 1.5% agarose
gels stained with SYBR Safe DNA Gel stain (Thermo Scientific, cat
[200]S33102), with all products showing a single amplicon at the
predicted size.
Western blot analysis
Primary human chondrocytes were treated with PBS or FN-f for 18 h as
described above. The cells were lysed in Cell Lysis Buffer (Cell
Signaling Technology, cat 98035), supplemented with Halt Protease and
Phosphatase Inhibitor Cocktail (ThermoFisher, cat 78440) and
phenylmethanesulfonyl fluoride solution (Cell Signaling cat 8553S).
Protein concentration was determined using the Micro BCA Protein Assay
Kit (Thermo Fisher, cat 23235), according to the manufacturer’s
instructions. Protein (12 µg per lane) was mixed with 4× Laemmli buffer
containing 2-mercaptoethanol, heated at 95 °C for 5 min, and separated
on polyacrylamide gels prepared with the TGX Stain-Free FastCast
Acrylamide Kit (Bio-Rad, cat 1610180), using 10% ammonium persulfate
and TEMED. The separated proteins were transferred to 0.2 μm PVDF
membranes using Trans-blot Turbo transfer buffer (Bio-Rad, cat
10026938). After blocking in 8% non-fat dry milk in Tris-buffered
saline with 0.1% Tween-20 (TBST) for 90 min at room temperature, the
membranes were incubated overnight (18 h) at 4 °C with SNRNP70
Recombinant Monoclonal Antibody (ThermoFisher, cat MA5-55557, 1:1000
dilution) and GAPDH (Trevigen, cat 2275-PC-100) as the loading control.
After washing with TBST, the membranes were incubated with secondary
antibody (Cell Signaling Technology, cat 7074 and cat 7076, 1:1000
dilution) for 90 min at room temperature. Protein bands were visualized
using Radiance Peroxide and ECL substrate (Bio-Rad, cat 1705062). Band
intensities were quantified using FIJI, and SNRNP70 expression was
normalized to GAPDH. A Student’s t-test was used to calculate the
p-value. Full western blots are provided in Supplementary Fig. [201]8.
sgRNA design and RNP complex formation
Two custom Alt-R crRNAs targeting SNRNP70 (5′-TTGGTTGAGGACAGCCCCCT-3′
and 5′-TATCACCCCACTGCCCGCTG-3′) were designed to target the plus and
minus strands, respectively. These crRNAs, synthesized by Integrated
DNA Technologies (IDT), exclude the PAM sequence. Guide RNA (gRNA)
duplexes were formed by combining 1.5 µL each of Alt-R tracrRNA (IDT,
cat 1073190) and crRNA in PCR tubes. The mixtures were annealed at
95 °C for 5 min and cooled to room temperature. Ribonucleoprotein (RNP)
complexes were assembled immediately before nucleofection by combining
1.2 µL of annealed gRNA with 2.6 µL of Alt-R Cas9 Nuclease V3 (IDT, cat
1081058).
Nucleofection of primary human chondrocytes
Primary human chondrocytes isolated from macroscopically normal
cartilage of cadaveric ankles were cultured in DMEM/F-12 medium
supplemented with 10% FBS and antibiotics. Approximately 3 × 105 cells
per condition underwent nucleofection using the Lonza Nucleofector
system (program ER100). Cells were resuspended in 16 µL nucleofector
solution, combined with 5 µL RNP complex (or PBS for mock control) and
1 µL Cas9 Electroporation enhancer (IDT, cat 1075916), then
electroporated. Post-nucleofection, cells were cultured in an
antibiotic-free medium with 20% FBS at 37 °C and 5% CO[2] for at least
one week.
Genome editing and colony screening
Genomic DNA was isolated using a column-based extraction method,
following the QIAGEN DNeasy Blood & Tissue Kit protocol (QIAGEN, Cat
69504). PCR screening targeted the SNRNP70 alternative exon 8 region
using specific primers (Forward: 5′-CTCCTCCCTCTGTTTCTGATG-3′, Reverse:
5′-CAGGAAAGGGGAGTCGTAGAG-3′). PCR reactions (20 µL total) contained
1–5 µL genomic DNA, 10 µL Platinum II hot start Green master mix
(Invitrogen, cat 14001012), and 0.5 µL each of forward and reverse
primers (10 µM stocks). PCR cycling conditions were 94 °C for 2 min; 40
cycles of 94 °C for 15 s, 57 °C for 30 s, 68 °C for 45 s; final
extension at 68 °C for 2 min. Products were analyzed by 1.5% agarose
gel electrophoresis and stained with SYBR Safe DNA Gel stain (Thermo
Scientific, cat [202]S33102). PCR final product gel electrophoresis is
provided in Supplementary Fig. [203]8.
Single-cell colony selection and expansion
Edited cells were seeded at ~200 cells per 6 cm^2 dish. After 20 days,
individual colonies were manually isolated. Selected colonies were
expanded from 24-well plates and plated as passage 1 cells, then
cultured in media supplemented with 1 ng/mL TGF-β1 (Life Technologies,
cat PHG9214) and 5 ng/mL bFGF (Life Technologies, cat PHG0264) to
expand. All 46 selected single colonies were confirmed to have edits by
PCR screening as described above. Further SNRNP70 alternative exon 8
deletion was confirmed by Sanger sequencing of PCR products from both
wild-type and potential heterozygous knockdown (Het-KD) clones. For
Het-KD clones, the 512 bp band was gel-extracted. PCR products were
purified using Exonuclease I (EXO1) and Shrimp Alkaline Phosphatase
(SAP) treatment, following Eurofins Genomics sample preparation
guidelines. Purified products were sequenced using SimpleSeq tubes
(Eurofins Genomics: [204]https://eurofinsgenomics.com/). To minimize
potential clone-specific effects and ensure that observed changes were
due to heterozygous deletion of SNRNP70 alternative exon 8 rather than
individual clonal variation, multiple independently derived clones for
each genotype were pooled. Five wild-type clones were combined into one
group, and four heterozygous knockout clones into another. Each group
was plated into six wells: two for sequencing, one as a backup, and two
for RNA-seq experiments.
RNA sequencing of Het-KD cells
Cells were expanded to 500,000 cells per well in six-well plates. Total
RNA was extracted from two sequencing replicates of each cell
population using the RNeasy kit (Qiagen #74104) with on-column DNase
digestion. RNA quality and quantity were assessed using the Qubit RNA
High sensitivity assay (Thermo Fisher Scientific #[205]Q32582) and the
Agilent TapeStation 4150 system. Library preparation was performed
using the KAPA RNA HyperPrep Kit with RiboErase (Roche #KK8560)
following the manufacturer’s protocol. Libraries were quantified using
the KAPA Library Quantification Kit (Roche #07960298001). Sequencing
was conducted on a NextSeq 2000 system using NextSeq™ 1000/2000 P1
Reagents (100 cycles, Illumina #20074933), generating an average of
24.5 million paired end reads (2 × 50 bp) per sample.
Edited cell quantification and differential expression
RNA-seq data from edited cells were processed using a pipeline similar
to the one described for the PBS-resting, FN-f-stimulated, and OA
chondrocyte RNA-seq data. Briefly, quality control with FastQC
(v0.11.9)^[206]69 and reads were aligned to the GENCODE release 45
(GRCh38.p14) reference genome using the STAR aligner (v2.7.7a)^[207]70.
Gene expression levels were quantified using Salmon (v1.9)^[208]71 with
–seqBias, –gcBias, and –validateMappings flags, utilizing GENCODE
version 45 transcript sequences. Gene-level scaled transcript analysis
was facilitated by the tximeta^[209]72 R package. Differential gene
expression analysis between wild-type and SNRNP70 alternative exon 8
heterozygous knockout cells was performed using DESeq2^[210]80. We
filtered out lowly expressed transcripts, retaining only those with at
least 10 counts in 10% or more of the samples. Differential expression
was defined by an FDR-adjusted p-value < 0.01 and absolute log[2] fold
change > 2, using the apeglm method for fold change shrinkage.
SNRNP70 Exon-level quantification
Exon-level counts for SNRNP70 were generated using featureCounts
(subread v2.0.6)^[211]81. The resulting counts were analyzed using the
EBSEA (v1.33.0)^[212]82 R package, which performs statistical testing
at the exon level before aggregating results to the gene level. This
approach increases statistical power while accounting for exon
dependence using empirical Brown’s method^[213]83. We considered exons
with p-value < 0.05 as significantly differentially used between
wild-type and heterozygous knockout samples.
Characterization of edited cells
Differentially expressed genes (log[2]FC > 2 and FDR-adjusted
p-value < 0.01) were analyzed for KEGG and Reactome pathway enrichment,
as well as GO terms. Genes were split based on up or downregulation in
heterozygous knockout cells compared to wild-type. Enriched pathways
and GO terms (p < 0.01) were identified using the findMotifs.pl script
from HOMER (v4.11)^[214]77. GO terms were further refined based on
semantic similarity using rrvgo (v1.14.2) with the following
parameters: simThresh = 0.7, simMethod = “Rel”^[215]78.
Comparative analysis with FN-f stimulation
We compared the differential genes between SNRNP70 alternative exon 8
heterozygous knockout cells and controls with those observed in
FN-f-stimulated chondrocytes (from our previous PBS-resting versus
FN-f-stimulated analysis). Log[2] fold changes were examined for genes
showing differential expression in both conditions. The significance of
these comparisons was evaluated using the Wilcoxon test, with p < 0.01
considered statistically significant.
Splicing quantitative trait loci (sQTL) analysis
We performed sQTL analysis to identify genetic variants associated with
alternative splicing events in chondrocytes under resting (PBS) and
osteoarthritis-like (FN-f stimulated) conditions. This analysis was
built upon the differential splicing using the same PSI values derived
from LeafCutter^[216]75. We tested genetic variants within ±100 kb of
each intron cluster for association with PSI values. This window size
was chosen based on previous studies that have demonstrated the
enrichment of sQTL SNPs in proximal regions of splice sites^[217]84.
The ±100 kb range allows for the capture of both proximal and distal
regulatory elements that may influence splicing^[218]85. We applied
stringent filtering criteria for variant inclusion using PLINK
(v1.90b3) to ensure robust associations^[219]64. For autosomal
chromosomes, we retained variants with either at least 5 heterozygous
donors or at least 2 homozygous minor allele donors. To account for
population structure, we calculated genotype principal components (PCs)
using bcftools (v1.10) and PLINK (v1.90b3)^[220]64,[221]86. Our models
included the first 10 principal components (PCs), explaining
approximately 60% of the genetic variance as covariates. We employed
QTLtools (v1.3), which extends FastQTL functionality, for sQTL
mapping^[222]30. We incorporated covariates to control for known
technical variables and hidden factors, optimizing selection through
iterative testing of global splicing principal components (PCs) 1–20.
The final sQTL models were:
For PBS samples:
[MATH: PSI~S
NP+10PCsofgenotype+5PCsofglobalsplicing+Fragm
entBatch+RNAextractionKitBatch+SequencingBatch :MATH]
2
For FN-f samples:
[MATH: PSI~S
NP+10PCsofgenotype+4PCsofglobalsplicing+FragmentBatch+RNAextractionKitBatch+SequencingBatch :MATH]
3
To determine significant associations, we implemented a
permutation-based approach using QTLtools^[223]30. We performed 1000
permutations for each phenotype to characterize the null distribution
of associations empirically. This permutation scheme approximates the
tail of the null distribution using a beta distribution. We used the
qtltools_runFDR_cis.R script to implement the Storey & Tibshirani False
Discovery Rate (FDR) procedure, identifying significant sQTLs at 5%
FDR^[224]87. The resulting nominal p-value threshold range was 9.99e-6
to 0.0001 for both PBS and FN-f conditions. To identify multiple
independent sQTLs per splice intron junction, we performed conditional
analysis using QTLtools^[225]30. This method uses a stepwise linear
regression approach to discover additional independent signals while
accounting for multiple testing. Secondary signals are variants
identified after conditioning on the primary signal (lead variant),
tertiary signals after conditioning on both primary and secondary
signals, and quaternary signals after conditioning on all previous
signals, with each level representing an additional independent genetic
effect. The process involves establishing a nominal p-value threshold
for each phenotype, iteratively identifying the most robust
associations, and assigning nearby variants to independent signals. For
genes with multiple sQTLs identified through conditional analysis, we
further assessed allelic heterogeneity by calculating pairwise LD using
PLINK (v1.90b3)^[226]64 and classified variants as independent if they
had R^2 < 0.7. This allowed us to systematically identify genes with
multiple independent genetic effects on splicing.
sQTLs and eQTLs overlap analysis
To investigate the relationship between genetic regulation of splicing
and gene expression, we analyzed the overlap between our sQTLs and our
previously published eQTLs derived from the same dataset^[227]9. We
identified variants in high LD (R^2 > 0.7) between lead sQTLs and
published eQTLs affecting the same gene. For overlapping variants, we
performed conditional analysis to determine if sQTL effects were
independent of eQTL effects. We constructed two linear models using the
lme4 package in R:
[MATH: Nullhypothesis=PSI~sQTLgenotype+covariates :MATH]
4
[MATH:
Alterna<
mi>tivehypothesis=PSI~sQTLgenotype+eQTLgenotype+covariates :MATH]
5
We used ANOVA to compare these models and applied the FDR correction to
determine significance at a 5% FDR.
Condition-specific and response sQTL analysis
To assess condition-specific effects and potential genetic
interactions, we applied inverse normal transformation of PSI ratios.
This transformation was necessary to normalize the distribution of PSI
values and reduce the impact of outliers, ensuring robust statistical
analysis. We then conducted interaction testing using linear
mixed-effects models, comparing reduced models (without interaction
terms) to full models (including interaction terms) for each
sQTL-splice intron junction pair. To investigate genetic variants
influencing splicing in a condition-specific manner, we constructed two
linear mixed models using the lme4 package^[228]31 in R:
[MATH: Nullhypothesis=PSI~genotype+co
variate<
mi>s+condition+(<
mn>1∣DonorID) :MATH]
6
[MATH:
Alterna<
mi>tivehypothesis=PSI~genotype+co
variate<
mi>s+condition+genotype:
conditio
n+(1∣
DonorID) :MATH]
7
Where covariates included the same technical factors and PCs used in
the standard sQTL mapping, the condition was coded as 0 (PBS) or 1
(FN-f), and (1|DonorID) accounted for donor-specific random effects. We
assessed the significance of the genotype:condition term using ANOVA,
with p < 0.05 considered significant for the interaction effect.
To identify high-confidence condition-specific sQTLs, we applied
additional filters: (1) sQTLs detected only in one condition (PBS or
FN-f), (2) at least 5 donors per minor allele to ensure robust
estimation, (3) Absolute beta difference ≥ 0.2 between conditions,
indicating a substantial change in effect size.
RNA binding protein motif analysis
We selected 140 RBPs from the ENCODE project, focusing on all eCLIP-seq
data generated for K562, HepG2, and SM-9MVZL cell lines^[229]88. The
eCLIP-seq data were processed, aligned to the GRCh38 genome assembly,
and provided in bed format. We used the Functional enrichment of
molecular QTLs (Fenrich) module from QTLtools to evaluate the overlap
between sQTLs and RBP binding sites^[230]30,[231]89. A permutation
scheme was employed to maintain the distribution of functional
annotations and molecular QTLs around molecular phenotypes. We
performed 1000 permutations to estimate the probability of the observed
overlaps of the lead sQTLs relative to expectation. For sQTLs
overlapping with RBP binding sites, we calculated the Pearson
correlation between the RBP gene expression (from the same donor
RNA-seq data) and the spliced intron junction usage pattern across all
samples. We considered correlations with R^2 > 0.15 and p-value < 0.05
as significant. We performed KEGG and Reactome pathway enrichment
analysis on the genes associated with sQTLs that overlap RBP binding
sites and show a significant correlation with RBP expression. Using the
findMotifs.pl script from the HOMER software suite (v4.11), the KEGG
and Reactome pathways (p < 0.01) were identified.
Colocalization between sQTL and OA GWAS
We conducted colocalization analysis between our sQTL results and OA
GWAS signals. Summary statistics for 11 OA phenotypes were obtained
from Boer et al., encompassing 826,690 individuals of European and East
Asian descent (Musculoskeletal Knowledge Portal,
[232]https://msk.hugeamp.org/). These statistics were converted to hg38
coordinates for consistency with our analysis. Overlapping signals were
identified by selecting lead sQTL variants in moderate LD (r^2 > 0.5)
with lead GWAS variants. LD was calculated using PLINK (v1.90b3.45)
with parameters –ld-window 200000 –ld-window-kb 1000, based on the 1000
Genomes European reference panel phase 3, chosen due to the East Asian
population sample size > 4000 individual^[233]4,[234]64,[235]65. We
defined a window encompassing the sQTL and GWAS lead variants for each
pair of overlapping signals, extending by 250 kb in each direction. To
ensure consistent identification, we mapped variants to rsIDs using
dbSNP 155, querying dbSNP VCF files with bcftools and custom scripts to
match variants based on chromosome, position, and alleles.
Colocalization analysis was performed using the Bayesian method coloc
package, implemented in R^[236]90. The coloc.abf function was used to
calculate posterior probabilities for various hypotheses, including H3
(association with sQTL only) and H4 (association with both sQTL and
GWAS traits). Signals with a posterior probability for H4 (PP H4) > 0.7
were considered as evidence of significant colocalization. We also
calculated the ratio (PP H4)/(PP H3 + PP H4) to quantify the strength
of colocalization evidence. This analysis was conducted separately for
sQTLs identified in PBS and FN-f conditions to examine
condition-specific effects on OA risk.
Visualization
Protein-protein interaction networks were visualized using Cytoscape
(v3.9)^[237]91. Heatmaps were generated using the ComplexHeatmap R
package (v3.19)^[238]92. Genomic signal tracks and locus zoom plots
were created with plotgardener (v1.9.4). All other plots were produced
using ggplot2 in R^[239]93. These visualizations were generated using R
version 4.3.
Reporting summary
Further information on research design is available in the [240]Nature
Portfolio Reporting Summary linked to this article.
Supplementary information
[241]Supplementary Information^ (2.4MB, pdf)
[242]41467_2025_63299_MOESM2_ESM.docx^ (20.5KB, docx)
Description of Additional Supplementary Files
[243]Supplementary Data 1^ (87.4KB, pdf)
[244]Supplementary Data 2^ (183.1KB, xlsx)
[245]Supplementary Data 3^ (52.4KB, xlsx)
[246]Supplementary Data 4^ (83.6KB, xlsx)
[247]Supplementary Data 5^ (2.3MB, xlsx)
[248]Supplementary Data 6^ (60.2KB, xlsx)
[249]Supplementary Data 7^ (19.6KB, xlsx)
[250]Supplementary Data 8^ (13.4KB, xlsx)
[251]Reporting Summary^ (1.6MB, pdf)
[252]Transparent Peer Review file^ (898.3KB, pdf)
Acknowledgements