RNA-seq datasets were retrieved from the Gene Expression Omnibus (GEO) database. GSE192391 included the scRNA-seq data for severe COVID-19; GSE211783 included the scRNA-seq data of six samples from three gout patients. GSE157103 consisted of the bulk RNA-seq data of 50 severe COVID-19 patients and 26 healthy individuals.

In this study, we performed MR Analyses by using the R software package TwoSampleMR to identify single nucleotide polymorphisms (SNPs) associated with severe COVID-19 and gout, which were defined as instrumental variables. The GWAS data were obtained from the website (http://gwas.mrcieu.ac.uk/datasets), including the severe COVID-19 dataset (ebi-a-GCST011075, n = 1,388,342) and the gout dataset (finn-b-M13_GOUT, n = 150,797). The expression quantitative trait loci (eQTL) data were also publicly available from the IEU OpenGWAS project accession code eQTL-a-ENSG00000162711 and eQTL-a-ENSG00000137331. Our study design was summarized in Figure 1.

Quality control of datasets was performed using the Seurat package (version 4.4.0) in R software (version 4.3.1). Samples were first converted to Seurat objects using the “CreateSeuratObject”. Cells with mitochondrial gene percentages < 25% and unique gene counts between 200 and 6,000 were used. The data were normalized using the “NormalizeData” function, scaled for all genes using the ScaleData function, and subjected to principal component analysis (PCA). Hypervariable genes were identified by the FindVariableFeatures function and used for downstream analysis. Since the data were obtained from different samples, batch correction was performed using the R package “Harmony” (1.0.3) to avoid any batch effects interfering with downstream analysis. Cell clustering and classification were performed by using the FindClusters function. The SingleR package (2.2.0) was then used to match the single-cell RNA-seq data to a known reference dataset and manually calibrated to improve the accuracy and reliability of cell type annotation. Marker genes are shown in Supplementary Figures S1A, S1B. Monocyte subpopulations were categorized as classical monocytes (CCR2, SELL, S100A8, S100A9, LYZ, SERPINB2, CD14), non-classical monocytes (NAP1L1, FCGR3A, FCGR3B, CSTA, CX3CR1, ITGAL), and intermediate monocytes (HLA-DRA, HLA-DPB1, EVL) based on markers genes listed in literature. Dot plots show the proportion and average expression of cell clusters expressing the marker gene in monocytes from COVID-19, gout flare and healthy samples (Supplementary Figure S1C), and the marker gene in monocytes from COVID-19, gout remission and healthy samples (Supplementary Figure S1D). The correspondence between cell clusters, cell marker and cell types from COVID-19, gout flare and healthy samples are shown in Supplementary Figure S2A, and cell types from COVID-19, gout remission and healthy samples are shown in Supplementary Figure S2B. The FindAllMarkers function was used to find differentially expressed genes between monocyte and other clusters. The communication between monocytes and other cells was determined by analyzing the ligand-receptor pairs using CellChat (version 1.6.1) R package, with CellChatDB.human as the reference database. The R package scMetabolism was used to quantify the metabolic activities of the different types of cells at the single-cell level.

To ensure the robustness and reliability of the findings, we followed the STROBE-MR guidelines and used a two-sample MR approach (Skrivankova et al., 2021). Two sample MR analysis was performed by using the TwoSampleMR package (0.5.7), with genes as the exposure and the disease as the outcome. Variants in multiple regions of the genome were selected as instrumental variables in this study. The significantly associated SNPs (P < 5 × 10⁻⁸) were screened, and the threshold for removing linkage disequilibrium (LD) was set to r² < 0.001. The instrumental strength for each SNP was assessed using the F statistic, with an F statistic >10 indicating a strong tool. The F statistic was calculated using the formula: F = R 2 × N − 2 1 − R 2 , where R 2 = 2 × β 2 × EAF × 1 − E A F 2 × β 2 × EAF × 1 − EAF + 2 × SE 2 × N × EAF × 1 − EAF (Papadimitriou et al., 2020; Levin et al., 2020; Gill et al., 2019; He et al., 2024). For genes with only one instrumental variable (IV), the Wald ratio method was used; otherwise, the inverse variance weighting method (MR-IVW) was used. In this study, P = 0.00052 (0.05/96) after correction of the Bonferroni method. Therefore, P < 0.00052 was considered significant.

We performed a Coloc test to examine the probability that SNPs associated with disease and gene expression (eQTL) are shared genetic causal variants. Colocalization analyses were performed respectively by using eqtl-a-ENSG00000162711 as exposure data, ebi-a-GCST011075 as outcome data, and eqtl-a-ENSG00000137331 as exposure data and finn-b-M13_GOUT as outcome data. The instrumental variables were determined based on 1 MB upstream and downstream of the SNP with the lowest p-value associated with gene. Coloc package is used in R to evaluate Bayesian factors under various colocalization assumptions. Based on Bayesian statistical modeling, five posterior probabilities were generated, corresponding to the five hypotheses described by Wallace et al. Two hypotheses were the focus of this study: 1) PPH3, association with the severe COVID-19 or gout risk and expression of the gene, with distinct causal variants; 2) PPH4, association with the severe COVID-19 or gout risk and expression of the gene, with a shared causal variant. Posterior probabilities (PP) were used to quantify support for all hypotheses, and co-localization analyses were restricted to genes that achieved PPH3+PPH4 ≥ 0.8 (Giambartolomei et al., 2014). Finally, we obtained the genotype data of the target genes including eqtl-a-ENSG00000162711 and eqtl-a-ENSG00000137331, which were used to drawing the regional associations plots.

Summary-data-based Mendelian randomization (SMR), based on Mendelian randomization, was proposed in 2016 by Zhu et al. (2016). The SMR analysis used single nucleotide variants of the top cis-eQTL as instrumental variables and combined GWAS and eQTL data to detect associations between gene expression and traits. In this study, we used the default settings in the SMR software. The P-value threshold for selecting the relevant eQTL for the SMR test is 5.0 × 10⁻⁸, and the window around the probe center for selecting cis-eQTL is 1 Mb. The eQTL data were obtained from two sources. Respectively, we used a study conducted by Westra et al. (Westra et al., 2013), which was an eQTL meta-analysis of peripheral blood samples from 5,311 healthy European individuals, and CAGE eQTL whole blood expression data, including 2,765 subjects of European ancestry (Lloyd-Jones et al., 2017). The criteria of significant results were as follows: P _SMR<0.05 and P _HEIDI > 0.05.

In this study, we focused on the relationship between severe COVID-19 and gout flare. The scRNA-seq data of PBMCs obtained from patients with severe COVID-19, gout flare, and healthy samples revealed significant differences in the cellular composition of different groups. Based on specific gene markers, PBMCs were clustered into 20 distinct subclusters. Six clusters (clusters 13, 14, 15, 17, 18, 19) were considered “mixed” as the marker genes were not associated with any common or specific cell type, and were subsequently excluded due to the small number of cells. The final cell clusters including natural killer (NK) cells, T cells, B cells, plasmacytoid dendritic cells (pDCs), classical DCs (cDCs), monocytes, neutrophils, and platelets, and the UMAP plot is shown in Figure 2A. The monocytes were analyzed further given their critical role in gout and severe COVID-19 progression. Based on the expression of marker genes, we identified seven distinct monocyte clusters, including classical monocytes (CMs; clusters 0, 1, 4, and 6), non-classical monocytes (NCMs; clusters 2 and 5), and intermediate monocytes (IMs; cluster 3). As shown in Figure 2B, the proportion of classical monocytes was higher in gout flare and severe COVID-19 samples compared to the healthy controls. The differentially expressed genes (Supplementary Table S1) in the monocytes included those related to cytokines (IL1B, NFKBIA, NLRP3), inflammation (CCR1), cell adhesion (JAML), and apoptosis (IER3, MCC1). Similar to the above results, scRNA-seq dataset from patients with severe COVID-19 and during gout remission were analyzed, and the distribution of identified cell populations is shown in Supplementary Figure S3A. The proportion of monocyte subsets is shown in Supplementary Figure S3B and differentially expressed genes are shown in Supplementary Table S2.

In the following analysis, we continued to focus on the association between severe COVID-19 and gout flares. The intercellular communication was analyzed by the CellChat algorithm. In the severe COVID-19 datasets, the CMs showed high interaction with IMs, NCMs, cDCs, and pDCs (Figure 3A). On the other hand, the major interaction partners of CMs in gout flare were IMs, NCMs, cDCs, and neutrophils (Figure 3C). We also identified several ligand-receptor pairs, and found that LGALS9 and CD45 were the predominant interacting pair in both diseases. RETN and CAP1 mediated the interactions between CM and the other cell subsets. In addition, MIF and (CD74+CXCR4) mediated significant interactions between CMs and B cells in the COVID-19 samples (Figure 3B). In the gout dataset, CMs and neutrophils showed strong interaction via ANXA1 and FPR1(Figure 3D).

The genes associated with severe COVID-19 infection and gout risk were identified by MR analysis as described in the methods. The volcano map highlights genes associated with COVID-19 that have significant p-values (Figure 4A). Forest plot shows the effect of NLRP3-associated SNPs on the risk of COVID-19 (Figure 4B). As shown in Figure 4C, higher expression levels of NLRP3 and CD93 were associated with an increased risk of severe COVID-19, whereas IL17RA was negatively correlated with severe COVID-19 risk. NLRP3 showed the association with the risk of severe COVID-19 (OR 1.12, 95% CI 1.01–1.25, P = 3.70 × 10⁻²). The funnel plot was drawn to assess heterogeneity (Figure 4D), and scatter plot of SNP effects on exposure and outcome is shown in Figure 4E. Furthermore, 4 genes were significantly associated with gout (Figure 5A). Forest plot shows the effect of IER3-associated SNPs on the risk of gout (Figure 5B). IER3 was found a positive correlation, whereas LRP1, MCL1 and RBP7 were negatively correlated with the risk of gout (Figure 5C). IER3 showed the association with gout risk (OR 1.2, 95% CI 1.04–1.39, P = 1.17 × 10⁻²). The scatter plot and funnel plot are shown in Figures 5D, E. MR-Egger test did not show evidence of multiplicity of effects in both diseases (Supplementary Table S3; Supplementary Table S4). Cochran^’s Q-test results showed no heterogeneity concerning severe COVID-19 but the analysis of gout showed heterogeneity. The results are detailed in Supplementary Tables S5, S6. We used an online tool (https://sb452.shinyapps.io/power) to calculate the minimum sample size required for this MR analysis. Calculations suggested the minimum sample size for COVID-19 and gout is 2,387,400 and 39,900, respectively, to reach 80% power.

We selected the SNP with most significant association with NLRP3 and IER3 as a representative for co-localization analysis. Colocalization analyses between NLRP3 and severe COVID-19 showed, PPH3 = 0.31, PPH4 = 0.11, and PPH3+PPH4 = 0.42. Colocalization analyses between IER3 and gout showed, PPH3 = 0.99, PPH4 = 0.00, and PPH3+PPH4 = 0.99. NLRP3 and severe COVID-19 were possibly associated with SNP loci in the genomic region, but were affected by different causal variants. On the other hand, IER3 and gout were significantly associated with SNP loci in the genomic region, and were affected by different causal variants. The genetic variants near the region of NLRP3 association with COVID-19 were shown in the regional association plot (Figure 6A), and the lead SNP was rs59215952. The regional association plot between IER3 and gout is shown in Figure 6B, and the lead SNP was rs35267732. We also performed SMR analysis to verify the causal relationship between NLRP3 and severe COVID-19, and the causal relationship between IER3 and gout. The results showed that the two probes of NLRP3 gene were ILMN_2310896 (P _SMR = 0.1324, P _HEIDI>0.05, β = 0.1325, OR = 1.1417), and the topSNP was rs12143966; and ILMN_1712026 (P _SMR<0.05, P _HEIDI>0.05, β = 0.1822, OR = 1.2000), and the topSNP was rs10925027. Meanwhile, we found that the expression of IER3 gene was associated with gout and there was no significant heterogeneity in the eQTL signaling (P _SMR<0.05, P _HEIDI>0.05, β = 0.0970, OR = 1.1018). The SMR results are shown in Table 1.

As shown in the UMAP plot (Figure 7A), the monocytes expressed higher levels of CD93, NLRP3, and IL17RA compared to the other cell populations. Furthermore, NLRP3 was highly expressed in the CMs, and showed lowest expression levels in the NCMs (Figure 7B). The network of interactions between the different subsets further showed strong interaction of the CMs with IMs, NCMs and cDCs (Figure 7C). The genes related with gout including IER3, LRP1, MCL1 and RBP7, were predominantly expressed in monocytes (Figure 8A), of which IER3 was expressed at higher levels in CMs compared to the other cell populations (Figure 8B). The network of interactions was shown in Figure 8C.

Results of the ligand-receptor interaction analysis showed that the NLRP3-positive CMs had more inter-cellular interactions compared to the NLRP3-negative CMs (Figure 9A). Furthermore, MIF and (CD74⁺CD44), MIF and (CD74+CXCR4) showed more cell-to-cell interactions. As shown in Figure 9B, the IER3-positive CMs showed more interactions than IER3-negative CMs involving MIF and (CD74+CD44), MIF and (CD74+CXCR4), RETN and CAP1, and RETN and TLR4 in gout flare. As a multifunctional cytokine, MIF mainly interacts with the CXC family of receptors (such as CXCR2, CXCR4, CXCR7, CXCR12), CD74, and CD44 to activate the downstream signaling pathways.

The metabolic activities in the cell populations were identified using scMetabolism. Notably, glycosphingolipid biosynthesis-related pathways, riboflavin metabolism and cytochrome P450-related metabolic pathways were enriched in the CMs. Sphingolipids are indispensable components of cell membranes, and not only maintain membrane stability but also regulate cell adhesion, proliferation, apoptosis, differentiation, and recognition. Sphingolipids can act as exogenous or endogenous ligands for some immune cells and stimulate the secretion various cytokines, including TNF, IL-6, and IL-8, and therefore play key roles in inflammation-related and autoimmune diseases. Riboflavin plays a critical role in mitochondrial energy metabolism, redox homeostasis, cell apoptosis and inflammatory response. Cytochrome P450 is widely involved in pathological processes such as apoptosis, inflammatory response, oxidative stress and endoplasmic reticulum stress. NLRP3-positive CMs showed higher porphyrin and chlorophyll metabolism in the COVID-19 patients (Figure 10A), while the glycosaminoglycan biosynthesis-keratin sulfate pathway was enriched in the IER3-positive CMs in gout patients (Figure 10B). In contrast, the IER3-negative and NLRP3-negative CMs showed increased biosynthesis of glycosylphosphatidylinositol-anchored proteins (GPI-Aps). These results are indicative of the metabolic differences between the NLRP3/IER3-positve and NLRP3/IER3-negative CMs. Furthermore, we created a heat map to illustrate the levels of NLRP3, CD93 and IL17RA genes expression by using bulkRNA-seq data from COVID-19 and healthy samples (Figure 10C). The bulk RNA-seq data confirmed the significantly higher expression of NLRP3 in the severe COVID-19 patients compared to the healthy samples (Figure 10D), indicating that NLRP3 plays a key role in disease progression.