Following a rigorous diagnostic and exclusion process ( Figure 1A ), individuals with HUA and non-HUA individuals were divided into a case group (n=30) and a control group (n=120), respectively. A total of 62 factors in 5 main categories were used for univariate analysis ( Supplementary Table S1 ). Significant differences (P < 0.05) were observed between the two groups in terms of drinking juice, keeping pets, exercise, abnormal stool, coarse grain and animal offal. Subsequently, 6 factors were used for multivariant logistic analysis via collinearity-test method for excluding confounders. Finally, drinking juice (OR = 2.271, 95%CI: 2.006~5.087, P = 0.012) and abnormal stool (OR = 4.061, 95%CI: 1.200~13.744, P = 0.024) were risk factors associated with HUA ( Figure 2A ). The adjusted model showed good fit through the Hosmer-Lemeshow test (χ² = 2.301, P = 0.593), and the receiver operator characteristic (ROC) curve was 0.935 ( Figure 2B ).

Lachnospiraceae, Bacteroidaceae, and Ruminococcaceae were the predominant bacterial families in both groups ( Figure 3A ). The relative abundance of Bifidobacteriaceae was higher in the NH group (7.72%) than in the H group (4.48%) (P < 0.001, Z-test). Enterobacteriaceae accounted for 5.04% in the NH group, while it only accounted for 1.83% in the H group (P < 0.001, Z-test). At the genus level, Prevotella constituted 14.50% in the NH group but was not listed among the top 15 genera in the H group. Conversely, Pseudomonas was ranked fourth in the H group but not in the NH group ( Figure 3B ). These data indicate that the bacterial composition at the genus level in the two groups was distinct.

Subsequently, the diversity and similarity of microbial communities at the genus level in the H group and NH group were analysed using α-diversity indices (Shannon, InvSimpson, Observed species) and β-diversity analyses (principal coordinate analysis [PCoA]). The α-diversity in the H group (Shannon: median 2.47, interquartile range [IQR]: 2.08–2.86; Observed species: median 65.00, IQR: 54–77.5) was significantly lower than that in the NH group (Shannon: median 2.87, IQR: 2.34~3.18, P = 0.020, U-test; Observed species: median 86.50, IQR: 69–93.5, P = 0.004, t-test, Figure 3C ). Moreover, the microbial community composition between the H and NH groups was significantly different (P < 0.001, permutational multivariate analysis of variance (PERMANOVA), Bray-Curtis; Figure 3D ).

To better explore the potential biomarkers for HUA, we employed Linear Discriminant Analysis Effect Size (LEfSe) to identify specific genera between the two groups. The results indicated 20 significantly different bacterial genera between the two groups ( Figure 3E ). Bacillus, Bacteroides, Escherichia, Klebsiella, Pseudomonas, and Streptococcus were enriched in the H group, whereas Catenibacterium, Clostridium, Collinsella, Coprococcus, Faecalibacterium, Gemmiger, Lactobacillus, Mitsuokella, Prevotella, Roseburia, Ruminococcus, Blautia, Lachnospira, and Proteus were enriched in the NH group.

The entire ARGs in gut was called resistome. A total of 18 ARG types (417 subtypes) and 19 ARG types (406 subtypes) were identified in the H and NH groups, respectively ( Figure 4A ). In general, 102 ARG subtypes were unique to the H group and 91 to the NH group, and 315 subtypes were shared between the two groups ( Supplementary Figure S1 ). Although there were certain commonalities, the distribution patterns of resistome across samples were not consistent ( Supplementary Figure S2 ). Total abundance of resistome in the H group (median: 1.10, IQR: 0.93–1.76) was significantly higher than that in the NH group (median: 0.76, IQR: 0.52–1.41, P = 0.039, U-test, Figure 4B ). Specifically, the abundance of tetracycline resistance genes in the H group (median: 0.46, IQR: 0.33–0.65) was significantly higher than in the NH group (median: 0.20, IQR: 0.15–0.37, P < 0.001, U-test). Similarly, vancomycin resistance genes in the H group (median: 0.03, IQR: 0.02–0.05) was significantly higher than in the NH group (median: 0.01, IQR: 0.004–0.02, P < 0.001, U-test, Figure 4B ). Furthermore, tetracycline resistance genes were the predominant ARGs, comprising 31.66% and 18.18% in the H and NH groups, respectively ( Supplementary Figure S3 ).

We further analysed diversity of resistome between two groups using α-diversity and β-diversity. The α-diversity was significantly greater in the H group (Shannon, median: 3.13, IQR: 2.83–3.62; InvSimpson, median: 13.47, IQR: 11.47–20.67) than that in NH group (Shannon, median: 2.80, IQR: 2.37–3.20, P = 0.030, t-test; InvSimpson, median: 8.37, IQR: 6.51–22.98, P = 0.040, U-test, Figure 4C ). Furthermore, different compositions of resistome were observed between the two groups (P = 0.004, PERMANOVA, permutations = 999, Figure 4D ).

LEfSe was used to explore the specific differences in ARG subtypes between the two groups. The results showed that 25 and 8 ARG subtypes were significantly enriched in the H and NH groups, respectively ( Figure 4E ). Specifically, tetracycline resistance genes (tetQ, tet40, tetO, tetM, tet32, tet37, and tetP), macrolide-lincosamide-streptogramin (MLS) resistance genes (ermF, ermC, ermG, ermX, and lsa), vancomycin resistance genes (vanR, vanU and vanS), and several multidrug resistance genes were significantly present in the H group ( Figure 4E ).

More ARGs associated with potential pathogens in hyperuriacaemic individuals.

A total of 9,265 assembled contigs carrying ARGs was obtained from all metagenomics data, in which the 3,957 and 3,835 ARG subtypes were identified in the H and NH group, respectively, after excluding 1,474 contigs without exactly genus information ( Supplementary Table S2 ). We further determined the ARGs subtypes in selected six potential pathogens and six non-pathogens that associated with HUA from LEfSe analysis. Total of non-duplicated ARG subtypes in potential pathogens in genus level including Klebsiella, Escherichia, Bacteroides, Streptococcus, Pseudomonas, and Bacillus possessed higher number in H group than those in NH group (175 vs. 147) ( Figure 5A ), and similar and opposite results were found in non-assembly data (128 in H group vs. 47 in NH group) ( Supplementary Figure S4 ) and non-pathogens including Prevotella, Ruminococcus, Blautia, Roseburia, Lachnospira,and Alistipe (35 in H group vs. 47 in NH group) ( Figure 5B ).

Furthermore, we also determined the number of ARG subtypes in pathogens and non-pathogens in H and NH group. We found that all pathogens in H group harbour significant higher number of ARG subtypes than those in NH group, especially in Escherichia genus (899 Vs. 550, P < 0.01, Z-test) ( Figure 5C ). However, the number of ARG subtypes in non-pathogens in H group was similar or lower than those in NH group ( Figure 5D ). In addition, multi-drug resistance genes were likely appeared in Escherichia, and the florfenicol resistance genes were found only in Bacillus in the H group. Above data indicates that potential pathogens harbour a greater diversity and number of ARGs in HUA individuals.

We further revealed that the ARG number (ρ = 0.358, P = 0.005), Shannon index (ρ = 0.24, P < 0.001) and InvSimpson index (ρ = 0.333, P = 0.010) were significantly correlated with the UA levels; however, they were negatively correlated with gut microbiota ( Figure 6A ). Procrustes analysis demonstrated a strong consistency of ARGs and gut microbiota in the H group (M² = 0.534, P = 0.001, permutations = 999), but not in the NH group (M² = 0.977, P = 0.001, permutations = 999, Figure 6B ).

Logistic regression analysis was performed to investigate the association between the ARGs and HUA. After excluding collinearity factors, age, sex, BMI, education, household income, exercise, basal metabolic rate (BMR), gut bacterial indicators (Shannon, number, and InvSimpson), and ARG indicators (Shannon, number, and InvSimpson) were included in our model using a forward selection method. We found that ARG number (OR = 1.090, 95%CI: 1.04–1.13, P < 0.001) and ARG InvSimpson (OR = 0.910, 95%CI: 0.836–0.992, P = 0.031)/gut bacterial number (OR = 0.941, 95%CI: 0.901–0.983, P = 0.006) were positively and negatively associated with HUA, respectively (iterated three times, Figure 6C ). The goodness and reliability of this model were confirmed using the Hosmer–Lemeshow test (χ² = 8.586, P = 0.378 > 0.05) and diagnostic receiver operating characteristic (ROC) curve area (0.900) to be satisfactory after adjustment, as evidenced by and the combined with 0.900 value ( Figure 6D ).

We included 5 demographic factors, 5 physiological factors, and 20 dietary intake levels to investigate the potential factors that may affect the differences in ARGs between the two groups. In total, 10 factors (8 of the 10 factors were associated with dietary consumption) were found to be significantly associated with the ARGs ( Figure 7A ). Pork consumption had the highest explanatory power (R² = 0.04, PERMANOVA, P = 0.014). We conducted a similar separate analysis for the H and NH groups. Our data showed more factors contributing to the effects on the ARGs in the H group (eight in the H group and three in the NH group) ( Supplementary Figure S5 ), especially diet-related factors, associated with the ARGs in HUA individuals, indicating a higher complexity.

Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathway analysis revealed differentially enriched genes between the two groups and indicated different microbial compositions and metabolic functions between the two groups (H vs. NH, Figure 7B ). The differentially genes between the H and NH groups were found to be enriched in pathways related to ARGs, amino acid biosynthesis and metabolism, phenylalanine, tyrosine, tryptophan biosynthesis, and energy metabolism. This emphasizes the differences in microbial composition between the H and NH groups, particularly in functions related to inflammation, immune response, bacterial secretion, and protein export. All these suggest that the microbiota in H group was prone to possess the higher ability of colonisation and pathogenesis in gut.

We further analysed virulence-associated genes in these two groups and identified 2682 and 2705 VGs in the H and NH groups, respectively ( Supplementary Figure S6 ). Genes related to adherence (P < 0.001, U-test), antimicrobial activity/competitive advantage (P < 0.001, U-test), and exoenzymes (P < 0.001, U-test) were significantly upregulated in the H group ( Supplementary Figure S7 ). This finding further confirmed our observations in KEGG analyses.

A flowchart of the study design is shown in Figure 1 . A total of 1,032 volunteers were recruited through advertising posters and face-to-face introductions in a large community in Shenzhen City from March 1^st 2018 to July 31^st 2022. Each volunteer was assigned a unique number to avoid duplicate enrolment. Subsequently, participants provided stool samples and completed a questionnaire. Individual consent forms were translated into Mandarin, and consent was obtained from all volunteers. We did not recruit volunteers for whom consent could not be obtained. All participants retained the right to withdraw from the study at any stage.

Given the study’s focus on elderly people; we applied the following exclusion criteria: age less than 65 years, absence of stool samples, and lack of questionnaire data. After exclusion, 504 participants were included in this study. Subsequently, subjects were included based on HUA diagnostic criteria: 1) a normal purine diet, 2) blood uric acid (UA) levels exceeding 360 μmol/L for females and 420 μmol/L for males, respectively (Luo et al., 2022). We also excluded individuals who were missing key covariate data such as age, gender, BMI, smoking status, alcohol consumption, education level, and income level, and also who had serious illnesses like other gastrointestinal disorders, HIV, or cancer and consumed antibiotics within the past month. Finally, 30 participants with HUA were included in the case group. A matching ratio of 1:4 was used to randomly select participants (n = 120, 1:4) without HUA as the control group. After 1:1 matching by age and sex, the 30 participants in the case group and the 120 participants in the control group were respectively divided into the H group and the NH group for subsequent analysis.

Demographics, lifestyle, medical history, and physical activity data were collected using a questionnaire. Participants’ habitual dietary intake was assessed using a food frequency questionnaire that recorded their food consumption over a week. Trained nurses from the community health centre measured the weight, height, waist circumference, hip circumference, and blood pressure. All data were entered using EpiData Version 3.1 with a double-entry method.

All stool samples were placed in a low-temperature transportation box and promptly sent to the laboratory (within 8 h) and were stored at -80°C until they were ready for genomic DNA extraction. The genomic DNA was extracted using Tiangen Magnetic Beads (Beijing, China). The purity and integrity of the extracted DNA were evaluated by 1% agarose gel electrophoresis. The DNA was quantified using the Qubit^® dsDNA Assay Kit by Qubit^® 2.0 Fluorometer (Thermo Fisher Technologies, CA, USA). Approximate 1 μg of DNA was used for library construction using the NEBNext^® Ultra DNA Library Prep Kit for Illumina (NEB, CA, USA). Preliminary quantification was performed using Qubit 2.0 after library construction, and subsequently, the insert size of the DNA library was checked using Agilent 2100 (Agilent, CA, USA). The effective concentration of the library was accurately quantified using the Q-PCR method (the effective concentration of the library > 3 nM). Paired-end sequencing reads were generated using an Illumina NovaSeq platform (Illumina, CA, USA) by Novogene Co., Ltd. All metagenomic sequence data had been uploaded to the public NCBI sequence read archive (SRA) database under the BioProject number PRJNA1059748.

Raw reads were quality-filtered using fastp (Chen et al., 2018) to remove sequences with a quality value > 38 or ambiguous nucleotides > 10. In cases where there was potential host contamination, default alignment was performed using Bowtie2 software (version 2.2.4, http://bowtie-bio.sourceforge.net/bowtie2/index.shtml) with the following parameters (Karlsson et al., 2013): –end-to-end, –sensitive, –I 200, –X 400. The reads with a match consistency ≥ 90% were filtered out. After quality filtering, approximately 10 GB of clean data were obtained for each sample. MetaPhlAn version 4.0.4 (17 Jan 2023) was used to analyse the relative abundance data at the family and genus levels (Truong et al., 2015). The abundances of ARGs, virulence genes (VGs), and MGEs were calculated using ARGs-SOP (Version 2.0) with default parameters using the SARG.2.2, VFDB, and Mobile Genetic Element Databases, respectively (Liu et al., 2022). To identify the hosts of antibiotic resistance genes, we used MEGAHIT version1.0.4-beta to assemble the data (Nielsen et al., 2014), and the contigs longer than 500 bp were kept for subsequent analyses. The filtered contigs were annotated for species using Kraken2 and we used DIAMOND to search the protein sequences against the SARG.2.2 database to identify ARGs containing contigs. The relative abundance of the KEGG pathway was determined using HUMANN3 (Franzosa et al., 2018).

For risk factors analyses, univariate analysis was conducted using the chi-square test or Fisher’s exact test (expected frequency < 5, use Fisher’s exact test; expected frequency > 5, use chi-square test). A collinearity-test was performed on the variables to select those with a variance inflation factor (VIF) < 3 for subsequent multivariate logistic regression analysis (VIF < 3 indicates no collinearity). Variables with P < 0.05 and VIF < 3 in previous analyses were further included in multivariable logistic regression model (enter, P < 0.05; excluded, P > 0.1, forward selection). We used the Wilcoxon rank-sum test (Mann–Whitney U-test) or t-test (depending on whether the data were normally distributed) to compare the number of observed features (i.e. richness), InvSimpson (i.e. richness), and Shannon index (i.e. diversity) (Anderson, 2001). The Z-test was used to compare the proportions of microbiota. We conducted PCoA (based on Bray–Curtis distances) to explore the similarities and differences in microbiota/ARGs/MGEs between the H and NH groups (Nocelli et al., 1999). LEfSe was used to distinguish differences in microbiota and ARGs between and within the H and NH groups (Segata et al., 2011). We also used Spearman’s correlation analysis to evaluate the linear associations between ARGs richness and diversity, and UA levels. All analyses were performed in R v3.6.1 (R Foundation for Statistical Computing, Vienna, Austria). P values < 0.05 were considered statistically significant. Heatmaps were generated using the heatmap package, while most other graphs such as box plot and bar chart were produced using the ggplot2 package (Chen et al., 2020). The α-diversity indices (Shannon, observed species, InvSimpson) of bacterial communities, ARGs, and MGEs for each sample were calculated using vegan package in R (Barberán et al., 2012). The β-diversity analysis was conducted using PCoA based on Bray–Curtis distances via vegan package (Cao et al., 2020), and PERMANOVA was used for statistical analysis (Anderson, 2001). We conducted Procrustes analysis based on Bray–Curtis dissimilarity to investigate the relationship between ARGs and gut microbiota. Network analyses between potential pathogens and ARGs were visualised using Gephi v0.9.2 (Chen et al., 2020).