Whole Exome Sequencing and Analysis of T Cell Receptor β & γ Repertoires in Rheumatoid Arthritis

1. Introduction

Rheumatoid arthritis (RA) is a major systemic autoimmune disease characterized by chronic inflammation of the synovial joints and bony destruction, resulting in a wide spectrum of symptoms such as polyarticular pain and swelling [1]. RA is associated with progressive damage, systemic complications, mortality, and morbidity, and imposes a significant socioeconomic burden [2]. Except for a few certain cases and populations, the prevalence of RA is generally constant at around 0.5 to 1.0% with tens of millions of patients worldwide [3].

RA has many complex etiologies involving genetic, immunological, and acquired factors [4]. Previous studies have reported that genetic background contributes to approximately 50% of patients with RA [5,6]. In addition, environmental factors, such as smoking, respiratory diseases, and changes in the gut microbiome, may increase the risk or cause epigenetic modifications in susceptible patients, which may lead to RA development [2,3]. As RA shows complex traits and rapid progression, identifying risk factors and preventing disease progression are important [1-4]. For these reasons, identifying genes for RA is very important because it can be used as basic data in medicine as well as in public healthcare, which can be applied to further understand the pathogenesis of RA and potential targets of RA therapy [7,8]. As part of the RA research, many studies have been conducted on human leukocyte antigen (HLA)-related alleles [1,2,9]. Many studies related to T cells have been conducted in RA. While existing studies have focused on the abnormal expression of T cells and the resulting autoantibody production, recent studies have focused on immunotolerance such as programmed cell death 1 (PD-1) [10,11,12]. Although many studies have been conducted on susceptibility to RA development using single nucleotide polymorphisms (SNPs) and have attempted to identify disease-specific genes or polymorphisms in RA, these SNPs alone cannot accurately explain the etiology because the RA-related genetic background is complex [13]. Genome-wide association studies (GWAS), suggested genetic loci for susceptibility to RA [14,15,16,17]. With the advent of next-generation sequencing (NGS), some potential genetic variants have been analyzed [18]. Abnormal and pathogenic T-cell responses that evade normal immune functions could be considered one of the mechanisms of RA development [2], and the T cell receptor (TCR) repertoire has been studied in patients with RA [19]. However, to date, there have been only a limited number of related studies, and there have been no published studies on Korean patients.

For these reasons, this study aimed to investigate the potential genetic variants for RA development in Korean patients using whole exome sequencing (WES) and the pathogenesis of disease course using T cell receptor (TCR) β & γ repertoire (TRB/TRG) analysis. This study aimed to provide basic data that could be applied to Korean patients by analyzing and evaluating results at the molecular genetic level.

2. Materials and Methods

2.1. Study design and data collection

From December 2020 to February 2022, study participants were recruited from patients who visited the outpatient department (OPD) of the Rheumatology of Internal Medicine at Wonju Severance Christian Hospital, a tertiary university-affiliated hospital located in Wonju, South Korea. The inclusion criterion of this study were naïve patients who were first diagnosed with RA in our hospital and those who have not yet started the disease-modifying antirheumatic drug (DMARD) treatment. A total of 17 patients were enrolled and scheduled to visit the OPD at various time points: the first time (before the initiation of DMARD treatment), and 6 and 12 months after treatment. However, three patients were lost to follow-up during the study period, and finally, 14 patients remained, and their blood samples were collected for up to 12 months and used for data analysis. Five healthy controls (HCs) without diagnosed RA were also enrolled for comparison, requiring only one visit. At each visit, blood samples were collected in two K₂-ethylenediaminetetraacetic acid (EDTA) tubes, one for measuring the erythrocyte sedimentation rate (ESR) and the other for genomic testing, and one serum separation tube (SST) for measuring other laboratory data.

Baseline characteristics were collected as follows: age, sex, medical history of hypertension and diabetes mellitus, height, body weight, and the duration before the first visit for clinical information; tender joint count of 44/28 joints (TJC44/28), swollen joint count of 44/28 joints (SJC44/28), disease activity score in 28 joints (DAS28), health assessment questionnaire (HAQ), visual analogue scale (VAS) score, simplified disease activity (SDAI), clinical disease activity index (CDAI), fulfillment of the 1987 American College of Rheumatology (ACR) criteria (whether ≥4 out of 7 cases are met) [20] and the 2010 ACR/European Alliance of Associations for Rheumatology (EULAR) criteria (whether a score of ≥6 out of 10 was met) [21] for clinical classification or scoring system.

Laboratory data were collected on the following: anti-nuclear antibody (ANA) determined by the QUANTA Lite^® ANA assay (Inova Diagnostics, Inc., San Diego, CA, USA), rheumatoid factor (RF) and C-reactive protein (CRP) were determined by the Cobas c 702 automated analyzer (Roche Diagnostics, Rotkreuz, Switzerland), anti-cyclic citrullinated peptide (ACCP) determined by the QUANTA^® Flash CCP3 assay (Inova Diagnostics), and ESR determined by the TEST-1 analyzer (SIRE Analytical Systems, Udine, Italy).

This study was approved by the Institutional Review Board (IRB) of Wonju Severance Christian Hospital (IRB No. CR320086). All the participants voluntarily participated in the study and provided written informed consent.

2.2. Genomic analysis

Peripheral blood mononuclear cells (PBMCs) were isolated from EDTA-treated whole blood by Ficoll-Hypaque density-gradient centrifugation. The detailed protocols for WES and TRB/TRG were as follows.

For WES, genomic deoxyribonucleic acid (DNA) was extracted by using the ChemagicTM Magnetic Separation Module I (MSM I) method (PerkinElmer Chemagen, Baesweiler, Germany) with the DNA blood 200 μL kit. The MGIEasy Exome Capture V5 Probe Set (MGI Tech Co., Ltd., Shenzhen, China) was used for library preparation, and sequencing was performed on the MGI DNBSEQ-G400 platform (MGI Tech Co., Ltd.) to generate 2 × 100 bp paired-end reads. DNA sequence reads were aligned to the reference sequence based on the public human genome build GRCh37/UCSC hg19. Alignments were performed with BWA-mem (version 0.7.17), duplicate reads were marked with biobambam2 and base quality recalibration, variant calling was performed with the Genome Analysis Tool kit (GATK, version 4.1.8), and annotation was performed with VEP101 (Variant Effect Predictor 101) and dbNSFP v4.1 (Liu et al, 2020).

For the TCR/TRB analysis, genomic DNA was extracted using the Chemagic^TM DNA Blood 200 kit (Chemagen, Baesweiler, Germany). For TCR repertoire sequencing, the Lymphotrack^® TRB and TRG assay (Cat No. 72270009 and 72270019, respectively) (Invivoscribe Inc., San Diego, CA, USA) for NGS were used. These assays amplify the DNA between primers that target the conserved V and J regions of the T-cell receptor genes, TRB and TRG. Polymerase chain reaction (PCR) amplicons were purified and quantified using a Qubit fluorometer (Thermo Fisher Scientific). Equimolar amounts of libraries were pooled and loaded onto a flow cell using a MiSeq Reagent kit v2 (500 cycles) (Illumina Inc., San Diego, CA, USA) and sequenced on a MiSeq instrument (Illumina).

2.3. Statistical analysis

The results of the numerical data were checked for normality using the Kolmogorov–Smirnov test, which determined a p-value greater than 0.05 had normality. For parametric data, results are presented as means ± standard deviation (SDs), and the Student’s t-test was used for comparison. For non-parametric data, results are presented as median and interquartile range (IQR), and the Mann–Whitney U test was used for comparison. The results of categorical data were presented as frequencies and percentages (%), and the χ² test was used for comparison.

For patients with RA, the results at baseline (naïve) were compared to those at 6 and 12 months of follow-up. To compare DAS28, DAS28-ESR (=0.70 × In [ESR]) and DAS28-CRP (=0.36 × In [CRP + 1] + 0.96) [22] were calculated. All three VAS scores (pain VAS, global VAS, and physician VAS) were evaluated, but all showed almost the same value; therefore, only the pain VAS scores were used for comparison.

All statistical analyses were performed using Analyse-it version 6.15.4 (Analyse-it Software, Ltd., Leeds, UK), an add-on program to Microsoft Excel 2019 (Microsoft Corp., Redmond, WA, USA), and SPSS (version 25.0; IBM Corp., Armonk, NY, USA). Statistical significance was set at P-value < 0.05.

3. Results

3.1. Baseline characteristics

The baseline characteristics of the 14 patients with RA and 5 HCs are presented in Table 1. The median age of the patients with RA was 57 years (IQR, 45–69), and 12 (85.7%) of were female. There was no statistically significant difference in the proportions of patients with hypertension, diabetes mellitus, or obesity between patients with RA and HCs. When comparing RA criteria, 10 patients (71.4%) satisfied the 1987 ACR criteria, whereas all 14 patients (100%) satisfied the 2010 ACR/EULAR criteria. Laboratory findings showed that the proportions of people positive for ANA, RF, and ACCP were significantly higher in patients with RA than in HCs, and the median ESR and CRP values were also higher in patients with RA than in HCs.

3.2. Genomic variants

Using WES, disease-related variants that have been reported to be associated with RA among several small nucleotide polymorphisms (SNPs) were detected in four (28.6%) patients with RA. Table 2 lists the variants detected using WES that are likely to be related to RA pathogenesis. Most were variants of uncertain significance (VUS), but there was one case of a pathogenic variant (PV).

3.3. Changes in clinical features after DMARD treatment

Table 3 displays the results of comparing the clinical features of naïve RA patients before starting DMARD treatment with the clinical features at 6 and 12 months after the initiation of DMARD treatment. Compared to baseline (before treatment), laboratory findings after 12 months of treatment decreased significantly. There were also changes in the number of joints; TJC44, SJC44, TJC28, and SJC28 counts at 6 months decreased significantly compared to baseline. These facts are also reflected in the disease measures as they are, and DAS28, HQA, VAS, SDAI, and CDAI scores at 6 months decreased significantly at 6 months, and even at 12 months they remained significantly decreased compared to the baseline.

3.4. Changes in TCR diversity after DMARD treatment

TRB and TRG diversities were assessed using the Shannon–Weiner diversity index [23,24]. The results of the comparison of TCR diversity between HCs and patients with RA at baseline and at 6 and 12 months after starting DMARD treatment are illustrated in Figure 1. For TRB diversity, the Shannon–Weiner diversity indices in HCs and patients with RA at baseline, and 6 months and 12 months after treatment were 0.69 (IQR, 0.63–0.74), 0.55 (IQR, 0.41–0.58), 0.56 (IQR, 0.43–0.73), and 0.96 (IQR, 0.66–1.52), respectively (Figure 1(A)). For TRG diversity, the Shannon–Weiner diversity indices in HCs and patients with RA at baseline, and 6 months and 12 months after treatment were 0.42 (IQR, 0.37–0.48), 0.29 (IQR, 0.26–0.31), 0.32 (IQR, 0.24–0.43), and 0.51 (IQR, 0.35–0.76), respectively. Naϊve RA patients showed a relatively lower TCR diversity than HCs, but after the initiation of DMARD treatment, TCR diversity has risen significantly according to the duration of the treatment.

3.5. Relationship between TCR diversity and RA-related factors

The relationship between TCR diversity and disease-related factors in RA was analyzed using Spearman's rank correlation analysis. For laboratory findings, data from patients with RA and HCs were included and analyzed, and only ESR showed a significant correlation with both TRB and TRG. Disease measures showed significant correlations in all items evaluated, especially in DAS28-ESR, SDAI, and CDAI, with correlation coefficient exceeding 0.5 (Table 4).

Additionally, the relationship between TCR diversity (Shannon–Weiner diversity) and the value of each RA-related factor was obtained through linear regression. Figure 2 shows a linear regression diagram between TCR diversity and laboratory values, with the thick line representing the regression line and the dotted line representing the 95% confidence interval. In Table 5, the slope and 95% confidence interval of the linear regression equation between the Shannon–Weiner diversity and RA-related factors were analyzed and presented using data from patients with RA and HCs. Because the analytical measuring interval (AMI) of each test item varies for laboratory items and the scale of each scoring varies for disease measures, it is difficult to compare uniformly using only the slope and intercept. However, it can be seen that, except for ACCP, compared to TRG diversity, it shows a negative slope; therefore, it can be seen that the improvement of laboratory findings and disease measures (decreases in numeric values) is related to the increase in TCR diversity.

4. Discussion

RA is a systemic autoimmune disease that causes chronic inflammatory synovitis, joint destruction, and multi-organ manifestations, resulting in morbidity and mortality [1,4]. RA is an autoimmune disease and many studies have been conducted on autoantibodies; however, a consensus on the crucial autoantigens or genetic and environmental factors that trigger the autoimmune process has not yet been established [25]. The cause of RA is not yet clearly known because of the complex traits of various factors; however, genetic factors account for approximately 50% of RA cases [5,6]. However, it is complicated because it is known that the causes of RA are related to type 1 diabetes, autoimmune thyroid disease, infection, smoking, obesity, and hormones in addition to genetic factors [3]. However, similar to other diseases, many studies have been conducted to identify genetic mutations in RA and identify treatment targets. Owing to the development of GWAS, many genetic mutations related to RA have been discovered [26]. Genetic factors are important because they are related not only to preventive activities for individuals with a high probability of RA but also to the treatment process, response, severity, and prognosis [9]. Current conventional and biological treatments sometimes fail or show a partial response. If the genetic marker is well identified, it is expected that the response to treatment and prognosis can be improved while minimizing toxicity through targeted treatment [25].

For these reasons, this study aimed to find possible variations for the onset of RA in Korean naϊve patients and to evaluate the usefulness of TCR repertoire analysis as a factor related to the disease course and the response to DMARD treatment. Of the 14 patients with RA finally enrolled in this study, all satisfied the 2010 ACR/EULAR criteria, but only 10 patients (71.4%) met the criteria for diagnosing RA according to the 1987 ACR criteria, which is known to have lower sensitivity for early RA diagnosis than the 2010 ACR/EULAR [21].

Based on the WES results, variants expected to be associated with RA were detected in 4 of 14 patients with RA. In the Janus kinase 3 (JAK3) gene, the variant c.1333C>T, p.Arg445Ter was detected, which has been reported as a pathogenic variant in severe combined immunodeficiency [27] but has not yet been reported in RA. However, its association with RA may be considered because JAK3 inhibitors may be used for RA treatment [28]. and follow-up studies are needed. Peptidyl arginine deiminase type IV (PAID4) is a gene that converts arginine residues into citrulline residues in the presence of calcium ions, and a mutation in this gene is known to overproduce citrulline and cause a loss of tolerance, thereby increasing vulnerability to RA [29]. Tumor necrosis factor ligand superfamily member 18 (TNFSF18) is known to activate T cells and B cells in connection with cell signaling [30], and tumor necrosis factor receptor-associated factor (TRAF) is known to be associated with TNF-α and interleukin (IL)-1 and cause inflammation [31]. NF-κB is known to be activated in the synoviums of patients with RA to regulate genes contributing to inflammation, such as TNF-a, IL-6, and IL-8, and treatment therapy targeting this gene is under study [7]. However, it does not match the SNP found in the European [15] and Japanese studies [32]; the clinical significance of the SNP shown in this study is somewhat limited, and it is necessary to clarify the clinical significance of this SNP or to further discover other SNPs through follow-up studies on Korean patients. Of course, it will be helpful to develop a drug or treatment agent that targets gene mutations as they affect the activity of certain proteins, but more GWAS studies will be needed as human genetic patterns are complex, making it difficult to set and establish targets and produce effective therapeutic effects [33].

Genetic and specific environmental factors combine to trigger RA, which activates the immune reaction and causes loss of immune tolerance and increases the number and binding affinity of autoantibodies [1]. The disease progression of RA occurs through epigenetic remodeling [1]. In the immune response, T cells and B cells are activated, resulting in an increase in inflammatory cytokines, production of matrix metalloproteinases, and induction and activation of osteoclasts, resulting in bony destruction [2]. Because the aberrant expression of CD4+ T cells plays an important role in the pathogenesis of RA, the identification and quantitative comparison of CD4+ T cells are being studied in relation to the etiology of RA for early diagnosis and as a possible indicator of the disease course or response to treatment [23]. A previous study reported that the TCR repertoire helps in the early diagnosis of RA [19]. It is known that the TCR repertoire in patients with RA is known to be less diverse than in healthy individuals and related to the disease activity [23,34].

In this study, a TCR repertoire analysis was performed. Naïve patients without prior DMARD treatment were enrolled, and the data of their baseline (before treatment) and 6 and 12 months after the initiation of DMARD treatment were compared. The baseline laboratory findings of patients with RA were significantly different from those of the HCs. The TRB/TRG diversities of patients with RA were lower than those of HCs. However, TRB/TRG diversity in RA patients have increased after the initiation of DMARD treatment. Laboratory values, affected joint counts, and disease measures were significantly decreased. As symptoms and disease-related factors improved according to treatment response, TRB/TRG diversity increased. According to previous studies, as the progression or severity of the disease increases, clinical expansion of T cells occurs, and their diversity of T cells decreases [35]. The diversity of the TCR repertoire may change according to treatment, and monitoring changes in the TCR repertoire is expected to play an important role in monitoring the disease course, evaluating the response to treatment, confirming recurrence, and predicting the prognosis of patients with RA [36]. Correlations between TCR repertoire diversity and disease-related factors were also evaluated. The factors were negatively correlated with TRB/TRG diversity, and this is in line with the result of the previous RA cohort study [23]. The DAS, SDAI, and CDAI scores were more strongly correlated with TRB/TRG diversity than the laboratory findings. The relationships between TCR/TRG diversity and RA-related factors were obtained using linear regression. However, because the reference value and AMI of each item are different, and the scale of disease measures is different, the correlation should not be determined using only the slope values of the equations.

This study has several limitations. First, the number of enrolled patients with RA was small. This is because only naïve patients without prior treatment who were diagnosed with RA at our hospital were enrolled. Secondly, using WES, the variants detected in this study had no definite associations with previous reports. Third, this study focused only on TCR diversity according to DMARD treatment. Follow-up studies with a larger number of patients with RA will be able to clarify the clinical significance of the variants detected by WES in this study, detect other variants, and increase the diagnostic and predictive value of TRB/TRG diversity as a marker for the diagnosis, treatment, and prognosis of RA.

5. Conclusions

This study investigated potential genetic variants for RA development in Korean patients using WES and evaluated the clinical significance of the TCR repertoire by analyzing TCR/TRB diversity. Some potential variants were detected using WES; however, further studies are required to confirm their clinical significance. This study is meaningful in that Korean patients showed similar results and these factors and TCR diversities showed negative correlations and linear relationships because the correlation coefficient was high and statistically significant. This study provides basic data that can be applied to Korean patients by analyzing and evaluating results at the molecular genetic level.