Analysis workflow

The entire analytical process is summarized in Fig. 1 to explore the association between plasma metabolites, gut microbes, transcript levels, non-synonymous single nucleotide variants (nsSNV), and RA disease activity. The contents include the discovery cohort population distribution (Fig. 1A) and multi-omics analysis of RA with different disease activity levels (Fig. 1B). The random forest model performs receiver operator characteristic curve (ROC) analysis in the external validation cohort (Fig. 1C).

Fig. 1

Flow chart of experiment. A Population distribution for subjects in the discovery cohort. B Multi-omics analysis of RA with different disease activity levels. C ROC analysis of the external validation cohort in a classification model constructed through a random forest. RA, rheumatoid arthritis; HC, healthy controls; DAS28, Disease Activity Score 28; ITS, Internally Transcribed Spacer; LC–MS, liquid chromatograph mass spectrometer; MADAs, metabolites associated with disease activity; RA-seq, RA sequencing; WES, whole exome sequencing; DEGs, differentially expressed genes; nsSNV, non-synonymous single nucleotide variation; ROC, receiver operating characteristic curves

Changes of characteristic plasma metabolites in RA with different disease activity levels

To identify the differential metabolic characteristics among the four groups, we performed a non-targeted metabolomics test on fasting plasma of all participants. A total of 281 metabolites (cationic modes: 133, anionic modes: 148) were annotated in HMDB and Metlin databases. The selection of significantly different metabolites was determined based on the variable importance for the projection (VIP) value of OPLS-DA model and P value by Wilcoxon test; VIP > 1, p < 0.05 were defined as significantly different metabolites. Compared with the HC group, there were 28, 38, and 49 different metabolites in the three disease groups (Supplementary Fig. 1A-C). Comparison of DAS28M vs. DAS28L, and DAS28H vs. DAS28M with 8 and 9 metabolites is shown in Supplementary Fig. 1D, E. The Venn diagram showed the same and specific differential metabolites between groups, suggesting that as the disease activity increases, the differential metabolites also increase (Supplementary Fig. 1F, G).

Next, based on Spearman correlation analysis, heat maps showed the relationship between clinical parameters and characteristic difference metabolites (Supplementary Table 3). L-threonine, linoleic acid, deoxycholic acid, docosahexaenoic acid, 1,4-dihydroxybenzene, and L-tryptophan (A1-A6) have a negative relationship with inflammatory biomarkers (CRP, ESR, and IL-6) and DAS28 scores in anion mode (Fig. 2A). Their relative expression abundance in three disease groups decreased (Fig. 2B). However, the reverse occurred for D-galacturonic acid (A7) (Fig. 2A, B). MG (18:2(9Z,12Z)/0:0/0:0), 1-oleoyl-sn-glycero-3-phosphocholine (OGPC), 1-stearoyl-2-hydroxy-sn-glycero-3-phosphocholine (18:0LYSO-PE), and glycerophosphocholine (GPC) (C5-C8) also have a negative relationship with inflammatory biomarkers and DAS28 score (Fig. 2C). Their relative expression abundance in three disease groups decreased (Fig. 2D). We performed KEGG pathway enrichment analysis with differential metabolites to explore the changes in metabolic pathways in RA progression. The metabolites in the disease group were mainly enriched in glycerophospholipid metabolism and glycine, serine, and threonine metabolism pathways compared with HC groups (Fig. 2E–G). Meanwhile, differential metabolites between the disease groups were enriched in arginine biosynthesis and linoleic acid metabolism pathways (Fig. 2H, I). These results suggest that metabolic pathways change during RA progression, especially glycerophospholipid and linoleic acid pathways, which were defined as lipid metabolic pathways.

Fig. 2

Analysis of characteristic plasma metabolites in RA with different disease activity levels. A, C Heat maps of correlation between characteristic metabolites in anionic (A) and cationic (C) modes and clinical indicators of RA. Spearman was used for correlation analysis. Red represents a positive correlation and blue represents a negative correlation. B, D The box plot showed that characteristic metabolites of anionic (B) and cationic (D) modules significantly changed between RA with different disease activity levels according to Wilcoxon rank-sum test. Boxes represent the inter-quartile ranges, and lines inside the boxes denote medians. E–I Differential metabolites of DAS28L vs. HC (E), DAS28M vs. HC (F), DAS28H vs. HC (G), DAS28M vs. DAS28L (F), and DAS28H vs. DAS28M (I) were enriched in KEGG pathway. The color bars indicate enrichment significance levels. P-value, *p < 0.05, **p < 0.01, ***p < 0.001

Changes in gut microbiota composition among RA with different disease activity levels

Intestinal microbiome composition of all subjects was detected by 16S rRNA and ITS sequencing; we found that, at the phylum level of gut bacteria, with the increase of disease activity, the relative abundance proportion of Firmicutes decreased gradually, while Proteobacteria was the opposite (Supplementary Fig. 2A). Gut fungi were mainly composed of Ascomycota and Basidiomycota. Ascomycota was enriched in the disease groups, while Basidiomycota was enriched in the HC groups (Supplementary Fig. 2B). We visualized the species composition of bacteria and fungi at the genus level and found that other genera also showed various degrees of change (Fig. 3A, B). Candida in fungi increased significantly in the disease group, but there was no significant difference among the disease groups (Fig. 3C). Penicillium was reduced in the disease group, especially in the das28H group (Supplementary Fig. 2C). Lactobacillus is considered an intestinal probiotic, and its relative abundance in the disease group was significantly higher than that in the HC group, but gradually decreased among the disease groups (Fig. 3D). Escherichia-Shigella was significantly higher than HC groups in the disease group, and there was a significant difference between DAS28M, DAS28H groups, and HC groups (Fig. 3E).

Fig. 3

Changes of gut bacteria and fungi between RA with different disease activity levels. A, B The distribution plot of relative abundance at the genus level of bacteria (A) and fungi (B). C–E Relative abundance of Candida (C), Lactobacillus (D), and Escherichia-Shigella (E) in the HC group and RA with different disease activity levels groups. F Intestinal bacterial PICRUSt2 function predicts changes in linoleic acid metabolism pathway in RA progression. G, H LefSe analysis was used to identify highly differentiated taxa for RA gut bacteria (G) and fungi (H) with different disease activity levels and to display LDA scores. The bar chart shows the mean value of each group. Error bars represent the standard error of the mean values. RA, rheumatoid arthritis; HC, healthy controls; PICRUSt2, Phylogenetic Investigation of Communities by Reconstruction of Unobserved States; LefSe, linear discriminant analysis (LDA) effect size; P-value, *p < 0.05, **p < 0.01, ***p < 0.001. NS, not significant

We also studied changes of gut microbial diversity. By the Venn analysis of ASV in the four groups, intestinal bacterial species diversity decreased in the DAS28L group and then gradually increased. Meanwhile, we also found that the number of specific bacteria increased with the increase of disease activity (Supplementary Fig. 2D). We also found the same results in gut fungi (Supplementary Fig. 2E). In our results, only the alpha diversity indices of intestinal fungi Shannon even, Simpson, and Simpson even showed significant differences between the disease group and HC, while there was no difference among the disease group (Supplementary Fig. 2F-H). Then, we performed PICRUSt2 function prediction on bacterial 16S ribosomal ribonucleic acid amplicon sequencing data, and the results showed that linoleic acid metabolic pathway was significantly enhanced in DAS28M and DAS28H groups (Fig. 3F). However, glycerophospholipid metabolic pathway and arachidonic acid (AA) metabolic pathway were changed in the disease group, but there was no significant difference (Supplementary Fig. 2I, J).

To further study the potential differences between the HC group and RA with different disease activity levels, the LEfSe method was used to analyze the composition of intestinal microbes. At the bacterial genus level, Lactobacillus, an unidentified Prevotella, and Eubacterium were the bacteria with the highest LDA scores in DAS28L vs. HC DAS28M vs. DAS28L and DAS28H vs. DAS28M groups, respectively (Fig. 3G). Similarly, the fungi with the highest LDA scores were Candida, an unidentified Phaeosphaeriaceae, and Schizosaccharomyces, respectively (Fig. 3H).

Transcriptome analysis of RA suggests changes in lipid metabolism pathways

We sequenced RNA-seq data from a large number of RA patients with different disease activities and used principal component analysis to demonstrate that these disease groups were similar (Supplementary Fig. 3A). Compared with the HC group, the differentially expressed genes (DEGs) in different disease groups increased with the increase in disease activity (Fig. 4A). Meanwhile, they had 162 shared DEGs, including 11 upregulated genes and 151 downregulated genes (Fig. 4B). A total of 15 genes were associated with DAS28 score (Spearman, |R|> 0.1, P < 0.05), and their fold change was altered in different disease groups (Fig. 4C, Supplementary Table 4). Genes related to disease activity were mainly enriched in biological processes related to phosphatidylinositol 3-kinase signaling (Fig. 4D). The most affected pathway was the chemokine signaling pathway (Fig. 4E). Similarly, genes related to chemokine signaling pathway were significantly upregulated in RA by GSEA analysis (Fig. 4F). In addition, we also found that linoleic acid metabolism genes were significantly downregulated in RA (Fig. 4G). The same results were also observed for glycine, serine, and threonine metabolism; glycerophospholipid metabolism; and arachidonic acid metabolism as linoleic acid metabolism (Supplementary Fig S3B-D).

Fig. 4

Transcriptome analysis of RA suggests changes in lipid metabolism pathways. A Volcano plots showing DEGs in RA with different disease activity, respectively. Red indicated upregulated genes, and green indicated downregulated genes. B Venn diagram showing 162 DEGs in common across RA with different disease activity. C Fold change of DAS28-related genes in RA with different disease activity. D, E Bar graphs showing the top 10 most significant pathways of DAS28-related genes in Gene ontology term (D) and Kyoto Encyclopedia of Genes and Genomes (E) enrichment analysis. F, G GSEA analysis showing DEGs of RA vs. HC were upregulated in the chemokine signaling pathway (F) and downregulated in the linoleic acid metabolism pathway (G). H, J Expression of PLA2G6 (H), PTGS1 (I), and GDE1 (J) gene in lipid metabolism pathway in RA progression. The bar chart shows the mean value of each group. Error bars represent standard error of the mean values. RA, rheumatoid arthritis; HC, healthy controls; DEGs, differentially expressed genes; P-value, *p < 0.05, **p < 0.01, ***p < 0.001

Then, to demonstrate the changes of lipid metabolism pathways in the progression of RA, we further analyzed 127 genes of three lipid metabolism pathways and found 61 differential lipid metabolism-related genes (LMRGs) (FDR < 0.05, Supplementary Fig. 3E). Among them, 9 differential LMRGs were associated with differential metabolites of lipid metabolism pathway (3 up-regulated and 6 down-regulated). Compared with the DAS28L group, PLA2G6Z was significantly decreased in DAS28M and DAS28H groups (Fig. 4H). However, PTGS1 increased gradually with the increase in disease activity (Fig. 4I). Another LMRG, GDE1, was significantly increased only in the DAS28H group (Fig. 4J). Overall, transcriptional data analysis indicated that lipid metabolic pathways were altered during RA progression.

nsSNV of the HLA-DRB1 and HLA-DRB5 gene locus were associated with the disease activity of RA

We collected PMBC of 20 healthy subjects and 33 RA patients for WES analysis. A total of 97 nsSNV genes with mutation rate greater than 10% in exons were selected for protein–protein interaction (PPI) analysis on STRING with key protease genes in three lipid metabolism pathways. Finally, a total of 74 genes participated in the protein interaction network of lipid metabolism (Fig. 5A). KEGG and GO enrichment analysis of these 74 genes showed that 3 genes were enriched in RA, namely HLA-DRB1, HLA-DRB5, and CCL3L3 (Fig. 5B). Meanwhile, leukotriene D4 metabolic process and leukotriene biosynthetic process are also significantly enriched (Supplementary Fig. 4A). Per-Arnt-Sim Kinase has the largest number of connections in the network, followed by neutrophil solute factor 1 (NCF1) and ANKRD36C. NCF1 was the gene with the largest number of nodes interacting with lipid metabolism genes (Supplementary Fig. 4B). According to the statistical analysis of mutation rates of the six genes mentioned above in disease activity degree, the mutation rate of HLA-DRB5 rs1071748 was the highest, and the mutation rates of the three groups were 40.0% (DAS28L), 53.8% (DAS28M), and 60.0% (DAS28H) respectively (Fig. 5C). Meanwhile, some gene loci of HLA-DRB1 and HLA-DRB5 were only mutated in the DAS28M groups and DAS28H groups. Similarly, we also found that the mutation rate of NCF1 rs201802880 in the DAS28L groups (20.0%) was lower than in the DAS28M groups (53.8%) and DAS28H groups (40.0%). Therefore, we hypothesized that nsSNV of the HLA-DRB1 and HLA-DRB5 gene locus may be associated with RA with the disease activity.

Fig. 5

nsSNV of the HLA-DRB1 and HLA-DRB5 gene locus were associated with the disease activity of RA. A PPI analysis of genes with nsSNV mutation rate greater than 10% and genes involved in lipid metabolism pathways. The blue dots represent genes with RA PMBC nsSNV mutation rate greater than 10%, and the yellow dots represent genes associated with lipid metabolism pathways. The size of the point represents the node size of the gene in PPI. B KEGG pathway enrichment analysis bars of related genes in PPI analysis. Color bars indicate enrichment significance. C Mutation rates of nsSNV in related genes in different disease activity groups. nsSNV, non-synonymous single nucleotide variation; RA, rheumatoid arthritis; PPI, protein–protein interaction; PMBC, peripheral mononuclear blood cell

Multi-omics analysis revealed the relationship between gut microbiota, plasma metabolites, LMRGs, and exons nsSNV and RA disease activity

In Fig. 2, metabolites enriched in glycerophospholipid and linoleic acid metabolic pathways were negatively correlated with RA disease activity. To further study the effect of gut microbe on metabolites associated with disease activity (MADAs) in lipid metabolism pathway, the bacteria and fungi with differences between the HC group and disease group and disease activity-related metabolites were used for correlation analysis. There were 28 fungal genera and 13 bacterial genera were significantly correlated with MADAs, and they were also correlated among themselves (Spearman’s correlation analysis, p < 0.05, Fig. 6A, Supplementary Fig. 5). We used RA transcription data to demonstrate that genes related to lipid metabolism were altered in RA progression and correlated with disease activity. Meanwhile, the genes of nsSNV participated in the interaction network of key proteins in lipid metabolism. Therefore, we reasoned that plasma metabolites, gut microbe, and nsSNV together constitute a regulatory network for RA disease activity (Fig. 6B).

Fig. 6

Multi-omics analysis revealed the relationship between gut microbiota, plasma metabolites, LMRGs, and exons nsSNV and RA disease activity. A Association analysis of plasma metabolites and gut microbiota with RA disease activity. Red dots represent DAS28-ESR scores, yellow diamonds represent metabolites, cyan squares represent intestinal bacteria, purple triangles represent intestinal fungi, solid yellow lines indicate positive associations, and gray dotted lines indicate negative associations. All lines indicate significant associations (p < 0.05). B RA pathway change pattern diagram. Plasma metabolites, gut microbe, and nsSNV together constitute a regulatory network for RA disease activity. LMRGs, lipid metabolism-related genes; nsSNV, non-synonymous single nucleotide variation; RA, rheumatoid arthritis; PGE, prostaglandin E; LT, leukotriene; HETE, hydroxyeicosatetraenoic acid

Identification and prediction of RA with different disease activities based on plasma metabolites and gut microbiota

The random forest algorithm was used to select characteristic parameters for the differential metabolites, bacteria, and fungi of HC groups and three disease groups for classification model construction, and the area under the curve (AUC) value determines the number of important characteristic parameters of the model. In DAS28L vs. HC, DAS28M vs. DAS28L, and DAS28H vs. DAS28M of the three models, the top 5, 4, and 5 characteristic parameters of importance were selected according to AUC values respectively to build the model (Supplementary Fig. 6A-F). ROC analysis was performed on the three models in the discovery cohort, and their AUC values were 0.987 (0.942,1.000), 0.769 (0.440,0.994), and 0.790 (0.700,0.880), respectively (Fig. 7A). Meanwhile, we also included an external validation cohort with the same criteria as the discovery cohort. Finally, the three models were verified by an external validation cohort, and their AUC values were 1.000 (1.000, 1.000), 0.689 (0.531, 0.847), and 0.682 (0.534, 0.829), respectively (Fig. 7B).

Fig. 7

Diagnostic outcomes of RA disease activity are shown by the receiver operating characteristic curve (ROC). A ROC of three random forest models (DAS28L vs. HC, DAS28M vs. DAS28L, and DAS28H vs. DAS28M) in the discovery cohort. B ROC of three random forest models in the external validation cohort