Background & Summary

Rheumatoid arthritis (RA) is a chronic, systemic autoimmune disease characterized by inflammation of the synovial joints, leading to joint destruction, disability, and decreased quality of life1,2,3,4. While the etiology of RA involves a complex interplay of genetic and environmental factors, recent studies suggest that dysbiosis of the gut microbiota plays a crucial role in the disease’s pathogenesis and progression. The gut microbiota, consisting of trillions of microorganisms, has emerged as a critical player in human health, influencing metabolic, immunologic, and even psychological processes. Alterations in the composition and function of this microbial community have been associated with a variety of diseases, including autoimmune disorders like RA1,5,6.

Studies on the relationship between the gut microbiome and RA using next-generation sequencing have been conducted to explore how alterations in the gut microbiota may influence the development and progression of RA1,5. However, the existing research has faced inconsistencies and controversies, mainly due to limited sample sizes, a lack of stratification by disease stage and treatment status, and insufficient depth of microbial profiling. Furthermore, the temporal dynamic changes in gut microbial communities related to RA onset, progression, and treatment remain poorly understood.

To address these knowledge gaps, we comprehensively analyzed the gut microbiota in a large, stratified cohort of 2,238 individuals, including 1,034 RA patients and 1,204 healthy controls. For the profiling of the bacterial communities, 16S ribosomal RNA (rRNA) V3-V4 amplicon sequencing was performed after fecal samples collection, DNA extraction, and library construction as described in the Methods section. The resulting sequencing dataset elucidated the role of gut microbiota in RA, spanning various stages of the disease, treatment regimens, and clinical outcomes.

The purpose of presenting our data is to communicate our research findings transparently, allowing for the validation of results and contributing to the collective knowledge in the field. We believe this approach leads to more comprehensive studies integrating diverse methodologies and perspectives, thus advancing our understanding of RA and its link with gut microbiota.

Methods

Study design

All participants were recruited from six centers, including Peking University People’s Hospital (Beijing 1) and Peking University Third Hospital (Beijing 2), Xinxiang Central Hospital in Henan, Shanxi Grand Hospital in Shanxi, Enshi Huiyi Hospital of Rheumatic Diseases in Hubei, and Peking University Shenzhen Hospital in Guangdong (Fig. 1A). It is important to note that certain facilities, such as those in Hubei and Beijing2, contributed only RA data and not HC data. This discrepancy is a potential source of bias that future researchers should be mindful of and account for in their analyses.

Fig. 1
figure 1

Overview of Cohort Composition and Data Processing in Gut Microbiota Study of RA patients (from APRAC) and HCs. (A) Pie charts represent the distribution of RA patients and HCs recruited from six centers. Each segment reflects the number of participants and their percentage from each location. (B) Bar graphs depict the cohort’s age distribution, categorized by gender and group (RA or HC). The x-axis indicates the number of individuals, while the y-axis lists the age ranges. (C) Pie chart displaying the age at disease onset for RA patients, with segments showing the number of patients and their corresponding age ranges. (D) Pie chart detailing the types and proportions of parameters documented by participants, categorized into demographic information, environmental factors, dietary habits, medication, clinical parameters, and extra-articular manifestations. (E) Pie chart illustrating the different extra-articular manifestations documented in the study, each segment indicates a system or organ involved and its percentage. (F) Flowchart summarizing the steps and tools used in data processing, from sampling and sequencing through bioinformatic analysis. Each icon represents a stage in the workflow, with a brief description underneath.

Full size image

The protocol was approved by the Peking University People’s Hospital Ethics Committee (Documentation ID: No. 2016PHB200), and written informed consent was obtained from all participating subjects for both the sharing of their data and their participation in the study. Participant information has been anonymised to protect their privacy. Personal identifiers have been removed, and data has been aggregated to prevent the identification of individual participants.

2,238 individuals, involving 1,034 RA patients and 1,204 healthy controls (HC), were enrolled in this study. Beijing cohort 1 registered the majority of RA (70%) and HC subjects (92%) (Fig. 1A). With a minority of subjects missing age information, data from 2,029 individuals showed that RA patients mainly consisted of female subjects aged between 40 to 69 (Fig. 1B). Of note, the female healthy controls enrolled in this cohort are young to middle-aged adults. We employed PERMANOVA to assess the relationships among age, gender, disease, and gut microbial composition. Our analysis revealed significant associations between these factors and microbial composition across various taxonomic levels. Importantly, even after adjusting for age and gender, we found that the disease remains significantly associated with the gut microbiota (Table 1).

Table 1 The effect size of metadata variables on gut microbial composition.
Full size table

Along with age, up to 91 demographic and clinical parameters were documented in this cohort (Fig. 1B–E). The metadata consists of seven demographic pieces of information (8%, such as gender, age, height, etc.), 22 clinical parameters (24%, including CRP, ESR, DAS28, etc.), 37 extra-articular manifestations involving six different systems or organs (41%, such as Oculopathy, thyroid disease, etc.), ten strategies of medication (11%, such as MTX, LEF, etc.), ten dietary habits (11%, involving milk, soybean, etc.) and five environmental factors (5%, including chemical reagent and radiation) (Fig. 1C–E). In the RA population, 85.6% of patients were positive for anti-CCP antibodies, and 70.3% were positive for rheumatoid factor.

In our study, 27% of patients were current smokers, while 4.5% had a history of smoking. Upon analyzing the association between smoking status and gut microbiota, we found no significant correlation between smoking and the bacterial Shannon index (r = -0.010, P = 0.770), observed OTUs (r = -0.024, P = 0.505), or microbial beta diversity (R² = 0.016, P = 0.875).

The datasets used in the current study are distinct and do not include any samples previously published.

Sample collection and DNA extraction

All fecal samples were collected using the same procedure in all centers. Briefly, fresh stool samples were frozen at -80 °C within 24 h of receipt. To avoid batch effects, all samples were transferred to a center (PKUPH) for the following process. Genomic bacterial DNA was extracted using the QIAamp DNA Stool Mini Kit (QIAGEN, Germany) according to the manufacturer’s recommendations. A stool mechanical disruption step with a bead-beater was performed per a previously described protocol7. DNA concentrations were determined using the Qubit® 3.0 fluorescent quantification kit (Thermo Fisher Scientific, Waltham, MA, USA). Extracted DNA was stored at -80 °C before sequencing.

16S rRNA gene amplification and sequencing

The V3-V4 hypervariable regions of the 16S rRNA gene were amplified. PCR reactions were performed using unique fusion primers designed based on the universal primer set, 338 F (5’-GTACTCCTACGGGAGGCAGCA-3’) and 806 R (5’-GTGGACTACHVGGGTWTCTAAT-3’), incorporating the Illumina adapters and a sample barcode sequence8. Triplicate PCR reactions were performed for each sample and visualized on 2% agarose gels (Thermo Fisher Scientific, Waltham, MA, USA). PCR amplicons were purified using the Agencourt AMPure XP kit (Beckman Coulter, Inc, Brea, CA, USA) and quantified by Qubit® 3.0 fluorescent quantification kit (Thermo Fisher Scientific, Waltham, MA, USA), then pooled at equimolar amounts as specified in the Illumina TruSeq Sample Preparation procedure. The amplicon library was constructed following the manufacturer’s instructions (Illumina, San Diego, CA, USA) and quantified using the KAPA Library Quantification Kit KK4824 (KAPA Biosystems, Woburn, MA, USA). The completed library was sequenced on an Illumina MiSeq (Illumina, San Diego, CA, USA) platform using a dual-index sequencing strategy according to the Illumina recommended protocol9.

The workflow for data processing is presented in Fig. 1F. In brief, 16S rRNA gene reads were compiled and processed using the Quantitative Insights Into Microbial Ecology (QIIME Version 2, http://www.qiime.org) pipeline10. Raw sequences were processed to concatenate reads into tags according to their overlapping relationship; then, reads belonging to each sample were separated with barcodes, and low-quality reads (Q-score < 20) were removed. After de-noising processing by DADA211, the clean tags were clustered de novo into amplicon sequence variants (ASVs) with a 100% similarity of merged reads. Instead of “feature”, the term “ASV” was used throughout the manuscript. ASVs were assigned to the different taxa by matching them to the Greengenes database (Release 13.8, https://greengenes.secondgenome.com)12 and chimeras were removed using the q2-feature-classifier plugin. The resulting rarefied ASV tables and their corresponding taxonomy profiles were used as input for downstream analyzes. All bioinformatic and statistical analyses were conducted in R software, version 4.1.3 (R Foundation for Statistical Computing, Vienna, Austria) unless otherwise stated.

Data Records

Raw fastq data

To share the original sequencing data for other researchers to verify, reprocess, and re-analyze according to their analytical needs and customized parameters, we deposited our raw 16S rRNA V3-V4 region amplicon data (FASTQ format) on the Genome Sequence Archive (GSA) platform13. This dataset is entitled “Fecal 16S rRNA genomics in patients with rheumatoid arthritis” and is available for download under GSA ID CRA003232 (https://ngdc.cncb.ac.cn/gsa/browse/CRA003232)14. The dataset contains 2,238 samples with 4,476 FASTQ files (two files per sample, representing forward and reverse reads). The pipeline used for the bioinformatic analysis can be downloaded through Figshare15. Other original data, such as metadata and taxonomic abundance datasets, are accessible through the R-data files tagged with the “.rds” extension.

Sample metadata

The demographic and clinical information of the study population is documented in the R-data file named “1.sample.info.original.rds”. As described previously, the metadata consists of seven demographic pieces of information, 22 clinical parameters, 37 extra-articular manifestations (EAMs) involving six different systems or organs, ten strategies of medication, ten dietary habits, and five environmental factors (Fig. 1C–E).

Taxonomic abundance files and quality report

To facilitate reanalysis for other researchers, we have shared the profiled ASV table, entitled “1.16S.ASV.profile.original.rds”, on the Figshare platform under the folder named “data”(https://figshare.com/articles/dataset/Data_for_publication_in_Scientific_Data/27603876/2)15. Additionally, we have provided the taxonomy data (“1.taxonomy.info.rds”) and the annotated bacterial ASV matrix with an abundance over 0.1% (“2-1a.tax_6Genus0.001_HC” and “2-1a.tax_6Genus0.001_RA”). The files named “stackplot” represent taxonomic data for the top 23 genera with highest abundance in both HC and RA groups and individuals. The files entitled “cuttree” were the clustering lists for two groups.

Technical Validation

Cohort control

Though individuals were enrolled from six different centers, operators in all centers followed parallel inclusion and exclusion criteria, strict documentation of epidemiological information, and fecal sample collection processes.

Sample and data

All fecal samples were collected and sent to the same center (PKUPH) for further processing. Bacterial DNA extraction was performed using the QIAamp DNA Stool Mini Kit (QIAGEN, Germany) according to the manufacturer’s recommendations, and the Thermo NanoDrop One instrument (Thermo Fisher Scientific, MA, USA) was used for controlling the quality and concentration of extracted DNA. Following amplification of 16S rRNA V3-V4 regions, PCR product quality was validated through 1% agarose gel electrophoresis and Qubit Fluorometer testing (Thermo Fisher Scientific, MA, USA). Raw sequencing data quality was checked using FastQC v0.11.916 and MultiQC v0.417 software, showing an average sequence count of 64,903 per sample, along with an average Q-score of 38, an average read length of 446 bp, and 52% GC content per sequence (https://figshare.com/articles/dataset/Data_for_publication_in_Scientific_Data/27603876/215, Fig. S2). QIIME (version 2) pipeline, q2-feature-classifier plugin, was used for assignment with referring to Greengenes database (Release 13.8). ASVs with a relative abundance of over 0.1% were kept. As a result, 4,567 ASVs involving 10 phyla, 16 classes, 26 orders, 57 families, and 204 genera were identified and included in further bioinformatic analysis.

Based on profiling, the “q2-vsearch” plugin in QIIME2 pipeline was used for de novo clustering, with a setting of 0.99 for the parameter of “p-perc-identity”. Consequently, 3 enterotypes of HC subjects and 5 RA patients were identified (Fig. 2A). In both HC and RA subjects, the genera Bacteroides, Roseburia, Escherichia-Shigella, and Prevotella_9 dominated the gut bacterial communities. Of note, though enterotype 1 in both HC and RA individuals was mainly constituted by Bacteroides, the abundance of Bacteroides was enhanced by approximately 6.5% in the RA group compared to the HC group, with an abundance of 27.0% versus 20.5%. Additionally, there was also a ~2.3% increase in the Ruminococcus gnavus group (3.76% versus 1.43% of abundance) in the RA group (34.6% versus 65.1% of individuals). In addition, the HC clusters were characterized by a lower abundance of specific enterotype drivers, such as Prevotella_9-enriched cluster 2 (15.8% of abundance and 20.1% of individuals) and Escherichia-Shigella-enriched 3 (16.7% of abundance and 14.6% of individuals). On the contrary, RA clusters were characterized by higher specific drivers, including Escherichia-Shigella-enriched cluster 3 (28.1% of abundance and 12.4% of individuals), Prevotella_9-enriched cluster 4 (18.8% of abundance and 12.8 of individuals), and Bifidobacterium-enriched cluster 5 (10.9% of abundance and 8.4% of individuals) (Fig. 2A). To infer the distribution of microbial data based on limited samples, we also performed a Kernel density estimation (KDE) with the “MASS” package18. The results also confirmed the data acquired from clustering analysis (Fig. 2B).

Fig. 2
figure 2

Community Structure of Gut Microbiota in RA and HC Groups. (A) The graphs display the relative abundance of the most prevalent genera in the gut microbiota of RA patients and HCs. Each bar represents an individual’s microbiota profile, with colors indicating different genera. Genera that do not fall into the top prevalent categories are grouped under ‘Others.’ The color-coded sidebar to the right of the plots indicates three major clusters in HC and five in RA, illustrating variations within the community composition. (B) The 3D PCoA plots illustrate the two groups’ overall gut microbiota composition and beta diversity based on genus-level profiles, using Bray-Curtis dissimilarity. The shape of the surfaces represents the density of points in the CPCoA space, with key genera labeled to show their relative influence on the separation between HCs and RA patients.

Full size image

Data sharing

The raw data concerning the 16s rRNA gene sequencing reported in this manuscript have been deposited in the Genome Sequence Archive at the National Genomics Data Center, Beijing Institute of Genomics (China National Center for Bioinformation), Chinese Academy of Sciences, under accession number CRA003232 that are publicly accessible at https://bigd.big.ac.cn/gsa.