Genetic alternative splicing regulation mapping of cartilage and synovium reveals tissue-specific mechanisms of joint-related traits

Domain: biology
Published: 2026-03-12
Editorial outcome: minor_revision

Abstract
Cartilage and synovium are essential tissues involved in joint-related diseases and traits, including osteoarthritis (OA), rheumatoid arthritis (RA), and human height. Although genome-wide association studies (GWAS) identify numerous risk loci for these traits, the molecular mechanisms underlying these associations, particularly those involving alternative splicing, remain poorly understood due to the lack of splicing-related genetic data in relevant tissues. To address this limitation, we generate a splicing quantitative trait loci (sQTL) resource for cartilage and synovium. We identify 2,796 independent cis-sQTLs and six trans-sGenes across the two tissues, including 179 tissue-specific cis-sQTLs. Fine-mapping analysis identifies 116 high-confidence functional sVariants predicted to affect splicing through splice site gain or loss, with approximately half located outside canonical splice site motifs. Integration of sQTL data with GWAS summary statistics reveals 12 osteoarthritis, 6 rheumatoid arthritis, and 183 height effector genes in joint tissues. Notably, seven of the 12 osteoarthritis effector genes show joint tissue-specific colocalization. Together with the finding that tissue-specific sGenes play crucial roles in tissue-related biological processes and diseases, our work highlights the significant impact of tissue-specific alternative splicing regulation on disease etiology and provides a valuable resource for further mechanistic validation.

Introduction
Cartilage and synovium are two essential joint tissues associated with multiple joint-related diseases and traits, including osteoarthritis (OA), rheumatoid arthritis (RA), and human height. OA is the most prevalent chronic joint disease, affecting over 500 million people worldwide and causing pain, disability, and loss of function1. Articular cartilage, the primary site of degeneration during OA progression, undergoes broad molecular disturbances that significantly drive OA pathogenesis2-4. Meanwhile, the synovium plays a crucial role in OA progression5,6, particularly through the production of inflammatory mediators that exacerbate cartilage damage7. RA, a systemic chronic inflammatory disease, primarily targets the synovium, which serves as the major site of inflammation and plays a critical role in disease progression8,9. Moreover, cartilage function is a main determinant of human height10,11. Thus, elucidating the molecular mechanisms governing cartilage and synovium function is essential for advancing our understanding of joint-related diseases and traits.

Genome-wide association studies (GWAS) have identified numerous genetic loci associated with joint-related diseases and traits12-14. However, uncovering the functional mechanisms of GWAS hits in relevant tissues remains challenging, as approximately 90% of trait-associated variants lie in non-coding regions15. Using expression quantitative trait loci (eQTL) data from joint tissues as a powerful tool, we and others have identified effector genes for dozens of GWAS loci associated with OA and RA16-19. Nevertheless, many GWAS loci remain functionally unresolved. A plausible hypothesis is that some of these unresolved loci might act by regulating other molecular processes, such as alternative splicing, which has been implicated in various pathological conditions and complex traits20-23. Therefore, the availability of splicing quantitative trait loci (sQTL) resources for joint tissues would be invaluable for resolving the remaining GWAS loci associated with joint-related diseases and phenotypes, offering further insight into the regulatory roles of genetic variants. However, no sQTL data from cartilage and synovium have been available until now.

In this work, we perform sQTL analyses in both cartilage and synovium tissues using samples from over 200 donors. We identify 2,340 genes with significant sQTLs. By complementary validation between deep learning-based predictions and sQTL effects, we pinpoint more than 100 high-confidence functional variants that influence splicing by inducing splice site gain or loss, including many located outside canonical splice site motifs. Additionally, we identify 179 tissue-specific sQTLs, and functional analyses suggest that tissue-specific sGenes play crucial roles in tissue-relevant biological processes and pathological mechanisms. Furthermore, we observe strong enrichment of joint-related trait heritability in the sQTLs from cartilage and synovium. Through colocalization analyses, we identify additional effector genes for OA, RA, and human height. These findings highlight the importance of integrating sQTL signals to illuminate the roles of genetically regulated alternative splicing in disease etiology and complex traits.

Results
To investigate the genetic regulation of alternative splicing in cartilage and synovium, we utilized our existing datasets from 238 donors who underwent total knee replacement surgery. These included 204 cartilage and 202 synovium RNA-seq datasets, as well as matched blood genotype data (4,293,131 common variants with MAF > 0.05)18,24 (Figure 1 and Supplementary Figure 1). We performed cis-sQTL mapping by testing associations between intron usage ratios (measured as PSI (percent spliced in) and calculated using LeafCutter25) and the genotypes of all variants located within ± 1 Mb of each intron cluster, followed by conditional analyses to identify independent signals. At a 5% false discovery rate (FDR), we identified 997 and 1,198 independent cis-sQTLs associated with alternative splicing events of 960 and 1,680 intron clusters in cartilage and synovium, respectively (Figure 2a, Supplementary Figure 2). Power analysis indicated that with our current sample size, the statistical power was already close to 1.0 for variants with MAF > 0.15, and approximately 0.8 for variants with MAF around 0.10 (Supplementary Figure 3).

We next annotated the introns with significant sQTLs using two annotation resources to assess whether they represented novel splicing events or had been previously annotated. First, we evaluated these introns against the GENCODE v19 annotation. Based on this comparison, 83% of significant introns in cartilage and 81% in synovium were fully annotated, with both splice junction boundaries matching GENCODE v19. An additional 11% of introns in cartilage and 13% in synovium had only one boundary annotated, whereas the remaining 6% of introns in both tissues had neither boundary annotated in GENCODE v19 (Supplementary Figure 4a). Second, we annotated the same set of introns using Snaptron26, a large-scale splice-junction database constructed from tens of thousands of public human RNA-seq samples. Using Snaptron, we found that only 11 introns (0.66%) in cartilage and 20 introns (0.64%) in synovium lacked supporting evidence in the database (Supplementary Figure 4b).

We aligned intron clusters with significant cis-sQTLs to 854 and 1,486 genes in cartilage and synovium, respectively (referred to hereafter as sGenes) (Supplementary Data 1,2). Over 97% of these genes were protein-coding (Supplementary Figure 5), and approximately 90% of variants were associated with a single sGene (Supplementary Figure 6). A total of 267 sGenes harbored more than one intron cluster with significant cis-sQTLs (Figure 2b and Supplementary Data 3,4). To assess the independence of regulatory signals within the same gene, we calculated the linkage disequilibrium (LD) R2 between the lead sVariants (the significant sVariants with the smallest P value in each sQTL signal) of different intron clusters. This analysis revealed a broad “U-shaped” distribution. Specifically, for 38.0% and 36.8% of these genes in cartilage and synovium, respectively, the lead-sVariant pairs showed LD R2 > 0.5, suggesting a potential coregulation model of different intron clusters for these genes (Supplementary Figure 7). In contrast, more than 60% of all lead-sVariant pairs showed LD R2 < 0.5 (Figure 2c), indicating a high proportion of independent genetic effects on alternative splicing of distinct intron clusters within the same gene. As an illustrative example of independent regulation, in synovium, rs2269219 and rs536765 were found to regulate the alternative splicing of intron cluster 1 and intron cluster 2 of the FECH gene, respectively (Figure 2d). The dosage of the C allele at rs536765 was associated with reduced usage of intron junction 3 (chr18:55222184-55230106) in intron cluster 2 (Figure 2d), leading to retention of exon 7 in the final transcript. Meanwhile, the A allele at rs2269219 was correlated with increased usage of intron junction 1 (chr18:55238835-55240478) in intron cluster 1, thereby promoting inclusion of the longer exon 4 in the final transcript (Figure 2d). These two variants were largely independent (LD R2 = 0.06). Consistent with these effects, individuals homozygous for C/C at rs536765 and A/A at rs2269219 showed the highest expression proportion of the isoform ENST00000585494, which contains both exon 7 and the longer exon 4 (Figure 2e). This finding illustrates how independent variants can act in combination on separate intron clusters to influence final transcript isoforms in a combinatorial manner. Additionally, we examined trans-effects on alternative splicing and identified six significant trans-sGenes across the two tissues (Supplementary Figure 8 and Supplementary Data 5). A larger sample size will be necessary to increase statistical power and identify additional trans-sGenes.

We next compared cis-sQTLs with cis-eQTLs from our previous analyses18,19. In cartilage and synovium, 60.3% and 53.2% of sGenes did not overlap with eGenes identified from cis-eQTL mapping, respectively (Figure 2f). Even among genes identified as both sGenes and eGenes, about 50% of sQTL-eQTL lead-sVariant pairs had LD R2 < 0.5 (Figure 2g), indicating largely independent genetic mechanisms underlying alternative splicing and gene-level expression. Colocalization analysis further confirmed that approximately 62% of sGenes in cartilage and 61% in synovium were either uniquely regulated by sQTLs or influenced by sQTL signals distinct from eQTLs (Supplementary Figure 9). Only about 25% of sGenes shared the same causal variant for both alternative splicing and gene expression. For example, rs149871778 was the lead variant for both the sQTL and eQTL affecting M6PR in both tissues, with a high posterior probability of sharing the same causal variant (PP4 ≈ 1) (Supplementary Figure 10). In contrast, H6PD exhibited both sQTL and eQTL signals, but these were driven by distinct causal variants in cartilage and synovium (LD R2 between lead variants < 0.1; Figure 2h). CPNE3, meanwhile, was regulated exclusively at the splicing level, showing no significant eQTL signals (Figure 2i). These results reinforce that splicing and gene-level expression can be regulated independently, underscoring the importance of investigating splicing QTLs to gain a more complete understanding of regulatory variation in joint tissues.

We next examined the genomic distance between lead sVariants and their target splice site. Approximately 75% of lead sVariants were located within 50 kb of the associated splice site, and more than 90% were within 200 kb (Supplementary Figure 11). Compared with eQTLs, we found lead sVariants were significantly closer to their associated splice sites than lead eVariants were to transcription start sites (TSS) (P < 2 × 10-16, Wilcoxon rank-sum test; Figure 2j), suggesting that genetic regulation of alternative splicing is, on average, more proximal than that of gene expression.

Epigenomic annotation and enrichment analyses further revealed that lead sVariants, similar to lead eVariants, were significantly enriched in open chromatin regions marked with ATAC-seq peaks and active histone modifications, including H3K4me3, H3K27ac, H3K4me1 and H3K36me3, while being depleted in repressive chromatin marks (H3K27me3 and H3K9me3) (Supplementary Figure 12). These results indicate that although alternative splicing and gene expression often involve distinct causal variants and regulatory architectures, they tend to localize to shared classes of active regulatory elements, reflecting partially overlapping chromatin contexts.

We conducted enrichment analysis of sVariants to characterize their genomic distribution based on SnpEff annotation27. Unlike eQTLs, sQTLs showed the strongest enrichment in splice donor and acceptor regions (Figure 3a), indicating that the perturbation of splice site sequences is a key regulatory mechanism for alternative splicing. To systematically identify functional sVariants and their effects on splice site gain or loss, we employed SpliceAI28, a deep learning- based model that predicts splice-altering consequences of genetic variants using a sequence-only approach and reports confidence scores (Δ scores) ranging from 0 to 1. SpliceAI predictions with Δ score ≥ 0.2 were compared with our sQTL results to validate both the genomic location and regulatory direction of the predicted splicing effects (Figure 3b and Supplementary Figure 13). Using this approach, we identified 116 functional sVariants that regulate alternative splicing by influencing splice site usage in cartilage and synovium (Figure 3b and Supplementary Data 6).

These functional sVariants were predominantly located near their target splice sites, implying direct disruption or creation of canonical splice site motifs (splice donor: 5’-G|GT-3’, splice acceptor: 5’-AG|-3’) (Figure 3c). For example, SpliceAI predicted that the G-to-A transition at rs7724759 disrupts the adjacent splice donor site (Δ score = 0.32). Consistent with this, our sQTL results from cartilage and synovium samples showed that this transition led to skipping of the exon, with its 3’ end defined by rs7724759 (Figure 3d and Supplementary Figure 14a). Among these nearby functional sVariants, we identified 15 and 22 sVariants that directly perturb the canonical GT sequence at splice donor sites or AG sequence at splice acceptor sites in cartilage and synovium, respectively (Supplementary Data 6). As expected, these sVariants exhibited clear genotype-dependent effects on splice-site usage, in which heterozygous individuals typically exhibited a substantial reduction in splicing, whereas individuals homozygous for the motif-disrupting allele exhibited almost complete loss of splicing at the affected site (Supplementary Figure 15).

Notably, 57 (49.1%) of the identified functional sVariants were located four base pairs (bp) away from their target splice sites, beyond the location of the canonical splice site motif, indicating the involvement of additional regulatory mechanisms (Figure 3c). One such example is rs4072037. SpliceAI predicted that it affects MUC1 splicing with high confidence (Δ score = 0.68), where the G-to-A transition leads to the loss of the splice acceptor site 34 bp upstream and the gain of another site 7 bp upstream, without disrupting the canonical splice site motif. This prediction was supported by our sQTL data from synovium and cartilage samples (Figure 3e and Supplementary Figure 14b).

To investigate the tissue specificity of sQTLs, we compared our results with data from the Genotype-Tissue Expression (GTEx)29 project using mash30 analysis. Consistent with previous studies20,29, we observed low tissue specificity by sign and significantly higher tissue specificity by effect size (magnitude) (Supplementary Figure 16). To provide a quantitative assessment of tissue specificity, we calculated a tissue specificity score as described previously31. The tissue-specificity scores for cartilage and synovium sQTLs were 0.27 and 0.26, respectively.

Focusing on our two joint tissues, we identified 51 cartilage-specific and 128 synovium-specific sQTL-sGene pairs (Figure 4a and Supplementary Data 7,8), respectively. Gene ontology (GO) enrichment analysis revealed that synovium-specific sGenes were significantly overrepresented in immune-related pathways (Figure 4b), consistent with the inflammatory processes known to occur in synovium during RA and OA32-35. Disease Ontology (DO) analysis further highlighted strong enrichment of these genes in immune disorders (Figure 4c). For example, IRF5, a sGene only detected in synovium encoding an interferon regulatory factor, plays a central role in synovial inflammation in both RA and OA33,36 (Figure 4d). For cartilage, although we did not detect significant pathway enrichment, we identified several cartilage-specific sGenes, such as ACAN37 and FOXO138, that are known to be involved in chondrocyte homeostasis and matrix formation (Figure 4e-f). These findings suggest that tissue-specific alternative splicing regulation plays a crucial role in tissue biology and may represent a potential mechanism underlying tissue-related diseases.

Additionally, we observed that synovium-specific sGenes tended to have higher pLI (Loss-of-function Intolerance probability) scores compared to shared sGenes (Figure 4g). A similar but statistically non-significant trend was observed for cartilage-specific sGenes, suggesting that tissue-specific sGenes may be subject to stronger purifying selection. Moreover, tissue-specific sQTL lead variants were located farther from their target splice sites than shared sQTLs (Figure 4h), potentially reflecting the involvement of tissue-specific distal regulatory elements in alternative splicing regulation, such as enhancers39,40. For instance, rs34816914 significantly regulated PAMR1 splicing in cartilage but not in synovium, and was located in a cartilage-specific enhancer marked by H3K27ac and ATAC-seq signals (Figure 4i).

OA is the most common chronic joint disorder. While GWAS have identified over 100 OA risk loci, the contribution of these loci to alternative splicing has remained unclear. By integrating our sQTL data with the OA GWAS dataset12 (Supplementary Data 9), we found significant heritability enrichment41 for multiple OA phenotypes in sQTLs from both cartilage and synovium (Figure 5a). Notably, the enrichment was stronger in sQTLs than in eQTLs, suggesting that alternative splicing plays a distinct and important role in OA pathogenesis. As negative controls, we evaluated heritability enrichment of several mental disorders in the same sQTL datasets, and as expected, observed no significant enrichment (Figure 5a).

To identify specific effector sGenes of GWAS loci, we performed colocalization analysis between our sQTL data and OA GWAS summary statistics. We identified a total of 12 effector genes across cartilage and synovium with high posterior probability of colocalization (PP4 ≥ 0.8) (Figure 5b, Supplementary Figures 17-28, and Supplementary Data 10). Notably, by further integrating sQTL data from non-joint tissues in the GTEx project, our sQTL data aided in identifying seven additional effector sGenes linked to OA risk loci, which could not be identified by sQTL data from any GTEx tissue. In addition, we found that nine effector sGenes were identified either in a cartilage-specific or synovium-specific manner, indicating potential tissue-dependent effects of the corresponding GWAS loci. Finally, for four effector sGenes, their eQTL signals showed no evidence of colocalization with the same GWAS loci in any of our joint-tissues or GTEx tissues (Figure 5c), providing strong support that these risk loci influence OA primarily through regulation of alternative splicing rather than gene-level expression.

COL2A1, which encodes type II collagen, plays a crucial role in cartilage biology42. Colocalization analysis nominated alternative splicing of COL2A1 as an effector event for OA in a cartilage-specific manner. COL2A1 has two annotated transcript isoforms, one containing exon 2 and the other lacking this exon. Previous studies have shown that exon 2-containing transcripts are predominantly expressed in chondroprogenitor cells and early chondroblasts, whereas transcripts lacking exon 2 are the major isoform in differentiated chondrocytes43,44 Our analysis showed that the G allele at rs1793955 was associated with both an increased risk of OA and the retention of exon 2 in the final transcript of COL2A1, but not with the gene-level expression of COL2A1 (Figure 5d-f). This suggests that the higher proportion of the isoform containing exon 2 contributes to OA progression. Further, using paired degenerated and intact cartilage samples from 31 donors (Supplementary Figure 1b), we found a significantly higher usage ratio of exon 2 in degenerated cartilage (Figure 5g). These findings collectively suggest that regulation of exon 2 retention in the final transcript of COL2A1 plays an important role in OA.

Cartilage is an important tissue involved in determination of human height10. By integrating GWAS summary statistics of human height13 with our cartilage sQTL data, we observed significant enrichment of height heritability in cartilage sQTLs (Figure 6a). Colocalization analysis between GWAS and sQTL signals identified 183 effector genes (Supplementary Data 11), of which 79.8% (146/183) could not be nominated using eQTL data (PP4 < 0.5, Supplementary Data 11). These effector genes were significantly enriched in biological processes related to cartilage extracellular matrix organization (Figure 6b), which is essential for cartilage homeostasis45,46 and exerts effects on human height47. For example, we identified PDGFRA (Figure 6c), for which rs6554162 was both a variant associated with human height and the lead variant regulating alternative splicing of PDGFRA without affecting gene-level expression. Previous studies have demonstrated that PDGFRA promotes chondrogenesis48,49. These findings support a role for cartilage alternative splicing in regulating human height and highlight the potential of sQTL analysis to reveal mechanisms underlying complex traits.

In parallel, we observed a significant enrichment of RA heritability in synovium sQTLs. Colocalization analysis identified six effector genes associated with RA risk (Figure 6d, Supplementary Figures 29-34, and Supplementary Data 12). Notably, all of these genes were uniquely nominated using synovium sQTL data rather than synovium eQTL data and four of them could not be nominated using sQTL data from any GTEx tissues (Figure 6d). TRAF1 (tumour necrosis factor receptor-associated factor 1), which encodes a negative regulator of Toll-like receptor signaling in monocytes, has reduced expression associated with increased inflammation and RA risk50. Our colocalization results suggest that alternative splicing of TRAF1 in synovium may contribute to RA susceptibility (Figure 6e). Additionally, FCGR2B, which encodes Fc receptor IIB for IgG, is involved in immune response and has been previously linked to several autoimmune diseases, especially systemic lupus erythematosus51,52. Our analysis indicates that alternative splicing of FCGR2B may contribute to RA susceptibility (Figure 6f).

Taken together, these findings provide a valuable set of candidate genes for future mechanistic exploration of splicing-mediated regulation in OA, RA, and human height, reinforcing the crucial role of alternative splicing in complex diseases and traits.

Methods
Ethics
This study was conducted in accordance with all relevant regulations regarding the use of human participants and follows the criteria set by the Declaration of Helsinki. The study protocol was approved by the Ethics Committee of Xi’an Jiaotong University Honghui Hospital (approval number: 202304005). All patients provided written informed consent before participating in the study.

Dataset and sample information
RNA-seq data of cartilage and synovium, together with genotype data from matched blood samples, were obtained from our previous studies. The raw data are available from Genome Sequence Archive (accession: HRA006224, HRA004624). The dataset comprises a total of 238 donors (163 females and 75 males) (Supplementary Figure 1a,d) who were clinically diagnosed with osteoarthritis (OA) and underwent total knee replacement surgery. All donors had matched blood-derived genotype data. Among these, 204 donors had cartilage RNA-seq data. For 173 of the 204 donors, only degenerated cartilage samples were available, whereas for 31 donors, both degenerated and intact cartilage samples were obtained. To ensure tissue-state consistency across all donors, RNA-seq data from degenerated cartilage samples were used for sQTL mapping. The 31 paired degenerated and intact cartilage samples were used for within-donor differential alternative splicing analysis. In addition, 202 donors had synovium RNA-seq data. Donors’ ages ranged from 46 to 84 years, with a median age of 66 years and a standard deviation of 7.44 years (Supplementary Figure 1c). Detailed quality control and data preprocessing procedures were described in our previous publications18,24.

Alternative splicing and transcript-level expression quantification
Clean RNA-seq reads were aligned to the human reference genome hg19 using STAR62 (v2.7.9) based on the GENCODE63 v19 annotation, with parameters set according to the GTEx v8 pipeline. WASP filtering64 was performed during read alignment by integrating genotype data from the same donor of RNA-seq data using the --varVCFfile parameter of STAR to correct for reference mapping bias. Only uniquely mapped reads that passed the WASP filtering were used for downstream analyses. Intron clustering and intron usage ratio quantification were performed using the leafcutter_cluster_regtools.py script from LeafCutter25 (v0.2.9), with intron junctions extracted from BAM files by RegTools65 (v0.5.2; parameters: -s 0 -a 8 -m 50 -M 500000) as input. The leafcutter_cluster_regtools.py script was run with parameters --minclureads 30 --maxintronlen 500000. Normalized intron usage levels were obtained using the prepare_phenotype_table.py script from LeafCutter (v0.2.9) with default settings and used as phenotypes for downstream sQTL mapping. Introns with zero read counts in more than 60% of samples or with insufficient variability (standard deviation < 0.005) were removed. Finally, introns were annotated using GENCODE v19 gene annotations. Intron clusters were retained if at least one intron within the cluster had a splice boundary matching a GENCODE v19-annotated intron, in which case all introns within the cluster were kept. In addition, only intron clusters that could be uniquely assigned to a single protein-coding gene or lincRNA were retained. Transcript-level expression was calculated using RSEM66 (v1.3.1) with BAM files aligned to GENCODE v19 transcriptome as input. The --fragment-length-max parameter was set to 1000, and all other parameters were left at their default values.

cis-sQTL mapping
cis-sQTL mapping was performed using tensorQTL67 (v1.0.10) to test the association between intron usage ratio of each intron that was calculated by LeafCutter and genotypes of variants located within a cis-window of ± 1 Mb upstream and downstream of the corresponding intron cluster. Age, sex, first two genotype principal components (PCs), and PEER factors68 were used as covariates. A range of 1 to 20 PEER factors was tested to maximize the number of significant intron clusters (Supplementary Figure 2). As a result, 11 PEER factors were used for cartilage and 5 for synovium. A grouped permutation procedure (by setting --phenotype_groups parameter for cis mode) was applied for each intron cluster to fit a beta distribution of the strongest association (the smallest P value of all variant-intron pairs tested), and a beta-approximated empirical P values for each intron cluster were obtained. Empirical P values of all intron clusters were subjected to multiple testing correction using the qvalue package69. The intron cluster-level significance threshold was set at FDR < 0.05. Variant-level significance threshold for each significant intron cluster was calculated using the fitted beta-distribution as described in previous studies20,29. Genes with at least one significant intron cluster were defined as sGenes.

To identify independent cis-sQTL signals for each intron cluster, a two-step conditional analysis was performed using cis_independent mode of tensorQTL67, which involved forward stepwise regression followed by backward selection. The lead variant for each independent sQTL signal was reported.

trans-sQTL mapping
Variants with 75-mer mappability scores below 0.9 were first removed. Trans-sQTL mapping was then performed using QTLtools65 (v1.3.1) following the two-step approximate permutation strategy implemented in QTLtools (https://qtltools.github.io/qtltools/). Briefly, a genome-wide null beta distribution of the strongest association P values was constructed based on 50,000 permutations. Empirical P values for each variant-intron pair were subsequently obtained by adjusting the nominal P values using the fitted beta distribution. Variants located at least 5 Mb away from the corresponding intron or on different chromosomes were considered for trans-sQTL analysis. Covariates were identical to those used in cis-sQTL mapping. To further reduce false-positive associations, cross-mapping filtering70 was applied as previously described29. Specifically, genes with mappability scores < 0.8 were excluded, and variants within ±1 Mb of genes that cross-mapped with candidate trans-sGenes were removed. Mappability scores were calculated using k-mer lengths of 75 for exonic regions and 36 for untranslated regions (UTRs). Gene-level P values were then computed to account for multiple introns tested per gene using the formula 1-(1-p))k, where p denotes the empirical P value of the strongest variant-intron association for a given gene and k represents the number of introns tested. Finally, false discovery rate (FDR) correction was applied to gene-level P values, and trans-sQTLs with FDR < 0.05 were considered.

Fine-mapping of sQTL signals
Fine-mapping was performed for each significant sQTL signal using the susieR package (v0.14.2) based on sQTL summary statistics. For each sQTL signal, linkage disequilibrium (LD) matrices were calculated using genotype data from the corresponding samples and provided as input to the susie_rss function. We set the maximum number of causal variants to L = 10, fixed the phenotype variance (var_y) to 1, and enabled residual variance estimation (estimate_residual_variance=TRUE), following recommended practice for summary-statistic-based fine-mapping. For each signal, 95% credible sets were derived when available and used for downstream analyses.

Epigenomic data processing
ATAC-seq data from cartilage and synovium tissues were processed using a unified analysis pipeline. For cartilage, ATAC-seq data from 16 cartilage samples derived from osteoarthritis donors were obtained from GEO (GSE108301)71. For synovium, ATAC-seq data from three synovium samples derived from osteoarthritis donors were obtained from our previous study18. Raw ATAC-seq reads were aligned to the human reference genome (hg19) using Bowtie2 (v2.3.5)72. Reads with a mapping quality score below 30 or marked as PCR duplicates were removed. Peak calling was performed using MACS2 (v2.2.7.1)73,74 following previously described procedures75. For each tissue, BAM files from all samples were merged and used to re-call peaks to define a unified peak set. Individual sample peak sets were then compared with the unified peak set, requiring a minimum overlap of 1 bp. Peaks detected in only a single sample were excluded from downstream analyses.

ChIP-seq peak data for six histone modifications (H3K27ac, H3K27me3, H3K36me3, H3K4me1, H3K4me3, and H3K9me3) were collected from publicly available datasets. For cartilage, processed ChIP-seq peak files generated from human mesenchymal stem cell–derived chondrocytes were obtained from the ENCODE project (GEO accession numbers: GSE188124, GSE187652, GSE187044, GSE187055, GSE187641, and GSE187471). For synovium, processed ChIP-seq peak files generated from synovium samples of 11 osteoarthritis donors were obtained from GEO (GSE112655)76. All ChIP-seq datasets were converted to the hg19 genome build using liftOver when necessary. For synovium ChIP-seq data, peaks were filtered to retain only those detected in more than one sample for each histone modification. The resulting peak sets were used for downstream analyses.

sVariants enrichment
All lead variants of independent sQTL signals were LD-pruned with a threshold of LD R2 < 0.8 to generate a query set. For each query variant, 50 random variants were extracted from the SNPsnap dataset77 as a background set, with matching minor allele frequency (MAF), distance to gene, gene density and LD buddy. We annotated both the query and background variant sets using epigenetic data and SnpEff27 (v5.2c), and performed enrichment analysis using the Fisher’s exact test.

SpliceAI prediction and cross-validation with sQTL effects
To investigate whether and which of the significant sVariants exert alternative splicing regulatory effects by mediating splice site gain and loss, we applied SpliceAI28 (v1.3.1) with default settings. SpliceAI is a deep residual neural network that predicts the impact of genetic variants on splice site gain and loss within genes by evaluating sequence context within ±50 bp of the variant position by default. The human hg19 genome sequence, GENCODE v19 gene annotation, and a VCF file containing all significant sVariants were used as inputs. SpliceAI reports predicted splice site gains and/or losses together with corresponding Δ scores for each variant. sVariants with Δ scores ≥ 0.2 that showed concordant effects supported by sQTL data were considered high-confidence functional sVariants.

Tissue-sharing and specificity of sQTLs
Multivariate adaptive shrinkage (mash)30 was applied using the mashr (v0.2.57) package to evaluate the tissue-sharing and specificity of sQTLs. We selected the variant-intron pairs detected in at least 40 tissues (including both our cartilage and synovium tissues, as well as at least 38 other GTEx tissues) and with the smallest P value in each gene as the strong set. The random set was generated by merging the 20 randomly selected variant-intron pairs from each gene. Missing beta value was set to 0, and the missing SE was set to 10. A mash model was then constructed. To estimate the tissue-sharing and specificity of sQTLs between our joint tissues and GTEx tissues, the fitted model was applied to the strong set to get the posterior means and LFSR. To identify tissue-specific sQTLs between cartilage and synovium, we applied the fitted model to all significant independent variant-intron pairs in cartilage and synovium, with the missing beta set to 0, and the missing SE set to 10. Variant-intron pairs with LFSR < 0.05 were considered significant. sQTL-sharing by sign means significant variant-intron pairs have the same direction of mashr-estimated effect size in two tissues. sQTL-sharing by magnitude means mash-estimated effect size is within a factor of 0.5 for significant variant-intron pairs in two tissues.

Joint-related traits associations
We obtained GWAS summary statistics of OA12, RA14 and human height13 from the studies with the largest sample size. The data details were provided in Supplementary Data 9. Additionally, we obtained GWAS summary statistics of multiple mental disorders to serve as the negative controls, including Attention deficit hyperactivity disorder (ADHD)78, Anorexia nervosa (AN)79, Autism spectrum disorder (ASD)80, Major depressive disorder (MDD)81, Obsessive-compulsive disorder (OCD)82, and Tourette syndrome (TS)83. Heritability enrichment of these traits in cartilage and synovium sQTLs was estimated using S-LDSC41 (v1.0.1) according to guidelines (https://github.com/bulik/ldsc/wiki/Partitioned-Heritability). To test the colocalization between GWAS and sQTL signals, we clumped the genome-wide significant variants (P < 5 × 10-8) using PLINK84 (v1.9) with parameters --clump-r2 0.1 --clump-kb 1000 to identify the index SNPs for each GWAS summary dataset. Beta values for all variants located within a ±1 Mb window of the index variant were used for colocalization analysis. For introns within intron clusters with beta-approximated P values less than 0.05 in cartilage or synovium, the beta values for all cis-variants (± 1 Mb) of each intron were extracted for colocalization analysis. For genes showing evidence of colocalization between their sQTL signal and GWAS, the corresponding eQTL signal (beta values for all cis-variants (± 1 Mb) around transcription start site) of the gene were also used for colocalization analysis. Further, to assess whether the colocalized GWAS signals exhibit regulatory effects on gene expression or alternative splicing in other tissues, we integrated sQTL and eQTL signals from 49 GTEx tissues for the corresponding introns and genes identified in cartilage and synovium, and performed additional colocalization analyses.

Colocalization analysis was performed using R package coloc85 (v5.1.0.1) with default settings, which implements a Bayesian test to assess whether two association signals share a consistent causal variant. The analysis reports posterior probabilities for the following five hypotheses: H0 (no association with either trait), H1 (association with trait 1, not with trait 2), H2 (association with trait 2, not with trait 1), H3 (association with both traits via two independent SNPs), H4 (association with both traits via one shared SNP). The results were denoted as PP0, PP1, PP2, PP3, PP4, with higher PP4 indicating greater probability that a single variant affects both traits. After preparation of the input files, we performed colocalization analysis using coloc.abf function with default settings. Genes of which the sQTL signal had PP4 ≥ 0.8 with GWAS signal were considered to be candidate effector genes of GWAS loci.

Differential splicing analysis
For 31 of the 204 donors, RNA-seq data were available for both degenerated and intact cartilage samples, which were used for within-donor differential splicing analysis. For each intron whose sQTL signal colocalized with an OA GWAS signal, differences in percent spliced in (PSI) values between degenerated and intact cartilage were assessed using a paired t-test. PSI values were calculated using LeafCutter (v0.2.9). Introns with a P value < 0.05 were considered differentially spliced, and no fold-change threshold was applied.

GO and DO enrichment analyses
GO and DO enrichment analyses were performed using the R package clusterProfiler86 (v4.16.0) with default settings. GO and DO terms with P value < 0.05 were considered statistically significant. For the tissue-specific sQTL gene enrichment, we used all genes included in the tissue-specific sQTL analysis as background set. For the height GWAS-colocalized sQTL gene enrichment, we used all genes that were eligible for the colocalization analysis as background set.

Data visualization
ChIP-seq and ATAC-seq signals were visualized using Gviz87 (v1.32.0). The Circos plot was generated using circlize88 (v0.4.16). The sashimi plot was generated using trackplot89 (v0.5.7). The gene model was constructed with plotgardener90 (v1.10.2). Other results were visualized using the R package ggplot2 (ref. 91) (v4.2.2).

Data availability
The raw RNA-seq data generated in this study have been deposited in the Genome Sequence Archive for Human (GSA-Human) at the National Genomics Data Center, China National Center for Bioinformation, under accession codes HRA006224 [https://ngdc.cncb.ac.cn/gsa-human/browse/HRA006224] and HRA004624 [https://ngdc.cncb.ac.cn/gsa-human/browse/HRA004624]. These data are available under restricted access due to ethical and legal restrictions regarding identifiable genetic information. Access can be obtained by submitting a data access request to the GSA-Human data access committee, which will review requests in accordance with ethical regulations and participant consent. The expected timeframe for review is approximately 4 weeks. Approved researchers will be required to comply with applicable data protection regulations and data use agreements, and data may be used for research purposes only and may not be redistributed.

The raw genotype data are protected and not publicly available due to data privacy laws. Access to the genotype data can be granted to qualified researchers upon approval by the relevant institutional ethics committee, subject to participant consent and data protection regulations. Requests for genotype data access should be directed to the corresponding author. Whole sQTL summary data have been deposited in the OMIX database (https://ngdc.cncb.ac.cn/omix/) under accession number OMIX015281. The GWAS-sQTL colocalization results and lead variants for cis-sQTL and trans-sQTL analyses are provided in Supplementary Data.

The GWAS data used in this study are provided in Supplementary Data 9. The ATAC-seq data of cartilage used in this study are available in the GEO database under accession code GSE108301 [https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE108301]. The ATAC-seq data for synovium used in this study are available in the GSA-human database under accession code HRA004624 [https://ngdc.cncb.ac.cn/gsa-human/browse/HRA004624]. The processed peak files for six histone modifications (H3K27ac, H3K27me3, H3K36me3, H3K4me1, H3K4me3, and H3K9me3) of cartilage are available in the GEO under accession code GSE188124 [https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE188124], GSE187652 [https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE187652], GSE187044 [https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE187044], GSE187055 [https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE187055], GSE187641 [https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE187641], and GSE187471 [https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE187471]. The processed peak files for six histone modifications (H3K27ac, H3K27me3, H3K36me3, H3K4me1, H3K4me3, and H3K9me3) of synovium are available in the GEO under accession code GSE112655 [https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE112655].

Code availability
All analyses were conducted using publicly available software tools, with no custom code or proprietary methods developed.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Key results: In this manuscript, the authors have identified genetic variants that impact alternative splicing of transcripts in cartilage and synovium tissue. They have identified novel sQTLs specific to the synovium and cartilage, and showed that many of the sQTLs colocalize with RA, OA and height GWAS hits, indicating splicing is likely the mechanism of these previously unexplained GWAS hits. The authors also showed that sVariants tended to be tissue-specific and sGenes enriched in trait-relevant pathways.\nValidity: The data interpretation and conclusions appear to be generally valid. Answering the questions posed below could help the reader further interpret these results.\nSignificance: This group reports previously unreported sQTLs for synovium and cartilage, of which some can help explain the mechanism of some GWAS hits. The results further strengthens the evidence that splicing constitutes a significant causal mechanism underlying disease process including genetic risk, in addition to gene expression regulation which historically received more attention from the genetics community.\nData and Methodology: The RNA-seq and genotyping data appear to be of high quality with a large number of samples and good QC methods as described in their previous papers. However, the following questions should be addressed:\n• One caveat with the dataset is the genotyping panel does not appear to be very dense, therefore could be lacking the true splice-altering SNPs. It would help if the authors could comment on this and provide a distribution of all the variants detected (regardless of if they affect splicing or expression); what percent of them are intronic, intergenic, exonic, splice site, 5’ UTR, 3’ UTR?\n• In comparing the potential role of sQTL vs. eQTL as causal mechanism underlying GWAS risk loci, the eQTL dataset used was previously generated also from joint tissues. This limited-scope comparison is not adequate for inferring causal biology as the biology behind any risk locus may very likely arise from other non-joint tissue and cell types. See suggestions below.\n• For sQTL mapping, the authors cite using LeafCutter, however the leafcutter manual recommends using FastQTL for mapping before performing the sQTL analysis in LeafCutter. This appears to be the methodology GTEx followed. Is this what the authors did? Could they expand their methodology accordingly?\nSuggested improvements:\n• Could the authors confirm the clinical phenotypes of the 238 donors? Specifically do they all have OA (it appears so based on the methods in the cited papers), or did some have RA? If there were non-OA comparators, would be very interesting to comment on how disease status impacts splicing, i.e. whether the sGenes overlapping with the GWAS hits were differentially spliced in OA vs non-OA.\n• Figure 2d: Could the authors also show an example of the same variant impacting two intron junctions, and explain how this can happen?\n• Figure 2j: Could the authors further elaborate on how close most of the sQTLs are to the target? (e.g. are 90% within 1kb? Etc). A zoomed-in plot like Figure 3c (not restricting to the SpliceAI-predicted variants) could help.\no For the variants further away (e.g. >50kb outside of the gene), is there a hypothesis on how those variants could be impacting splicing? Are they in LD with a closer variant that is not imputed which is the functional one?\no How does the distance of sQTLs from synovium and cartilage compare to the sQTLs from different tissues from GTEx? The GTEx dataset has a whole genome sequencing data & dense genotyping panel, therefore perhaps this comparison could inform if the density of the genotyping panel is affecting these results.\no From what can be gathered in the main text and methods, the eQTLs used for comparison with sQTLs for explaining disease biology and genetic risk are from the authors’ previous eQTL analyses in the same tissues. One cannot assume that RA/OA pathogenesis is concentrated in the cartilage/synovium, therefore it is not adequate to restrict the search space to joint tissues (or any tissue) when making statements about sQTL vs. eQTL underlying genetic associations. The authors should consider either repeating/expanding parts of the analysis related to GWAS loci interpretation by including larger, multi-tissue eQTL datasets such as GTEx, and/or discuss the caveats around findings of sQTL vs. eQTL in explaining genetic associations.\n• What percentage of eQTLs and sQTLs impacted multiple genes vs just a single gene?\n• What percentage of the alternative splicing events were novel vs previously annotated?\n• Figure 3A: Please define ‘splice region’. Are there also otherwise ‘exonic’ variants? The confidence that the functional variant is the causal variant is likely highest for the splice acceptor, donor and splice region variants. For the remaining regions, do you suspect the actual functional/causal variant was not genotyped/imputed? Or could it be true biology that intergenic variants are affecting splicing.\n• It is nice to see that the variants that impact the canonical AG or GT splice motifs have higher odds ratios than variants at other locations, as expected (Figure 3A). It would be helpful to also show the percent change splicing they result in? (average and range would be helpful). They should result in 100% missplicing for homozygotes and 50% for heterozygotes, is this true?\n• Figure 3B: Please add the total number of significant sVariants to this diagram, to show what % of these are predicted to be functional by SpliceAI.\n• Could the authors show the distribution of the SpliceAI scores for all sQTLs? Would be interesting to learn more about what percentage SpliceAI misses. Of those, is the SpliceAI algorithm less-well trained on synovium and cartilage-specific data? Or do you think the lead sQTL is not the functional variant; the actual functional variant is not typed?\n• “Only about 27% of sGenes shared the same causal variant for both alternative splicing and gene expression regulation” -> In these examples, what is the suspected mechanism of a variant resulting in missplicing and differential expression? Are the genes differentially expressed due to splicing causing NMD? Or could transcript annotations be biasing the expression levels? If the authors could go through a few examples to potentially illustrate any patterns of their mechanisms.\n• Regarding the specificity of synovium and cartilage sQTLs: How many were shared between the synovium and cartilage? Perhaps these tissues can be added to Supplementary Figure 8, to also show if synovium and cartilage are particularly closer. Also, how many were shared with the other GTEx sQTLs? Supplementary Figure 8 implies that most are shared with other GTEx sQTLs, but it would be interesting to understand the # of non-specific sQTLs in comparison to the # of specific in Figure 4A. In addition, could the authors calculate a ‘tissue specificity score’ similar to what is in Figure 2d of Garrido-Martín, et al. (2021) Nat Commun. doi: 10.1038/s41467-020-20578-2.? This could give us a sense of how specific the synovium and cartilage are relative to the other tissues.\n• The authors mentioned that tissue-specific sQTL lead variants were located 195 farther from their target splice sites than shared sQTLs. Please clarify the term ‘sQTL lead variants’; were these the variants deemed functional through SpliceAI? If so, SpliceAI could be introducing a bias in this finding.\n• The colocalization of sQTLs with RA, OA and height GWASs are a very exciting finding. It would be very helpful to expand more on this. Could the authors show the specific splice events (similar to what was shown on COL2A1) in the supplement? The intron_id doesn’t seem to be sufficient to understand exactly what is happening in each of these splice events. Particularly for the OA and RA hits, if plots similar to what is in Figure 3d could be added for each of these events in the supplement (14 plots for OA + 11 for RA), that would make these findings even more actionable.\n• Could the authors share the coordinates of exon 2 to ensure the reader understands the correct exon (due to differences in annotations, isoforms, etc “exon 2” could mean different things). Is the function of exon 2 known, is this an annotated exon? Similarly I would be curious if anything was previously known about any of the other 7 OA-splice-specific loci.\n• Figure 6A: Could the authors add eQTL values to this plot as well?\n• Of the synovium and cartilage sQTLs that colocalize with GWAS hits for OA, RA and height, how many were also sQTLs in GTEx? How many were novel?\n• Line 251 and Figure 6f: The RA GWAS signal at GSDMB had also been previously demonstrated to colocalize with an eQTL in (peripheral) T cells https://pubmed.ncbi.nlm.nih.gov/24968232/#&gid=article-figures&pid=figure-3-uid-2\n• It would help if the authors could provide a bit more information on individual sQTL splice events for the reader to understand the specific splice changes, and whether each of these spliceQTLs are novel, especially for the sQTLs that are GWAS hits.\n• There are a few minor typos/grammatical errors throughout the manuscript and figure captions:\no Line 62: “-related diseases and traits”\no Line 72: “for cartilage or synovium”\no Line 194: missing a period after “selection”\no Line 223: “COLA2A1” should say “COL2A1”.\no Line 280: “in a tissue-specific manner”\no Line 346: “inferred to be sGenes”?\no Figure 2g: typo in x-axis title", "reviewer_2": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_3": "This manuscript presents a high-quality and comprehensive analysis of splicing quantitative trait loci (sQTLs) in cartilage and synovium—two key joint tissues implicated in osteoarthritis (OA), rheumatoid arthritis (RA), and human height. Using RNA-seq and genotype data from 238 donors, the authors identified 4,133 independent cis-sQTLs and three trans-sGenes, including 309 tissue-specific sQTLs. Integration with GWAS data revealed novel effector genes for OA (14 genes), RA (11 genes), and height (276 genes), many of which were not detectable via eQTL analysis. The study also employed SpliceAI to fine-map 207 functional splice-altering variants, highlighting both canonical and distal regulatory mechanisms. These findings provide a valuable resource and offer new insights into the tissue-specific genetic regulation of alternative splicing and its role in complex trait etiology.\nOverall, the manuscript is well-conceived, methodologically rigorous, and biologically meaningful. It represents a significant advance in the field of post-transcriptional regulation and joint disease genetics. The integration of sQTLs with GWAS data is particularly impactful, revealing regulatory mechanisms that are not captured by gene expression alone.\nMajor Strengths\n• Novel Resource: First large-scale sQTL map for cartilage and synovium, filling a critical gap in joint tissue genomics.\n• Methodological Rigor: Use of advanced computational tools (LeafCutter, tensorQTL, SpliceAI, mash, coloc) and robust statistical frameworks.\n• Biological Insight: Identification of tissue-specific splicing events and their links to disease-associated loci provides mechanistic depth.\n• Data Integration: Effective combination of sQTL and GWAS data enables discovery of effector genes not captured by eQTL analysis.\nMinor Points\n• SpliceAI Scope: SpliceAI predictions are limited to variants within annotated genes, potentially missing regulatory variants in intergenic or enhancer regions.\n• Trans-sQTL Power: Only three trans-sGenes were identified, likely due to limited sample size. This aspect could be expanded in future studies with larger cohorts.\n• Visualization Density: Some figures are dense and could benefit from clearer annotations and legends to improve accessibility for readers.\n• Terminology Consistency: Ensure consistent use of terms such as “sVariants,” “sGenes,” and “effector genes” throughout the manuscript.", "reviewer_4": "This manuscript is a complementary study to a recent publication by the same group, which used what appear to be largely the same samples for an eQTL analysis. In the present work, the authors perform a comprehensive splicing QTL (sQTL) analysis on more than 200 samples from knee synovium and cartilage obtained from patients undergoing knee replacement surgery.\nFrom the standpoint of general novelty, there is limited conceptual innovation, as sQTLs have been recognized since 2013 as a major mechanism through which human genetic variants and particularly those identified in GWAS affect gene regulation. That said, the quality of this study is superb. The tools used are standard but well-implemented, the analyses are rigorous, and the overall presentation is of high quality. This paper will be of interest to researchers studying rheumatoid arthritis, inflammatory joint disorders, and other autoimmune conditions. The data quality appears excellent, the illustrative examples are compelling, and the lists of significant sQTLs provided in the supplementary tables will be valuable to the community.\nWe have no majo comments, mostly specific comments and clarifications:\n• It is unclear why the number of samples used in this study (238) differs from the previous eQTL study (243), perhaps explaining the sample retention criteria would help especially if it differed from their previous study.\n• Figure 2e has a typo in the x axis label. Please change “ sVraint” to “sVariant”.\n• Figure 3a: Please clarify what the blue and red distributions represent. I assume that red corresponds to sVariants and blue to eVariants, but this should be explicitly stated.They are specified in figure 6.\n• SpliceAI threshold: The rationale for using a 0.3 cutoff for SpliceAI is unclear. Was this based on prior benchmarking or empirically determined in this dataset? Further, calling SpliceAI overlap “cross-validation” is misleading, cross validation has a very specific meaning in machine learning.\n• Unannotated splice events: One potential missed opportunity was the identification of de novo or previously unannotated splice junctions arising from genetic variants. Did the analysis pipeline allow for detection of unannotated intron-exon junctions, or was it restricted to reference-annotated events?\n• Enrichment analysis: Please clarify how enrichment was computed. Was the background a random set of all genes, or a random subset of expressed genes matched for tissue expression level? This distinction is important for interpreting enrichment magnitudes.\n• Terminology (“synovium-specific”): Consider avoiding the term synovium-specific for genes detected only in synovium but not in cartilage. For instance, IRF5 was only observed in synovium here, yet it is a well-known transcription factor induced in multiple inflammatory contexts. “Detected only in synovium” or “preferentially expressed in synovium” might be less misleading than “specific.”\n• Figure 4i and distal sVariants: I assume that the example shown in Figure 4i is an illustrative case rather than a general trend, and that most distal sVariants were not located within clearly defined or shared DNA-accessible or transcriptionally active regions. Please clarify whether this was generally observed. Which fraction of lead or fine mapped sQTLs reside within accessible regions.\n• Negative control (“mental disorders”): Is “mental disorders” an appropriate negative control category, given the vastly different expression profiles between synovium/cartilage and brain? If the goal was to include a trait set with minimal overlap in expression context.\n• Power and saturation: Do the authors have an estimate of how far from saturation the current sQTL discovery is? In other words, how many additional samples would be needed to approach saturation in sQTL detection for these tissues? Such estimates would help contextualize the study’s power.\n• Data reuse and new datasets: The HRA004624 dataset was already used in the authors’ previous eQTL study (Nat. Commun., 2024; [DOI: 10.1038/s41467-024-45652-x]*). However, dataset HRA006224 was not available from GSA at the time of review, making it unclear whether it constitutes new data. The authors should clarify whether this study includes any non-overlapping samples relative to the prior eQTL analysis and what is the HRA006224 contributing to this study since it seems that all data was already available from the previous study.\n• While new methods may not have been developed. We still think the methods reported lack in detail, for example the exact pipeline ran (with software versions and parameters) would help readers be able to replicate their study.", "reviewer_5": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "Response to reviewers’ comments:\nReviewer #1:\nKey results: In this manuscript, the authors have identified genetic variants that impact alternative splicing of transcripts in cartilage and synovium tissue. They have identified novel sQTLs specific to the synovium and cartilage, and showed that many of the sQTLs colocalize with RA, OA and height GWAS hits, indicating splicing is likely the mechanism of these previously unexplained GWAS hits. The authors also showed that sVariants tended to be tissue-specific and sGenes enriched in trait-relevant pathways.\nValidity: The data interpretation and conclusions appear to be generally valid. Answering the questions posed below could help the reader further interpret these results.\nResponse: Thank you for your positive assessment. We appreciate your comments and have addressed the specific questions raised below to further clarify our interpretation of the results.\nSignificance: This group reports previously unreported sQTLs for synovium and cartilage, of which some can help explain the mechanism of some GWAS hits. The results further strengthens the evidence that splicing constitutes a significant causal mechanism underlying disease process including genetic risk, in addition to gene expression regulation which historically received more attention from the genetics community.\nResponse: Thank you for recognizing the significance of our work. We are pleased that you find the newly identified sQTLs and their relevance to GWAS loci valuable, and we agree that our findings further support the important role of splicing in disease mechanisms.\nData and Methodology: The RNA-seq and genotyping data appear to be of high quality with a large number of samples and good QC methods as described in their previous papers. However, the following questions should be addressed:\n1. Comment: One caveat with the dataset is the genotyping panel does not appear to be very dense, therefore could be lacking the true splice-altering SNPs. It would help if the authors could comment on this and provide a distribution of all the variants detected (regardless of if they affect splicing or expression); what percent of them are intronic, intergenic, exonic, splice site, 5’ UTR, 3’ UTR?\nResponse: Thank you for your thoughtful comment. We acknowledge that, as with any genotyping-based study, limited variant density may result in some true splice-altering variants not being captured by the genotyping panel. We have now explicitly discussed this limitation in the revised Discussion (Main text, Lines 368-370). We would also like to note that this limitation is not unique to array-based genotyping. Even in sQTL analyses based on whole-genome sequencing (WGS) datasets, variants with low or rare minor allele frequency (MAF < 0.01) are often excluded due to the limited statistical power arising from the small number of carriers in typical cohort sizes. As a result, some true functional variants may remain undetected even in WGS-based studies.\nFurthermore, in response to your suggestion, we annotated the genomic locations of all variants included in our analysis and present their distribution across functional categories in Figure 1. To evaluate whether our genotyping panel exhibits any systematic bias in genomic coverage, we compared the distributions with those from the GTEx WGS datasets used for QTL analyses. The distributions were highly similar, indicating that our genotype data do not show an obvious bias towards or against specific genomic regions. These results have been included in the revised manuscript as a supplementary figure (Supplementary Figure 1e in revised manuscript).\nFigure 1. The (a) number and (b) percentage of variants annotated to different genomic regions in the GTEx dataset and in our study. Note: The percentages sum to more than 100% because some variants map to multiple genomic regions owing to transcript-specific annotations.\n2. Comment: In comparing the potential role of sQTL vs. eQTL as causal mechanism underlying GWAS risk loci, the eQTL dataset used was previously generated also from joint tissues. This limited-scope comparison is not adequate for inferring causal biology as the biology behind any risk locus may very likely arise from other non-joint tissue and cell types. See suggestions below.\nResponse: Thank you for this important comment. We agree that restricting comparisons to sQTL and eQTL datasets derived solely from cartilage and synovium may be insufficient to fully infer the regulatory mechanisms underlying GWAS loci, as some loci may exert their effects through gene regulation in other tissues or cell types. To address this concern, we expanded our colocalization analyses to include sQTL and eQTL datasets from 49 GTEx tissues, thereby enabling a multi-tissue assessment of regulatory mechanisms.\nThese extended analyses revealed that a subset of identified candidate disease-associated genes, such as SERPINF1 and INTS2 for OA (Figure 2a), HIC2 and GATS for RA (Figure 2b), colocalized exclusively with sQTL signals identified in our joint tissue dataset and were not detected in any GTEx tissue. These observations strengthen the evidence that the regulatory effects at these loci are likely mediated through splicing regulation in cartilage or synovium, highlighting the importance of constructing splicing maps directly from disease-relevant joint tissues. We have integrated these results into the revised manuscript (Main text, Lines 248-250, Lines 289-290, Lines 533-537; Figure 5b,c, Figure 6d; Supplementary Tables 10-12).\nAt the same time, we acknowledge that colocalization analyses alone cannot definitively establish causal genes or regulatory mechanisms, nor can they exclude potential contributions from other tissues. Accordingly, we have revised the main text to avoid definitive language and now describe these loci as potentially exerting their effects through splicing regulation in cartilage or synovium. Further functional validation will be required to determine the precise molecular mechanisms involved.\nFigure 2. Colocalization of OA and RA GWAS risk loci with sQTL and eQTL signals across our joint tissues and GTEx tissues. (a) Scatter pie plots summarizing colocalization results between osteoarthritis (OA) GWAS loci and QTL signals derived from cartilage, synovium, and 49 GTEx tissues. Left panel, colocalization with sQTLs. Right panel, colocalization with eQTLs. Each row represents a colocalized gene, and each column corresponds to an OA GWAS phenotype. For each gene–phenotype pair, coloured sectors within a circle indicate the tissues in which a colocalization signal was detected: dark blue, colocalization with cartilage QTLs; dark red, colocalization with synovium QTLs; grey, colocalization with GTEx QTLs. (b) Bar plots summarizing the colocalization results between rheumatoid arthritis (RA) GWAS loci and QTL signals derived from synovium and GTEx tissues. Each row represents a colocalized gene. Left panels, colocalization with synovium sQTLs and eQTLs, respectively. Right panels, colocalization with sQTLs and eQTLs from GTEx tissues. Bars show the posterior probability of colocalization (PP4) for each gene–tissue pair.\n3. Comment: For sQTL mapping, the authors cite using LeafCutter, however the leafcutter manual recommends using FastQTL for mapping before performing the sQTL analysis in LeafCutter. This appears to be the methodology GTEx followed. Is this what the authors did? Could they expand their methodology accordingly?\nResponse: Thank you for pointing this out. We realize that our original wording may have caused some confusion. To clarify, LeafCutter is a tool developed for alternative splicing quantification by generating intron cluster-based splicing phenotypes (PSI values), which represent the relative usage of different splicing events. These PSI values serve as the phenotypic data for subsequent sQTL analysis.\nFor the sQTL mapping itself, we used tensorQTL, which integrates genotype, phenotype, and covariates to perform cis-sQTL mapping. TensorQTL implements the same statistical framework as FastQTL but is significantly faster. This approach aligns with the general methodology used in the GTEx project. We have updated the Results and Methods section to provide a more accurate description of the roles of LeafCutter and tensorQTL in our analysis pipeline (Main text, Lines 89-90, Lines 403-404, Lines 415-417).\nSuggested improvements:\n4. Comment: Could the authors confirm the clinical phenotypes of the 238 donors? Specifically do they all have OA (it appears so based on the methods in the cited papers), or did some have RA? If there were non-OA comparators, would be very interesting to comment on how disease status impacts splicing, i.e. whether the sGenes overlapping with the GWAS hits were differentially spliced in OA vs non-OA.\nResponse: Thank you for this insightful comment. We confirm that all 238 donors included in this study were clinically diagnosed with osteoarthritis (OA), and none had rheumatoid arthritis (RA). We agree that, if non-OA comparators were available, a direct comparison of splicing profiles between OA and non-OA tissues would have provided a more explicit assessment of disease-related splicing changes, particularly for sGenes colocalizing with OA GWAS loci. Unfortunately, such comparator samples were not available in our cohort.\nTo address this question to the extent permitted by our dataset, we examined whether within-donor comparisons could offer informative evidence. For cartilage, paired degenerated and macroscopically intact samples were available from 31 donors (Figure 3a). Although only degenerated cartilage samples were used for sQTL and eQTL mapping to maintain consistency in genetic association analyses, these paired samples nonetheless provided an opportunity to explore disease-associated splicing alterations. Notably, COL2A1, a sGene that colocalizes with OA GWAS signals, exhibited consistent and significant differential splicing between degenerated and intact cartilage (Figure 3b-e). Degenerated tissue showed markedly higher retention of exon 2. Together with our GWAS-sQTL colocalization results and the observation that the GWAS risk allele G promotes exon 2 inclusion, this finding further supports a model in which increased exon 2 retention may contribute to OA pathogenesis. These data, analyses and interpretations have been incorporated into the revised manuscript (Main text, Lines 265-269, Lines 384-389, Lines 549-555; Figure 5g; Supplementary Figure 1b).\nFigure 3. Colocalization and differential splicing analysis reveal that increased exon 2 retention of COL2A1 may contribute to OA pathogenesis. (a) The sample information for cartilage. (b) LocusZoom plot showing a KneeOA GWAS signal colocalized with COL2A1 sQTL signal in cartilage but not with its eQTL signal. (c) Sashimi plot showing association between rs1793955 genotypes and the exon 2 usage ratio of COL2A1 in cartilage. (d) Boxplot showing association between rs1793955 genotypes and usage ratio of introns of COL2A1 as shown in panel c. The P values were reported by tensorQTL. (e) Intron usage ratio comparison of COL2A1 between degenerated and intact cartilage samples. Data are shown for 31 paired samples. The P values were calculated using a paired Student’s t-test.\nNo paired samples were available for synovium, precluding similar analyses in that tissue. We also acknowledge that the absence of non-OA comparator tissues and the limited number of paired cartilage samples restrict our ability to systematically evaluate disease-associated splicing changes across all colocalized sGenes. These limitations and their implications have been discussed in the revised Discussion section (Main text, Lines 364-370).\n5. Comment: Figure 2d: Could the authors also show an example of the same variant impacting two intron junctions, and explain how this can happen?\nResponse: Thank you for this comment. In the revised manuscript (Main text, Lines 117-119; Supplementary Figure 7), we have added an example illustrating a case in which a single variant is associated with two distinct intron clusters, as requested (Figure 4).\nFigure 4. Example illustrating that rs56167434 regulates alternative splicing of two intron clusters in SH3YL1. (a) Transcript model of SH3YL1 showing the two intron clusters (highlighted in red and blue, respectively) and the genomic location of rs56167434. (b) Boxplots showing the association between rs56167434 genotypes and normalized intron usage for each of the two intron clusters.\nRegarding the underlying regulatory mechanism, we conducted an extensive literature review but did not identify prior studies that directly address how a single genetic variant can simultaneously regulate multiple intron clusters. To further explore this question, we performed comparative enrichment analyses between variants associated with multiple intron clusters and those associated with a single intron cluster, examining genomic context and epigenetic features. However, we did not observe clear or consistent differences between these two groups.\nBased on these analyses, we must acknowledge that we are currently unable to propose a definitive mechanistic explanation for this phenomenon. We believe that the coordinated regulation of multiple intron clusters by a single variant represents an interesting and potentially important biological question, but a detailed investigation of the underlying mechanisms would require additional experimental and analytical approaches and therefore lies beyond the scope of the present study.\n6. Comment: Figure 2j: Could the authors further elaborate on how close most of the sQTLs are to the target? (e.g. are 90% within 1kb? Etc). A zoomed-in plot like Figure 3c (not restricting to the SpliceAI-predicted variants) could help.\nResponse: Thank you for this helpful suggestion. As recommended, we generated histograms showing the distance between lead sVariants (the sVariants with smallest P value for each significant sCluster) and their corresponding target splice sites (Figure 5). Within the 1-Mb cis window used for sQTL mapping, we observed that ~75% of lead sVariants lie within 50 kb of the target splice site, and more than 90% fall within 200 kb. These results have been added to the revised manuscript (Main text, Lines 155-157; Supplementary Figure 11).\nFigure 5. Distribution of distances between lead sVariants and their target splice sites in (a) cartilage and (b) synovium. Histograms show the number of sVariants falling into each 10-kb bin, and the red curve indicates the cumulative fraction of sVariants within increasing distances. The lower panels provide a zoomed-in view of the 0-200 kb range.\n7. Comment: For the variants further away (e.g. >50kb outside of the gene), is there a hypothesis on how those variants could be impacting splicing? Are they in LD with a closer variant that is not imputed which is the functional one?\nResponse: We greatly appreciate your interest in the regulatory mechanisms underlying distally located sVariants (e.g., >50kb).\nWe would like to clarify that the variants shown in Figure 2j in the manuscript used for distance analysis are the lead sVariants, which are defined as the variant with the smallest P value within each sQTL signal. However, it is important to note that lead sVariants do not necessarily represent the true functional variants, and may be in linkage disequilibrium (LD) with more proximal variants that may have functional effects.\nIn response to this, we performed an additional analysis in which, for each sQTL signal, we identified both the lead sVariant and the variant within the 95% fine-mapping credible set (which represents the region containing the true functional variant with 95% confidence) that is closest to the target gene (Figure 6a). We then categorized the distances between these variants and the target genes, with variants within the gene assigned a distance of zero (Figure 6b). As illustrated in Figures 6c and 6d, we observed that for 60% sQTLs with lead sVariants outside gene, another variant within the 95% fine-mapping credible set resided within the gene itself, suggesting that the true functional variant may be intragenic, potentially influencing splicing by altering canonical splice site motifs, splicing regulatory element sequences (including exon/intron splicing enhancers and silencers) 1, and mRNA secondary structure (including G-quadruplex) 2-4. These changes could, in turn, affect the binding of trans-acting splicing factors, including the spliceosome and other auxiliary RNA-binding proteins (RBPs) 5-8, ultimately modulating splicing events.\nInterestingly, for about 20% of genes in cartilage and synovium, all the variants within the credible set were located outside the target gene (Figure 6d, f). While this may be partly due to the incomplete variant coverage, it is also consistent with distal regulatory mechanisms. Previous studies have shown that enhancers9,10, which often reside far from their target genes, can influence alternative splicing. To explore this possibility, we performed epigenomic annotations on the variants in the credible sets for these genes. We found that over half of these genes harbored at least one credible-set variant overlapping enhancer elements (identified by ATAC-peak, H3K27ac-peak, or H3K4me1-peaks). We hypothesize that these variants may modulate enhancer activity, thereby influencing splicing of the target gene.\nIn the revised manuscript, we discussed these candidate alternative regulatory mechanisms and provided two systematic annotation tables for all sVariants in fine-mapped 95% credible set to show whether they impact the RBP binding or located in candidate enhancer regions (Main text, Lines 328-340; Supplementary Tables 13,14). However, we also acknowledge that these distal regulatory mechanisms remain speculative, and require further experimental validation.\nFigure 6. Analysis of sVariants and their relationship with target genes and splicing regions. (a) Schematic diagram illustrating how gene regions, splice regions, and the distance to the target gene are defined. (b) Classification of lead sVariants and their closest fine-mapped sVariants partner from the 95% credible set based on their locations relative to the target gene. The variants are categorized into gene regions, splice regions, and distal regions, with color coding to represent different spatial relationships. (c, e) Bar plot showing the distribution of distances between lead sVariants and target genes of cartilage (c) and synovium (e), categorized by bin width (10 kb). The lead sVariants are classified based on their distance from the target gene, with the closest fine-mapped sVariants shown in relation to the lead sVariants. (d, f) Percentage of genes with no fine-mapped sVariant located within the gene (left bar plot) and the percentage of genes with a fine-mapped sVariant within the gene and located in enhancer regions (right bar plot).\n8. Comment: How does the distance of sQTLs from synovium and cartilage compare to the sQTLs from different tissues from GTEx? The GTEx dataset has a whole genome sequencing data & dense genotyping panel, therefore perhaps this comparison could inform if the density of the genotyping panel is affecting these results.\nResponse: Thank you for this insightful comment. To address this question, we compared the median distance between sVariants and their target splice sites in two tissues of our study with those from the GTEx project. As shown in Figure 7, the distances between sVariants and their target splice sites in cartilage and synovium were consistent with those observed in other tissues. This suggests that our results are not substantially influenced by genotyping panel density.\nFigure 7. The median distance between sVariants and their target splice site across our joint tissues and GTEx tissues.\n9. Comment: From what can be gathered in the main text and methods, the eQTLs used for comparison with sQTLs for explaining disease biology and genetic risk are from the authors’ previous eQTL analyses in the same tissues. One cannot assume that RA/OA pathogenesis is concentrated in the cartilage/synovium, therefore it is not adequate to restrict the search space to joint tissues (or any tissue) when making statements about sQTL vs. eQTL underlying genetic associations. The authors should consider either repeating/expanding parts of the analysis related to GWAS loci interpretation by including larger, multi-tissue eQTL datasets such as GTEx, and/or discuss the caveats around findings of sQTL vs. eQTL in explaining genetic associations.\nResponse: Thank you for raising this important point. As noted in our response to the previous comment, we agree that relying solely on eQTL and sQTL maps derived from cartilage and synovium would not be sufficient to capture the full range of biological mechanisms through which GWAS loci may act. Regulatory variants influencing OA or RA risk could, in principle, exert their effects in tissues outside the joint, and limiting the analysis to joint-derived QTLs could therefore bias interpretation.\nTo address this concern, we broadened our colocalization analyses by incorporating eQTL and sQTL datasets from 49 GTEx tissues, thereby enabling a cross-tissue evaluation of regulatory mechanisms at disease/trait associated loci. This expanded analysis revealed that several loci, including SERPINF1 and INTS2, showed convincing colocalization only with the sQTL signals detected in our joint-tissue dataset and not with QTLs from any GTEx tissue (Figure 8). These findings suggest that, for at least a subset of OA and RA risk loci, splicing regulation in cartilage or synovium may play a particularly prominent role. We have integrated these results in our revised manuscript (Main text, Lines 248-250, Lines 289-290, Lines 533-537; Figure 5b,c, Figure 6d; Supplementary Tables 10-12).\nNevertheless, we acknowledge that colocalization alone does not establish causality, and we cannot exclude the possibility that additional regulatory effects may arise in tissues not assayed in our dataset. We have revised the manuscript to explicitly highlight these caveats and to adopt more cautious phrasing regarding the regulatory mechanisms that may underlie individual GWAS loci. Additional functional studies will be required to delineate the tissue- and mechanism-specific roles of these variants.\nFigure 8. Colocalization of OA and RA GWAS risk loci with sQTL and eQTL signals across our joint tissues and GTEx tissues. (a) Scatter pie plots summarizing colocalization results between osteoarthritis (OA) GWAS loci and QTL signals derived from cartilage, synovium, and 49 GTEx tissues. Left panel, colocalization with sQTLs. Right panel, colocalization with eQTLs. Each row represents a colocalized gene, and each column corresponds to an OA GWAS phenotype. For each gene–phenotype pair, coloured sectors within a circle indicate the tissues in which a colocalization signal was detected: dark blue, colocalization with cartilage QTLs; dark red, colocalization with synovium QTLs; grey, colocalization with GTEx QTLs. (b) Bar plots summarizing the colocalization results between rheumatoid arthritis (RA) GWAS loci and QTL signals derived from synovium and GTEx tissues. Each row represents a colocalized gene. Left panels, colocalization with synovium sQTLs and eQTLs, respectively. Right panels, colocalization with sQTLs and eQTLs from GTEx tissues. Bars show the posterior probability of colocalization (PP4) for each gene–tissue pair.\n10. Comment: What percentage of eQTLs and sQTLs impacted multiple genes vs just a single gene?\nResponse: Thank you for your comment. As shown in Figure 9, the vast majority of variants regulate only a single target gene. Specifically, 92.2% of significant sVariants in cartilage and 89.9% in synovium were associated with only one gene. For eQTLs, the corresponding proportions were 79.2% in cartilage and 74.6% in synovium. These results indicate that most regulatory variants in our dataset impact a single gene rather than multiple genes. We have added these results in our revised manuscript (Main text, Lines 111-112; Supplementary Figure 6).\nFigure 9. Distribution of the number of target genes regulated by each significant sVariant and eVariant in cartilage and synovium. Histograms (grey bars) show the count of variants associated with specific number of target genes, and the red line represents the cumulative fraction of variants.\n11. Comment: What percentage of the alternative splicing events were novel vs previously annotated?\nResponse: Thank you for this important question. We would like to clarify that the classification of alternative splicing events as “novel” or “previously annotated” depends strongly on the annotation reference used. To provide a comprehensive assessment, we annotated all significant sQTL-associated introns using two independent resources.\nFirst, we used the GENCODE v19 annotation, which matches the hg19 reference genome used for read alignment in our study. Based on this annotation, 83% of significant introns in cartilage and 81% in synovium were fully annotated (i.e., both splice boundaries matched GENCODE v19). An additional 11% (cartilage) and 13% (synovium) had only one annotated boundary, while the remaining 6% of introns in both tissues had neither boundary annotated in GENCODE v19 (Figure 10a).\nSecond, we annotated the same set of introns using Snaptron11, a large-scale splice-junction database constructed from tens of thousands of public human RNA-seq samples. Snaptron aggregates junction-level evidence from widely used datasets, including SRA, GTEx, TCGA, as well as large uniformly processed resources such as recount2 and recount3, thus providing population-scale support for both annotated and unannotated splice junctions. Using this database, we found that only 11 introns (0.66%) in cartilage and 20 introns (0.64%) in synovium lacked support in Snaptron (Figure 10b).\nThese results have been added to the revised manuscript, where we now provide a detailed description of the annotation strategy and its implications for determining whether sQTL-associated splicing events are novel (Main text, Lines 98-107; Supplementary Figure 4).\nFigure 10. Novelty annotation of significant sQTL-associated introns in cartilage and synovium. (a) Classification of introns according to whether both intron boundaries, one intron boundary, or neither intron boundary are annotated in GENCODE v19. (b) Annotation of the same introns using the Snaptron splice-junction database, which aggregates junction evidence from SRA, GTEx, TCGA, and large-scale uniformly processed datasets such as recount2 and recount3.\n12. Comment: Figure 3A: Please define ‘splice region’. Are there also otherwise ‘exonic’ variants? The confidence that the functional variant is the causal variant is likely highest for the splice acceptor, donor and splice region variants. For the remaining regions, do you suspect the actual functional/causal variant was not genotyped/imputed? Or could it be true biology that intergenic variants are affecting splicing.\nResponse: Thank you for this helpful comment. Figure 3A presents the enrichment of lead sVariants across different genomic regions. As you noted, this includes a category labelled as the “splice site”. We would like to clarify that this functional annotations in our study were generated using SnpEff, and the definitions therefore follow SnpEff’s terminology. Specifically, the splice site region is defined as a sequence variant occurring within 1–3 bases inside the exon or 3-8 bases inside the intron adjacent to an exon–intron boundary. The splice donor region corresponds to the two bases immediately downstream of the coding exon end (except for the last exon), and the splice acceptor region corresponds to the two bases immediately upstream of the exon start (except for the first exon). These definitions have now been explicitly added to the revised figure legend of Figure 3A.\nWith respect to the reviewer’s second question regarding variants located outside the canonical splice acceptor, donor, or splice-site regions, we agree that confidence in direct functional effects on splicing is generally highest for variants residing within these canonical motifs. For sVariants located in other exonic, intronic, or intergenic regions, we acknowledge that we cannot fully exclude the possibility that the true causal variant was not directly genotyped.\nNevertheless, prior studies have demonstrated that variants located within genes but outside canonical splice-site regions can still influence alternative splicing by altering RNA-binding protein (RBP) interactions5-8. Here we give an illustrative example: sQTL analysis in cartilage samples revealed that rs6720 regulates the usage of a downstream intron (chr7:75677544-75677894), located approximately 40 bp from the splice site. Specifically, the T allele is associated with increased intron usage, whereas the C allele is associated with reduced usage (Figure 11a).\nBased on the allele-specific RBP binding data from K562 cells6, the T allele at this locus is associated with reduced binding of DROSHA, a protein that has been reported to participate in the regulation of alternative splicing12,13. This observation suggests a negative regulatory role of DROSHA binding on intron usage at this site. Consistent with this model, we further observed a negative correlation between DROSHA expression levels and intron usage (Figure 11b), providing additional support for the proposed regulatory mechanism.\nFigure 11. In silico analysis suggest that rs6720 regulate alternative splicing by disrupting RNA-binding of DROSHA. (a) Sashimi plot showing the association between rs6720 genotypes and usage ratio of intron located 40 bp downstream of it in cartilage samples. (b) Scatter plot showing the correlation between DROSHA expression and the usage ratio of its potential target intron.\nIn addition, for sVariants located in intergenic regions, there is evidence that enhancers, which are frequently located outside gene bodies, can result in exon skipping by promoting transcription elongation rate. Genetic variants affecting enhancer activity may therefore influence splicing through long-range regulatory mechanisms. As an example, in our results, we found that rs79194220, which simultaneously regulates both splicing and expression of the CYP3A5. This variant is located in an enhancer region downstream of CYP3A5, and the G allele promotes both gene expression and an intron skipping event (highlighted in red in Figure 12c).\nTaken together, these observations suggest that both technical limitations and genuine biological mechanisms may contribute to the presence of sVariants outside canonical splice-site regions. In revised manuscript, we discussed these candidate alternative regulatory mechanisms and provided two systematic annotation tables for all sVariants in fine-mapped 95% credible set to show whether they impact the RBP binding or located candidate enhancer regions (Main text, Lines 328-340; Supplementary Tables 13,14).\nFigure 12. Epigenomic annotation suggested that rs79194220 might regulate CYP3A5 alternative splicing through enhancer mechanism. (a) Gviz plot showing the QTL signal of CYP3A5 and the enhancer annotation of rs79194220. (b) Boxplot showing the association between rs79194220 genotypes and the normalized expression of CYP3A5. (c) Sashimi plot showing the association between rs79194220 genotypes and the splicing of CYP3A5.\n13. Comment: It is nice to see that the variants that impact the canonical AG or GT splice motifs have higher odds ratios than variants at other locations, as expected (Figure 3A). It would be helpful to also show the percent change splicing they result in? (average and range would be helpful). They should result in 100% missplicing for homozygotes and 50% for heterozygotes, is this true?\nResponse: Thank you for this suggestion. As expected, our analysis of variants disrupting the canonical AG/GT splice site motifs showed clear genotype-dependent effects on splice-site usage. Heterozygous individuals typically exhibited a substantial reduction in splicing (ranging from 20% to 70% of the normal level), whereas individuals homozygous for the motif-disrupting allele showed almost complete loss of splicing at the site (Figure 13). These results have been added to the revised manuscript (Main text, Lines 188-194; Supplementary Figure 15).\nFigure 13. Splicing ratios at splice sites disrupted by canonical AG or GT splice-site variants across different genotypes. For the homozygous genotype carrying the canonical AG or GT splice motif, the splicing ratio at the affected site was used as the baseline and set to 100%. Splicing ratios for the other two genotypes were calculated relative to this baseline.\n14. Comment: Figure 3B: Please add the total number of significant sVariants to this diagram, to show what % of these are predicted to be functional by SpliceAI.\nResponse: Thank you for the suggestion. We have added a stepwise summary showing the number of sVariants input and retained at each stage of the SpliceAI prediction and subsequent validation (Supplementary Figure 13), which allows the proportion predicted to be functional by SpliceAI to be easily assessed.\n15. Comment: Could the authors show the distribution of the SpliceAI scores for all sQTLs? Would be interesting to learn more about what percentage SpliceAI misses. Of those, is the SpliceAI algorithm less-well trained on synovium and cartilage-specific data? Or do you think the lead sQTL is not the functional variant; the actual functional variant is not typed?\nResponse: Thank you for this insightful suggestion. In the revised manuscript, we now show the full distribution of SpliceAI Δ scores for all significant sVariants in each tissue and explicitly quantify the proportion of sVariants that are not prioritized by SpliceAI under commonly used thresholds. We also added a stepwise summary showing how many sVariants are retained at each stage of the SpliceAI prediction and subsequent validation (Supplementary Figure 13).\nUsing cartilage as an illustrative example (Figure 14), we started from 80,754 significant sVariants. Because SpliceAI requires variants to be mapped to gene models, 50,903 sVariants (63.0%) fell within gene regions were analyzed by SpliceAI. The resulting Δ score distribution was strongly enriched near zero, with only a small subset (358/50903, 0.7%) exceeding the prioritization threshold used in our study (Δ score ≥ 0.2), and therefore carried forward for validation against the observed sQTL signals. Among these, 62 sVariants (17.3%) showed concordant splicing effects between SpliceAI predictions and our sQTL results, and were considered high-confidence causal candidates within their corresponding sQTL signals.\nFigure 14. Diagram showing the detailed steps of SpliceAI prediction procedure used in this study. The number of sVariants that were used as input and retained after each step was shown. The distribution of the SpliceAI scores for all sVariants was also shown.\nThese results clarify what “SpliceAI misses” in this context. Importantly, this does not necessarily reflect a limitation of the algorithm, but rather its design scope: SpliceAI is a sequence- and gene annotation-based model that primarily captures local sequence changes affecting splice sites or nearby splicing motifs (typically within a limited window around splice junctions). In contrast, a substantial fraction of sQTLs, particularly tissue-specific sQTLs, may act through mechanisms that are not well captured by local sequence features alone, such as cell-type-specific RNA-binding protein availability, chromatin accessibility, distal regulatory elements, or long-range coupling between enhancers and splicing. Because the SpliceAI model is not trained on any tissue-specific features, it is not expected to systematically capture these tissue-dependent regulatory contexts.\nWe also agree with the reviewer’s alternative explanation that the lead sVariant is not always the functional variant. For this reason, we did not restrict the analysis to lead sVariants but instead included all significant sVariants within each signal, acknowledging that the lead variant may be an LD proxy. Moreover, we now explicitly note in the discussion section of revised manuscript that the true causal variant may be untyped or excluded (e.g., due to genotyping/imputation limitations or allele frequency filtering) (Main text, Lines 362-364), which would further reduce the apparent overlap between SpliceAI-prioritized variants and sQTL signals.\nOverall, the added Δ score distributions and stepwise counts provide a transparent view of SpliceAI prioritization in our datasets.\n16. Comment: “Only about 27% of sGenes shared the same causal variant for both alternative splicing and gene expression regulation” -> In these examples, what is the suspected mechanism of a variant resulting in missplicing and differential expression? Are the genes differentially expressed due to splicing causing NMD? Or could transcript annotations be biasing the expression levels? If the authors could go through a few examples to potentially illustrate any patterns of their mechanisms.\nResponse: Thank you for raising this question. When an sGene shows evidence that the same causal variant drives both alternative splicing and gene expression regulation, several non-mutually exclusive mechanisms could explain the coupling.\nOne plausible mechanism is nonsense-mediated decay (NMD), whereby a variant increases the production of unproductive isoforms, leading to reduced steady-state mRNA abundance. In principle, NMD-mediated expression changes can be evaluated by comparing variant effects derived from nascent RNA (pre-decay) versus steady-state RNA (post-decay). As proposed in a recent work14, if a shared sQTL/eQTL signal is driven by NMD, the genotype-dependent splicing differences should be most evident in nascent RNA, whereas gene-level expression differences would be absent or substantially attenuated at this stage (i.e., weak or non-significant eQTL effects pre-decay) and become apparent only in steady-state RNA due to selective degradation of NMD-prone isoforms. Under this model, an inverse relationship between splicing changes that increase the production of NMD-sensitive isoforms and reduce steady-state gene abundance is often observed.\nHowever, our study relies on steady-state bulk RNA-seq only, therefore, the estimated sQTL effects reflect splicing outcomes after RNA surveillance and decay, and we cannot directly test the nascent-versus-steady-state predictions required to conclusively implicate NMD. Moreover, to attribute the expression change specifically to NMD, it is also necessary to exclude the possibility that the eQTL effect arises from transcriptional regulation (e.g., promoter/enhancer-driven effects on transcription initiation or elongation), which cannot be disentangled with the current dataset alone. Given these considerations, we believe it would be inappropriate to make definitive NMD-mediated mechanistic claims based solely on steady-state RNA-seq.\nRegarding the reviewer’s concern that transcript annotations could bias expression estimates, we consider this is unlikely to explain our findings. First, alternative splicing quantification was performed using LeafCutter, which quantifies intron excision/junction usage without relying on transcript isoform annotations, thereby minimizing annotation-driven bias in sQTL discovery. Second, gene expression quantification was conducted at the gene level using a consistent reference annotation across all samples, meaning that any residual annotation limitations would be expected to affect samples similarly and should not generate systematic genotype-dependent artefacts.\nOn the other hand, enhancer regulation has been established as a mechanism that can affect both splicing and gene expression, as shown in a previous study9. Enhancers can increase gene expression while simultaneously reducing the inclusion of cassette exon (an exon that can be either included in or skipped from the mature transcript) in the mature transcript. This occurs because enhanced transcription reduces the dwell time of the spliceosome at certain splice sites, increasing the likelihood of exon skipping.\nIn our study, we hypothesize that some of the variants that simultaneously affect gene expression and splicing may be acting through enhancers. To explore this hypothesis, we performed fine-mapping of our sQTL data and then annotated the variants in 95% credible set with multiple enhancer epigenetic marks. We identified that one third of the genes have at least one 95% credible set fine-mapping sVariants located within enhancer regions (marked by ATAC and H3K27ac), suggesting that these variants may regulate both expression and splicing via enhancer activity.\nOne such example is rs79194220, which simultaneously regulates both splicing and expression of the CYP3A5 gene (Figure 15). This variant is located in an enhancer region downstream of CYP3A5, and the G allele promotes both gene expression and an exon skipping event (highlighted in red in Figure 15c). This mechanism aligns with the model proposed in the previous study, where enhancer activity can influence both splicing and expression in a coordinated manner. In revised manuscript, we have provided two systematic annotation tables for all sVariants in fine-mapped 95% credible set to show whether they are located candidate enhancer regions (Supplementary Tables 13,14).\nHowever, we must acknowledge that the regulatory mechanism we propose here, including the example of rs79194220, still requires further experimental validation to confirm the functional roles of these variants in gene expression and splicing.\nFigure 15. Epigenomic annotation suggested that rs79194200 might regulate both expression and splicing of CYP3A5 through enhancer-mediated mechanism. (a) Gviz plot showing the sQTL and eQTL signal of CYP3A5 and the enhancer annotation of rs79194220. (b) Boxplot showing the association between rs79194220 genotypes and the normalized expression of CYP3A5. (c) Sashimi plot showing the association between rs79194220 genotypes and the splicing of CYP3A5.\n17. Comment: Regarding the specificity of synovium and cartilage sQTLs: How many were shared between the synovium and cartilage? Perhaps these tissues can be added to Supplementary Figure 8, to also show if synovium and cartilage are particularly closer. Also, how many were shared with the other GTEx sQTLs? Supplementary Figure 8 implies that most are shared with other GTEx sQTLs, but it would be interesting to understand the # of non-specific sQTLs in comparison to the # of specific in Figure 4A. In addition, could the authors calculate a ‘tissue specificity score’ similar to what is in Figure 2d of Garrido-Martín, et al. (2021) Nat Commun. doi: 10.1038/s41467-020-20578-2.? This could give us a sense of how specific the synovium and cartilage are relative to the other tissues.\nResponse: Thank you for this suggestion. In the revised manuscript, we have updated the original Supplementary Figure 8 (Supplementary Figure 16 in revised manuscript) to explicitly include the sharing of sQTLs between synovium and cartilage, allowing their relative proximity to be directly assessed. We also added two bar plots showing the absolute numbers of sQTLs shared between synovium/cartilage and each GTEx tissue. In addition, as described by Garrido-Martín et al. (2021), we calculated a tissue specificity score for sQTLs of cartilage and synovium, providing a quantitative measure of synovium and cartilage specificity relative to other GTEx tissues (Main text, Lines 208-210).\n18. Comment: The authors mentioned that tissue-specific sQTL lead variants were located 195 farther from their target splice sites than shared sQTLs. Please clarify the term ‘sQTL lead variants’; were these the variants deemed functional through SpliceAI? If so, SpliceAI could be introducing a bias in this finding.\nResponse: Thank you for this helpful comment. We are happy to clarify the terminology. In our analysis, sQTL lead variants refer to the variants with the smallest 𝑃 value within each significant sQTL signal, as identified by tensorQTL. Thus, lead variants are defined purely by statistical association strength and not by SpliceAI predictions. SpliceAI was used only for downstream functional annotation and was not involved in determining the set of lead sQTL variants or the distance metrics.\n19. Comment: The colocalization of sQTLs with RA, OA and height GWASs are a very exciting finding. It would be very helpful to expand more on this. Could the authors show the specific splice events (similar to what was shown on COL2A1) in the supplement? The intron_id doesn’t seem to be sufficient to understand exactly what is happening in each of these splice events. Particularly for the OA and RA hits, if plots similar to what is in Figure 3d could be added for each of these events in the supplement (14 plots for OA + 11 for RA), that would make these findings even more actionable.\nResponse: Thank you for this insightful suggestion. We agree that visualizing the underlying splice events would greatly enhance the interpretability of the colocalization results. In response, we have expanded the supplementary materials to include a comprehensive set of figures for all OA- and RA-associated colocalized loci (Supplementary Figures 17,18).\nTaking SPATA33 as an example (Figure 16), for each OA and RA colocalized gene, we now provide a three-part visualization consisting of:\na) A GWAS-sQTL colocalization LocusZoom plot, illustrating the shared association signal between the GWAS phenotype and splicing event.\nb) A Sashimi plot, showing the specific splice junction(s) affected, the direction of effect, and the transcript structures involved. This enables readers to directly observe the splicing alteration corresponding to each intron cluster.\nc) An sQTL effect boxplot, presenting the association between genotypes and normalized intron usage, which offers an orthogonal and often more quantitative view of the splicing effect, complementing the Sashimi plot.\nFigure 16. Example of supplementary figures illustrating detailed information for a GWAS-colocalized splicing event. (a) LocusZoom plot showing colocalization between the GWAS association signal and the sQTL signal at the locus. (b) Sashimi plot illustrating the effect of the GWAS- and sQTL-associated variant on alternative splicing. (c) Boxplot showing the association between variant genotypes and normalized intron usage, providing quantitative support for the splicing pattern observed in the Sashimi plot.\n20. Comment: Could the authors share the coordinates of exon 2 to ensure the reader understands the correct exon (due to differences in annotations, isoforms, etc “exon 2” could mean different things). Is the function of exon 2 known, is this an annotated exon? Similarly I would be curious if anything was previously known about any of the other 7 OA-splice-specific loci.\nResponse: We thank the reviewer for this helpful comment. In the revised version, we now explicitly provide the genomic coordinates of the COL2A1 exon 2 in Figure 5e (based on hg19). This exon is an annotated exon of the canonical COL2A1 transcript in GENCODE v19 genome annotation file.\nCOL2A1 encodes type II collagen, the major fibrillar collagen of hyaline cartilage, which plays a crucial role in cartilage extracellular matrix structure, development and homeostasis. Exon 2 of COL2A1 is a well-characterized developmentally regulated alternative exon that encodes a cysteine-rich domain in the N-propeptide of type II procollagen. Alternative inclusion or skipping of exon 2 gives rise to the IIA (+exon 2) and IIB (−exon 2) isoforms of type II procollagen. Multiple studies have shown that exon 2–containing IIA transcripts are predominantly expressed in chondroprogenitor cells and early chondroblasts, whereas IIB transcripts lacking exon 2 are the major isoform in differentiated chondrocytes15,16, indicating a marked shift in exon 2 usage during chondrogenesis. We have added this well-established biological background to the revised manuscript (Main text, Lines 259-262)\nTo our knowledge, there is currently no definitive in vivo evidence directly demonstrating that changes in COL2A1 exon 2 usage per se cause osteoarthritis onset. However, our data provide several lines of support suggesting that dysregulation of this developmentally controlled exon is related to OA pathogenesis. First, at this locus, the OA risk allele G of rs1793955 significantly increases exon 2 inclusion in cartilage, and this genetically driven splicing effect colocalizes with the OA GWAS association at the COL2A1 locus. Second, when comparing paired degraded and intact cartilage samples from the same individuals, we also observed significantly higher exon 2 inclusion in degraded cartilage. These independent and concordant observations indicate that elevated exon 2 usage is consistently associated with OA-related molecular and tissue-level changes. Together, these findings support the hypothesis that aberrant regulation of COL2A1 exon 2 splicing may contribute to OA susceptibility, although further functional studies will be required to establish causality.\nFor the remaining seven OA-splice-specific loci that colocalize with OA GWAS signals, we systematically searched the literature for prior functional evidence linking the corresponding splicing events to OA. We did not find studies that directly implicate these specific splice events in OA risk. Thus, our work appears to provide the first evidence connecting these splicing changes to OA susceptibility, and we agree that further functional studies will be required to elucidate their precise roles in disease.\n21. Comment: Figure 6A: Could the authors add eQTL values to this plot as well?\nResponse: Thank you for this helpful suggestion. As shown in Figure 17 below, we have now added the heritability enrichment estimates for both RA and height based on cartilage and synovium eQTL annotations (open circles), alongside the corresponding sQTL enrichments (filled circles). This update allows direct comparison between eQTL- and sQTL-based annotations within each tissue and trait.\nThe same updates have been incorporated into Figure 6A in the revised manuscript.\nFigure 17. Heritability enrichment (proportion h2/proportion variant) of human height and RA GWAS in sQTLs and eQTLs of cartilage and synovium, respectively. The x-axis represents GWAS traits, and the y-axis represents enrichment. Error bars present the 95% CIs, calculated as enrichment ± 1.96 × s.e.\n22. Comment: Of the synovium and cartilage sQTLs that colocalize with GWAS hits for OA, RA and height, how many were also sQTLs in GTEx? How many were novel?\nResponse: Thank you for this helpful question. To address it, we compared all synovium and cartilage sQTLs that colocalize with GWAS signals for OA, RA, and height against the GTEx v8 sQTL data.\nFor OA, among the 12 colocalized sQTLs, 5 are also detected as sQTLs in GTEx, whereas 7 appear to be novel and tissue-specific to cartilage or synovium. For RA, among the 6 colocalized sQTLs, 2 are also GTEx sQTLs and 4 are novel. For height, among the 183 colocalized sQTLs, 80 are present in GTEx, while 103 represent novel sQTLs not detected by GTEx.\nThese results indicate that a substantial proportion of disease-relevant sQTLs identified in joint tissues are not captured in GTEx, highlighting the importance of analyzing primary cartilage and synovium to uncover potential tissue-specific regulatory mechanisms underlying OA, RA and height. The analysis results have been added to the revised manuscript (Main text, Lines 248-250, Lines 289-290; Figure 5b, Figure 6d; Supplementary Tables 10-12).\n23. Comment: It would help if the authors could provide a bit more information on individual sQTL splice events for the reader to understand the specific splice changes, and whether each of these spliceQTLs are novel, especially for the sQTLs that are GWAS hits.\nResponse: Thank you for this helpful suggestion. As noted in our response to the comment #19, we agree that providing clearer information on the underlying splice events would substantially improve the interpretability of the GWAS-sQTL colocalizations. Accordingly, we have expanded the supplementary materials to include detailed visualizations for all OA- and RA-associated colocalized splice events as shown in figure 16 above.\n24. Comment: There are a few minor typos/grammatical errors throughout the manuscript and figures captions:\no Line 62: “-related diseases and traits”\no Line 72: “for cartilage or synovium”\no Line 194: missing a period after “selection”\no Line 223: “COLA2A1” should say “COL2A1”.\no Line 280: “in a tissue-specific manner”\no Line 346: “inferred to be sGenes”?\no Figure 2g: typo in x-axis title\nResponse: Thank you for pointing these out. We apologize for the errors in the manuscript and any inconvenience caused. All errors have been corrected, and the manuscript has been carefully revised to improve clarity and accuracy. We appreciate the reviewer’s careful reading and helpful suggestions.\nReviewer #2 (Remarks to the Author):\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts\nResponse: We thank the co-reviewer for their time and thoughtful evaluation of our manuscript. We appreciate the constructive comments and suggestions provided, which have helped us improve the clarity and quality of the work.\nReferences:\nReviewer #3:\nThis manuscript presents a high-quality and comprehensive analysis of splicing quantitative trait loci (sQTLs) in cartilage and synovium—two key joint tissues implicated in osteoarthritis (OA), rheumatoid arthritis (RA), and human height. Using RNA-seq and genotype data from 238 donors, the authors identified 4,133 independent cis-sQTLs and three trans-sGenes, including 309 tissue-specific sQTLs. Integration with GWAS data revealed novel effector genes for OA (14 genes), RA (11 genes), and height (276 genes), many of which were not detectable via eQTL analysis. The study also employed SpliceAI to fine-map 207 functional splice-altering variants, highlighting both canonical and distal regulatory mechanisms. These findings provide a valuable resource and offer new insights into the tissue-specific genetic regulation of alternative splicing and its role in complex trait etiology.\nOverall, the manuscript is well-conceived, methodologically rigorous, and biologically meaningful. It represents a significant advance in the field of post-transcriptional regulation and joint disease genetics. The integration of sQTLs with GWAS data is particularly impactful, revealing regulatory mechanisms that are not captured by gene expression alone.\nMajor Strengths\n• Novel Resource: First large-scale sQTL map for cartilage and synovium, filling a critical gap in joint tissue genomics.\n• Methodological Rigor: Use of advanced computational tools (LeafCutter, tensorQTL, SpliceAI, mash, coloc) and robust statistical frameworks.\n• Biological Insight: Identification of tissue-specific splicing events and their links to disease-associated loci provides mechanistic depth.\n• Data Integration: Effective combination of sQTL and GWAS data enables discovery of effector genes not captured by eQTL analysis.\nResponse: We sincerely thank the reviewer for the positive evaluation and recognition of our work. We are pleased that the reviewer found the construction of the large-scale sQTL map for cartilage and synovium, the methodological rigor of our analyses, and the biological insights gained from integrating sQTLs with GWAS data to be valuable and impactful. We greatly appreciate these encouraging comments, which reinforce the significance of our study and motivate our continued efforts to advance the understanding of splicing regulation in joint diseases.\nMinor Points\n1. Comment: SpliceAI Scope: SpliceAI predictions are limited to variants within annotated genes, potentially missing regulatory variants in intergenic or enhancer regions.\nResponse: Thank you for pointing this out. We agree that SpliceAI predictions are restricted to variants located within annotated genes, and therefore may miss functional splice-altering variants in intergenic or enhancer regions. In the revised manuscript, we have discussed candidate alternative regulatory mechanisms for functional sVariants located in other genomic regions, such as disruption of RNA-binding protein (RBP) binding and perturbation of enhancer activity. We also provided two systematic annotation tables for all sVariants within fine-mapped 95% credible sets to show whether they are predicted to affect RBP binding or overlap candidate enhancer regions. In addition, we have added a statement in the Discussion noting that future work incorporating advanced computational algorithms and high-throughput experimental approaches may provide a more systematic framework for deciphering the functional mechanisms of sVariants whose regulatory roles were not fully resolved in this study (Main text, Lines 328-340; Supplementary Tables 13,14).\n2. Comment: Trans-sQTL Power: Only three trans-sGenes were identified, likely due to limited sample size. This aspect could be expanded in future studies with larger cohorts.\nResponse: We appreciate the reviewer’s observation. Indeed, the identification of only a limited number of trans-sGenes is likely a consequence of limited sample size and the inherently modest effect sizes of trans-regulatory signals. We now explicitly discuss this limitation in the revised Discussion and emphasize that larger cohorts will be needed to more comprehensively characterize trans-acting splicing regulation in joint tissues (Main text, 360-362).\n3. Comment: Visualization Density: Some figures are dense and could benefit from clearer annotations and legends to improve accessibility for readers.\nResponse: Thank you for bringing this to our attention. We have carefully reviewed all figures and legends and revised them for clarity, including improving annotation density and adding explanatory labels where needed. We hope these refinements enhance readability and accessibility for a broader audience.\n4. Comment: Terminology Consistency: Ensure consistent use of terms such as “sVariants,” “sGenes,” and “effector genes” throughout the manuscript.\nResponse: We appreciate this helpful suggestion. We conducted a thorough review of the manuscript to ensure consistent use of terms such as “sVariants,” “sGenes,” “effector genes,” and related nomenclature across all sections, figures, and supplementary materials. These terms have now been standardized throughout the revised version.\nReviewer #4:\nThis manuscript is a complementary study to a recent publication by the same group, which used what appear to be largely the same samples for an eQTL analysis. In the present work, the authors perform a comprehensive splicing QTL (sQTL) analysis on more than 200 samples from knee synovium and cartilage obtained from patients undergoing knee replacement surgery.\nFrom the standpoint of general novelty, there is limited conceptual innovation, as sQTLs have been recognized since 2013 as a major mechanism through which human genetic variants and particularly those identified in GWAS affect gene regulation. That said, the quality of this study is superb. The tools used are standard but well-implemented, the analyses are rigorous, and the overall presentation is of high quality. This paper will be of interest to researchers studying rheumatoid arthritis, inflammatory joint disorders, and other autoimmune conditions. The data quality appears excellent, the illustrative examples are compelling, and the lists of significant sQTLs provided in the supplementary tables will be valuable to the community.\nWe have no major comments, mostly specific comments and clarifications:\n1. Comment: It is unclear why the number of samples used in this study (238) differs from the previous eQTL study (243), perhaps explaining the sample retention criteria would help especially if it differed from their previous study.\nResponse: Thank you for this important feedback, and we apologize for the lack of clarity in the original description of sample numbers. We would like to clarify that the 238 samples refer to all qualified donors included in the present study after quality control, each with matched RNA-seq and genotype data. Among these, 204 donors had high-quality, matched cartilage RNA-seq and genotype data, consistent with the sample size reported in our previous cartilage eQTL study, and 202 donors had high-quality, matched synovium RNA-seq and genotype data. 168 donors have both qualified cartilage and synovium RNA-seq data (as shown in Figure 18a below).\nThe number 243 refers to the number of donors included in our previous synovium eQTL study after genotype data quality control but prior to RNA-seq data quality control (Figure 18b). Of these, 41 samples were further excluded due to either missing matched synovium RNA-seq data or RNA-seq data failed quality control, and 202 samples were retained for downstream analyses, which corresponds exactly to the synovium samples analyzed in the current study. Therefore, the samples used here are fully consistent with those in our previous studies, and no additional quality control was applied beyond what has already been described.\nTo minimize potential confusion regarding sample numbers, we have added the Figure 18a below as a new supplementary figure in the revised manuscript and added more detailed sample donor description in Methods section to provide a clear overview of sample information (Main text, Lines 384-389; Supplementary Figure 1a).\nFigure 18. Clarification of sample sizes. (a) Venn diagram showing the number of donors with high-quality, matched cartilage and synovium RNA-seq and genotype data included in this study. (b) Overview of donor inclusion and exclusion at each quality control step in the previous synovium eQTL study, indicating the number of donors retained or removed due to missing synovium RNA-seq data or failing RNA-seq or genotype quality control, resulting in 202 donors with matched genotype and synovium RNA-seq data.\n2. Comment: Figure 2e has a typo in the x axis label. Please change “sVraint” to “sVariant”.\nResponse: Thank you for pointing this out. We apologize for the erroors in the manuscript and any inconvenience caused. All errors have been corrected, and the manuscript has been carefully revised to improve clarity and accuracy. We appreciate the reviewer’s careful reading and helpful suggestions.\n3. Comment: Figure 3a: Please clarify what the blue and red distributions represent. I assume that red corresponds to sVariants and blue to eVariants, but this should be explicitly stated. They are specified in figure 6.\nResponse: Thank you for your thoughtful comment. We apologize for the missing figure legend and the inconvenience caused. We have added the legend for Figure 3a in the revised version.\n4. Comment: SpliceAI threshold: The rationale for using a 0.3 cutoff for SpliceAI is unclear. Was this based on prior benchmarking or empirically determined in this dataset? Further, calling SpliceAI overlap “cross-validation” is misleading, cross validation has a very specific meaning in machine learning.\nResponse: Thank you for this insightful comment. We apologize for the earlier lack of clarity. In our analysis, we used a SpliceAI cutoff of 0.2, which was selected based on the benchmarking reported in the original SpliceAI publication1. In that study, Jaganathan et al. showed that a threshold of 0.2 achieves a high sensitivity of approximately 80%, although the validation rate at this cutoff is relatively modest (~30%). Because our primary objective was to comprehensively capture potential splice-altering variants, we prioritized high sensitivity to avoid missing functionally relevant candidates. Importantly, the relatively lower validation rate associated with this threshold can be effectively compensated by integrating our sQTL data, which provides orthogonal genetic evidence for splice-altering effects. By leveraging the complementarity between computational predictions (high sensitivity) and statistical association signals (high specificity), this combined strategy enables us to identify as many bona fide functional variants as possible, while reducing false positives that may arise from SpliceAI alone.\nWe also agree with the reviewer that referring to this overlap as “cross-validation” was inappropriate, as the term has a specific technical meaning in machine learning. We have revised the manuscript accordingly and now describe the integration as a complementary validation strategy rather than cross-validation.\n5. Comment: Unannotated splice events: One potential missed opportunity was the identification of de novo or previously unannotated splice junctions arising from genetic variants. Did the analysis pipeline allow for detection of unannotated intron-exon junctions, or was it restricted to reference-annotated events?\nResponse: Thank you for this valuable comment. Our analysis pipeline does allow for the discovery of previously unannotated splice junctions. Specifically, we used LeafCutter for splicing quantification, which is annotation-independent and identifies splice events directly from the genomic coordinates of junction-spanning reads. As a result, our workflow is capable of detecting de novo or unannotated splice junctions that are not present in existing gene models.\nTo assess this explicitly, we compared all significant introns identified in cartilage and synovium with the GENCODE v19 annotation. We found that 17% of significant introns in cartilage and 19% in synovium correspond to unannotated splice junctions (Figure 19). These results highlight the ability of our approach to uncover novel splicing events that would be missed by annotation-restricted methods. These results have been added to the revised manuscript (Main text, Lines 98-104; Supplementary Figure 4).\nFigure 19. Classification of introns according to whether both intron boundaries, one intron boundary, or neither intron boundary are annotated in GENCODE v19.\n6. Comment: Enrichment analysis: Please clarify how enrichment was computed. Was the background a random set of all genes, or a random subset of expressed genes matched for tissue expression level? This distinction is important for interpreting enrichment magnitudes.\nResponse: Thank you for raising this important point. We have clarified the definition of the background sets used in each enrichment analysis. For the tissue-specific sQTL gene enrichment, the background consisted of all genes included in the tissue-specific sQTL analysis, ensuring that both the test set and background were drawn from the same pool of genes tested for tissue specificity.\nFor the height GWAS-colocalized sQTL gene enrichment, the background set included all cartilage genes that were eligible for the colocalization analysis, thereby matching the background to the set of genes with sufficient expression and genetic signal for colocalization testing.\nThese clarifications have been added to the Methods section to make the enrichment framework fully transparent and interpretable (Main text, Lines 559-562).\n7. Comment: Terminology (“synovium-specific”): Consider avoiding the term synovium-specific for genes detected only in synovium but not in cartilage. For instance, IRF5 was only observed in synovium here, yet it is a well-known transcription factor induced in multiple inflammatory contexts. “Detected only in synovium” or “preferentially expressed in synovium” might be less misleading than “specific.”\nResponse: Thank you for your helpful suggestions. As suggested, we have revised all expressions like “synovium-specific sGene” to “sGene only detected in synovium” to make them less misleading.\n8. Comment: Figure 4i and distal sVariants: I assume that the example shown in Figure 4i is an illustrative case rather than a general trend, and that most distal sVariants were not located within clearly defined or shared DNA-accessible or transcriptionally active regions. Please clarify whether this was generally observed. Which fraction of lead or fine mapped sQTLs reside within accessible regions.\nResponse: Thank you for this thoughtful comment. We agree that the locus shown in Fig. 4i is intended as an illustrative example; however, the observation that distal sVariants overlap accessible/active regulatory regions is not limited to a few isolated cases. To clarify this point, we performed a systematic epigenomic annotation analysis and added the results as new supplementary materials (Main text, Lines 162-168; Supplementary Figure 12; Supplementary Tables 13,14).\nFirst, using lead sVariants, we observed significant enrichment in open and active chromatin marks, including ATAC-seq peaks and H3K27ac peaks (Fig. 20a), supporting the idea that many lead sQTL variants fall within regulatory regions that are accessible and/or transcriptionally active.\nSecond, we fine-mapped all sQTL signals and annotated variants in the 95% credible sets. We found that more than two-thirds of sQTL signals had at least one credible-set variant overlapping ATAC-seq and/or H3K27ac peaks (Fig. 20b). Moreover, approximately one-third of sQTL signals had at least one credible-set variant supported by both ATAC-seq and H3K27ac annotations (Fig. 20b). Together, these results indicate that overlap with accessible/active regulatory regions is a general pattern, rather than being driven by a small number of loci.\nFinally, to address the reviewer’s point on “distal” sVariants, we examined the distance between annotated sVariants and their target splice sites. We observed that >20% of annotated sVariants are located >50 kb from the target splice site (Fig. 20c), indicating that distal regulatory sVariants are also a frequent feature of our data rather than rare exceptions.\nFigure 20. Epigenetic annotation of sVariants. (a) The enrichment of lead sVariants across different epigenetic markers. Dots represent odds ratios, and the whiskers represent 95% CIs. (b) Bar plot showing epigenetic annotation results of fine-mapped sVariants. (c) Histogram showing the distribution of distances between the enhancer marker annotated credible set sVariants and their target splice sites.\n9. Comment: Negative control (“mental disorders”): Is “mental disorders” an appropriate negative control category, given the vastly different expression profiles between synovium/cartilage and brain? If the goal was to include a trait set with minimal overlap in expression context.\nResponse: Thank you for raising this point. We would like to clarify that stratified LD score regression (S-LDSC) is used to identify functional genomic regions in which GWAS heritability is enriched, by testing whether disease heritability is preferentially associated with specific molecular annotations. When S-LDSC is applied across tissue- or cell-type–specific regulatory features, the results show that disease heritability tends to be enriched in functional elements derived from established disease-relevant tissues or cell types, while little or no enrichment is observed in tissue-specific regulatory annotations from etiologically unrelated tissues2-8. Notably, multiple studies have consistently reported that the heritability of mental disorders is predominantly enriched in molecular features derived from neural or brain-related cell types, with minimal enrichment in non-neural tissues2,5,8.\nBased on this rationale, we selected mental disorders as a negative control under the assumption that cartilage and synovium are not causally involved in the pathogenesis of mental disorders, rather than on differences in gene expression profiles between joint tissues and brain, similar to the negative control selection strategy by other studies5,7. Thus, we would not expect the heritability of psychiatric traits to be enriched in cartilage- or synovium-derived eQTL or sQTL annotations, making mental disorder GWAS an appropriate biological negative control in this context.\n10. Comment: Power and saturation: Do the authors have an estimate of how far from saturation the current sQTL discovery is? In other words, how many additional samples would be needed to approach saturation in sQTL detection for these tissues? Such estimates would help contextualize the study’s power.\nResponse: Thank you for this helpful suggestion. In the revised manuscript, we have performed a formal power analysis to assess the current degree of sQTL discovery saturation and to estimate the sample sizes required to detect additional signals.\nBased on this analysis, with our current sample size (~200 donors), the statistical power is already close to 1.0 for variants with MAF > 0.15, and approximately 0.8 for variants with MAF around 0.10. In contrast, power is substantially lower for variants with MAF < 0.05, indicating that sQTL discovery for low-frequency variants is not yet saturated. The analysis further shows that increasing the sample size to ~400 donors would raise power for MAF around 0.05 variants to approximately 0.8, whereas detecting sQTLs driven by rarer variants (MAF < 0.01) would require >1,000 samples to achieve comparable power.\nFigure 21. Power analysis of cis-sQTL showing statistical power to detect cis-sQTL across different MAF and sample sizes. We assume α = 0.05/2,500,000, reflecting Bonferroni correction for 2,500,000 hypotheses based on 50,000 intron clusters and an average of 50 non-redundant SNPs in the region ±1 Mb of each intron cluster.\nThese results indicate that our current dataset is well powered for common sQTLs, but that additional samples would primarily increase discovery for low-frequency and rare variants, rather than substantially altering conclusions for common variants. We have added these power estimates to the revised manuscript to better contextualize the study’s power and to provide guidance for future sQTL studies in synovium and cartilage (Main text, Lines 94-96; Supplementary Figure 3).\n11. Comment: Data reuse and new datasets: The HRA004624 dataset was already used in the authors’ previous eQTL study (Nat. Commun., 2024; [DOI: 10.1038/s41467-024-45652-x]*). However, dataset HRA006224 was not available from GSA at the time of review, making it unclear whether it constitutes new data. The authors should clarify whether this study includes any non-overlapping samples relative to the prior eQTL analysis and what is the HRA006224 contributing to this study since it seems that all data was already available from the previous study.\nResponse: We thank the reviewer for raising this important point. We confirm that the current study uses exactly the same set of samples as our previous eQTL publication, and no additional or non-overlapping samples were included. The dataset HRA006224 corresponds to the same sequencing data and was generated as part of the same cohort; it has now been released and is publicly accessible through the GSA repository at the time of revision (https://ngdc.cncb.ac.cn/gsa-human/browse/HRA006224). We have clarified this in the revised manuscript to avoid any misunderstanding regarding data reuse and availability.\n12. Comment: While new methods may not have been developed. We still think the methods reported lack in detail, for example the exact pipeline ran (with software versions and parameters) would help readers be able to replicate their study.\nResponse: Thank you for this helpful comment. We agree that providing sufficient methodological detail is essential for reproducibility. In the revised manuscript, In the revised manuscript, we have substantially expanded the Methods section to include a clearer and more comprehensive description of the analysis pipeline, including the software tools used, their versions, key parameters, and data-processing steps. These additions ensure that readers can accurately replicate our workflow and analytical results.\nReviewer #5:\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.\nResponse: We thank the co-reviewer for their time and thoughtful evaluation of our manuscript. We appreciate the constructive comments and suggestions provided, which have helped us improve the clarity and quality of the work.\nReferences:\n\nResponse to reviewers’ comments:\nReviewer #1:\nWe have reviewed the revised manuscript and response to our comments, and are impressed with the authors' thorough and thoughtful responses to all reviewer comments. We are impressed with the revisions and clarifications provided and have no further questions or requests.\n1. Comment: We recommend (if possible, not required) that the authors inquire whether the results from this study could be made available through the GTEx portal.\nResponse: We sincerely thank the Reviewer for the positive evaluation of our revised manuscript and for the helpful suggestion regarding data sharing.\nWe have inquired about the possibility of depositing our results in the GTEx portal; however, the GTEx portal does not accept externally generated datasets. Accordingly, we have deposited the summary-level results from this study in the OMIX database under accession number OMIX015281, as indicated in the Data Availability section of the revised manuscript.\nReviewer #4:\nThe authors have addressed most of my initial minor comments and in response to reviewer’s 1 points have improved their presentation of their analysis. I have 2 small, suggestions\n2. Comment: Although the authors do not explicitly quantify the number of sGenes harboring multiple independent sQTLs, Fig. 2a suggests a markedly higher proportion of multi-sQTL genes in cartilage compared to synovium. A rough estimation from the figure indicates that approximately 25–30% of cartilage sGenes harbor two or more independent sQTLs, versus roughly 10–15% in synovium (an odds ratio of 2!). A simple contingency analysis based on these approximate counts yields a highly significant difference (χ² test, p < 1015), consistent with the similar total number of sQTL signals despite substantially fewer sGenes in cartilage. The authors do not comment on this, is it a result of some technical artefact or do the authors think there is a greater number of genes with a more complex splice regulation in the cartilage?\nResponse: We thank the Reviewer for this careful observation. We agree that the proportion of genes harboring multiple independent sQTL signals is an important feature of splicing regulation.\nUpon explicit quantification at the intron-cluster level, we found that among the 960 significant sQTL intron clusters identified in cartilage, 33 harbored more than one independent sQTL signal (3.4%). In contrast, among the 1,680 significant sQTL intron clusters identified in synovium, 102 harbored multiple independent signals (6.1%). Thus, the proportion of multi-signal intron clusters is approximately 1.8-fold higher in synovium than in cartilage. This pattern is consistent with previous large-scale studies, including GTEx and other recent analyses1,2, which have shown that as the total number of detected sQTL genes increases, the fraction of genes harboring multiple independent sQTL signals also rises. This reflects increased statistical power to resolve secondary signals rather than intrinsic differences in regulatory complexity.\nTo improve clarity for readers, we have updated Figure 2a such that the left y-axis now represents the proportion of intron clusters harboring multiple independent sQTL signals rather than the absolute counts.\n3. Comment: Is the case displayed in figure 2h,i general, meaning that lead sQTL SNPs frequently differ between synovium and cartilage for the same regulatory region suggests that the lead SNP is not, in general, the causal variant but rather a tissue-specific proxy within an LD block. This presents an opportunity to further refine causal inference by leveraging cross-tissue differences. For example, in cases where the synovium lead SNP lies closer to a splice junction while the cartilage lead SNP is more distal (or vice versa), functional proximity could help prioritize the more likely causal variant within the shared LD region. Given the concerns raised by Reviewer 1 regarding the distance of causal sQTL variants from splice sites, such an analysis could help re-estimate the fraction of sQTLs whose likely causal variants lie near splice regions versus more distal regulatory elements. This may either strengthen the conclusion that most causal variants act locally on splicing, or alternatively provide additional evidence that a substantial fraction act through more distal mechanisms.\nResponse: We thank the Reviewer for this thoughtful suggestion. We agree that, for shared sQTL events, tissue-specific lead variants may often represent proxies within the same LD block rather than the true causal variant, and that cross-tissue differences could, in principle, provide additional resolution for causal inference.\nTo further evaluate this point, we performed a splice event-level assessment across n = 1141 shared significant sQTL events between cartilage and synovium. As shown in Figure 1, the lead sVariant was identical in 26.5% (302/1141) of events, whereas 73.5% (839/1141) showed different lead variants across tissues. Among these discordant lead pairs, 78.3% (657/839) were in strong LD (R² > 0.8), consistent with tissue-specific proxy variants within a shared LD block rather than independent signals. For the strong-LD discordant events, we compared the distance of each tissue’s lead variant to the corresponding target splice site. We observed frequent differences in proximity, with a median absolute distance difference of 13,207 bp; 53.0% (348/657) of events had the synovium lead variant closer to the target splice site, whereas 47.0% (309/657) had the cartilage lead variant closer. These findings indicate that cross-tissue lead discordance is common and may provide additional information within shared LD regions.\nHowever, because proximity to the target splice site alone does not establish causality, we did not further prioritize variants solely based on distance. Instead, we have incorporated additional discussion in the revised manuscript highlighting the potential of cross-tissue fine-mapping approaches to refine causal variant inference in future studies (Main text, Lines 333-335).\nFigure 1. Cross-tissue comparison of lead sVariants for shared sQTL events. a, Distribution of linkage disequilibrium (LD R²) between lead sVariants of shared significant splice events identified in cartilage and synovium. Each bar represents the number of shared splice events within a given LD bin. b, Scatter plot showing the genomic distance between the lead sVariant and the target splice site for each shared splice event in cartilage (x-axis) and synovium (y-axis). Distances are shown in kilobases (kb). Each point represents one shared splice event. Blue points indicate events in which the cartilage lead sVariant is closer to the target splice site, whereas red points indicate events in which the synovium lead sVariant is closer. The dashed diagonal line represents equal distance in the two tissues.\nReferences:\n1. Consortium, G.T. The GTEx Consortium atlas of genetic regulatory effects across human tissues. Science 369, 1318-1330 (2020).\n2. Tian, C. et al. Single-cell RNA sequencing of peripheral blood links cell-type-specific regulation of splicing to autoimmune and inflammatory diseases. Nature Genetics (2024)."}]

Decision letter