Multi-omics Analysis of the Biomarkers and Molecular Mechanism of Rheumatoid Arthritis With Bone Destruction

Objectives Our study aimed to elucidate the role of metabolites, bacteria, and fungi in rheumatoid arthritis (RA) patients with bone destruction (BD(+)) and nd some biomarkers to predicate bone progression of RA. Methods We conducted plasma metabolites of the 127 RA patients and 69 healthy control by using nontargeted liquid chromatography-mass spectrometry (LC-MS), and the gut bacteria and fungi were assessed by 16S rRNA and internal transcribed spacer (ITS). Results Compared with RA patients without bone destruction (BD(-)), some metabolites, bacteria, and fungi altered in BD(+). Sever metabolites were selected as key metabolites for classifying the BD(+) and BD(-) groups with moderate accuracy (AUC=0.71). Metabolites-groups, metabolites-metabolites, and metabolites-clinical factors had a certain correlation, and 7 metabolites were enriched in glycerophospholipid metabolism and L-arginine and proline metabolism pathways. The bacteria and fungi of the BD(+) group showed signicant differences in composition and function compared with BD(-) group. The changed 4 bacteria and 12 fungi yielded accuracy (AUC=0.74 and AUC=0.87, respectively) for the two groups. Taken 7 metabolites, 4 bacteria and 12 fungi as a panel for AUC analysis, an improved AUC of 0.99 signicantly discriminated the two groups. The changed metabolites, gut bacteria, and fungi may affected the pathway related to L-arginine.


Introduction
Rheumatoid arthritis (RA) is a systemic disease mediated by autoimmunity and characterized by in ammatory manifestations of the joints and synovitis, as well as bone destruction caused by joint injury, cartilage destruction and osteoclast activation, which eventually leads to the destruction and even deformity of bone, cartilage and tendon [1][2][3]. In recent years, the research of metabolites and intestinal microbes has surpassed genetic factors and the environmental triggers in RA bone destruction (BD(+)), but the involved mechanisms of among metabolites, bacteria, and fungi is stilling unknown.
Metabolism is a general term for a series of ordered chemical reactions that take place in living organisms to maintain life [4] and is indispensable for the maintenance of life [5]. At present, metabolomics is an effective tool to identify biomarkers of RA and other diseases, Hiroshi Furukawa et al. have proved that metabolomic pro ling will be useful for discovering candidate screening biomarkers for interstitial lung disease in RA [6], Lun Zhang et al. successfully identi ed and validated a simple, highperforming, metabolite-based test for detecting early stage (I/II) non-small cell lung cancer patients in plasma [7].
In recent years, balancing the human gut microbes has been used as a powerful tool in the treatment of a variety of diseases [8][9][10]. Increasing evidence has shown that the composition and function of intestinal bacteria are closely related to autoimmune diseases [11]. Studies have demonstrated that intestinal ora imbalances can lead to the occurrence and deterioration of RA and a series of rheumatic diseases, Toshihiro Kishikawa et al. have proved that microbiome plays a important role in RA aetiology. Deshire Alpizar-Rodriguez also certi cated that Prevotella spp. enrichment in individuals in pre-clinical stages of RA, before the onset of RA, suggests a role of intestinal dysbiosis in the development of RA [12].
Fungi in the gastrointestinal tract are now recognized as a signi cant part of the gut microbes, and they may play an important role in human health [13]. Although there are many studies on fungi in humans [14][15][16][17], in contrast to bacteria, the characteristics of the fungi in the gastrointestinal tract of RA patients has never been described, and there is little knowledge about the fungus in the gastrointestinal tract of RA patients with bone destruction.
As far as we know, there is no any study reported a mechanism driving the imbalance of bone metabolism involving differences in metabolites, bacteria and fungi in RA patients. In our study, a systematic study was used to expound the role of the changed metabolites, bacteria, and fungi in bone destruction and the related mechanisms during the pathogenesis of BD(+).

Participant recruitment
Our study was approved by the Dazhou Central Hospital (Ethical review number: IRB-022). The participants in this study included BD(+), BD(-), and HC. The patients with RA in this study were diagnosed according to the 2010 American College of Rheumatology (ACR)-European League Against Rheumatism (EULAR) classi cation criteria. They were recruited from the Department of Rheumatology and Immunology, Dazhou Central Hospital, from November, 2017 to July, 2020. Two senior radiologists reviewed the X-ray lms of both the hands and wrists of the same patient to determine whether there was bone destruction. If the judgments were inconsistent, another senior radiologist made the nal decision after reviewing the lm again. Sixty-nine HC were recruited from the Department of Physical Examination at Dazhou Central Hospital in July, 2020. Finally, a total of 127 RA patients, including 91 BD(+), 36 BD(-), and 69 HC were considered for this analysis and had sufficient plasma for metabolomics analysis. Fecal samples of 50 BD(+) and 23 BD(-) were used to analyze gut bacteria and fungi by 16S rRNA and ITS sequencing (Fig. S1), and the participants' characteristics are shown in Table S1.
Plasma and fecal sample collection Each participant's whole blood was centrifuged at 3,500 rpm/min for 10 minutes rst, and the supernatant was moved to an enzyme-free 1.5-ml EP tube and centrifuged in a high-speed centrifuge at 12,000 rpm/min for 10 minutes. Then, the supernatants (plasma) was moved to another enzyme-free 1.5ml EP tube, frozen rapidly in liquid nitrogen, and stored at -80°C until extraction. The fecal samples were divided into 200 mg sections, frozen rapidly in liquid nitrogen, and then stored at -80°C until extraction.

Fecal sample DNA extraction and Illumina MiSeq sequencing
According to the manufacturer's protocol, total DNA was extracted from fecal samples. All DNA samples were quality checked, and the concentration was quanti ed by NanoDrop 2000 spectrophotometers (Thermo Fisher Scienti c, Wilmington, DE, USA). Bacterial 16S rRNA gene fragments (V3-V4) were ampli ed from the extracted DNA using the primers 338F ACTCCTACGGGAGGCAGCAG and 806R GGACTACHVGGGTWTCTAAT, and fungal ITS gene fragments (V3-V4) were ampli ed from the extracted DNA using the primers ITS1F CTTGGTCATTTAGAGGAAGTAA and ITS2R GCTGCGTTCTTCATCGATGC.

Statistical analysis
The processed data were analyzed by Pareto-scaled principal component analysis (PCA) and orthogonal partial least-squares discriminant analysis (OPLS-DA). The variable importance in the projection (VIP) value of each variable in the OPLS-DA model was calculated to indicate its contribution to the classi cation. Metabolites with a VIP value >1 was further included in a Student's t-test at the univariate level to measure the signi cance of each metabolite. The differentiation performance of selected classi ers to distinguish the BD(+) and BD(-) groups was quanti ed using receiver operating characteristic (ROC) curve analysis. DADA2-denoised sequences are usually called amplicon sequence variants (ASVs). The taxonomic assignment of ASVs was performed using the Naive Bayes consensus taxonomy classi er implemented in QIIME 2 and the SILVA 16S rRNA database (v138). ASV analysis, community diversity analysis, species diversity analysis, model predictive analysis and Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt2) function prediction analysis for the 16S rRNA bacteria and the ITS fungi sequencing data were performed using the free online platform microbiomeanalyst (https://www.microbiomeanalyst.ca/). All tests were performed using GraphPad Prism (v5.0) (GraphPad Software, Inc., CA, USA), SPSS Statistics (V.24.0.0.0) (SPSS Inc., Chicago, USA) or R software (Version 3.4.4).

Nontargeted LC-MS/MS and correction analysis
In this study, a total of 404 metabolites were identi ed among the BD(+), BD(-), and HC groups. Unsupervised PCA revealed some degree of separation between the metabolites of RA and HC, indicating that there was a difference in metabolic pro les between RA and HC. However, the PCA results revealed that it was di cult to separate the BD(+) and BD(-) groups (Fig. 1A, B). Seven metabolites were selected based on the VIP (VIP >1) and P-value (P < 0.1) as key metabolites for classifying the BD(+) and BD(-) groups ( Table 1). The heatmap analysis results of the relative intensity of metabolite showed that 6 downregulated metabolites negatively and 1 upregulated metabolite positively correlated with BD(+) group (Fig. 1C). The results of Spearman correlation analysis, which was used to identify potential metabolite-metabolite correlations, showed that 6 metabolites had a correlation (Fig. 1D). The correlations between the relative content of metabolites and clinical factors were analyzed, and the results shown that there were negative correlations among the three metabolites, L-arginine, 1-Oleoyl-snglycero-3-phosphocholine, and glycerophosphocholine, and Disease Activity Score 28-joint count erythrocyte sedimentation rate (DAS28-ESR), interleukin 6 (IL6), erythrocyte sedimentation rate (ESR) and C-reactive protein (CRP) (Fig. 1E) and the statistical analysis, among 7 metabolites, IL6, ESR, and CRP, was shown in Table S2.

Enrichment of altered metabolites and functional analysis
Pathway analysis of the 5 identi ed metabolites revealed signi cant enrichment in two pathways, including glycerophospholipid metabolism and L-arginine and proline metabolism ( Fig. 2A). A schematic diagram of the L-arginine and proline metabolism signaling pathway is shown in Fig. 2B, this pathway is related to L-arginine and creatine. Different levels of L-arginine and creatine were observed among the BD(+), BD(-) and HC groups (Fig. 2C-D). A schematic diagram of the glycerophospholipid metabolism signaling pathway is shown in Fig. 3E; this pathway is related to 1-Oleoyl-sn-glycero-3-phosphocholine, choline and glycerophosphocholine. Different levels of 1-Oleoyl-sn-glycero-3-phosphocholine, choline and glycerophosphocholine were observed among the three groups ( Fig. 2F-H).
Alterations in bacterial composition in RA patients based on 16S rRNA data and bacteria-based prediction In order to visualize the overall community structure of gut bacteria among fecal samples, Principal coordinates analysis (PCoA) was applied to visualize the overall community structure of gut micro ora among all fecal samples on the ASV level. There was no signi cant difference between the BD(+) and BD(-) groups based on PCoA at the ASV level (Fig. S2A, R 2 =0.011, P=0.788). Moreover, the rank abundance curve was analyzed, which re ected the richness and evenness of the species. These results suggested that although some genera differed in their abundance rank, the community composition was largely concordant between the two groups (Fig. S2B). A Venn diagram displayed 62 unique bacteria in the BD(+) group, 17 unique bacteria in the BD(-) group and 177 bacteria shared by both groups based on the genus level (Fig. 3A). Community diversity analysis was used to analyze the percent of community abundance of the two groups at the genus level. The results showed that these bacteria have subtle differences between the two groups (Fig. 3B). To explore the speci c bacterial taxa characterized in the two groups, linear discriminant analysis (LDA, score >2) effect size (LEfSe) analysis was applied.
The results showed that 10 bacteria were signi cantly changed, and Morganella levels was more abundant in the BD(-) group, while other bacteria were more abundant in the BD(+) group (Fig. 3C,  D). Then, the result of the Wilcoxon rank-sum test bar plot at the genus level showed that 8 bacteria were signi cantly different between the two groups (Fig. 3E). To test whether potential diagnostic biomarkers can be used to predict BD(+), we developed a random forest model based on the bacteria. Finally, the optimized model utilized 12 bacteria, which provided the best discriminatory power (AUC=0.74) ( Fig. 3F and Fig. S2C). Thus, the prediction model showed a high discriminatory power to predict the RA status. We predicted corresponding changes in modules and pathways using PICRUSt2 between the two groups. PICRUSt2 analysis suggested that these changes in the relative ASV abundance may be associated with the upregulation of pathways involving tetracycline biosynthesis; the mTOR signaling pathway; caffeine metabolism (Fig. 3G); and the regulation of modules involving polyamine biosynthesis, arginine => ornithine =>..., F420 biosynthesis, cholesterol biosynthesis, squalene 2,3-epoxide =>..., cationic antimicrobial peptide (CAMP) resistance and the NADH dehydrogenase (ubiquinone) 1 alpha subcomplex (Fig. 4H).
Analysis of the fungi in the BD(+) and BD(-) groups ITS was used to analyze the differences between the BD(+) and BD(-) groups. First, alpha diversity was used to obtain information about the abundance and diversity of species in the fecal samples, and the results showed that Shannon, Shannoneven, Simpson, and Simpsoneven indexes were signi cant difference between two groups (P < 0.05), while the Ace, Chao, and Sobs indexes were not different between two groups (P > 0.05) (Fig. 4A). Rank abundance curves evaluated relative fungi evenness , the result showed that all samples were similar patterns (Fig. S3A). To measure the extent of the similarity of the fungal macrobiotics, non-metric multidimensional scaling (NMDS) at the genus level was used to test the homogeneity of dispersion among different groups. Our results demonstrated that there was a signi cant difference between the two groups (Fig. 4B, R 2 =0.036, P=0.027). The composition of these fungi has subtle differences between the two groups, and largely made up of Candida, Aspergillus, and Debaryomyces (Fig. 4C). Moreover, the Venn diagram showed that 212 fungi were shared between the two groups, and 86 fungi were unique for the BD(-) group and 165 fungi were unique for the BD(+) group on the genus level (Fig. S3B). To identify fungal taxa that differed in relative abundances between the two groups, LEfSe analysis was performed. The results showed that there were some signi cantly differential fungi according to the following criteria: LDA score > 2 and P-value cutoff < 0.05 ( Fig. 4D and Fig. S3C ). The result of a Wilcoxon rank-sum test bar plot on the genus level also showed that 25 fungi had signi cant differences between the BD(+) and BD(-) groups, the top 15 fungi are shown in Fig. 4E. Then, random forest was used to select the important characteristic fungi. Finally, the optimal model utilized 20 fungi, which provided the best discriminatory power (AUC=0.83) (Fig. 4F and Fig. S3D). To assess the potential functional changes in the fungi, we predicted the abundance of pathways and enzymes using PICRUSt2 and compared the differences between the BD(+) and BD(+) groups. We found that the BD(+) group was associated with some pathways, including signi cant upregulation of L-proline biosynthesis II (from arginine) (Fig. 4G), and 2 enzymes (NAD(+)-protein-arginine ADP-ribosyltransferase and arginine kinase) were related to arginine (Fig. 5H).
The correction of metabolites, bacteria and fungi between the BD(+) and BD(-) groups and AUC analysis To explore the potential relationships among the metabolites, gut bacteria and fungi, the Spearman rank correlation coe cient was used to evaluate the correlation among the metabolites, gut bacteria and fungi. Combined with the previous analysis, we focused on the relationship between L-arginine and others. The results showed a strong positive correlation between L-arginine and glycerophosphocholine, and a strong negative correlation between L-arginine and 1oleoyl−sn−glycero−3−phosphocholine, choline, Roussoella, Hannaella, Debaryomyces (P < 0.05) (Fig. 5A). To further observe the predictive effect of metabolites, bacteria, and fungi on BD(+), we used the 7 identi ed metabolites as a panel, and the AUC was 0.71. Based on LEfSe, Wilcoxon rank sum test and random forest model analysis, we respectively selected 4 bacteria or 12 fungi as the research objects for AUC analysis, and the AUC were 0.74 and 0.85. Finally, 7 metabolites, 4 bacteria, and 12 fungi were combined as a panel for analysis, and an improved AUC of 0.99 signi cantly enhanced the power to discriminate between the two groups (Fig. 5B). Together, these results show that the changed of 7 metabolites, 4 bacteria, and 12 fungi has a certain ability to predict bone destruction in RA patients.

Discussion
With the rapid development of technology, microbiology and metabolomics have progressed greatly, and they have been widely used in various studies [18][19][20][21]. Microbial gene sequencing methods or metabolomics cannot directly clarify the mechanism of disease occurrence and development, and they also cannot be used alone to identify which members of the gut microbiota affect the host and other key problems. The disadvantages of single-omics studies have emerged, while the advantages of multiomics studies have gradually become clear. The purpose of this study was to elucidate the role of metabolites, bacteria, and fungi in BD(+) pathology and nd some biomarkers to predicate bone progression of RA. At present, although there are many studies on patients with rheumatoid bone destruction [22][23][24], only few studies have examined changes in metabolism or changes in intestinal bacteria in patients with bone destruction.
In our study, three novel metabolites were used to describe the metabolic pro le, gut bacteria and fungi change spectrum in the BD(+) group. Although nontargeted metabolomics and bacterial analyses did not clearly discriminate the BD(+) and BD(-) groups, some metabolites and gut microbiota changed in patients, which suggested that disrupts biochemical homeostasis may lead bone destruction. In our study, pathway analysis revealed that 7 dysregulated metabolites were potentially related to the metabolic pathways of arginine and proline metabolism and glycerophospholipid metabolism, at the same time,the changed gut bacteria, and fungi may affected the pathway related to L-arginine. Some studies have indicated that L-arginine and proline have effects on in ammation [25][26][27], and these research results are consistent with our experimental results.
After analyzing the relationships between the identi ed metabolites and clinical parameters, we concluded that three metabolites, 1-Oleoyl-sn-glycero-3-phosphocholine, glycerophosphocholine and Larginine, related to the severity of RA patients might represent a panel of potential small molecule biomarkers for assessing the severity of RA patient. Julia S. Brunner found that L-arginine plays an extremely important role in the process of bone destruction in RA and concluded that a high content of Larginine promotes bone destruction in patients with RA based on animal experiments [28]. Interestingly, in our study, we found that RA patients had a higher content of L-arginine than HC, but when we compared the content of L-arginine between the BD(+) and BD(-) groups, we found that the BD(+) group had a lower content of L-arginine. This may be because L-arginine was consumed excessively in the process of bone destruction, so the content of L-arginine was reduced.
Glycerophospholipids can be hydrolyzed to produce lysophospholipids (LPs). LP was initially regarded as a common intermediate in the synthesis of phospholipids. However, later studies showed that LP exhibits biological characteristics similar to extracellular growth factors or signaling molecules [29][30][31]. Some studies have shown that important LPs, including lysophosphatidylcholine (LPC), lysophosphatidic acid (LPA) and certain sphingomyelins, can participate in disease processes, such as atherosclerosis, vascular dementia, and spinal cord injury, by activating the PPARγ pathway [32][33][34]. At present, a large number of studies have proven that there is a clear correlation between glycerophospholipids and in ammation [35][36][37]. In our study, the glycerophospholipid signaling pathway was changed in BD(+) patients, which showed that the glycerophospholipid signaling pathway plays an important role in the process of bone destruction. The results of this experiment are consistent with those of other previous studies, which proves that our experimental results have a certain degree of reliability.
LEfSe bar analysis, the Wilcoxon rank-sum test and random forest analysis were used to compare the changes in intestinal bacteria and fungi between the BD(+) and BD(-) groups. Finally, 4 species of bacteria and 12 fungi were selected as markers to predict the progress of bone destruction. AUC analysis was used to analyze the metabolites, bacteria and fungi of single or combined changes to distinguish the BD(+) and BD(-) groups. These results showed that the AUC of 7 metabolites was 0.71, the AUC of 4 bacteria was 0.74, the AUC of 12 fungi was 0.87, and the combined AUC was 0.99, which suggested that the combination of metabolites, bacteria and fungi was more effective for the evaluation of bone destruction of RA.
When studying the potential relationships among metabolites, bacteria and fungi, we found that there were positive and negative correlations among metabolites, bacteria and fungi. After reviewing previous studies, we focused on the relationship between arginine and metabolites, microbes and fungi. The results proved that the metabolites choline and 1-Oleoyl-sn-glycero-3-phosphocholine were negatively correlated and glycerophosphocholine was positively correlated with L-arginine, and the fungi Debaryomyces, Hannaella, and Roussoella were negatively correlated with L-arginine. However, when examining the relationship between L-arginine and bacteria, we found no obvious correlation.
Combined analysis of the changed pathways among the three omics' analyses revealed a common pathway involving L-arginine. These results indicate that L-arginine plays an important role in the BD(+) group, but the speci c mechanism by which metabolites, bacteria and fungi affect the L-arginine pathway is not clear. We should perform further research on the mechanism of L-arginine-related pathways in subsequent experiments.

Conclusions
Collectively, our nontargeted LC-MS/MS, 16S rRNA and ITS analysis study of RA patients suggests some disrupted metabolites, bacteria and fungi form a distinct functional pro le in RA patients with bone damage. This study is the rst step in assessing the bone progression of RA patients by using metabolic and microbiological methods. It also has some certain limitations. First, because our study cohort did not include a validation cohort, it was impossible to accurately determine the prognostic value of the identi ed potential metabolic biomarkers for bone progression in RA patients. Subsequently, since there were no plasma or fecal samples from RA patients before bone destruction in this study, we were unable to examine the changes in the concentrations of metabolites, bacteria and fungi over time, which may be predictive of disease progression, treatment response or clinical outcome. In addition, we did not perform cell-based experiments or animal experiments to clarify the speci c role of the identi ed metabolites, bacteria and fungi in the BD(+) group. Therefore, we need to conduct a large-scale external cohort study to verify the utility of the identi ed potential biomarkers and basic experiments to clarify the speci c role of the identi ed metabolites, bacteria and fungi in the BD(+) group. The objective is to perform multiomics analysis and study the mechanism of the occurrence and development of BD(+) and provide new treatment ideas and methods for the clinical treatment of RA patients.   Figure 1 Metabolite analysis of serum sample among HC, BD(+) and BD(-) groups. A. PCA analysis of anion metabolites. B. PCA analysis of cation metabolites. C. Expression of differential metabolites in BD(+) group and BD(-) group. D. The correlation analysis of differential metabolites. E. The correlation analysis between metabolites and clinical factors. * 0.01 < P ≤ 0.05, ** 0.001 < P ≤ 0.01. HC: healthy control. BD(+): RA patients with bone destruction. BD(-): RA patients without bone destruction. PCA: Pareto-scaled principal component analysis.

Figure 2
Various methods were used to analyze the bacterial on genus level between BD(+) and BD(-) groups. A.