Background: Rheumatoid arthritis (RA) is a clinically heterogeneous autoimmune disease. Our recent work in early, treatment-naïve RA identified joint involvement patterns (JIP) as key discriminators, defining four subsets based on joint counts and anatomical regions (JIP-feet, JIP-oligo, JIP-hand, JIP-poly) [1]. These subsets differed in treatment response and synovial histology: JIP-hand responded best to treatment and showed classical synovial inflammation, whereas JIP-oligo responded less well to methotrexate. However, genetic variants characterizing these clinically defined JIP subclusters have not yet been established. In early RA, tender–swollen joint difference (TSJD) has been reported as a clinical marker associated with prognosis [2]. Given that the subclusters are defined by joint involvement patterns, genetic variants distinguishing these subclusters may also relate to TSJD, contributing to a better understanding of heterogeneity in RA.

Objectives: To investigate genetic associations between clinically defined early RA subclusters using a penalized regression framework, and to assess whether associated variants are related to TSJD.

Methods: We analyzed two different treatment-naïve RA cohorts from Leiden University Medical Center. The discovery cohort was collected up to 2020, and the replication cohort was collected between 2021 and 2024. The discovery (n = 771) and replication (n = 239) cohorts included all four JIP subsets (foot/oligo/hand/poly: 166/226/188/191 and 70/53/62/54, respectively). Genotype data were imputed using the 1000 Genomes Phase 3 reference panel. Pairwise genetic analyses were performed across all subcluster combinations. We applied binomial logistic regression with elastic net penalization using RA GWAS–derived lead SNPs [3] as predictors. Model sparsity was determined using deviance-based 10-fold cross-validation (CV). The number of SNPs included in the model was defined using a CV–based criterion (λ.1se) favoring a parsimonious solution [4]. Stability selection with iterative subsampling (80% of samples, 1,000 iterations) was used to rank SNPs. The top-ranked SNPs corresponding to the number defined by λ.1se were retained for refitting unpenalized logistic regression models. Fixed-effects inverse-variance weighted meta-analysis was performed to obtain pooled effect estimates and assess between-cohort heterogeneity. For variants associated with the comparison between two JIP subclusters, genotype–TSJD trends were assessed using the Jonckheere–Terpstra trend test.

Results: The strongest genetic association was observed for the JIP-oligo versus JIP-hand comparison, captured by three SNPs selected through cross-validation (Figure 1A). Accordingly, subsequent analyses focused on the JIP-oligo versus JIP-hand comparison. Stability selection identified three SNPs for refitting. (Figure 1B).The refitted model including these three SNPs was significant in both cohorts (discovery likelihood ratio test p = 0.0028 < 0.0083, Bonferroni-corrected for six pairwise tests; replication p = 0.035). Among the three variants, rs4409785 showed a significant association with consistent direction of effect in both the discovery (OR = 1.49, p = 0.048) and replication cohorts (OR = 2.11, p = 0.044). In contrast, rs6495979 and rs2300373 were associated in the discovery cohort (OR = 0.71/1.66, p = 0.034/0.030, respectively) but not replicated (OR = 1.57/0.44, p = 0.185/0.109), with inconsistent effect directions. Meta-analysis showed a significant pooled effect for rs4409785 (OR = 1.61, p = 0.007) with consistent effects across cohorts and no evidence of heterogeneity (I ² = 0%) (Figure 2A). No significant genotype–TSJD trend was observed when JIP-oligo and JIP-hand were analyzed together (p = 0.631; Figure 2B, top). In contrast, subcluster-stratified analyses revealed opposite trends, with rs4409785 associated with a significant decrease in TSJD>0 in JIP-oligo (p = 0.044) and a consistent increase in JIP-hand (p = 0.093) (Figure 2B, middle and bottom).

Conclusions: Pairwise genetic analyses identified a reproducible association of rs4409785 with JIP-oligo versus JIP-hand, driven primarily by rs4409785.This regulatory variant has been implicated in chromatin architecture and SESN3 regulation in memory T cells [5]. Opposite TSJD associations across JIP subclusters suggest endotype-dependent clinical translation of genetic risk. Notably, divergent TSJD trends were observed in the stromal-rich JIP-hand subset compared with JIP-oligo, suggesting that tissue context may influence downstream clinical effects of genetic variation. Our current findings merit a larger replication in independent cohorts of patients with early RA.

Figure 1.

SNP selection for pairwise genetic analysis (JIP-oligo vs JIP-hand)

(A) Deviance-based cross-validation curve from penalized logistic regression comparing JIP-oligo and JIP-hand in the discovery cohort. Vertical dashed lines indicate the value of λ that minimizes the deviance (λ.min; left) and the value within one standard error of the minimum (λ.1se; right). (B) Stability selection results for the JIP-oligo vs JIP-hand comparison.

Figure 2.

Genetic association of rs4409785 with JIP-oligo versus JIP-hand and subcluster-specific associations with TSJD

(A) Forest plot showing cohort-specific odds ratios (ORs) and 95% confidence intervals for rs4409785, and the pooled fixed-effects meta-analysis estimate. Heterogeneity is quantified using I ² . (B) Distribution of TSJD categories, defined as previously reported [2], across rs4409785 genotype dosage within JIP subclusters. Stacked bars show the proportion of TSJD categories within each genotype group for JIP-oligo + JIP-hand (top), JIP-oligo (middle), and JIP-hand (bottom).

REFERENCES: [1] NPJ Digit Med. 2025 23;8:623.

[2] ACR Open Rheumatol. 2024 6;6:347–355.

[3] Nat Genet. 2022;54:1640-1651

[4] J Cheminform. 2014 29;6:10.

[5] Genome Biol. 202510;26:26.

Acknowledgments: NIL.

Disclosure of Interests: None declared.