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

This study investigated the potential genetic variants of rheumatoid arthritis (RA) using whole-exome sequencing (WES) and evaluated the disease course using T cell receptor (TCR) repertoire analysis. Fourteen patients with RA and five healthy controls (HCs) were enrolled. For the RA patient group, only treatment-naïve patients were recruited, and data were collected at baseline as well as at 6 and 12 months following the initiation of the disease-modifying antirheumatic drug (DMARD) treatment. Laboratory data and disease parameters were also collected. Genetic variants were detected using WES, and the diversity of the TCR repertoire was assessed using the Shannon–Wiener diversity index. While some variants were detected by WES, their clinical significance should be confirmed by further studies. The diversity of the TCR repertoire in the RA group was lower than that in the HCs; however, after DMARD treatment, it increased significantly. The diversity was negatively correlated with the laboratory findings and disease measures with statistical significance. Variants with a potential for RA pathogenesis were identified, and the clinical significance of the TCR repertoire was evaluated in Korean patients with RA. Further studies are required to confirm the findings of the present study.


Introduction
Rheumatoid arthritis (RA) is a systemic autoimmune disease characterized by chronic inflammatory synovitis, joint and bony destruction, and multi-organ manifestations, resulting in morbidity and mortality [1,2].RA has many complex etiologies involving genetic, immunological, and acquired factors [2].In addition, environmental factors, such as smoking, respiratory diseases, and changes in the gut microbiome, may increase the risk of or cause epigenetic modifications in susceptible patients, which may lead to RA development [3,4].For these reasons, identifying target genes for RA is very important for preventing and/or treating patients with RA, which can also be applied to further understand the pathogenesis of RA and potential targets of therapy [5,6].
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 [7].The cause of RA is not yet clearly known because of the complex traits of various factors; however, it is widely established that genetic factors account for approximately 50% of RA cases [8,9].For many decades, many studies have been conducted to identify genetic variants in RA and identify treatment targets.Owing to the development of genome-wide association studies (GWAS), many genetic variants related to RA have been discovered [10].Genetic factors are important not only for preventive activities in individuals with a high probability of RA but also for the treatment process, response, severity, and prognosis [11].Current conventional and biological treatments sometimes fail or show only 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 [7].Although many studies have been conducted on the susceptibility to RA development using single nucleotide polymorphisms (SNPs), which have attempted to identify disease-specific genes or polymorphisms, these SNPs alone cannot accurately explain the etiology because the RA-related genetic background is complex [12].GWAS have suggested genetic loci for susceptibility to RA [13][14][15][16].With the advent of next-generation sequencing (NGS), some potential genetic variants have been analyzed [17].Meanwhile, abnormal and pathogenic T cell responses that evade normal immune functions could be considered one of the mechanisms of RA development [3], and the T cell receptor (TCR) repertoire has been studied in patients with RA [18].However, to date, there have been only a limited number of related studies, with no published research 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 TCR β and γ repertoire (TRB and TRG) analysis.This study aimed to provide basic data that can be applied to Korean patients by analyzing and evaluating the results at a molecular genetic level.

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 was treatment-naïve patients who were first diagnosed with RA in our hospital and who had not yet started the diseasemodifying 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 initiating DMARD treatment, and then 6 and 12 months after treatment.However, 3 patients were lost to follow-up during the study period, and finally 14 patients remained; 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 purposes, requiring only one visit.At each visit, blood samples were collected in two K 2 -ethylenediaminetetraacetic acid (EDTA) tubes, one for measuring the erythrocyte sedimentation rate (ESR) and the other for genomic testing, as well as one serum separating tube 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 analog scale (VAS) score; simplified disease activity index (SDAI); clinical disease activity index (CDAI); and fulfillment of the 1987 American College of Rheumatology (ACR) criteria (whether ≥4 out of 7 cases are met) [19] and the 2010 ACR/European Alliance of Associations for Rheumatology (EULAR) criteria (whether a score of ≥6 out of 10 was met) [20] for clinical classification or scoring system.
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.

Genomic Analysis
Peripheral blood mononuclear cells were isolated from EDTA-treated whole blood by Ficoll-Hypaque density gradient centrifugation.The protocols for WES and TRB/TRG were as follows.
For WES, genomic deoxyribonucleic acid (DNA) was extracted using the Chemagic TM Magnetic Separation Module I method (PerkinElmer Chemagen, Baesweiler, Germany) with a 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 pairedend reads.DNA sequence reads were aligned to the reference sequence based on the public human genome build GRCh37/UCSC hg19.Alignments were performed with BWAmem (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 [21].
For the TRB and TRG 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 assays (Cat No. 72270009 and 72270019, respectively) (Invivoscribe Inc., San Diego, CA, USA) were used for NGS.These assays amplify the DNA between the primers that target the conserved V and J regions of the T cell receptor genes, TRB and TRG [22,23].
Polymerase chain reaction amplicons were purified and quantified using a Qubit fluorometer (Thermo Fisher Scientific, Waltham, MA, USA).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).

Statistical Analysis
The numerical data results were checked for normality using the Kolmogorov-Smirnov test, which determined that a p-value greater than 0.05 had normality.For the parametric data, the results are presented as the means ± standard deviations (SDs), and the Student's t-test was used for comparison.For the non-parametric data, there results are presented as the median and interquartile range (IQR), and the Mann-Whitney U test was used for comparison.The results of the categorical data are presented as frequencies and percentages (%), and the χ 2 test was used for comparison.
For patients with RA, the results at baseline (treatment-naïve) were compared to those at the 6 and 12 month follow-up.To compare the DAS28, DAS28-ESR (=0.70 × ln [ESR]) and DAS28-CRP (=0.36 × ln [CRP + 1] + 0.96) [24] 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 the comparisons, except for the analysis of variance, were made individually for each subgroup.
All the 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 a p-value of < 0.05.

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 them were female.There was no statistically significant difference in the proportions of patients with hypertension, diabetes mellitus, or obesity between the patients with RA and the HCs.When comparing the 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 the patients with RA than in the HCs, and the median ESR and CRP values were also higher in the patients with RA than in the HCs.

Genomic Variants
Using WES, disease-related variants that have been reported to be associated with RA among several SNPs were detected in four (28.6%) of the patients with RA.Table 2 lists the variants detected using WES that are likely related to RA pathogenesis.Some of the patients with RA had these variants, and one patient even had multiple variants.However, these SNPs were not found in the HCs, leading us to believe that these variants may be associated with RA pathogenesis.Most of the variants were of uncertain significance (VUS), but there was one case of a pathogenic variant (PV).

Changes in Clinical Features after DMARD Treatment
Table 3 displays the results of comparing the clinical features of the treatment-naïve patients with RA before starting DMARD treatment with their clinical features at 6 and 12 months following the initiation of the DMARD treatment.For the DMARD treatment, mainly 10 to 15 mg of methotrexate was used per week, and 1 g of sulfasalazine per day was added for some patients.Notably, there were no non-responders or poor responders to the DMARD treatment.Compared to the baseline (before treatment), the laboratory findings after 12 months of treatment decreased significantly.There were also changes in the number of joints; the TJC44, SJC44, TJC28, and SJC28 counts at 6 months decreased significantly compared to the baseline.These facts are also reflected in the disease measures; the DAS28, HQA, VAS, SDAI, and CDAI scores decreased significantly at 6 months, and even at 12 months they remained significantly decreased compared to the baseline.

Changes in TCR Diversity after DMARD Treatment
The TRB and TRG diversities were assessed using the Shannon-Wiener diversity index.The formula is as follows: where R is the number of TCR CDR3 amino acid clones, and p i is the frequency of the ith clonotype [25].In general, the values of the Shannon-Weiner diversity index are usually between 1.5 and 3.5 and are known to be rare in excess of 4.5 [26].However, in previous RA studies, the values of the Shannon-Wiener diversity indices were very high and exceeded 4.5 [25,27].
The results of the comparison of TCR diversity between the HCs and patients with RA at baseline as well as at 6 and 12 months after starting DMARD treatment are illustrated in Figure 1.For TRB diversity, the Shannon-Wiener diversity indices in the HCs and patients with RA at baseline as well as at 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 1A).For TRG diversity, the Shannon-Wiener diversity indices in the HCs and patients with RA at baseline as well as at 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.The TRB and TRG diversities of the patients with RA after 12 months of treatment were much higher than when they were treatment-naïve as well as after 6 months of treatment.
1.52), respectively (Figure 1A).For TRG diversity, the Shannon-Wiener diversity indices in the HCs and patients with RA at baseline as well as at 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.The TRB and TRG diversities of the patients with RA after 12 months of treatment were much higher than when they were treatment-naϊve as well as after 6 months of treatment.

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.The numerical relationships between the TCR diversity and RA-related factors were considered, but the disease status (patients with RA or HCs) was not considered.Therefore, the data of the patients with RA and HCs were combined and analyzed together.For the laboratory findings, data from the patients with RA and HCs were included and analyzed, and only the ESR showed a significant correlation with both TRB and TRG.Disease measures showed significant correlations with all the items evaluated, especially with the DAS28-ESR, SDAI, and CDAI, with the correlation coefficient exceeding 0.5 (Table 4).

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.The numerical relationships between the TCR diversity and RA-related factors were considered, but the disease status (patients with RA or HCs) was not considered.Therefore, the data of the patients with RA and HCs were combined and analyzed together.For the laboratory findings, data from the patients with RA and HCs were included and analyzed, and only the ESR showed a significant correlation with both TRB and TRG.Disease measures showed significant correlations with all the items evaluated, especially with the DAS28-ESR, SDAI, and CDAI, with the correlation coefficient exceeding 0.5 (Table 4).Additionally, the relationship between the TCR (TRB and TRG) diversity (Shannon-Wiener diversity indices) and the value of each RA-related factor was obtained through linear regression.Figure 2 shows a linear regression diagram between the 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-Wiener diversity and RArelated factors are analyzed and presented using data from the 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 the disease measures, it is difficult to uniformly compare using only the slope and intercept.However, it can be seen that, except for ACCP, compared to the TCR diversity, the slope is negative; therefore, it can be concluded that the improvement in the laboratory findings and disease measures (decreases in numeric values) is related to the increase in the TCR diversity.■, healthy controls; ▲, treatment-naïve patients; •, patients with RA after 6 months of treatment; ♦, patients with RA after 12 months of treatment.

Discussion
This study aimed to find possible variations for the onset of RA in Korean treatmentnaï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 who were 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 [20].
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 [28] but has not yet been reported in RA.However, its association with RA should be considered because JAK3 inhibitors may be used for RA treatment [29]; follow-up studies on this are required.Peptidyl arginine deiminase type IV (PADI4) is a gene that converts arginine residues into citrulline residues in the presence of calcium ions; mutations in this gene are known to overproduce citrulline and cause a loss of tolerance, thereby increasing vulnerability to RA [30].Tumor necrosis factor ligand superfamily member 18 (TNFSF18) is known to activate T cells and B cells in connection with cell signaling [31], and tumor necrosis factor receptor-associated factor (TRAF) is known to be associated with TNF and interleukin (IL)-1 and to cause inflammation [32].NF-κB is known to be activated in the synovium of patients with RA to regulate the genes contributing to inflammation, such as TNF, IL-6, and IL-8, and treatment therapy targeting this gene is under study [5].However, it does not match the SNP found in the European [14] and Japanese studies [33]; the clinical significance of the SNP shown in this study is somewhat limited, and it is necessary to clarify this 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; however, more GWAS will be required as human genetic patterns are complex, making it difficult to set and establish targets and produce effective therapeutic effects [34].
Genetic and specific environmental factors combine to trigger RA, which activates an immune reaction and causes loss of immune tolerance; this also 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; this results in bony destruction [3].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 [27].A previous study reported that the TCR repertoire helps in the early diagnosis of RA [18].The TCR repertoire in patients with RA is known to be less diverse than in healthy individuals; this is related to the disease activity [25,27].
In this study, a TCR repertoire analysis was performed.Treatment-naïve patients without prior DMARD treatment were enrolled, and the data at baseline (before treatment) as well as at 6 and 12 months after the initiation of DMARD treatment were compared.The baseline laboratory findings of the patients with RA were significantly different from those of the HCs.The TCR (TRB and TRG) diversities of the patients with RA were lower than those of the HCs.However, the TCR diversity in the patients with RA increased after the initiation of DMARD treatment.Laboratory values, affected joint counts, and disease measures were significantly decreased.As the symptoms and disease-related factors improved according to the treatment response, TCR diversity increased.According to previous studies, as the progression or severity of a disease increases, the clinical expansion of T cells occurs, and their diversity decreases [35].The diversity of the TCR repertoire may change according to treatment; 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 the TCR repertoire diversity and disease-related factors were also evaluated.The factors were negatively correlated with TCR diversity; this is in line with the result of a previous RA cohort study [27].The DAS, SDAI, and CDAI scores were more strongly correlated with TCR diversity than the laboratory findings.The relationships between TCR 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, correlations 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 treatment-naïve patients without prior treatment who were diagnosed with RA at our hospital were enrolled.Furthermore, age-or gender-matched controls were not collected.We only recruited those who voluntarily visited the OPD through a notice of recruitment.Further studies with larger sample numbers and age-or gender-matched controls are required.Secondly, by using WES, the variants detected in this study have 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 and TRG diversity as a marker for the diagnosis, treatment, and prognosis of RA.

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 the Korean patients showed similar results; the disease-related 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 the results at a molecular genetic level.

Figure 1 .
Figure 1.Comparison of the (A) TRB and (B) TRG diversity between the healthy controls and patients with RA at various time points following the initiation of DMARD treatment.TRB, T cell receptor β repertoires; TRG, T cell receptor γ repertoires; DMARD, disease-modifying antirheumatic drug.* p < 0.05 represents statistical significance.Black solid lines represent the p values between patients after DMARD treatment and treatment-naϊve patients.And, brown dotted lines represent the p values between patients after DMARD treatment and healthy controls.

Figure 1 .
Figure 1.Comparison of the (A) TRB and (B) TRG diversity between the healthy controls and patients with RA at various time points following the initiation of DMARD treatment.TRB, T cell receptor β repertoires; TRG, T cell receptor γ repertoires; DMARD, disease-modifying antirheumatic drug.* p < 0.05 represents statistical significance.Black solid lines represent the p values between patients after DMARD treatment and treatment-naïve patients.And, brown dotted lines represent the p values between patients after DMARD treatment and healthy controls.

Figure 2 .
Figure 2. Linear regression between the (A) TRB and (B) TRG diversity and laboratory findings.TRB, T cell receptor β repertoires; TRG, T cell receptor γ repertoires; ESR, erythrocyte sedimentation rate; CRP, C-reactive protein; RF, rheumatoid factor; DAS28, disease activity score in 28 joints;ACCP, anti-cyclic citrullinated peptide.A total of n = 47 data (data of 14 patients with RA before and after DMARD treatment, as well as at 6 and 12 months, and 5 healthy controls' data) were used for analysis.■, healthy controls; ▲, treatment-naϊve patients; •, patients with RA after 6 months of treatment; ◆, patients with RA after 12 months of treatment.

Figure 2 .
Figure 2. Linear regression between the (A) TRB and (B) TRG diversity and laboratory findings.TRB, T cell receptor β repertoires; TRG, T cell receptor γ repertoires; ESR, erythrocyte sedimentation rate; CRP, C-reactive protein; RF, rheumatoid factor; DAS28, disease activity score in 28 joints; ACCP, anti-cyclic citrullinated peptide.A total of n = 47 data (data of 14 patients with RA before and after DMARD treatment, as well as at 6 and 12 months, and 5 healthy controls' data) were used for analysis.

Table 1 .
Baseline characteristics of the study participants.

Table 2 .
Variants likely to be associated with rheumatoid arthritis detected using whole-exome sequencing.

Table 3 .
Comparison of the clinical features of patients with rheumatoid arthritis before and after DMARD treatment.

Table 4 .
Correlations between TRB/TRG diversity and disease-related factors.