The data for this controlled cross-sectional study came from a prospective cohort of incident cases recruited between 2007 and 2011 in the Department of Rheumatology, Hospital Universitario Regional de Málaga, Málaga, Spain (9). All the patients were aged ≥18 years, fulfilled the 2010 criteria of the American College of Rheumatology/European League Against Rheumatism for RA (14), and had been diagnosed and treated for the first time during the 12 months after onset of their disease. All participants provided their written informed consent before inclusion. The study was conducted according to the principles of the Declaration of Helsinki, and the study protocol was approved by the Ethics Committee of Málaga (Project identification code 4/2016, P19).

Since the creation of the prospective cohort (2007–2011), all patients have been followed up at the outpatient clinic every 3–6 months by a rheumatologist using a systematic clinical data collection protocol. The data collected included inflammatory activity and physical function throughout follow-up.

During 2016–2018, all the patients in the cohort consented to a relevant modification to the protocol aimed at new, broader study objectives that enabled, among other options, the present study. Peripheral venous blood samples were collected after overnight fast, and fecal samples were refrigerated immediately and transported to the laboratory, where they were stored at −80°C for subsequent analysis. The date of this visit was the index date.

At the last visit in 2022 (final visit in the present study), we identified 2 groups of patients with extreme RA phenotypes: (1) a group comprising patients with difficult-to-treat RA (D2T RA); and (2) a group of easy-to-treat RA (E2T RA) patients at a 2:1 ratio, matched by age, sex, and time since diagnosis.

D2T RA was defined according to the EULAR criteria for D2T RA, as follows: failure of ≥2 biological or targeted synthetic disease-modifying antirheumatic drugs (b/tsDMARDs) (with different mechanisms of action) after failure of a csDMARD (unless contraindicated); signs suggestive of active/progressive disease; and management of signs and/or symptoms perceived as problematic by the rheumatologist and/or the patient (3). E2T RA was defined as treatment throughout follow-up with a single csDMARD and neither active nor progressive disease. Following the EULAR recommendations (3), active/progressive disease was defined as ≥1 of the following: at least moderate disease activity (according to validated composite measures including 28-joint Disease Activity Score for Rheumatoid Arthritis with erythrocyte sedimentation rate [DAS28-ESR] >3.2), signs and/or symptoms suggestive of active disease, inability to taper glucocorticoid treatment (below 7.5 mg/day prednisone or equivalent), rapid radiographic progression, and symptoms that diminish quality of life.

On the index date, we collected demographic, clinical, and treatment-related data, including comorbid conditions associated with traditional cardiovascular risk factors (smoking, obesity, arterial hypertension, diabetes mellitus, history of cardiovascular disease, and sedentary lifestyle).

Similarly, on the index date, we evaluated inflammatory activity at the visit and calculated the cumulative activity. Activity was estimated as an arithmetic mean of all the values collected regularly since diagnosis (time-averaged disease activity). Inflammatory activity was measured using the DAS28-ESR (range, 0–9.4) (15). A DAS28-ESR value >3.2 was considered high and ≤3.2 was considered low. The erythrocyte sedimentation rate (ESR; mm/h) was measured. High-sensitivity C-reactive protein (hsCRP; mg/L) was measured for all participants using nephelometry (MMAGE-Immunochemistry Systems, Beckman Coulter, Brea, CA, United States). Physical function on the index date was assessed (average value) using the Health Assessment Questionnaire (HAQ).

We recorded the presence of autoantibodies and their titers, as follows: rheumatoid factor (positive if >10 IU/mL) and anti–citrullinated peptide antibody (positive if >20 IU/mL), as well as the presence or absence of radiologic erosions. We recorded all csDMARDs, such as methotrexate, leflunomide, and sulfasalazine, and bDMARDs, such as tumor necrosis factor inhibitors (anti-TNF), tocilizumab, abatacept, and rituximab. We also included tsDMARDs, such as the Janus kinase inhibitors (JAKi) tofacitinib and baricitinib, and glucocorticoids. Adherence to the Mediterranean diet was evaluated using a validated questionnaire in which adherence was defined as a score of >9 out of 14 (16).

Inflammatory mediators such as TNF-α, IL-1β, and IL-6 in plasma were quantified using enzyme-linked immunosorbent assay (ELISA) Quantiglo kits (R&D Systems Inc., Minneapolis, United States) according to the manufacturer’s instructions. Plasma levels of insulin-like growth factor I were analyzed using ELISA (Mediagnost GmbH., Tuebingen, Germany). Malondialdehyde-oxidized low-density lipoprotein (LDL) was also measured in plasma using an ELISA kit (Biomedica GmbH., Vienna, Austria).

DNA was extracted using the QIAamp DNA Stool Mini Kit (Qiagen, Hilden, Germany) according to the manufacturer’s instructions. The concentration and purity of DNA were determined using a Nanodrop spectrophotometer (Nanodrop Technologies, Wilmington, DE, United States). Ribosomal 16S RNA (16S rRNA) gene sequences were amplified from DNA using the Ion 16S Metagenomics Kit (Thermo Fisher Scientific, Waltham, MA, United States). The kit includes 2 primer sets (V2-4-8 and V3–6, 7–9) that selectively amplify the corresponding hypervariable regions of the 16S region in bacteria. Libraries were built with the Ion Plus Fragment Library kit (Thermo Fisher Scientific). Barcodes were added to each sample using the Ion Xpress^™ Barcode Adapters kit (Thermo Fisher Scientific). Emulsion PCR and sequencing of the amplicon libraries were performed on an Ion 530 chip (Ion 530^™ Chip Kit) using the Ion Chef System and Ion Torrent S5^™ system (Ion 510^™/520^™/530^™ Kit-Chef, Thermo Fisher Scientific), respectively, according to the manufacturer’s instructions.

Base calling and run demultiplexing were performed using Torrent Suite™ Server software (Thermo Fisher Scientific), version 5.4.0, with default parameters for targeted sequencing of 16S (bead loading ≤30, key signal ≤30, and usable sequences ≤30). Quality sequences were further translated into amplicon sequence variants (ASVs) using DADA2 with adapted parameters for Ion Torrent data within the open-source Quantitative Insights into Microbial Ecology software (QIIME2, version 2022.2) (17). Taxonomic assignment was performed through clustering with VSEARCH and the reference base Greengenes version 13_8 at 97% identity. Before analysis, samples with fewer than 1,500 sequences, features with a count sum of less than 10 across all samples, those presented only in one sample, mitochondrial features, and features unidentified at the phylum level were removed in the pre-processing step. Alpha and beta diversities were evaluated using the core-metrics-phylogenetic plugin in QIIME2 after rarefaction to the minimum number of sequences per sample. Alpha diversity was assessed through different indexes, including Pielou-evenness, Faith-PD, observed-features, and Shannon index. Statistically significant differences in alpha diversity between groups were determined by using the Kruskal-Wallis test with a significance level of p = 0.05. The beta diversity was assessed using non-phylogenetic Bray–Curtis dissimilarity index. The structure of the gut microbial community of the different groups was examined using PCoA plots for Bray-Curtis distances. The significance of variations among groups was determined using permutational multivariate analysis of variance (PERMANOVA) with a p-value of 0.05. Linear discriminant analysis (LDA) Effect Size (LEfSe) was used to identify potential biomarker taxa. LEfSe was performed on the online Galaxy web application¹ (18), where data were normalized using counts per million (CPM). The nonparametric factorial Kruskal–Wallis sum rank test was used to detect differential abundant features (at genera, family, class, order, and phylum level) within the two groups. As a last step in the LEfSe analysis, LDA was used to determine the effect size of each differentially abundant feature. The differences in microbial relative abundance were considered significant when the LDA < 2 and p < 0.05 (19). The cladogram from the LEfSe method indicates the phylogenetic distribution representing differentially abundant taxonomic groups. The size of each node represents its relative abundance. The Phylogenetic Investigation of Communities by Reconstruction of Unobserved States plugin (PICRUSt2) (20) was used to predict metagenome function within QIIME2 with the DADA2 output. MetaCyc pathways (21) were normalized within QIIME2 and further analyzed using the open-source software STAMP (Statistical Analysis of Metagenomics Profiles) with Welch’s t test option (22).

A descriptive analysis of the main outcome measures was performed. Values are expressed as frequencies and percentages or as mean (standard deviation [SD]) or median (interquartile range [IQR]), as applicable. Normality was assessed using the Kolmogorov–Smirnov test. We compared clinical and laboratory characteristics and inflammatory activity between patients with D2T RA and patients with E2T RA using the Pearson χ² test or the t test, as applicable. A multivariate analysis was performed, as was a Cox regression analysis, to identify factors associated with D2T RA adjusted for disease duration. The variables entered into the models were those that proved to be significant in the bivariate analysis and of clinical interest. Statistical significance was set at p < 0.05. The statistical analyzes were performed using IBM SPSS Statistics for Windows, Version 28 (IBM Corp., Armonk, NY, United States).

Table 1 shows the clinical and epidemiological characteristics of D2T RA patients and E2T RA patients on the index date. Data from both groups were consistent for most of the clinical-epidemiological characteristics and comorbid conditions, except that erosions were more frequent in D2T RA patients than in E2T RA patients (p = 0.021). There were no significant differences in the studied comorbidities related to cardiovascular risk, osteoporosis, fibromyalgia, or depression between the two groups. As for treatment, most D2T RA patients were receiving treatment with an anti-TNF agent on the index date (69.2%) or an IL-6 inhibitor (23.1%) (Table 1).

However, at the final visit (see Table 2), most of the 13 patients with D2T RA were receiving rituximab (30.8%) or tofacitinib (23.1%) after a mean (SD) of 2.6 (1.3) switches of biologics, with a mean (SD) retention period of 95.8 (56.3) months. In patients with D2T RA, the main reason for switching biological or targeted synthetic disease-modifying antirheumatic drugs (b/tsDMARDs) was loss of efficacy (18/34 treatments used [52%]), followed by insufficient response (11/34 [32%]) and nonserious adverse events (5/34 [14%]).

In contrast, as expected, the 26 E2T RA patients maintained the same csDMARDs throughout follow-up. Similarly, a larger number of patients with D2T RA were receiving glucocorticoids on the index date (p = 0.018) and at a higher median dose than the E2T RA patients (p = 0.040).

Of the 39 patients included, 28 (71.8%) were in remission or with low disease activity at the index date according to their DAS28-ESR values, and 30 (76.9) had maintained an average DAS28-ESR value indicating low activity. D2T RA patients had higher average DAS28-ESR values at the cut-off than the E2T RA patients (Table 3). The same was true of physical functioning according to the HAQ. However, laboratory values were generally similar for both groups, except for some notable differences, such as homocysteine (p = 0.010) and CRP (p = 0.029), which were superior in D2T RA.

Regarding proinflammatory cytokines, D2T RA patients had higher levels of IL-6 (p = 0.031), and numerically higher levels of TNF-α (p = 0.085). In the case of lipoproteins and human growth factors, the values remained similar in both groups. Additionally, D2T RA patients were similar to those with E2T in their adherence to the Mediterranean diet (69.2% vs. 73.0%; p = 0.901).

There were no differences in any of the alpha diversity measured indexes (Pielou’s evenness, p = 0.121; Faith’s phylogenetic diversity, p = 0.121; observed features, p = 0.318; and Shannon’s diversity, p = 0.101) between E2T RA and D2T RA patients (Figure 1A). While microbiota populations tended to differ based on beta diversity (Bray–Curtis dissimilarity; PERMANOVA, p = 0.089; Figure 1B).

At the phylum level, the dominant bacterial phyla were Bacteroidetes, Firmicutes and Proteobacteria, while Actinobacteria, Synergistetes, Lentisphaerae, and Verrucomicrobia accounted for smaller proportions, between 1 and 10%, in both RA groups (Figure 2A). Analysis of the abundance of Bacteroidetes and Firmicutes did not reveal significant differences between E2T RA and D2T RA (p = 0.251 and p = 0.486, respectively; Figure 2B). However, the Firmicutes/Bacteroidetes ratio was lower in the D2T RA patients than in the E2T RA patients (p = 0.011; Figure 2C).

Additionally, we performed LEfSe to identify changes in the gut microbiota between E2T RA and D2T RA patients. This analysis revealed an enrichment in the phylum Bacteroidetes and in the genus Megasphaera in the D2T RA group. In contrast, the E2T RA group exhibited a greater abundance of the phylum Proteobacteria, the family Pasteurellaceae and its genus Haemophilus, and the family Lachnospiraceae and its genus Coprococcus (LDA > 2; p < 0.05; Figure 3).

Metacyc pathway analysis was performed using PiCRUSt2 to increase our understanding of the role of gut microbiota in each of the groups studied. A Kruskal-Wallis analysis showed that 19 pathways differed between the groups (p < 0.05). The main pathways affected belonged to Degradation/Utilization/Assimilation pathways and Biosynthesis. In D2T RA patients, increased values were reported for the thiamine diphosphate salvage II pathway (PWY-6897) and the carbohydrate degradation pathway (PWY-6507). Increased values were also reported for E2T RA patients in the remaining pathways, some of which are involved in peptidoglycan biosynthesis and generation of precursor metabolites and energy (PWY0-1586 and REDCITCYC, respectively), and most of which are involved in degradation, specifically, aromatic compound (PWY-6185, PWY-5417, PWY-5431, PWY-7002, PWY5F9-12), toluene degradation (PWY-5178 and PWY-5181), fatty acid and lipid degradation (PWY-6946 and LIPASYN-PWY), carbohydrate degradation (PWY-7644), secondary metabolite degradation (PWY-6507), and nucleoside and nucleotide degradation (PWY-6353) (Figure 4).

Table 4 shows the results of the Cox multivariate analysis (DV: D2T RA), in which 39 patients with RA were included over a mean (SD) follow-up of 103.8 (37.8) months. A total of 13/39 patients had D2T RA. The multivariate analysis showed that the Firmicutes/Bacteroidetes ratio was associated with a reduced risk of D2T RA (HR, 0.288; 95% CI, 0.092–0.907; p = 0.033), whereas the variables associated with a greater probability of D2T RA were greater average inflammatory activity according to the DAS28-ESR (HR, 2.649; 95% CI, 1.225–5.732; p = 0.013) and treatment with prednisone on the index date (HR, 3.794; 95% CI, 1.098–10.990; p = 0.008). Thus, for each 0.1-point increase in the Firmicutes/Bacteroidetes ratio, the risk of D2T RA decreased by approximately 71%.