INTRODUCTION
Rheumatoid arthritis (RA) is a chronic inflammatory autoimmune disease that primarily affects the joints [1]. However, extra-articular involvement, including involvement of the lungs, is not rare [2]. The prevalence of RA-associated interstitial lung disease (RA-ILD) has been reported to be approximately 1–10% [2-8] and its presence is linked to reduced survival compared with RA patients without ILD (RA-nonILD) [4,8,9].
A telomere is a DNA-binding protein complex composed of TTAGGG repeats located at chromosome ends, protecting these ends during replication [10]. Shortened telomere length (TL) is associated with age-related diseases and cancers. In a UK Biobank study, shorter TL was associated with a higher risk of idiopathic pulmonary fibrosis (IPF) [11]. Compared with ILD in other connective tissue diseases, patients with RA-ILD frequently exhibit the usual interstitial pneumonia (UIP) pattern, which shares characteristics with IPF [12]. The cumulative evidence from two meta-analyses and a Mendelian randomization analysis consistently demonstrates that TL is shorter in patients with RA compared to healthy controls [13,14]. A previous cross-sectional study showed that RA-ILD patients had shorter TLs than RA-nonILD patients [15]. Yet, only one study has evaluated the association between TL and RA-ILD and compared the results with TL in RA-nonILD. No study has examined the relationship between TL and RA-ILD patterns or risk factors for shortened TL in patients with RA-ILD. In this study, we investigated whether TL differs between RA and controls, between RAILD and RA-nonILD, and between UIP and nonspecific interstitial pneumonia (NSIP) subtypes among RA-ILD patients. Furthermore, we investigated the factors influencing TL in RA-ILD, with a particular focus on the UIP pattern.
METHODS
Study population
This cross-sectional study included RA patients with or without ILD. The study participants were part of the RA cohort at Hanyang University Hospital for Rheumatic Diseases and the RA-ILD cohort at Hanyang University Guri Hospital. Patients were registered between October 2001 and July 2022 and between November 2021 and July 2022, respectively [16,17]. The study protocol was approved by the Institutional Ethics Review Boards (IRBs) of Hanyang University Hospital (IRB No. 2021-06-011) and Hanyang University Guri Hospital (IRB No. 2021-08-031). Written informed consent was obtained from all participants. The study was conducted in accordance with the World Medical Association’s Declaration of Helsinki.
The study included RA patients who underwent chest computed tomography (CT) scans after cohort enrollment and had peripheral blood samples collected at the time of enrollment. Based on CT scan results, study participants were categorized into RA-ILD and RA-nonILD groups. The RA-nonILD group comprised patients whose CT scans revealed no clinically significant abnormalities. The RA-ILD group consisted of patients exhibiting UIP or NSIP subtypes on their CT scans. Patients with CT findings other than UIP or NSIP, such as organizing pneumonia, indeterminate for UIP, or unclassifiable ILD, were excluded from the RA-ILD group. Patients with other lung abnormalities, including bronchiectasis, pneumonia, emphysema, bronchiolitis, pulmonary edema, pleural effusion, malignancy, nontuberculous mycobacterium, or mycobacterium tuberculosis, were also excluded from the RA-ILD and RA-nonILD group. RAILD patients with more than a 1-year interval between the chest CT scan and blood sampling were excluded. Two experienced radiologists reviewed and confirmed all CT scan findings.
Data from healthy controls, including age, sex, and TL, were obtained from MEDIZEN HUMANCARE. These data were then matched with RA patients (with or without ILD) in a 1:5 ratio based on age and sex using propensity score matching (PSM) to minimize bias from confounding variables. RA-ILD and RA-nonILD patients were also matched 1:1 by age and sex using propensity score matching (PMS). Among patients with RA-ILD, those with NSIP or UIP patterns were matched in a 1:4 ratio. This study calculated and selected the optimal matching ratios for each group (control vs. RA, RA-ILD vs. RA-nonILD, and NSIP vs. UIP) to maximize statistical power, while effectively utilizing the limited sample size available for analysis. All RA patients fulfilled either the 2010 European League Against Rheumatism-American College of Rheumatology criteria or the 1987 American College of Rheumatology revised criteria for RA [18,19].
Clinical data collection
At the time of enrollment, demographic data (age at sampling, age at RA diagnosis, sex, and disease duration) and smoking status (never, ever, or current) were collected. Rheumatoid factor (RF) and its titer, anti-cyclic citrullinated peptide (CCP) antibody and its titer, disease activity score 28 (DAS28)–erythrocyte sedimentation rate (ESR), DAS28–C-reactive protein (CRP), RA medication, and lung function data were obtained within 6 months of the chest CT scan.
TL measurement and genotyping of the TERC and TERT variants
Leukocytes were extracted using the Gentra Puregene kit (Qiagen, Germantown, MD, USA). The process involved combining 3 mL of blood with 9 mL of RBC Lysis Solution, followed by centrifugation to obtain a white pellet. Subsequent treatment with Cell Lysis Solution, Protein Precipitation Solution, and isopropanol, along with intermediate centrifugation steps, produced a final DNA pellet, which was washed with 70% ethanol and air-dried. TL was determined from these leukocytes using a standard quantitative polymerase chain reaction (qPCR) assay, following the protocol described by O’Callaghan and colleagues [20]. The normalizing control gene used was 36B4, which encodes the acidic ribosomal phosphoprotein. Oligomer sequences are shown in Supplementary Table 1. In each reaction with a total volume of 20 μL, 20 nanograms per microliter of target DNA, along with telomere or 36B4 primers and SYBR Green mixed with low ROX reference dye (Cat No: DQ384-40h; BioFACT, Daejeon, Korea), were added following the manufacturer’s instructions. PCR commenced with 1 cycle of 95°C for 10 minutes, followed by 40 cycles of 95°C for 15 seconds and 60°C for 1 minute using a QuantStudioTM 6 Flex Real-Time PCR System (Thermo Fisher, Waltham, MA, USA), and analysis was conducted with QuantStudio™ Real-Time PCR Software (Thermo Fisher) in accordance with the manufacturer’s guidelines. Samples deviating more than 0.25 Ct in standard deviation were excluded from the analysis, in accordance with Thermo Fisher’s guide. Absolute TL was calculated as previously described using the kb/reaction value to determine total TL in kb per human diploid genome [20]. The telomere kb per reaction value was divided by diploid genome copy number to yield the total TL in kb per human diploid genome. The results were presented as the length in kb of individual telomeres, calculated by dividing the telomeric sequence per diploid genome in kb by 92 (the number of individual telomeres in the diploid genome). Additionally, we genotyped two polymorphisms of TL-related genes―the telomerase ribonucleic acid component (TERC) rs12696301 (C/G) and telomerase reverse transcriptase (TERT) rs2736100 (A/C)―using the QuantStudio™ 6 Flex Real-Time PCR System (Thermo Fisher) with a 96-well block. PCR conditions included 95°C for 10 minutes, followed by 40 cycles of 95°C for 15 seconds and 60°C for 1 minute.
Quantitative CT texture analysis
We assessed the extent of honeycombing using commercially available artificial intelligence (AI) lung texture software (AVIEW Lung Texture, v48.8; Coreline Soft, Co. Ltd, Seoul, Korea). All high-resolution CT (HRCT) images were obtained using standard HRCT protocol parameters, and the lung parenchyma of axial HRCT images was reconstructed. Among the reconstructed images, honeycombing was characterized by extensive fibrosis with lung destruction, leading to a cystic reticular appearance [21]. The software automatically measured the honeycombing extent using a two-dimensional-based deep learning algorithm in each section of the HRCT images [21,22] and expressed this as a percentage of the total lung volume. An experienced radiologist (S.J.H.) supervised these processes.
Statistical analysis
We compared baseline characteristics between patients with RA-ILD and RA-nonILD using an independent Student’s t-test for continuous variables and the chi-square test for categorical variables. The TL between RA-ILD patients with UIP and NSIP was compared with a Mann–Whitney U test. In the PSM conducted to compare TL between RA-ILD patients with UIP and NSIP, only age and sex variables were included. Propensity scores for the study participants were calculated using the “glm()” function in R statistics, based on the logistic regression model logit(p) = β0 + β1X1 + β2X2, where p represents the propensity score, and X1 and X2 are the independent variables age and gender, respectively. The matching process employed the nearest neighbor method with a caliper of 0.2, ensuring that control group subjects (e.g., RA-nonILD) were matched with treatment group subjects (e.g., RA-ILD) only if their propensity scores were within 0.2 of each other. After matching, the standardized mean difference was calculated for each variable to confirm the balance between the two groups. Following this, to further assess other potential factors influencing TL in RA-ILD patients beyond just the RA-ILD patterns, linear regression analyses were performed. This additional analysis aimed to identify other relevant variables that may have an impact on TL in these patients. Absolute TL was set as the dependent variable in univariate analysis, while age at sampling, age at RA diagnosis, sex, smoking status, RF positivity, RF titer, anti-CCP antibody positivity, anti-CCP antibody titer, DAS28-ESR, DAS28-CRP, diffusing capacity of the lung for carbon monoxide (DLCO), and the presence of UIP were considered as independent variables. Variables with a significance level of p < 0.05 in the univariate analysis were chosen for inclusion in subsequent multivariable analyses.
We used Spearman correlation analysis to evaluate the relationship between TL and honeycombing extent in patients with RA-ILD. All tests were two-sided, and p values less than 0.05 were considered statistically significant. We performed all statistical analyses using R (version 4.1.0) and SPSS Statistics for Windows, version 17 (SPSS Inc., Chicago, IL, USA).
RESULTS
Of the 2,536 patients in the RA cohort, 749 underwent a chest CT scan. From these, 161 patients with RA-ILD and 588 patients with RA-nonILD were initially selected. Subsequently, out of the 161 RA-ILD patients, 70 patients were chosen after the exclusion of those with lung parenchyma abnormalities other than ILD findings and those with ILD subtypes other than UIP or NSIP patterns. Additionally, from the 588 RA-nonILD patients, a total of 146 patients were selected after excluding those with lung parenchyma abnormalities (Fig. 1). From the 216 selected patients with RA (70 with RA-ILD and 146 with RA-nonILD), specific individuals were chosen for analysis to ensure age and gender matching for each analysis.
TL comparison: control participants versus RA patients
We began by comparing the TL of control participants (n = 14) with that of age- and sex-matched RA patients (n = 70) (Fig. 2). The mean ages of the control participants and RA patients were 55.0 ± 3.4 years and 57.0 ± 6.1 years, respectively (p = 0.100) (Supplementary Table 2). The proportions of women in the control and RA groups were 64.3% (n = 9) and 70.0% (n = 49), respectively (p = 0.754). Half of the RA patients (n = 35) were RA-ILD. The TL of the RA patients was 8.3 ± 2.5 kb, significantly shorter than 9.5 ± 0.8 kb TL of control group (p = 0.002; Fig. 2A).
TL comparison: RA-ILD versus RA-nonILD patients
Next, we assessed whether the presence of ILD in RA patients was associated with TL. After 1:1 age- and sex-matching among the 216 selected patients with RA (70 with RA-ILD and 146 with RA-nonILD), we compared TL between 53 patients with RA-ILD and 53 patients with RA-nonILD. Baseline demographics and clinical characteristics, including autoantibody positivity, disease activity, and medication, did not differ significantly between the two groups (Table 1). The mean absolute TL for RA patients (n = 106) was 8.0 ± 2.6 kb, with no significant difference observed between RA-ILD and RA-nonILD patients (8.2 ± 2.8 kb and 7.7 ± 2.4 kb, respectively, p = 0.271; Table 2). There were also no differences in TERC and TERT polymorphisms between the two groups.
Factors associated with TL in patients with RA-ILD
While RA-ILD itself did not exhibit a difference in TL compared to RA-nonILD, we explored whether specific conditions or subgroups within patients with RA-ILD might be associated with TL. This analysis is grounded in the heterogeneous characteristics of ILD, encompassing various pathologic patterns with different clinical outcomes. Specifically, our focus was on the UIP pattern, which shares characteristics with IPF, a condition known to be associated with shortened TL. To explore independently associated factors with shortened TL, we performed univariate and multiple linear regression analyses in RA-ILD patients, incorporating 1:4 age- and sexmatched patients with NSIP (n = 8) and UIP (n = 32). In a univariate linear regression model, significant associations of TL were observed with age at sampling (beta coefficient = -0.14, p = 0.0004), DAS28-CRP (beta coefficient = -0.71, p = 0.017), DLCO (beta coefficient = -0.06, p = 0.004), and UIP (beta coefficient = -2.57, p = 0.004). On multiple linear regression analysis, UIP (beta coefficient = -1.38, p = 0.038) remained independently associated with TL in RAILD patients with age at sampling (beta coefficient = -0.08, p = 0.011) and DAS28-CRP (beta coefficient = -0.60, p = 0.018) (Table 3). The mean absolute TL of UIP patients was significantly shorter than that of NSIP patients (7.4 ± 1.9 kb and 10.0 ± 3.1 kb, respectively, p = 0.026; Fig. 2B). Clinical characteristics of patients with UIP and NSIP, except TL, did not exhibit significant differences (Supplementary Table 3).
Association of honeycombing extent with TL
Lastly, we investigated whether the extent of honeycombing, a hallmark of UIP, correlated quantitatively with absolute TL. Among the total 45 RA-ILD patients with UIP, honeycombing extent ranged from 0.1 to 11.2%. The correlation analysis between TL and honeycombing extent showed a trend of a negative correlation of TL with the percentage of tissue exhibiting a honeycombing appearance within the entire lung (Rho = -0.131, p = 0.387; Fig. 3A), but there was no significant association between them. However, a strong negative correlation was observed between honeycombing extent and TL in the 11 RA-ILD patients with UIP belonging to the highest quartile for honeycombing extent (Rho = -0.745, p = 0.008; Fig. 3B).
DISCUSSION
Our study validated earlier findings indicating that RA patients exhibit shorter TL than healthy individuals, but our results did not establish a significant association between the presence of RA-ILD and TL in patients with RA. However, we did observe differences in TL based on the specific patterns of RA-ILD. The presence of the UIP subtype was associated with shortened TL when compared with that of age- and sex-matched patients with NSIP. Furthermore, although statistical significance was not reached, honeycombing extent within RA-ILD patients with UIP showed a tendency toward a negative correlation with TL. To our knowledge, our study is the first to investigate TL in relation to the patterns of RA-ILD and explore potential correlations between honeycombing extent and TL.
Our analysis contrasted with a study by Natalini et al. [15], which reported significantly shorter TL associated with RA-ILD than with RA-nonILD. As the sample size in our study (n = 53) was comparable to that in the prior study (n = 54), the observed discrepancy may be attributed to differences in the characteristics of the patient cohorts. Natalini et al. [15]’s study predominantly involved males (94%), participants with White or Caucasian ethnicity (75%), and current or former smokers (91%), whereas our study primarily included females (76%), those with Asian ethnicity (100%), and non-smokers (79%). Since research on the TL between RA-ILD and RA-nonILD is currently very limited, further large-scale studies with diverse patient characteristics are needed to validate these findings.
We demonstrated an independent association of UIP subtype with shortened TL compared with that in NSIP in a multiple linear regression analysis. Previous reports have shown that IPF patients had significantly shorter TL in both peripheral leukocytes and lung tissue than age-matched healthy controls [11,23,24]. Snetselaar et al. [25] compared TL across a spectrum of ILDs, including IPF, connective tissue disease-related ILD, sarcoidosis, idiopathic NSIP, cryptogenic organizing pneumonia, etc., consistently finding TL in IPF to be shorter than that in idiopathic NSIP and controls. One study, comparing TL between the UIP pattern and non-UIP patterns such as NSIP, organizing pneumonia, and lymphocytic interstitial pneumonia in CTD-ILD, found that patients with the UIP pattern have significantly shorter telomeres [26]. To date, however, no study has directly compared TL between RA-ILD patients with the UIP pattern and those with the NSIP pattern. In our study, we observed that RA patients with the UIP pattern exhibited shorter TL, in line with previous reports from clinically and pathogenetically similar IPF studies and consistent with results from UIP studies in CTD-ILD.
This finding suggests that RA-ILD with a UIP pattern shares a pathophysiology similar to that of IPF, involving telomere dysfunction. The telomere system is associated with biological aging-like disease processes, which also contribute to the development of IPF [27]. It is therefore believed that, in RA-ILD with a UIP pattern, aging-like disease pathologies, in addition to autoimmunity, may similarly play a role in lung fibrosis. On the other hand, in NSIP, where inflammation rather than fibrosis is the predominant pathological mechanism, TL does not appear to be as shortened, suggesting that autoimmunity and inflammation, rather than aging-like disease processes, are more likely the primary contributors to pathology. Interestingly, even in RA patients, where synovial inflammation is the dominant pathological feature rather than fibrosis, TL is also found to be shorter compared to control groups. Although the precise mechanism underlying this observation is not fully understood, some potential links can be inferred from shared non-genetic risk factors between RA, IPF, and shortened TL. Smoking is a common risk factor in all three conditions [28], and it has been reported that smokers have shorter TL compared to non-smokers, with a negative correlation between the number of cigarettes smoked and TL [29,30]. This suggests that smoking may contribute to the process of telomere shortening in RA. However, further studies are required to clarify the potential mechanisms behind these suggestions.
The presence of honeycombing in HRCT images is a crucial prognostic factor in IPF, with the extent of honeycombing serving as a predictor of mortality in these patients [31]. In the context of RA-ILD, honeycombing may be more indicative of prognosis than the presence of the UIP pattern alone. Radiological honeycombing and a high fibrosis score have both been independently associated with mortality among patients with RA-ILD [32]. Ito and colleagues found that fibrosis score (the combined extent of reticulation and honeycombing) was strongly associated with worse survival in RA-ILD. Patients with fibrosis scores > 20% had an approximately nine-fold increased risk of mortality [33]. In our study, honeycombing extent within RA-ILD with UIP showed a tendency toward a negative correlation with TL, although statistical significance was not reached. However, in RA-ILD patients with UIP belonging to the highest quartile for honeycombing extent, a strong correlation between honeycombing extent and TL was observed. These observations suggest that, while a subtle change in honeycombing appearance may not have a significant impact, a substantial extent of honeycombing appearance is significantly negatively correlated with TL. Shortened TL has previously been independently associated with worse survival among patients with IPF [34]. These results suggest that shortened TL is associated with the pathogenesis of RA-ILD with a UIP pattern and the progression of honeycombing and that the combined effects of these factors might contribute to poor survival. Given that no studies have directly examined the association between TL, honeycombing extent, and survival, especially among patients with RA-ILD with a UIP pattern, further research is warranted.
To quantify honeycombing extent in this study, AI analysis was utilized instead of human analysis. This AI analysis was conducted on patients displaying characteristic CT findings of UIP or NSIP, as identified by experienced radiologists. The AI analysis has been validated in prior studies [21,22]. Remarkably, the AI demonstrated high levels of similarity in image analysis compared gold standard by manual correction, distinguishing ILD patterns with a similarity of 98.84 ± 0.55% [22]. Considering that interobserver agreement among radiologists in identifying UIP patterns on HRCT in RA-ILD patients can be relatively modest [35], the use of AI analysis tools is anticipated to be advantageous in future research on RA-ILD, particularly in relation to UIP patterns.
The limitation in the study is that our sample size was relatively small, which may have limited the power of the study in comparing participants. Particularly, there were few patients in the control groups (fourteen) due to the limited availability of control subjects and the NSIP (only eight) due to the exclusion of RA-ILD patients with potential confounding factors such as ILD with patterns other than UIP or NSIP (including indeterminate UIP and OP); patients diagnosed with BE, lung cancer, or non-tuberculosis mycobacteria infection; and patients who had more than a 1-year interval between the time of chest CT scan and blood sample collection. Moreover, during the age- and sex-matching process, the number of subjects available for analysis further decreased, leaving us with a reduced sample size for the final analysis. The effect size calculated using Cohen’s d when comparing RA patients to control groups was 0.52, with a statistical power of 0.42 (42%) at α = 0.05, as determined by G*Power. Despite the lower power, our results align with previously published studies. Furthermore, for the interest of our study―the TL difference between UIP and NSIP―we observed a large effect size (1.34) and high statistical power (0.91, or 91%) despite the small sample size, indicating robust findings in this comparison. Our strict selection strategy, while mitigating bias, combined with a small sample size necessitate future studies with larger cohorts to confirm our findings.
In conclusion, our study confirmed the previous finding that TL in patients with RA was significantly shorter than in control participants. Although we found no significant difference in TL when comparing RA patients with or without ILD, we found that TL varies depending on the subtype in patients with RA-ILD, suggesting a possible relationship between shortened TL and the pathogenesis of UIP, especially regarding the progression of changes in honeycombing.