Gut microbiome-immune interactions and their role in rheumatoid arthritis development

A study population

In this observational study, all participants were offered informed consent to join the research after they had an explanatory discussion about the study’s objectives and structure (they received information about the study, and a copy of the informed consent was provided). All participants who were part of the study underwent a registration process and were assigned an individual code to maintain confidentiality principles. The recruitments were carried out at City polyclinics No.10, No.12 and the Treatment and Diagnostic Center based at City Multidisciplinary Hospital No. 2 Astana, Kazakhstan. All RA patients were examined by the rheumatologist.

The research received approval from the ethics committee at the National Laboratory Astana, Nazarbayev University, with Protocol No. 03-2019, in Astana, Kazakhstan. Participants were informed about the study’s objectives and provided their consent by signing an informed consent document. We confirm that the informed consent was received from all participants of our study in the article form. Interactions with the participants adhered to the established guidelines and regulations.

The main group of the study were Kazakh female patients with the established RA diagnosis according to American College of Rheumatology 2010 classification criteria for RA (Radner et al., 2014), and the disease duration more than 1 year. Disease activity was determined by Disease Activity Score Calculator (DAS28), X-ray stage was established according to X-ray examination result conducted within last current year, and anamnestic data were collected. The control group consisted of healthy female individuals, matched by age (older than 18 years) and ethnicity with the RA patients, ensuring they had no family history of RA or other autoimmune diseases. Healthy controls underwent the identical laboratory examination as RA patient did, including specific autoimmune disease testing such as antinuclear antibody (ANA) tests, rheumatoid factor (RF), and anti-citrullinated protein/peptide antibody (ACPA) assays, and the acute phase reactants (C reactive protein and erythrocyte sedimentation rate). Those who tested positive were excluded from the study.

Exclusion criteria for study participants included the following: age less than 18 years, pregnancy, and lactation at the time of inclusion in the study, cancer at the study time and in anamnesis, severe concomitant somatic pathology in the decompensation stage, use of antibiotics, probiotics, antiviral medication in the last 3 months prior study, parasitosis at the study time or in anamnesis, acute respiratory disease at the study time, resection of the stomach and/or intestines, and inflammatory bowel disease in anamnesis. Additional exclusion criteria for control group were the presence of heredity for RA and other autoimmune diseases. The research study design is depicted using a STARD diagram (Fig. 1). This study employs data from the identical set of patients as presented in our earlier publication (Kozhakhmetov et al., 2023).

Fecal and blood sample collection, DNA extraction and sequencing

All patients included in the study were sampled in the form of feces and blood in compliance with infectious safety measures. Feces were collected using the DNA/RNA Shield Fecal Collection Tube (Zymo Research, R1101; Tustin, CA, USA). Blood sampling, for the subsequent determination of the levels of cytokines, chemokines, and immunoglobulins, was carried out from a vein into vacutainers with ethylenediaminetetraacetic acid. Biomaterial labeling was carried out according to the coding.

All activities for the study of the gut microbiome were carried out subject to the requirements of good laboratory practice (GLP). Samples were stored as recommended by Zymo Research, Tustin, CA, USA. All genomic DNA extraction steps from fecal samples were performed using the ZymoBIOMICS DNA Miniprep Kit (Zymo Research, D4300; Tustin, CA, USA) according to the manufacturer’s protocol (Zymo Research, Tustin, CA, USA). The DNA isolation underwent a qualitative check through electrophoresis in a 1% agarose gel. To assess the concentration and purity of each DNA sample, an Invitrogen Qubit 3.0 Fluorometer (Invitrogen, Carlsbad, CA, USA) was employed. A negative control was maintained using sterile water. The analysis of gut microbiome composition involved the sequencing of the V1-V3 region of the 16S ribosomal RNA gene. This sequencing was performed using the Illumina NovaSeq 6000 platform (Illumina, San Diego, CA, USA) at the Novogene laboratory in Beijing, China, in accordance with established Illumina protocols.

Analyzing sequencing data

The 16S amplicon sequencing data, starting from raw reads, underwent processing into taxon density tables through the utilization of the LotuS2 pipeline (Less OTU Scripts 2). Key steps in this process included demultiplexing, quality filtering, and dereplication of reads, which were carried out with the assistance of a straightforward demultiplexer (sdm). Additionally, chimeric sequences were identified and eliminated using UCHIME algorithms. Taxonomic postprocessing and sequence clustering, employing the SILVA database, were performed using LCA and UPARSE, respectively. Total number of reads accounted for 16,796,296, employing a similarity threshold of 97% for distance comparison. After filtering, the sequences were categorized into 3,681 OTUs and were attributed to the bacterial domain,15,176,348 reads remained in the matrix, an average sequencing depth was 41,990 reads per sample.

Analysis of metabolic pathways

The analysis of metabolic pathways was conducted using PICRUSt2 (Phylogenetic Investigation of Communities by Reconstruction of Unobserved States) version 2.5.0, following its default settings, to predict functional metagenomic profiles from 16S rRNA sequencing outcomes. To determine the gene family copy numbers for each operational taxonomic unit (OTU), a reference tree was used with a Nearest Sequenced Taxon Index (NSTI) threshold of 2. The abundance of metabolic pathways in bacteria was predicted using the MetaCyc Metabolic Pathway database.

Measurement of the immunological profile

A volume of 20 μL from the serum sample was utilized for immunological analysis. The measurement of cytokine levels using the MILLIPLEX MAP #1 Human Cytokine/Chemokine Panel 1 kit (Merck, Rahway, NJ, USA) for multiplex analysis via the Luminex xMAP method in accordance with the manufacturer’s standard protocol. The levels of inflammatory cytokines, chemokines, and immunoglobulins in fecal and plasma samples were measured using the Multiplex immunoassay. Multiple cytokine and chemokine assays and immunoglobulin isotyping were performed using the Bio-Plex 200 system (Bio-Rad, Hercules, CA, USA).

Statistical analysis

All statistical evaluation was done in Python 3.9.16. Operational taxonomic unit (OTU) tables underwent rarefaction, retaining only OTUs detected in a minimum of 25% of samples per group. Two-group comparisons utilized the independent T-test, Welch’s t-test, or Mann-Whitney U-test as considered fitting, using the “SciPy 1.10.1” library. For analyses involving multiple groups, the Kruskal-Wallis test followed by Dunn’s post-hoc test was employed, with the Benjamini-Hochberg FDR correction for multiple testing applied via “statsmodels 0.13.5”, maintaining a significance threshold of q ≤ 0.05. Alpha diversity metrics, including Shannon, Chao1, and Faith’s indices, and beta diversity assessments utilizing Bray-Curtis distance through ANOSIM and PERMANOVA with 9,999 permutations were executed as per “scikit-bio 0.5.6”. Principal coordinate analysis (PCoA) facilitated ordination of between-sample compositions. Spearman’s r correlation coefficient, adjusted for FDR, identified taxonomic features significantly related to ordination results. Differential analysis of taxonomic features was determined using “LEfSe 1.0.8”, applying a significance criterion of LDA > 2, p ≤ 0.05, consistent with methods reported by Issilbayeva et al. (2024). Differential analysis of functional features was based on differences in medians with 95 CI estimation using the Hodlges-Leghman estimator, considering only those with non-overlapping intervals as significant, following the methods described by Kozhakhmetov et al. (2023). Correlation analyses was performed using Spearman’s r, and feature importance analysis for more robust estimation was performed using a combination of gradient boosting decision trees (GBDT) and binomial logistic regression (LR) with leave-one-out cross-validation using “scikit-learn 1.2.2” and “LightGBM 3.3.5” libraries.

The demographic information and clinical serum biochemical data of study responders

Our study included 190 female participants, all Kazakh ethnicity, with 77 diagnosed with RA and 113 healthy individuals (as shown in Fig. 1). The average age of RA patients was 44 years, with the onset of the disease typically occurring in the third decade of life. The average disease duration amounted to 9 years. Most of the RA patients (81.8%) demonstrated from moderate to high disease activity.

Out of the 77 RA patients in the study cohort, 77.9% had RF positive and 22.1% RF negative RA form. While ACPA positive RA form accounted for 54.5% and ACPA negative RA form was 45.5% (Table 1).

The results of gut microbiome sequencing

The results revealed significant differences in both alpha and beta diversity between the control and RA groups. No significant difference was shown by the Shannon index (p = 0.4) between the two groups, while the α-diversity of microbial taxa, as assessed by Chao 1 index (p = 0.019) and Faith’s PD (p = 0.006) of the gut microbiota, was observed to be somewhat lower in patients with RA when compared to healthy controls, and these differences were deemed statistically significant (Fig. 2A). The assessment of beta diversity at the group level showed a significant dissimilarity in community comparisons according to PERMANOVA (F = 2.66, p = 0.0009), whereas ANOSIM (p > 0.05) did not reveal any significant large-scale distinctions (Fig. 2B).

The findings regarding the most abundant phyla indicated shifts in relative abundance. The analysis of the relative abundance showed decreased levels of Proteobacteria and Actinobacteriota in the RA group compared to the control group (Fig. 2C). A histogram depicting the structure of the gut microbiota is shown in Fig. 2D. The most abundant taxa on the family level among RA group were Lachnospiraceae, Clostridiaceae, Bacteroidaceae, Ruminococcaceae and Veillonellaceae. While in control group the most abundant taxa were Staphylococcaceae, Coriobacteriaceae, Rhizobiaceae, Sphinogomonadaceae, Puniceicoccaceae and Eubacterium Coprostanoligenes group. The most abundant genera among RA group were Dorea, Tyzzerella, Sellimonas which belongs to Lachnospiraceae family (Fig. 2E).

The connection between gut microbiome, RF and ACPA

We conducted an analysis to examine the correlation between microbial taxa and the presence of autoantibodies, such as ACPA and RF, as well as the levels of inflammation in the blood, quantified by measuring CRP, in these patients. The examination of alpha (Shannon/Chao1, Faith PD) and beta (Bray-Curtis, PERMANOVA) diversity, based on the presence of RF and ACPA antibodies, did not reveal any significant differences within the RA group (Fig. S1).

Based on LEfSe analysis focusing on RF status, it was observed that Eubacterium_ruminantium, and Monoglobus were linked with the RF-positive form of the disease, whereas Prevotella, Colidextribacter, and Butyricimonas were associated with the RF-negative form (Fig. 3A).

In relation to ACPA status, taxa like g_UCG_005, Odoribacter, and Ligilactobacillus were linked to ACPA-positive cases, while Faecalibacterium and Haemophilus were associated with the ACPA-negative disease form (Fig. 3B).

The study of immunological profile of the study groups

We conducted analysis to examine the relationship between the immunological profile and clinical parameters in RA. The findings revealed that RA patients with high disease activity exhibited elevated levels of MIG, IP10, FLT3L, IL6, TNFa, VEGF_A, IL17E_IL25, IL1a, GROa, and FGF2. The most substantial differences were observed in the chemokines MIG and IP10 (p < 0.01) (Fig. 4A). Increased levels of IL5, FGF2, and TNFb were noted in RA patients with RF positivity (Fig. 4B), while no significant differences were found in the analysis of the immunological profile and ACPA status (Fig. S2).

The analysis of correlations between taxa and cytokines-chemokines revealed various associations. For instance, Sellimonas showed a positive correlation with TGFa, Fractalkine, MIP1a, IL17a, IFNy, IL7F, FGF2, MCSF, and IL1b, but had negative correlations with IgM, IgG3, IgG4, IgA, and IgG2 (Fig. 5A).

Similarly, Clostridium_sensu_stricto_1 was positively correlated with GROa, MIP1a, IL17A, IL2, FGF2, MCSF and negatively correlated with IgG2. Ruminococcus displayed a positive correlation with MCSF, TNFb, IL5, and negative correlations with IgM, IgG3, IgG4, IgA, and IgG2. Additionally, g_Dorea had positive correlations with MCSF and IL5, while being negatively correlated with IgM and IgG3. Tyzzerella was positively correlated with MCSF and negatively correlated with IgA. Bacteroides had positive correlations with IL10 and MCSF.

Both Lachnospiraceae_NK4A136_group and Eubacterium_siraeum_group were negatively correlated with PDGF_AB_BB. Collinsella exhibited a positive correlation with IgG2, while Staphylococcus was positively correlated with IL13.

We also identified differentially abundant taxa in fact follow the same pattern as cytokines/chemokines levels, demonstrating a correlation with RA/RA forms. Notably, the strongest correlation is observed in taxonomic features associated with RA, while features associated with Control show a weaker association, although consistent across all such features (Fig. S3). Additional information on differentially abundant taxa distribution in RA/RA forms is provided in Fig. S3. In study group with RA, there is a notable prevalence of the sulfur oxidation super pathway, which demonstrates a positive correlation with proinflammatory cytokines such as IL17E, IL25, IL15, IL6, and MCSF. Additionally, the biosynthesis pathways for L-glutamate and L-glutamine exhibit positive correlations with IL17A, TNFb, IL5, and MCSF. Conversely, these metabolic pathways exhibit a negative correlation with IgG3 and PDGF_AB_BB (Fig. 5B). Furthermore, Mevalonate pathway I and geranylgeranyldiphosphate biosynthesis I display positive correlations with IgG2 and IgA. The flavin biosynthesis I pathway is positively associated with IL6 (Figs. S4, S5).

Analysis of lachnospiraceae and oscillospiraceae families

The family-level relative median abundance revealed an increased level of Lachnospiraceae and a decreased level of Oscillospiraceae in the RA group (Figs. 6A–6C). Key taxa within the Lachnospiraceae family included Dorea, Sellimonas, Lach._NK4A136. and Tyzzerella. Similarly, key taxon within the Oscillospiraceae family was identified as NK4A214 (Fig. 6E). Among these families, two specific species, Dorea and NK4A214_group, exhibited significant differences in abundance (Fig. 6D) and emerged as the most crucial predictors in the dataset (Fig. 6E).

When examining significant correlations between the relative abundance of Lachnospiraceae and Oscillospiraceae (LO) families and both immunological markers and clinical parameters, significant associations emerge (Figs. 6F–6J). Notably, Lachnospiraceae correlates with interleukin-5 (IL5), which, in turn, is associated with various clinical indicators. Oscillospiraceae, on the other hand, exhibited a negative correlation with interleukin-22 (IL22) (Fig. 6H), key taxa of both families Dorea and NK4A214 also demonstrated specific correlations with certain immunological parameters (Fig. 6G).

LO families display divergent correlations with various metabolic pathways. Specifically, two pathways demonstrate inverse associations over LO: lactose degradation I (Lactosecat-PWY) and mixed acid fermentation (Fermentation-PWY) pathways (Fig. 6H). Additionally, a strong positive correlation is observed between the Lachnospiraceae family and L-glutamate and L-glutamine biosynthesis (PWY-5505), super pathway of sulfur oxidation (PWY-5304), sucrose degradation III pathway (PWY-621) (Fig. 6H). The LO families also had the highest impact on prediction (Fig. 6F), and when used alone independently from other taxonomic features exhibited AUC values of 0.57 and 0.59 (Fig. 6J). Additional information on feature importance analysis for all genera is provided in Fig. S6.

Strengths and limitations of the study

The strength of the study lies in its thorough examination of the complex interplay among the microbiome, metabolic pathways, and immunological markers. To the best of our knowledge, this study signifies pioneering research on the gut microbiome in individuals with RA in Kazakhstan, and potentially Central Asia at large. It constitutes a significant contribution to the existing knowledge base regarding the research landscape of RA in this region. Future studies undertaken by researchers worldwide with a broader range of demographic groups are required. This approach will, probably, ensure our findings are applicable across different cultures and genetic profiles, thereby broadening their relevance and impact on a global scale. Our findings offer insights into the RA pathogenetic aspects and their potential diagnostic and therapeutic implications, however further research is needed to unravel the underlying mechanisms and relevance of these differences in RA. Future research should also widen its focus to include long-term studies on individuals with RA from the debut of their condition, looking at how different treatments affect the microbiome over time. This method will help us better understand the microbiome’s role in RA development and responds to treatment.

Certainly, there are several limitations in our study. The study emphasizes the necessity for larger-scale research to enhance statistical power and increase the robustness of the findings. The current sample size may not adequately represent the diversity within the RA patient population. The study’s findings are limited by the focus on patients with RA with long-term disease duration, potentially limiting the generalizability of the results to individuals in the early stages of the disease. All RA patients received the treatment before and at the time of the study, which included disease-modifying and symptomatic treatment, while the control group included healthy individuals without getting any drugs. The functional compartment of the microbiome was assessed via 16S rRNA sequence profile, and we assume that further metagenome sequencing is required, which facilitates a more comprehensive metabolic profiling of microbial communities by capturing the diversity of all related genes, enabling predictions about their behavior under different conditions, including RA, offering insights into how the microbiome influences the disease’s development and progression.