Structural Characteristics and Functional Analysis of Gut Microbiome in Patients with Osteoarthritis

Objective Recent studies illustrated that changes of gut microbiome in human body may be related to the pathogenesis of osteoarthritis (OA). We aimed to determine the difference of gut microbiome between patients with osteoarthritis and normal people. Methods The samples were collected from 56 patients with osteoarthritis (34 females) and 52 healthy subjects (28 females). High-throughput sequencing of 16S rRNA gene was used to analyze the classication and composition of gut microbiome in OA patients, and then we used Linear discriminant analysis Effect Size (LEfSe) to nd out the main differential bacteria and biomarkers in OA patients. At last, PICRUSt was applied to predict the genomic function of gut microbiome in OA patients, and we explored its mechanism in the development of the OA. Results


Introduction
Osteoarthritis (OA) is a common musculoskeletal disease. Recently, it has been reported that the prevalence of osteoarthritis is increasing, both in Europe and China, and the prevalence rate of osteoarthritis in China can reach 8.2% [1,2]. The progress of the disease is slow, usually accompanied with pain, joint swelling, deformity and other symptoms. In the global burden of disease study (GBD), osteoarthritis is listed as the second fastest rising disease related to disability, second to diabetes [3]. At present, the treatment of osteoarthritis is to relieve pain and delay the progress of the disease in clinical work, while for patients with severe osteoarthritis, surgical treatment can be considered when drug therapy can not fully alleviate the symptoms. But so far, the etiology and pathogenesis of osteoarthritis has not been clari ed in detail. It is generally believed that osteoarthritis is the result of multiple factors, such as age, gender, joint injury, heredity, or life factors(obesity, diet, work, smoking) [4,5].
In 2019, an European expert consensus mentioned that the increase in the prevalence of osteoarthritis is not only due to the extension of human life expectancy, but the changes in modern lifestyle, especially the lack of physical exercise and taking in too much sugar or saturated fatty acids [6]. In addition, they suggested that there is a potential link between osteoarthritis and gut microbiome. This provides us with a new entry point for the study of osteoarthritis. It is reported that the number of gut microbiome can reach 100 trillion, which is 10 times the total number of normal human somatic cells and germ cells, and far exceeds the number of other microorganisms on the surface of the human body [7]. In general, intestinal microorganisms maintain a balance between the human body and the external environment.
However, once this balance is broken, it will lead to the imbalance of community structure of gut mircobiome, resulting in metabolic disorders and low-grade in ammatory response in the body, both of which are important components of the pathogenesis of OA [6,8]. Studies have con rmed that gut microbiome can produce a variety of bioactive products, such as short-chain fatty acids, proteins and some enzymes, some of which are secreted in vesicles, and these vesicles are not sensitive to proteases, besides, some of these products can affect the permeability of mucosal epithelium, such as lipopolysaccharide (LPS). As a result, the components of Gram-negative bacteria (LPS) may enter the human blood circulation, increasing the concentration of lipopolysaccharide (LPS) in the blood circulation [9]. In addition, speci c dietary interventions may be associated with systemic, low-grade in ammation [10]. Dennis found that a controlled diet can reduce the incidence of osteoarthritis in a 15years animal experiment [11]. Besides, K.H.Collins found that the abundance of Lactobacillus species and Methanobrevibacter in the feces of the experimental group with high glucose and fat was strongly correlated with the improved Mankin score of articular cartilage, and the score could be predicted by the change of bacterial abundance [12]. Lee also reported that the ratio of Phylum Firmicutes in the stool of patients with osteoarthritis was signi cantly higher than that in the rheumatoid joint group [13]. Research has demonstrated that abundance of Streptococcus species is associated with osteoarthritis-related knee pain [14]. These ndings suggested that there is a potential link between gut microbiome and osteoarthritis.
In this study, we recruited knee osteoarthritis patients and healthy subjects, performed 16S rDNA sequencing on their gut ora, and compared the gut ora of the two groups through bioinformatics methods. In addition, we also explored the relationship between changes in the gut ora and knee osteoarthritis, and predicted the function of each ora in patients with knee osteoarthritis, hoping to discover the role of gut ora in knee osteoarthritis.

Study Design and Participants
This study recruited patients with knee osteoarthritis and healthy subjects from 2018 to 2019. All participants were from the Orthopedic Department of the First A liated Hospital of Zhengzhou University. At the same time, the body mass index (BMI) of the two groups was calculated, the BMI of OA group was (17)(18)(19)(20)(21)(22)(23)(24)(25)(26)(27)(28), and that of healthy control group was (16)(17)(18)(19)(20)(21)(22)(23)(24)(25)(26)(27)(28). The diagnostic criteria of knee osteoarthritis are detailed in the table (Table 1 ). There is a bone friction sound (sense) during the activity * All of the two major criteria and at least one of the three minor criteria are required for diagnosis.

Inclusion and Exclusion Criteria
The clinical data, use of medication and habits of all participants were from the hospital's electronic medical records and questionnaires. All participants should meet the following criteria: (1) No drug or alcohol abuse; (2) No antibiotics, prebiotics, probiotics or hormones were taken in the one months before stool samples were collected; (3) No special diet. In addition, the following situations should be excluded in both groups: (1) History of digestive tract diseases or gastrointestinal surgery; (2) Metabolic diseases (such as diabetes, hyperthyroidism, hyperlipidemia, etc.), infectious diseases, immune system diseases or severe systemic diseases, such as cardiac insu ciency, renal insu ciency; (3) Regularly take nonsteroidal anti-in ammatory drugs, opioids or proton pump inhibitors within 3 months (Among the subjects we recruited, the use of PPI and NSAIDs was generally intermittent, so the observation time was extended. And if the subjects took the drugs several times within three months, they should also be excluded). Besides, all subjects were told not to smoke, drink alcohol or take any medication within 1 week before collecting stool samples. In the end, 56 patients with osteoarthritis and 53 healthy subjects were included in the study (in the follow-up experiment, one case was removed because of pollution, 52 cases were nally included in the control group). Each participant signed an informed consent form that was approved by the Institutional Review Board of the First a liated Hospital of Zhengzhou University.
Collection of fecal samples, extraction of sample DNA and PCR Each subject was required to use a speci c fecal collection device to collect fresh fecal samples in the morning. At rst, it was stored in -20℃ in the hospital, and then it was stored in -80℃ in the laboratory within two hours.
According to the manufacturer's instructions of QIAamp 96 PowerFecal QIAcube HT Kit (QIAGEN), the DNA of the samples were extracted. Quality and quantity of DNA was veri ed with ultraviolet spectrophotometer (NanoDrop) and agarose gel. Extracted DNA was diluted to a concentration of 1 ng/μl and stored at -20°C until further processing. The diluted DNA was used as template for PCR ampli cation of bacterial 16S rRNA genes with the barcoded primers and Takara Ex Taq (Takara). The V3-V4 region primers:343F, 5'-TACGGGAGGCAGCAG-3'; 798R, 5'. AGGGTATCTAATCCT. 3'. Amplicon quality was visualized using gel electrophoresis, puri ed with AMPure XP beads (Agencourt), and ampli ed for another round of PCR. After puri ed with the AMPure XP beads again, the nal amplicon was quanti ed using Qubit dsDNA assay kit (Life Technologies). Equal amounts of puri ed amplicon were pooled for subsequent sequencing.The library was constructed using the lllumina platform and sequenced by an Illumina MiSeq PE300 sequencer.

OTU Clustering and Species Annotation
Illumina MiSeq sequencing generates the original double-terminal sequence, called raw data. First, the low-quality sequence was removed by Trimmomatic software [15], and then we used Flash software to assemble quali ed double-ended raw data in the previous step [16]. After that, split_libraries software in QIIME was used to remove the sequences containing N bases and sequences with a length less than 200bp [17]. Finally, the sequence was obtained after removing the chimera by UCHIME software and eventually got the data, which was used to divide the OTUs [18].
OTU, Operational Taxonomic Unit, is to divide the sequences into many groups according to their similarity greater than or equal to 97%. Vsearch software was used for OTU classi cation, and then we selected the most abundant sequence in each OTU as the representative sequence. In order to obtain the species information of each OTU, we compared and annotated the representative sequence with the Silva database through RDP classi er Naive Bayesian algorithm [19,20].

Comparison of Gut Microbiome between two groups
First, we constructed a Rarefaction Curves to analyze whether the sequencing depth of the experimental samples was reasonable. Then, the abundance and evenness of OTUs in different groups were compared by Alpha diversity analysis, including Shannon-Wiener Index, Simpson's diversity Index, Chao1, Specaccum species accumulation curve and Rank Abundance curve. In addition, Beta diversity analysis (PCoA analysis, PCA analysis and NMDS analysis) were applied to explore whether there were signi cant differences in ora composition between the two groups.
According to the Wilcoxon rank-sum test, species with signi cant differences between the two groups at different levels were obtained, and then the biomarkers which could effectively distinguish the osteoarthritis group from the normal group were found according to the Linear discriminant analysis Effect Size LEfSe , and the validity of the biomarkers was veri ed by Indicator analysis, Random Forest model analysis and 10-fold cross validation.

Function Prediction
Data obtained by sequencing was analyzed by PICRUSt, and then the gene functional spectrum was inferred based on the Greengenes database. Finally, compared with the KEEG database and COG database, the function of gut microbiota was predicted.

Statistical Analysis
Comparison of the clinical characteristics between the OA group and the healthy control group (continuous data that meet the normal distribution and homogeneity of variance, use two-tailed Student's t test, or use the Mann-Whitney U test when the conditions are not met). Using RDP classi er Naive Bayesian algorithm to complete the alignment and annotation between the representative sequence and the database. Statistics of OTUs and differences at different levels (phylum, class, order, family, genus, species) were performed by Wilcoxon rank-sum test. Using Adonis and Anosim analysis to test whether there were signi cant differences between two groups. LEfSe was used to nd out the main differential biomarkers. Indicator analysis, Random Forest model analysis and 10-fold cross validation were used to test the effectiveness of biomarkers. All the statistical analysis of data was carried out with the R packages, and p values less than 0.05 were considered signi cant. Spearman's correlations were used to associate the main differential bacteria  peculiar OTUs in the OA group (361) was greater than that in the HC group (284). (Fig. 1a).
In the experiment, the Rarefaction Curves of the samples are gradually at, indicating that the sequencing depth of the experimental samples was reasonable. In addition, the Good's Coverage index of each sample is more than 97%, which also shows that the sequencing has covered almost all the species in the sample. (Fig. 1b, Table 3) According to the composition of microbes at different levels (phylum, class and family), we draw histograms about the composition of the gut microbiome in two groups (Fig. 2a, 2b, 2c). And use box plots to show microorganisms with signi cant differences at each level (Fig. 2d, 2e, 2f).

Differences in abundance and composition
Alpha diversity analysis can re ect the abundance and evenness of microbial community. groups. (Fig. 1c, Table 3) Beta diversity analysis is to compare the composition of microbial communities in different samples. In this study, based on Weighted Unifrac and Unweighted Unifrac distance algorithms, PCoA analysis and NMDS analysis were used to compare the differences of microbial community composition between OA group and HC group. Through PCoA analysis and NMDS analysis, it can be seen that the OA group and the HC group have different aggregation tendencies, suggesting that there are differences in the composition of community between two groups. Besides, through Anosim analysis and Adonis analysis, based on four different distance algorithms, the difference of community composition between two groups was statistically analyzed. The results showed that there was signi cant difference in community composition between two groups. (Fig. 3a, 3b, Table 4) In this study, Linear discriminant analysis Effect Size (LEfSe) was used to nd the Biomarker with statistical difference between the two groups. As can be seen from the gures, Actinobacteria, Bi dobacteriaceae, Bi dobacterium and Alistipes were signi cantly enriched in osteoarthritis patients. While the abundance of Prevotellaceae and Faecalibacterium was lower than that of the healthy group. Combined the main difference bacteria at the genus level obtained by the Wilcoxon rank-sum test, it was found that Bi dobacterium, Alistipes and Faecalibacterium could be used as biomarkers for osteoarthritis (Fig. 4a, 4b, Fig.S1). In addition, Index analysis and random forest analysis con rmed the effectiveness of biomarkers (Fig. S2)

Functional prediction
Combined with the KEGG database, the 16S sequencing data were analyzed, and the results showed that the number of genes which were related to lipid metabolism, glycan biosynthesis and metabolism involved in gut microbiome in OA group was signi cantly lower than that in HC group at level 2 KEGG signaling pathway (p = 0.049, p = 0.005), and in the level 3 KEGG signal transduction, the enrichment of OA group in signal pathways such as glycan biosynthesis and metabolism, lipopolysaccharide biosynthesis and Adipocytokine signaling pathway were also signi cantly lower than that of the control group (p = 0.009, p = 0.001, p = 0.012), but in signal pathways related to apoptosis, the degree of enrichment was signi cantly higher than that of healthy people (p = 0.006).
In addition, the sequencing data were also analyzed by combining with COG database. The results showed that the proteins related to Na + -transporting NADH: ubiquinone oxidoreductase ( subunit NqrA, subunit NqrB and subunit NqrF ) of OA group were signi cantly lower than that of the control group (p<0.001, p<0.001, p<0.001).

Discussion
With the exploration of gut microbiome, it has been found that the richness and community structure of gut microbiome are associated with many diseases [21,22,23]. Clinical studies have con rmed that the gut microbiome can participate in bone metabolism by regulating the immune system and endocrine system [24,25]. Alessandro reported that gut microbiome may regulate human joints through the "Gut-Joint Axis" [26]. Boer has demonstrated that abundance of Streptococcus species is associated with osteoarthritis-related knee pain14. However, the effects of gut microbiome on the development of host osteoarthritis have not been discussed in detail. Therefore, in this study, bacterial 16SrRNA gene sequencing was used to analyze the diversity of gut microbiome of OA group, in order to explore the potential relationship between gut microbiome and osteoarthritis.
In this experiment, Alpha diversity associated indexes (including Chao1, Shannon-Wiener Index, Simpson's diversity Index) and Rank Abundance curve showed that Chao1, Shannon index and Simpson index in OA group were higher than those in HC group showed that there were differences in community composition between the two groups, combined with Anosim analysis and Adonis analysis, we concluded that the difference is statistically signi cant. Through the method of bioinformatics, we demonstrated the differential bacteria between patients with osteoarthritis and normal people at each level (phylum, class, order, family, genus), and then we found the main differential bacteria in the gut microbiome of patients with osteoarthritis by LEfSe. Combined with the signi cantly different bacteria from the Wilcoxon rank-sum test, we considered that Bi dobacterium, Faecalibacterium and Alistipes may be used as biomarkers. In addition, we con rmed the possibility of biomarkers by Indicator analysis and Randomforest analysis.
With the development of epidemiology and basic research, studies have reported that the low-level in ammatory and metabolic disorders may play important roles in the progression of osteoarthritis [27]. And the imbalance of gastrointestinal micro ora may be one of the main causes of obesity-related, lowlevel in ammation [28]. In our study, we also found that the enrichment of lipid metabolism, glycan biosynthesis and metabolism related genes in the OA group was signi cantly lower than that in HC group. In addition, some studies have found that the abundance of Bi dobacterium, Faecalibacterium were related to different markers of low-grade in ammation [29]. And Secretion of Faecalibacterium participates in the activation of NF-κB signal pathway [30]. As for NF-κB signal pathway, it is related to the in ammation that occurs when host cells come into contact with microorganisms, and more studies have reported that the gut microbiome can use a variety of strategies to regulate NF-κB signaling pathway [31,32].
In our experiments, we found that genera Lactobacillus abundance and family Lactobacillaceae abundance were signi cantly higher than those of healthy people. Similarly, Imhann found that in different cohorts, the use of PPI was associated with the increase of order Lactobacillales, family Streptococcaceae, Lactobacillaceae and genera Lactobacillus, Streptococcus [33]. In a large healthy twin cohort, Jackson con rmed that there was an increase in Lactobacillales, particularly Streptococcaceae, in PPI users [34]. However, in our study, the abundance of family Streptococcaceae and genera Streptococcus was not signi cantly different from that of healthy people, and the increase in the family Lactobacillaceae abundance was driven by Lactobacillus. Besides, the a-diversity was not signi cantly different between the two groups in our study, which was different from the use of PPI.
Although the users of PPI have been screened out, the time frame of our study essentially limits our ability to fully describe the connection between the observed osteoarthritis group and its related treatments. And non-time-matched questionnaires also brought additional interference to the dataset. We cannot rule out whether the frequency and dosage of PPI use will affect the distribution of the intestinal ora of patients in a longer period of time. Therefore, more rigorous experimental design is needed to prove the in uence of genera Lactobacillus on OA.
In summary, this study found that there are some differences in the distribution of gut microbiome between OA patients and normal people. In additon, Bi dobacterium, Faecalibacterium and Alistipes can be used as biomarkers in patients with osteoarthritis. Besides, gut microbiome may induce the development of osteoarthritis through metabolic pathways and lipopolysaccharide biosynthesis. These ndings provide a new theoretical basis for the study of osteoarthritis, and provide a new target for the prevention and treatment of osteoarthritis.

ETHICS STATEMENT
The studies were approved by Ethics Committee of the First A liated Hospital of Zhengzhou University and all patients/participants signed an informed consent. show the difference in species richness between samples, and can also be used to assess whether the sequencing depth of the samples is reasonable. (c) Rank Abundance is used to explain the richness and uniformity of the species contained in the sample. The richness of species is re ected by the length of the curve on the horizontal axis, the wider the curve (the larger the span of the horizontal axis), the richer the composition of species; the uniformity of species composition is re ected by the shape of the curve, the atter the curve (the smaller the span of the vertical axis), the higher the uniformity of species composition.

Figure 2
Histogram of community composition. Differential species between the two groups at the Phylum, Class, and Family. The left side is the main bacterial species at each level in the two groups a, b, c , and the right side is the main differential bacteria between two groups (d, e, f).

Figure 3
Differences in abundance and composition between OA group n=56 and HC group n=52 (a) PCoA Principal Co-ordinates Analysis. Based on Weighted Unifrac distance and Unweighted Unifrac distance, we performed PCoA analysis to extract the most important elements and structures from multidimensional data, each point in the graph represents a sample, and the samples of the same group are represented by the same color. (b) NMDS, Non-Metric Multi-Dimensional Scaling. NMDS is a nonlinear model, which is designed to overcome the shortcomings of linear models (including PCA and PCoA) and better re ect the nonlinear structure of ecological data. Each point in the graph represents a sample, and the distance between points indicates the degree of difference. The samples of the same group are represented by the same color, which can re ect the differences between groups and within groups through the distance between points.

Figure 4
The main differential bacteria (a) In the cladogram, the circle radiating from the inside to the outside represents the classi cation level from the phylum to the genus. Each small circle represents a classi cation at that level, and the diameter of the small circle is proportional to the relative abundance size. Different colors indicate different groups, red nodes indicate microbial groups that play an important role in the OA group, and blue nodes indicate microbial groups that play an important role in the HC group. (b) The LDA value distribution histogram shows the species with signi cant differences in abundance between two groups. The length of the histogram represents the size of the in uence of the different species (LDA Score).

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.