Metagenomic Study Revealed the Potential Role of the Gut Microbiome in Gout

2.1 Intestinal microbial dysbiosis in gout patients

To investigate the gut microbiota profiles in gout patients, we performed metagenomic shotgun sequencing of 307 fecal samples from 188 individuals. There were 140 samples in the discovery cohort, including 77 gout patients and 63 healthy controls, and 48 samples in the validation cohort, including 25 gout patients and 23 healthy controls. In addition, 70, 40 and 9 fecal samples were collected from gout patients in the discovery cohort longitudinally at three time points (weeks 2, 4 and 24), to assess the effect of therapeutic intervention on the gut microbiome in gout. All patients’ information is summarized in Table S1 (Supporting Information). Importantly, statistical analysis of clinical parameters (Table 1) showed that differences mainly focused on the inflammation, kidney and gout-associated indices, such as C-reactive protein (CRP), erythrocyte sedimentation rate (ESR), serum creatinine (SCr), estimated glomerular filtration rate (eGFR), serum uric acid (SUA), visual analog scale (VAS), arthralgia, joint swelling scores (JSS) and joint tenderness scores (JTS). Moreover, we observed an elevated BMI in gout patients, consistent with the view that gout was associated with OB^[15]. Furthermore, gout-associated indices declined after receiving medication, suggesting the relief of gout symptoms. Gout patients and healthy controls showed no significant difference on specific dietary habits, including drinking, probiotics/prebiotics and completely vegetable-based diet.

We carried out rarefaction analysis by random sampling 100 times with replacement to assess the metagenomic gene counts in discovery samples (healthy control, n = 63; gout patient, n = 77). As the sample size increased, both curves tended to be flat with few new genes detected, suggesting that the sequencing data were sufficient (Figure 1A). We found that gout patients had lower microbial gene richness and diversity than healthy controls (gene number, P = 0.0024; Shannon index, P = 0.0016; Wilcoxon rank-sum test; Figure 1B, C). Beta diversity analysis indicated a significantly greater dissimilarity in microbial community between samples in disease and healthy control groups than those within each group (P = 0.0069, Wilcoxon rank-sum test; Figure 1D).

• Download figure
• Open in new tab

Figure 1.

Gut microbial alterations in gout patients. (A) The gene rarefaction curves based on the Chao2 values in healthy controls (n = 63) and gout patients (n = 77). (B) Box and whisker plots of gene count in the healthy controls and gout patients. (C, D) Box and whisker plots of alpha diversity (Shannon index) and beta diversity (Bray Curtis distance) at the gene level. Wilcoxon rank-sum test was used to determine significance. ‘**’ P < 0.01. To exclude the influence of the various data sizes among the samples, Figures A, B, C and D were based on 11 M matched reads per individual. (E) Principle component analysis (PCA) based on the genus abundance profile. Considering PC1 significantly related with enterotype, PC2 and PC3 were display. (F) The Bacteroidetes/Firmicutes ratio (Wilcoxon rank-sum test; ‘**’, P < 0.01). (G, H, I) The relative abundance of differential phyla (G, top 3), genera (H, top 30) and species (I, top 30) between gout patient and healthy control groups (FDR < 0.05, Wilcoxon rank-sum test). The color bar above genera or species names were colored according to the phylum.

To better understand the gut microbial differences between gout patients and healthy controls and to identify microbial taxa associated with gout, all genes were annotated, and microbial relative abundances at different taxonomic levels were obtained. In discovery cohort, a total of 18 phyla, 400 genera and 861 species were identified, respectively. The 5 predominant phyla were Bacteroidetes (52.65%), Firmicutes (15.24%), Proteobacteria (3.40%) Fusobacteria (0.68%) and Actinobacteria (0.63%) (Figure S1A, Supporting Information). At the genus level, genera with an average relative abundance more 1% were Bacteroides (34.89%), Prevotella (9.68%), Faecalibacterium (3.07%), Roseburia (1.89%), Alistipes (1.86%), Parabacteroides (1.67%), Blautia (1.20%), Escherichia (1.18%) and Megamonas (1.07%) (Figure S1B, Supporting Information). These 9 genera accounted for more than 50% of gut microbial composition. Considering that enterotype has been reported to be related to diseases^[10], we investigated enterotype of each sample in discovery cohort. We identified two distinct enterotypes in 140 samples (enterotype1, 99 samples; enterotype2, 41 samples) (Figure S2A, Supporting Information) at the genus level, each containing healthy controls and gout patients. Predominant genera in enterotype1 and enterotype2 were Bacteroides (P = 1.96e-18, Wilcoxon rank-sum test) and Prevotella (P = 2.67e-20, Wilcoxon rank-sum test) (Figure S2B, Supporting Information), respectively. Moreover, no significant difference of the two enterotypes distribution between gout patients and healthy controls was observed (P = 0.8541, Fisher’s exact test). And all clinical indices except SCr (P = 0.0118, Wilcoxon rank-sum test) also showed no significant difference between individuals of two enterotypes (Figure S2C, Supporting Information). These findings suggested that there was no significant association between enterotype and gout. By examining the top 5 principal components (PC), the first PC was significantly correlated with enterotype (false discovery rate (FDR) = 1.53e-20, Wilcoxon rank-sum test), whereas the second was significantly correlated with gout (FDR = 2.25e-05, Wilcoxon rank-sum test; Figure 1E), which indicated that gout status was one of determining factors of gut microbial composition in our study. Similar results were obtained at functional gene level (PC3, FDR = 1.01e-06; Figure S2D, E, Supporting Information). Though BMI was significantly different between two groups, none of the top 5 PCs were significantly correlated with it (Table S3, Supporting Information). Furthermore, no genus was detected to be significantly correlated with BMI using MaAsLin analysis (q value < 0.25). These results indicated that BMI had limited influence on gut microbiota in our study.

Wilcoxon rank-sum test was used to detect the significantly different bacteria between two groups. At the phylum level, Bacteroidetes (FDR = 0.022, Wilcoxon rank-sum test) and Fusobacteria (FDR = 0.0183, Wilcoxon rank-sum test) were enriched in gout patients, whereas Proteobacteria (FDR = 0.0183, Wilcoxon rank-sum test) were enriched in healthy controls (Figure 1G). Furthermore, we observed that the ratio of Bacteroidetes to Firmicutes (two main phyla of the human gut microbes) was higher in gout patients than in healthy controls (P = 0.0098, Wilcoxon rank-sum test; Figure 1F), which suggested that the gut microbial composition of gout patients had changed. A total of 110 genera and 223 species were significantly different in abundance between the two groups (FDR < 0.05, Wilcoxon rank-sum test; Figure 1H, I and Tables S4, 5, Supporting Information). All differential Bacteroides spp. (3 species), Prevotella spp. (13 species) and Fusobacterium spp. (4 species) were enriched in gout patients. Bacteroides has been reported to be significantly higher in the gut microbiota of gout patients and may be associated with inflammation^[13]. Similarly, increased Prevotella-mediated gut mucosal inflammation leads to systemic diseases^[16]. Prevotella copri (FDR = 0.0169, Wilcoxon rank-sum test) was observed to be increased in untreated rheumatoid arthritis^[17] and displayed the ability to increase inflammation in a mouse model of colitis^[17]. Fusobacterium spp., especially Fusobacterium nucleatum (FDR = 0.0035, Wilcoxon rank-sum test), have been demonstrated to correlate with colorectal cancer (CRC) and were often considered a pathogen. They can generate a proinflammatory microenvironment that is conducive to the development of CRC^[18]. Lactobacillus spp. (15 species) and Streptococcus spp. (11 species) were also higher in gout patients. By contrast, butyrate-producing species (P < 0.05, FDR < 0.1, Wilcoxon rank-sum test), such as Roseburia spp. (3 species), Coprococcus spp. (3 species), Eubacterium spp. (3 species), Faecalibacterium prausnitzii and butyrate-producing bacterium SS3/4 were enriched in healthy controls. These species were reported to have anti-inflammatory effects^[19,20]. Some species belonging to Enterobacteriaceae, including Escherichia spp. (2 species), Klebsiella spp. (7 species), Enterobacter spp. (9 species) and Citrobacter spp. (7 species), were reported to degrade uric acid^[21]. These species were higher in healthy controls, possibly helping reduce uric acid accumulation in the healthy controls.

2.2 Functional characterization of the gut microbiome in gout patients

To investigate the functional differences in the microbiome between gout patients and healthy controls in discovery cohort, a KEGG Orthology (KO) profile was constructed according to the reference gene catalogue annotation. In total, 7,289 KOs were identified in all samples (n = 259), of which 2,666 KOs (2,260 and 406 KOs enriched in healthy controls and gout patients, respectively) were differentially enriched between the two groups (FDR < 0.05, Wilcoxon rank-sum test; Table S6, Supporting Information). KOs profile was used to perform PCA, and the top two principal components were significantly different between gout patients and healthy controls (PC1, P = 0.0008; PC2, P = 0.0043, Wilcoxon rank-sum test; Figure S3A, B, Supporting Information), which indicated that the microbial function of gout patients was altered.

To further explore the functional changes, the different KEGG pathways and modules between the two groups were identified (Figure S3C, D and Table S7, Supporting Information). The gut microbial function of gout patients displayed higher potential for carbohydrate metabolism, including fructose and mannose metabolism, starch and sucrose metabolism, amino sugar and nucleotide sugar metabolism, the citrate cycle and 3 modules of glycolysis/gluconeogenesis. Low carbohydrates was showed to reduce SUA and gout attacks^[22]. The increase in fructose was reported to be associated with the risk of gout^[23], and a previous study showed that a rapid infusion of fructose could raise the level of SUA^[24]. Concurrently, we observed increased potential for lipopolysaccharide (LPS) biosynthesis and lipid A biosynthesis (Figure S4A, B, Supporting Information). LPS was reported to be a common component of gram-negative bacteria and could be an endotoxin associated with inflammation and disease^[25]. In the two groups, Bacteroides and Prevotella were the major contributors to lipid A biosynthesis. We note that lipid A biosynthesis genes from Proteobacteria (FDR = 0.017, Wilcoxon rank-sum test) species were enriched in healthy controls (Figure S4C, Supporting Information), including Escherichia coli (FDR = 0.005, Wilcoxon rank-sum test) and Klebsiella pneumonia (FDR = 0.014, Wilcoxon rank-sum test). In contrast, the samples of healthy controls showed higher potential for cell motility, containing bacterial chemotaxis and flagellar assembly, similar to the previous results observed in AS^[11] and type 2 diabetes (T2D)^[26], and such mechanisms may be beneficial for health. Previous studies indicated that the gut microbiota may take part in purine metabolism^[13,14]. Consistently, we observed that healthy controls had a higher potential for urate degradation (FDR < 0.05, Wilcoxon rank-sum test; Figure 2A, B) and the major contributor was Enterobacteriaceae spp. (FDR < 0.05, Spearman’s rank correlation; Figure 2C and Table S8, Supporting Information) in which Klebsiella showed the highest Spearman’s rank correlation (mean Spearman’s rank correlation of 5 urate degradation KOs was 0.53). Meanwhile, SUA was negatively correlated with Enterobacteriaceae spp., especially Klebsiella (r = −0.23, P = 0.0057, Spearman’ rank correlation; Figure 2D and Figure S5, Supporting Information), which indicated that the gut microbiome had the ability to degrade uric acid. To verify this result, we investigated the distribution of KOs involved in urate degradation among 1,135 species (Figure S6 and Table S9, Supporting Information). Proteobacteria (706 species) was the predominant contributor, and the largest family was Enterobacteriaceae (119 species), which is consistent with our results. As SCFAs have been reported to be linked to healthy gut microbiota, we analyzed the abundance of genes encoding key enzymes involving SCFAs production. We found that the genes responsible for propionate and butyrate were significantly different between gout patients and healthy controls, whereas genes for acetate did not show any difference (Figure S7, Supporting Information).

• Download figure
• Open in new tab

Figure 2.

Gout-associated microbial gene functions related to urate degradation. (A) Module for urate degradation. (B) Relative abundance of KOs involved in urate degradation. Significantly enriched KOs were identified by a Wilcoxon rank-sum test, and the boxes or KO names were colored according to the direction of enrichment. Green, enriched in healthy controls (FDR < 0.05). Boxes with no color or KO names with black, no difference; boxes with gray, not detected in samples. (C) Correlations between gout-associated genera and urate degradation-associated KOs (red and green for positive and negative correlation, respectively). ‘+’ denotes FDR < 0.05; ‘*’ denotes FDR < 0.01; ‘#’ denotes FDR < 0.001. (D) The associations between SUA and healthy control-enriched species. Spearman’s rank correlation was calculated by taking the species relative abundance and SUA content. An inverse correlation was observed between SUA and Proteobacteria, Enterobacteriaceae and Klebsiella. Abbreviation: hpxO, FAD-dependent urate hydroxylase; uraH, 5-hydroxyisourate hydrolase; hpxQ, 2-oxo-4-hydroxy-4-carboxy-5-ureidoimidazoline decarboxylase; allB and hpxB, allantoinase.

2.3 Metagenomic species (MGS) associated with gout

To further explore the features of the gout-associated microbiome and the correlations between clinical indices and the microbial alteration in gout patients, 178,233 differentially enriched genes (60,511 and 117,722 enriched in healthy controls (n = 63) and gout patients (n = 77), respectively) were clustered into 3,831 co-abundance gene groups (CAGs) based on the co-variation of relative abundance across all samples using the Canopy-based algorithm. CAGs with fewer than 50 genes were filtered, and the remaining 134 were termed metagenomic species (MGSs) for follow-up analysis (Table S10, Supporting Information). All MGSs were significantly different between the two groups, where 73 MGSs were more abundant in gout patients, and the remaining 61 MGSs were more abundant in healthy controls (FDR < 0.05, Wilcoxon rank-sum test). Notably, only 24 of the healthy control-enriched MGSs could be annotated as specific species compared with 53 of gout patients. Subsequently, we carried out a co-occurrence network analysis for all MGSs based on SparCC correlation coefficients (pseudo P < 0.05; Figure S8, Supporting Information). We observed more intensive and sophisticated interactions between gout-enriched and gout-depleted MGSs than within each group.

Permutational analysis of variance (PERMANOVA) was performed to assess the influence of clinical indices on MGSs and we observed that the disease status was the strongest factor (lowest P-values in PERMANOVA; Table S11, Supporting Information). Other clinical indices such as SUA, CRP, eGFR, ESR and SCr also showed correlation with microbial alteration. Canonical coordinate analysis (CCA) showed that gout-enriched and control-enriched MGSs were clearly separated (Figure S9, Supporting Information). All clinical indices were associated with gout-enriched MGSs except for eGFR, and gout status was the main factor affecting the composition of the gut microbiome, followed by SUA. Spearman’s rank correlation coefficient analysis was conducted for physiological indices and MGSs (FDR < 0.05; Figure S10 and Table S12, Supporting Information). A high level of SUA is a major feature of gout patients. We found SUA was positively correlated with many gout-enriched MGSs but negatively correlated with many control-enriched MGSs. Interestingly, we observed that almost all Bacteroides spp. enriched in gout patients were correlated with SUA, whereas no association was found between Bacteroides spp. enriched in healthy controls and SUA. The indicator of inflammation CRP showed a similar result. SCr displayed positive correlation with several gout-enriched MGSs. An increase in the level of SCr suggests impaired renal function. BMI, age, eGFR, ESR and urea displayed a correlation with a small number of MGSs. Furthermore, we carried out PCoA and PERMANOVA on the MGS level for gout related phenotypes, including gout duration, VAS, JSS, JST and tophi of gout patient samples (Figure S11A-E and Table S11, Supporting Information), only tophi significantly separated the gout patient samples. In a subsequent step, significantly different MGSs (P < 0.05, Wilcoxon rank-sum test) were detected between gout patients with and without tophi. Control-enriched MGSs, such as F. prausnitzii, were higher in patients without tophi, whereas gout-enriched MGSs, especially Bacteroides spp., were enriched in patients with tophi (Figure S11F, Supporting Information).

2.4 Gut microbiota-based gout classification

The distinct microbial differences between gout patients and healthy controls prompted us to investigate if the gut microbiome had the potential to distinguish gout patients from healthy controls. A disease classifier was constructed based on the MGS profile using a random forest model in the discovery cohort (healthy control, n = 63; gout patient, n = 77), and the corresponding receiver operating characteristic (ROC) curve was constructed. All MGSs were entered into the disease classifier, and 16 MGSs were selected as biomarkers to establish a classification model with a fivefold cross-validation approach (Figure 3A-C and Table S13, Supporting Information). Notably, the area under the receiver operating curve (AUC) reached 0.93. We found that only three MGS biomarkers were enriched in healthy controls, whereas others were enriched in gout patients. Five of the top ten most important MGSs for classification, including Gout-107, Bacteroides sp. (100), Megamonas sp. (102), Blautia sp. (92), Parabacteroides distasonis (74), had the lowest P value. To evaluate the accuracy of this classification model, an independent cohort consisting of 25 gout patients and 23 healthy controls was used for validation, and the AUC value was 0.89 (Figure 3D, E and Table S14, Supporting Information). In addition, 16 MGS biomarkers were used to distinguish patients from healthy controls in the AS (n = 211), RA (n = 169), T2D (n = 268) and OB (n = 200) cohorts (Figure 3F). The AUC was lower than in the gout cohort, suggesting that these 16 MGS biomarkers were gout-specific. These results indicated that gout-associated microbial biomarkers have the potential to be used as a novel and noninvasive diagnostic strategy for gout.

• Download figure
• Open in new tab

Figure 3. Classification of gout status by relative abundance of MGS.

(A) Distribution of 5 trials of 10-fold cross-validation error in random forest classification of gout as the number of MGS increases. The model was trained using the relative abundance of the MGS (> 50 genes) in the controls and gout samples. The black curve indicates the average of the five trials (gray lines). The red line marks the number of MGSs in the optimal set. (B) The probability of gout in the cross-validation discovery samples according to the 16 MGS biomarkers. (C) Receiver operating curve (ROC) for the discovery samples. (D) Classification of the validation samples (healthy control, n = 63; gout patient, n = 77). (E) ROC for the validation samples (healthy control, n = 23; gout patient, n = 25). (F) ROCs were obtained for five cohorts using 16 gout-associated MGS biomarkers. The AUC for each disease is displayed in brackets.

2.5 Effect of uric-acid-lowering and anti-inflammatory drugs intervention on the gut microbiome

Previous studies of metagenomics have demonstrated the potential influence of therapeutic drugs on the microbiota^[27]. Here, we assessed the effect of drug intervention on intestinal microbiota in this gout cohort. We conducted a 2-week (2W, n = 61), 4-week (4W, n = 38), and 24-week (24W, n = 7) treatment intervention for gout patients. Most of the treated gout patients received uric-acid-lowering and anti-inflammatory drugs intervention (Figure S12A and Table S1, Supporting Information). We excluded 9, 2, and 2 fecal samples from the 2W, 4W and 24W groups, respectively, from this analysis, because these patients did not use any drugs or used other non-conventional drugs (the sample size was too small), such as NaHCO3 and traditional Chinese medicine. Subsequently, we compared their microbial composition with untreated gout patients (n = 77) and healthy controls (n = 63). The gene numbers, Shannon index and PCoA showed that the gut microbial composition was significantly different between healthy controls and gout patients. (P < 0.05, Wilcoxon rank-sum test; Figure S12B, C, Supporting Information). However, no significant difference was found between before and after treatment in gout patients, and PERMANOVA analysis also supported this result (Table S15, Supporting Information). Considering that the unmatched sample size at different time points may influence the results, we selected five gout patients whose fecal samples were collected at all the time points for further analysis. Interestingly, we observed samples of 24W treatment and found that the gut microbial composition was further away from other time points and more similar to samples of healthy controls (Figure S12D, Supporting Information). This may indicate that the effect of uric-acid-lowering and anti-inflammatory drugs on the gut microbiome was modest and partially restored the gut microbiome over time, similar to disease-modifying anti-rheumatic drugs (DMARD)^[12]. Allopurinol, colchicine, celecoxib and etoricoxib had no impact on human gut bacteria in vitro in a previous report^[28]. Benzbromarone has been shown to inhibit gut bacterial growth in vitro^[28] Our results suggested a more complicated relationship between these drugs and gut microbiota in vivo. However, due to the insufficient data in our study, the mechanisms of specific drugs on the gut microbiome needs to be further studied in intervention trials of larger cohorts.

2.6 Comparison of the gut microbiome between gout and other diseases

To better understand the potential role of intestinal microbes in gout, we compared the gut microbial composition of gout with AS, RA, T2D and OB, of which the first two are autoimmune diseases, and the last two are metabolic diseases. However, as the comorbidity data is not available, we cannot exclude that there might be individuals with a combination of multiple diseases. From the P value distribution, we inferred that the degree of gut microbial dysbiosis was intermediate in gout patients when compared with the other four diseases, with OB being more severe and RA being least severe (Figure 4A). The number of significant microbiota markers also supported this conclusion (Figure 4B). Very few significant microbiota markers (Wilcoxon rank-sum test, FDR < 0.05, only 6 genes FDR less than 0.05, Figure 4B) were detected in RA patients. Surprisingly, the gut microbial composition of gout patients was closer to AS patients (Figure 4C and Figure S13A, Supporting Information), and they shared 40 different species. A significant enrichment of Acidaminococcus fermentans, Streptococcus anginosus, Megasphaera elsdenii and Megasphaera micronuciformis was observed in all diseases except RA. In addition, Streptococcus spp. also had a higher abundance in patients. For gut microbial functions, gout patients were also more similar to autoimmune diseases RA and AS (Figure 4D and Figure S13B, Supporting Information). The functions of oxidative phosphorylation; alanine, aspartate and glutamate metabolism; carbon fixation pathways in prokaryotes; and carbapenem biosynthesis were enriched in gout, AS and RA patients, whereas the functions of bacterial chemotaxis and flagellar assembly were depleted in these three diseases and T2D. Except for fructose and mannose metabolism and S. aureus infection, no differential function was shared among gout, OB or T2D. Collectively, these results indicated that the gut microbiome of gout patients was more similar to autoimmune diseases than metabolic diseases.

• Download figure
• Open in new tab

Figure 4.

Comparison of gut microbiome and functions in different diseases. (A) The distribution of P-values in AS (n = 211), gout (n = 140), OB (n = 200), RA (n = 169) and T2D (n = 268). (B) The number of gene markers for each disease. Statistical differences between control and case of each disease was evaluated using Wilcoxon rank-sum test and adjusted by the Benjamin-Hochberg method. (C) Comparison of differential species in AS, gout, OB, RA and T2D. Green, enriched in healthy controls; red, enriched in patients; ‘*’ denotes FDR < 0.05. (D) Comparison of microbial gene functions in AS, gout, OB, RA and T2D. Green, enriched in healthy controls; red, enriched in patients. ‘*’ denotes reporter score of maps > 1.65 or < −1.65.

DNA extraction and library construction

The collected fecal specimens were centrifuged at 12000 x g at room temperature for 5 min and the supernatant was discarded. 200 mg pellet was weighted from each sample and used for total bacterial DNA extraction with the E.Z.N.A Stool DNA Kit (OMEGA Bio-tek, USA) according the manufacturer’s instructions. The quality of DNA was analyzed using Qubit (Invitrogen, USA) and 1% agarose gel electrophoresis. The detail of DNA library construction was described in Supporting Information file. The final DNA library was determined the average insert size using the Agilent 2100 Bioanalyzer (Agilent Technologies, USA) and quantified by ABI StepOnePlus Real-Time PCR system (Applied Biosystems, USA).

Metagenomic sequencing

Paired-end metagenomic sequencing was performed on the Illumina HiSeq 4000 platform with an insert size of 350 bp and paired-end (PE) reads of 150 bp for each sample. After removing adaptors and low quality and ambiguous bases from the raw reads, the remaining reads were aligned to human genome reference (hg19) by SOAPaligner (v2.22, parameters: -m 280 -x 420 -r 1 -l 32 -s 75 -c 0.9) to remove human host DNA contamination. The average rate of host contamination was 0.52 ± 2.06%. Finally, 2,768.76 Gb of high-quality PE reads for the 307 samples were acquired with an average of 9.02 Gb per sample (Table S2, Supporting Information). The metagenomic shotgun sequencing data for all samples have been deposited in the CNGB Nucleotide Sequence Archive (CNSA) under accession code CNP0000284.

Construction of gene, phylum, genus, species and KO profiles

The clean reads were aligned to the 11,446,577 genes in the reference gene catalogue^[52] using SOAPaligner (v2.22, parameters: -m 200 -x 1000 -r 2 -v 13 -l 32 -s 75 -c 0.95), and 72.79 ± 2.89% reads (n = 307, data are mean ± s.d.) were mapped. The gene abundance profile was calculated as previously described^[11]. After removing genes detected in less than 10% of the discovery samples (n = 259), 1,564,977 genes remained. To improve the gene taxonomy annotation result, genes were aligned to the National Center for Biotechnology Information (NCBI) microbial reference genomes (including 5,847 microbial genomes, v20171114) and NT database (v20170924) using BLAT (v.36) and Megablast (v2.2.26) with default parameters, respectively. The cutoff was at least 70% overlap of query and 65% identity for phylum, 85% identity for genus, and 95% identity for species^[53]. The relative abundance of phyla, genera, species and KOs were calculated from the relative abundance of their respective genes using previously published methods^[53].

Rarefaction curve, gene counts and biodiversity analysis

Rarefaction analysis was performed to assess the gene richness in the healthy controls and gout patients. For a given number of samples, we performed random sampling 100 times in the cohort with replacement and estimated the total number of genes that could be identified from these samples by the Chao2 richness estimator^[54]. The total gene counts in each sample were calculated as previously described^[55]. The alpha diversity and beta diversity were estimated by the Shannon index and Bray-Curtis distance, respectively. To adjust for the effect of the various data sizes among the samples, these analyses were based on the gene profile that was randomly sampled to 11 million (M) matched reads.

Identification of enterotypes

The enterotypes of each sample in discovery cohort were identified by the PAM clustering algorithm using the relative abundance of genus. Calinski-Harabaze (CH) index was used to assess the optimal number of clusters^[56]. Finally, PCA was performed using Jensen-Shannon distance to cluster the samples.

Gene functional analysis

Differentially enriched KEGG pathways and modules were identified according to reporter score^[57] from the Z-scores of individual KEGG Orthologies (KOs) (KEGG database release 79). Pathways or modules were considered significantly different if the |reporter score| > 1.65, corresponding to 95% confidence to a normal distribution.

For analysis of urate degradation genes, species that contain at least one urate degradation gene in the KEGG database (v79) were collected, and the corresponding 16S sequences were downloaded from the Ribosomal Database Project (RDP, set14). Then, 16S sequences were aligned by PyNAST (v1.2.2), and trees were constructed by fasttree (v2.1.7). The tree and heat map were visualized in iTOL^[58].

The protein sequences of SCFA-producing enzymes were obtained from the NCBI database^[59]. Genes in the reference gut microbiome gene catalogue were aligned to these sequences using BLASTP (v2.2.26, best hit with e-value < 1e-5, identity > 70% and coverage > 70%), and their relative abundance were calculated as KOs.

Co-abundance gene groups (CAGs)

A total of 178,233 genes differed in relative abundance between healthy controls and gout patients in the discovery cohort (Wilcoxon rank-sum test, P < 0.01, FDR < 0.05) from 1,564,977 genes (present at least 10% of the samples in discovery cohort). These marker genes were clustered into CAGs based on their relative abundance variation across all samples using the Canopy-based algorithm with default parameters^[60]. In brief, a gene was randomly selected as a ‘seed’ and this and other genes with similar abundance profiles were clustered as a CAG if Pearson correlation coefficient of those genes was more than 0.9. In order to avoid the profile being too sparse, relative abundance of gene in the third quartiles of each group was selected as a cluster relative abundance. Finally, 3,831 CAGs (contain at least two genes) were obtained. Clusters including at least 50 genes remained for further analysis and were called metagenomic species (MGSs). MGSs were assigned to a species, genus or phylum if more than 50% of their genes had the same level of assignment^[60]. The relationship between the abundance of each MGS was according to correlation calculated by SparCC^[61], and the co-occurrence network was visualized with Cytoscape (v3.4.0).

Canonical correspondence analysis (CCA)

CCA was performed on the MGS abundance profile to assess the effect of each of the indices listed. The plot was generated by R (v3.4.0, vegan package).

PERMANOVA on the influence of phenotype

Permutation multivariate analysis of variance (PERMANOVA) was performed on the MGS abundance profile to assess the impact from each of the clinical indices listed^[62]. We used Bray-Curtis distance and 9,999 permutations in R (3.4.0, vegan package).

MGS-based classifier

We constructed a classifier to discriminate healthy controls and gout patients based on a random forest model (randomForest 4.6-14 package) using the profile of MGS relative abundance^[62]. A fivefold cross-validation approach was used to evaluate the predictive model. We calculated the average cross-validation error curves and the minimum error in the averaged curve, plus the standard deviation at that point was used as the cutoff to filter the predictive model. The ones containing the smallest number of MGSs among all sets with an error below the cutoff were chosen as a diagnosis classifier. Then, the probability of gout was calculated based on this optimal set, and the receiver operating characteristic (ROC) for both the discovery cohort and validation cohort were drawn using the pROC package^[63].

Comparison of the gut microbiota between gout and other diseases

Published metagenomic sequencing data of ankylosing spondylitis (AS) (SRP100575 and ERP005860; 114 healthy controls and 97 patients with AS)^[11], rheumatoid arthritis (RA) (ERP006678; 74 healthy controls and 95 patients with RA)^[12], type 2 diabetes (T2D) (SRP008047 and SRP011011; 185 healthy controls and 183 patients with T2D)^[26] and obesity (OB) (ERP013562; 105 healthy controls and 95 patients with OB)^[9] were downloaded from NCBI or the European Bioinformatics Institute (EBI). All sequencing reads were mapped to the same reference gut gene catalogue and analyzed for gene function by the same pipeline as the gout cohort. To avoid methodological differences, taxonomic analyses of genes for all cohorts was based on the reference gene catalogue taxa result, including the gout cohort.

Statistical analysis

All statistical analyses were performed by R (v3.4.0). Differential relative abundance of genes, taxa and KOs were detected by Wilcoxon rank-sum with an adjusted P value (corrected by the Benjamini-Hochberg) < 0.05. Enrichment in healthy controls or gout patients was determined according to the higher mean rank-sum. The correlations between the relative abundance of differentially genera and KOs involved to urate degradation KOs, the relative abundance of MGSs and clinical indices were calculated by Spearman’s rank correlation coefficient and visualized by heatmap in R using the ‘pheatmap’ package. The association between species and BMI was calculated by MaAsLin with default parameters^[64].