Introduction

Cholecystitis refers to a group of inflammatory diseases with different etiologies and clinical courses¹. While the corresponding pathological mechanisms of cholecystitis are complicated, systemic inflammation is widely recognized as a primary catalyst of cholecystitis pathology^2–6. Cholecystectomy is typically the preferred treatment option for patients with severe cholecystitis and gallstones, whereas medication is prescribed for those with mildly symptomatic cholecystitis or calculous cholecystitis. However, a substantial percentage of patients face the possibility of complications or recurrence^7,8. Compelling evidence supports that the proportion, gene expression profile, and functional status of immune cells hold significant prognostic value^9–14. Therefore, it is critical to comprehensively comprehend the immune cell dysfunction linked with cholecystitis, its corresponding altered characteristics, and their prognostic implications.

Natural killer (NK) cells are the major innate lymphocyte subset and are found in most organs^15,16. They produce lytic granules that eliminate infected or cancerous cells, coordinate the inflammation process, form immunological memory, and modulate antigen-presenting cell function in inflammatory diseases^17–19. In the context of tumors and chronic infections, NK cells exhibit an exhausted phenotype, termed NK cell exhaustion, which is characterized by downregulated expression of activating receptors (e.g., NKG2D, FCGR3A), upregulated expression of inhibitory receptors (e.g., KLRC1, PD-1, TIGIT, Tim-3), decreased production of effector cytokines (e.g., IFNγ, GZMB), and impaired cytolytic activity^20,21. Recently, blockade of KLRC1 and TIGIT has been shown to reverse NK cell exhaustion and enhance the therapeutic effects of anti-tumor and anti-infective immunotherapies^21–23. In addition, NK cell exhaustion has been reported to predict poor prognosis in hepatocellular carcinoma and acute myeloid leukemia; however, the correlation of NK cell exhaustion with prognosis in cholecystitis and other inflammatory diseases remains largely unexplored^24–26.

The liver plays a critical role in immune defense. It contains a substantial population of diverse immune cells, particularly NK cells^27–29. The liver modulates these cells to mount immune responses upon encountering specific inflammatory triggers^28,30. Furthermore, liver pathologies affect mesenchymal immune cell function and influence systemic inflammatory responses, impacting the development and outcome of inflammation-associated diseases such as cholecystitis^31–33. Wilson’s disease (WD) is a rare genetic disorder caused by mutations in the gene encoding the transmembrane copper-transporting P-type ATPase protein³⁴. It leads to deficits in mitochondrial structure and function, as well as specific alterations in liver metabolism, resulting in a variety of liver diseases. Additionally, WD is commonly associated with inflammatory conditions in the clinic^35–38. Thus, WD, with its unique etiology, serves as an excellent model for exploring the regulation of immune cell function by hepatic metabolic abnormalities.

In this retrospective clinical study of over six hundred WD patients, it was found that they had a significantly higher cholecystitis incidence than the general population. This, in turn, is associated with a poor prognosis for cholecystitis in these patients. A thorough examination of the immune cell landscape in the hepatic mesenchymal stromal microenvironment of WD patients was performed using single-cell RNA sequencing (scRNA-seq). Specifically, there were significant changes in the composition and function of the innate immune system, including an increase of mononuclear phagocytes (MPs) and NK cells, as well as enhanced pro-inflammatory characteristics. Additionally, there was an enhancement in the biological processes of antigen presentation, activation of immune response, and activation of lymphocytes. One key event among these changes was the exhaustion of NK cells, as demonstrated by an increase in the expression of inhibitory receptors KLRC1 and TIGIT and a decrease in the expression of cytotoxic molecules. Multiplex immunohistochemistry and flow cytometry techniques confirmed the exhaustion of NK cells in both clinical liver tissue and peripheral blood samples. Finally, our bioinformatics analysis has confirmed a positive correlation between NK cell exhaustion and poor prognosis in cholecystitis and other inflammatory diseases. Our study has demonstrated abnormal liver mesenchymal immune cell function triggered by specific metabolic dysfunction in WD, focusing on NK cell exhaustion and its correlation with poor healing of cholecystitis. These findings provide foundation for further studies of effect of hepatocytes metabolic changes on immune cell function and insights into improvement of inflammatory diseases by assessing immune cell function.

Results

WD causes poor prognosis of cholecystitis

WD is a genetic disorder resulting from mutations in the ATP7B gene, exhibiting significant genetic heterogeneity with over 600 mutations or polymorphisms³⁹. Through sequencing DNA samples from 20 patients with WD, we found that c. 2333G>T (p. R778L) accounted for the majority (65%) of the three most common mutations in Chinese, while c. 2804C>T (p. T935M) and c. 2975C>T (p. P992L) were not detected (Figure 1a)⁴⁰. In hepatocytes, ATP7B deficiency in the trans-Golgi network (TGN) leads to abnormalities in mitochondrial function and copper metabolism, and ultimately to a wide range of liver diseases (Figure 1b)^41–44. Typical clinical manifestations are liver fibrosis and Kayser-Fleischer Ring in the cornea that caused by Cu accumulation in the cornea (Figure 1c, d)⁴¹. Additionally, a lower total bilirubin and metabolic symptoms of WD patients, a lower risk of hyperglycemia (Cohort 1 and Cohort 2), are observed (Figure 1e, Figure S1a). Mitochondrial stress and glycolysis assays revealed that knockout of the ATP7B (ATP7B-KO) in human hepatocyte cell lines WRL 68 significantly impaired the basal oxygen consumption rate (OCR) and capacity of ATP production, while the extracellular acidification rate (ECAR) of glycolysis and glycolysis capacity remained unchanged (Figure S1b, c). These results demonstrate that WD is a disease characterized by abnormal hepatocyte metabolism.

Figure 1.

Interestingly, chronic cholecystitis emerged as the most common complication, with a notably higher incidence in WD patients (49.51%) compared to HBV patients (9.46%) in Cohort 2 (Figure 1e). Despite the comparable incidence of gallbladder in the two groups and the significantly reduced plasma concentration of total bilirubin in WD patients (Cohort 1 and Cohort 2), this finding holds (Figure 1e). Moreover, in cholecystitis patients, the proportion of cases with abnormal index was significantly higher when combined with WD, which included white blood cell (WBC), neutrophil (NEU), and lymphocyte (LYM) counts (Cohort 3, Figure 1e). Overall, these findings suggest that cholecystitis outcomes are negatively affected by WD.

Single-cell profiling of non-parenchymal cells in the liver of WD patients

The hepatic mesenchymal cells from three WD patients and three liver hemangioma patients were analyzed using the 10X Genomics platform for scRNA-seq (Figure 2a) to better understand the complicated immune cell abnormalities in WD patients. 57,384 high-quality single transcriptomes were generated, with 30,068 cells in the group of WD patients (CASE) and 27,316 in the group of liver hemangioma patients (CON). We validated twenty-six principal components via PCA and substantiated the grouping through UMAP dimensionality reduction cluster graph (Figure 2b, c, Figure S2a). We then utilized individual cell expression profiles and reference to cell marker gene expression to implement an automated annotation method, which was further complemented by manual review and fine-tuning. This initial annotation process resulted in the classification and display of 12 significant cell clusters. The UMAP algorithm was used to segregate and display the clusters. The study identified various cell types, including hepatocytes^45,46, cholangiocytes⁴⁵, endothelial cells (ECs)^47–49, hepatic stellate cells (HepSCs)^46,50, proliferating cells^48,51,52, B cells^53,54, plasma cells^52,53, T and NK cells^53,55,56, neutrophils⁵⁷, mast cells^48,55,57, mononuclear phagocytes (MPs)⁵⁸, and plasmacytoid dendritic cells (pDCs)^48,59 (Figure 2d, e and Figure S2b, c). Figure S2d exhibited the top 10 differently expressed genes (DEGs) for each cell type.

Figure 2.

Based on the initial annotation, there were notable variations in the quantities and ratios of different cell types between CASE and CON. Our study observed that T and NK lineages were the primary immune cell population among all the examined samples, accounting for 56.54% in CASE and 73.55% in CON, followed by MPs which constituted 17.33% in CASE and 9.97% in CON (Figure 2f, Figure S2e, and Supplementary Table S1). The significant increase in the proportion of MPs and decrease in T and NK cells suggest that MPs likely dominate the immune microenvironment remodeling in WD patients, while lymphocytes play a minor role. B cell proportions (4.53% in CASE and 2.27% in CON), plasma cell proportions (2.36% in CASE and 0.42% in CON), and pDC proportions (0.42% in CASE and 0.16% in CON) were elevated, while neutrophil proportions (3.53% in CASE and 6.45% in CON) declined. The proportion of remaining cells is shown in Supplementary Table S1.

An overview of the immune microenvironment in WD patients

To demonstrate dissimilarities in the composition of the immune microenvironment between the two groups of individuals in a more detailed manner, a secondary subtype clustering and annotation process was applied to MPs and T and NK cells (Figure 3a, b, and Figure S3a, b). MPs are myeloid immune cells that serve phagocytic and antigen-presenting functions. These cells, which include monocytes, macrophages, and dendritic cells, play an essential role in determining the nature of the immune response⁶⁰. A total of 6,962 MPs were identified (4,439 in CASE and 2,523 in CON), and they were subdivided into 7 clusters through annotation (Figure 3a, Figure S3a, and Supplementary Table S1). Cell marker genes for annotating cell subtypes were depicted in Figure S3c and d, and Figure S3e demonstrated the top 10 DEGs between the two groups^{46,48,49,51,52,55–59,61–66}. The DEGs were related to various cell types, including proliferating cells, macrophages, monocytes, mature dendritic cells (MatureDCs), conventional type 1 dendritic cells (cDC1), conventional type 2 dendritic cells (cDC2), and Kupffer cells (KCs). In CASE, there was a significant increase in both macrophages (25.03% in CASE versus 14.33% in CON) and KCs, the resident macrophages in the liver (26.03% in CASE versus 22.96% in CON), suggesting their involvement in promoting hepatic inflammation in metabolic liver disease (Figure 3a)⁶⁷. By contrast, the percentages of monocytes (35.45% in CASE and 41.78% in CON), cDC1 (2.14% in CASE and 5.99% in CON), and cDC2 (10.39% in CASE and 14.13% in CON) were decreased in CASE (Figure 3a).

Figure 3.

Seventeen clusters were obtained through reclustering for the largest population of T and NK cells, totaling 34,339 cells (15,403 cells in CASE and 18,936 in CON)^{53,57,68–83}. The clusters included ILC3_IL1R1, NKT_GNLY, NKT_IFNG, NKT_KLRG1, NK_FCER1G, NK_KLRC1, NKT_XCL2, NK_FCGR3A, CD4NaiveT_CCR7, CD4Tmem_CD40LG, CD4Treg_CTLA4, CD8MAIT_KLRB1, CD8MAIT_SLC4A10, CD8Teff_GZMH, CD8Teff_GZMK, CD8Teff_ISG15, and CD8Teff_MT1X (Figure 3b, Figure S3b, Supplementary Table S1). The proportion of NK cells (40.97% and 17.07% in CASE and CON, respectively) and CD4Treg (0.92% and 0.31% in CASE and CON, respectively) was significantly increased in CASE. On the other hand, a decrease of NKT cells (28.72% in CASE and 35.26% in CON), CD8MAIT (7.76% in CASE and 22.14% in CON) and CD8Teff (14.95% in CASE and 17.88% in CON) was observed in CASE. CD4 regulatory T cells are thought to inhibit T cell activation in autoimmune and virus-related liver diseases, indicating a correlation between elevated CD4 regulatory T cells and reduced CD8 T cells in WD patients⁸⁴. The distribution of other types of cells between the two groups is shown in Figure 3b and Supplementary Table S1.

To systematically evaluate the key functional alteration, we analyzed the top 20 differently expressed genes in all immune cell populations. We found that the expression of genes encoding complement C1q (C1QA, C1QB, and C1QC) and HLA-DRA was significantly upregulated in CASE. Conversely, the expression of genes encoding heat shock protein (HSPA1A, HSPA1B, HSPM1, HSPA6, HSP90AB1, HSPE1, and HSPD1) was markedly downregulated in CASE (Figure 3c). The Gene Ontology (GO) analysis results demonstrate that WD enhances the abilities of antigen-presenting cells to activate immune response and related lymphocytes, through upregulating antigen processing and presentation, activation of immune response, and positive regulation of lymphocyte activation. However, heat shock proteins delivering antigens may be responsible for the downregulation of lymphocytes differentiation and T cell activation⁸⁵. Next, we analyzed the cellular communication network using Cellchat. The primary signal senders and receivers were shown, with their respective strength. Our findings reveal that CD8Teff, KCs, NKT, CD8MAIT, monocytes, cDC2, and macrophage were the primary participants in intercellular communication, with the high value in both incoming and outgoing interaction strength (Figure 3e). Furthermore, the incoming interaction strength of CD8Teff and NK cells was higher in CASE, suggesting the regulation of biological process by extracellular signals, while neutrophils and monocytes exhibited the lower interaction strength in both outgoing and incoming (Figure 3e). Collectively, these results demonstrate changes in the composition and effects of innate immunity in the liver of WD patients.

The dysfunction of main immune cells in WD patients

To obtain a detail analysis of function comparison in immune cells between two groups, we scored and evaluated the core functional signature of major cell types. We further clustered macrophages, which accounted for 14.33% of the MPs in CON and increased to 25.03% in CASE, resulting in the identification of four distinct subtypes (Figure S4a). Interestingly, these four subtypes in the WD patients exhibited distinct gene expression profiles. Macrophages_3 in CASE displayed elevated expression levels of S100A family genes (S100A8, S100A9, S100A12, S100A6, VCAN, FCN1, and S100A4), indicating a potential contribution to MDSC-like macrophages. In contrast, Macrophages_4 in CASE were characterized by upregulation of lipid metabolism genes (FABP5, APOC1, CTSD, and PLA2G7) (Figure S4b)^86,87. By scoring various gene signatures, we compared the M1 polarization, M2 polarization, pro-inflammatory signature, and interferon response signature of four subtypes within two groups (Figure 4a). These four subtypes displayed consistently overlapped M1/M2 signatures. Significantly, the Macrophages_3 subtype exhibited a more pronounced pro-inflammatory signal in CASE, while the other subtypes showed a slight decrease. In addition, the interferon response signature showed slight elevation in four subtypes in CASE. Kupffer cells, which are liver-resident macrophages with unique functions in maintaining homeostasis, represent the most abundant population of tissue-resident macrophages in the human body⁶⁷. KCs were also further clustered and three subtypes were obtained (Figure S4c). Significant differential expression of genes related to the inflammatory response, including IL1B, CXCL2, CCL3C1, CXCL8, and CCL2, were identified between the two groups (Figure S4d). Similarly, KCs subtypes in two groups also showed overlapping M1/M2 signatures, but all subtypes in CASE showed varying degrees of increase in pro-inflammatory and interferon-responsive signatures (Figure 4b).

Figure 4.

Neutrophils, the most abundant leukocytes in human blood, respond to inflammatory stimuli by migrating to the site of inflammation. Its immunological functions, such as phagocytosis, reactive oxygen species generation, degranulation, and production of cytokines and chemokines, play a crucial role in the resolution of inflammation and tissue repair^88,89. In CASE, a lower percentage of neutrophils (3.53% in CASE and 6.45% in CON) was observed, and the functional characteristics were subsequently evaluated. The results showed that the maturation, chemotaxis, phagocytosis, and chemokine activity scores were decreased in CASE, while the activation of the type I interferon signaling pathway was enhanced (Figure 4c).

Since the biological process of antigen processing and presentation was upregulated, we calculated the relative expression intensity of genes encoding MHC molecules and complement molecules in cDC1, cDC2, and mature DCs (Figure 4d). Notably, genes encoding MHC class D molecules (HLA-A, -B, -C, -E, and -F) were significantly upregulated in the mature DCs in CASE, while genes encoding MHC class D molecules (HLA-DQA1, -DQA2, -DQB1, -DQB2, -DPA1, -DPB1, -DRA, -DRB1, and -DRB5) were all downregulated in the mature DCs. Besides, genes encoding MHC molecules were not the top significant DEGs in B cells, while the immunoglobulin genes (IGHD and IGHM) and the naïve B cell marker TCL1A were significantly upregulated in CASE (Figure S4e). Then the immunoglobulin genes (IGHA1, IGHA2, IGHG1, IGHG2, IGHG3, IGHG4, IGHD, and IGHM) expression was compared in B cells and plasma cells (Figure 4e). IGHA1, IGHA2, IGHG1, IGHG2, IGHG3, and IGHG4 expression were unchanged in B cells between the two group and all of examined genes were downregulated except IGHG3 in plasma cells in CASE.

Monocytes accounted for the highest proportion (35.45% in CASE and 41.78% in CON), and most of them were CD14⁺ classical monocytes (Figure S4f). Compared to CON, S100A8, S100A9, and S100A12 were the top 3 upregulated genes in classical monocytes of CASE (Figure S4g). S100A8 and S100A9 are confirmed to play vital roles in modulating the inflammatory responses by recruiting leukocytes and inducing cytokine secretion during inflammation, and as an alarmin, their elevation indicates inflammation in the local environment^90,91. GO analysis of DEGs in monocytes showed a decrease in biological processes involving cytokine production, response to lipopolysaccharide and bacteria, as well as lymphocyte activation (Figure 4f).

For T cells, we performed GO analysis on DEGs in CD8 T cells between two groups. Notably, the biological processes linked to energy and lipid metabolism were upregulated in several subtypes, including CD4Tmem_CD40LG, CD4Treg_CTLA4, CD8MAIT_KLRB1, CD8MAIT_SLC4A10, CD8Teff_GZMH, and CD8Teff_GZMK. In addition, biological processes associated with response to copper ion were upregulated in the CD8Teff_MT1X subtype (Figure S4h). Meanwhile, T cell subtypes showed comparable cytotoxicity scores, and there were slight changes between two groups (Figure 4g). However, CD8 T cells in CASE exhibited lower levels of exhaustion, except for CD8Teff_GZMH and CD8Teff_ISG15 (Figure 4g). It was also discovered that while the proportion of CD4Treg_CTAL4 was increased in CASE, the Treg signature was slightly lower (Figure 4h). Therefore, these results suggest that WD patients exhibit variations in immune cell function and certain biological processes.

NK cell exhaustion in WD patients

The KEGG analysis of DEGs in NK cells between CASE and CON revealed down-regulation of the natural killer cell-mediated cytotoxicity, along with the NF-kappa B signaling pathway and chemokine signaling pathway, in CASE (Figure S5a). Interestingly, GO analysis showed upregulated biological processes of oxidative phosphorylation and ATP metabolic process, implying the altered metabolic program in NK cells and immune microenvironment (Figure S5b). Since the primary effect of NK cells is cytotoxicity on target cells, we compared the gene expression levels of activating receptors, inhibitory receptors, and effector molecules between two groups. The results showed significantly upregulated expression in KLRC1 and TIGIT, and significantly downregulated expression in FCGR3A, TNF, IFNG, GZMB, GNLY, and CCL4 (Figure 5a). NK cells were clustered into three NK subtypes, NK_FCER1G, NK_KLRC1, and NK_FCGR3A, by gene expression profile. Analysis of the expression of these genes in three NK subtypes yielded similar results, except for increased expression of GZMB and GNLY (Figure S5c). The results were consistent across the three subtypes, suggesting that the results in total NK population were contributed by all three subtypes and not affected by a single composition. Additionally, the proportion of NK cell was higher in CASE, while the proportion of NK_FCER1G and NK_KLRC1 was higher and the proportion of NK_FCGR3A was lower in CASE (Figure 3b and Figure S5d). The respective transcriptome profile was compared in CASE and CON. Consistently, KLRC1 was top-ranking in all three subtypes among the significantly upregulated DEGs, while CCL4 and genes encoding heat shock protein were top-ranking in all three subtypes among significantly downregulated DEGs (Figure S5e). KLRC1 and TIGIT are inhibitory receptors expressed on NK cell surface, and their activation inhibits NK cell functions^21,92. It was collectively indicated that WD patients have NK cells with exhausted status, which results in poor effector function.

Figure 5.

In addition, trajectory analysis showed that NK cells in two groups exhibited a distinct differentiation trajectory. Most of the NK cells in CON differentiated towards cell fate 1, while those in CASE displayed a higher tendency for differentiation towards cell fate 2, with cell grouping occurring at the final stage (Figure S5f). To ascertain the gene expression profile and the respective cell statuses of the two outcomes, we conducted pseudotime analysis of gene expression in two differentiating directions. The results showed that the gene expression of KLRC1 (included in gene cluster 1) was gradually increased along the way to cell fate 2, while the gene expression of FCGR3A, TNF, and GZMB (included in gene cluster 3) was weakened in the same direction, suggesting that cell fate 2 is a more exhausted state of NK cells (Figure S5g, see full gene list in Supplementary Table S2). Consequently, these results demonstrate that NK cells differentiate towards an exhausted status in WD patients. The development, maturation, and function of NK cells are strictly regulated by transcription factors (TFs)⁹³. Next, to assess the transcription factors (TFs) underlying differences in NK cells between CASE and CON, we applied PYSCENIC analysis in NK cells. We found that the gene expression of T-bet (encoded by TBX21) was downregulated and Eomes (encoded by EOMES) was upregulated in CASE (Figure S5h). These two constitute a pair of TFs that play significantly complementary roles in directing the transcriptional program of NK cells. Consistent with previous studies, T-bet expression was downregulated in exhausted NK cells, whereas Eomes expression was increased^24,94,95.

To validate NK cell exhaustion in WD patients, we performed multiplex immunohistochemistry (mIHC) staining on liver tissue sections and flow cytometry on peripheral blood samples. The mIHC results indicated that KLRC1⁺ cells were more frequently observed in NK cells marked by CD56 in the liver tissue sections of WD patients (Figure 5b). Moreover, flow cytometry was also performed to test the number and indicators expression of NK cell from the peripheral blood. After isolation from peripheral blood, NK cells were determined by CD3^-CD56⁺ and NK cells with positive indicators (KLRC1, TIGIT, and IFNγ) were recorded (Figure 5c). The results showed a significant increase in the percentage of KLRC1⁺ and TIGIT⁺ NK cells and a significant decrease in the percentage of IFNγ⁺ NK cells in peripheral blood lymphocytes from WD patients (Figure 5d). Meanwhile, the proportion of NK cells in lymphocytes from peripheral blood and the proportion of CD56^Bright NK cells in total NK cells were also calculated. The significantly higher proportion of NK cells in WD patients was consistent with the scRNA sequencing data, and more CD56^Bright NK cells suggested a relatively weaker effector function in WD patients (Figure S5i). Therefore, these results indicated that, despite the increased proportion and the upregulated oxidative phosphorylation, NK cell exhibited exhausted status in WD patients.

NK cell exhaustion predicts poor prognosis of cholecystitis

Cholecystitis patients with WD exhibit poorer clinical outcomes. There is evidence demonstrating that NK cell exhaustion is correlated with poor liver cancer prognosis^24,26,96. Thus, we hypothesize that NK cell exhaustion also contributes to the worse prognosis of cholecystitis patients with WD.

Prognostic models for inflammatory related diseases in GEO were developed using machine learning methods, including random survival forest, ROC curve, and Kaplan-Meier survival curves. A logistic regression algorithm was used to construct a multigene prediction model based on GSE166915 (cholecystitis) and other inflammatory diseases, namely GSE91035 (pancreatitis), GSE75037 (lung adenocarcinoma), and GSE41804 (hepatocellular carcinoma). Using stepwise regression analysis, 12 genes, KLRC1, TIGIT, HAVCR2, PDCD1, CD96, LAG3, FCGR3A, NCR3, NCR2, NCR1, KLRK1, and PRF1, were selected to obtain the best model. The results showed that the predictive model constructed from these genes had good prognostic performance, with AUC of 0.77, 0.80, 0.97, 0.90 (Figure S6a-d). The AUCs of the models in GSE75037 and GSE41084 were relatively high, 0.97 and 0.90, respectively.

Univariate Cox regression analysis revealed that the TIGIT (in GSE166915) and PDCD1 (in GSE75037) were risk genes with hazard ratio (HR) > 1, while NCR2 (in GSE91035), HAVCR2, KLRC1, NCR3 (in GSE75037), NCR3, LAG3 (in GSE41804) were protective gene with HR < 1 among 12 differently expressed markers of NK cell exhaustion (Figure S6a-d). Consistently, the Kaplan-Meier survival curves suggested a shorter survival time in the group with high markers of NK cell exhaustion than in the group with low markers in the GSE166915 (P= 6.1e-6), GSE91035 (P= 1.9e-5), GSE75037 (P= 3.8e-16), and GSE41804 (P= 3.0e-3) datasets (Figure 6a-d). Patients with cholecystitis, pancreatitis, lung adenocarcinoma, and hepatocellular carcinoma who show high markers of NK cell exhaustion are associated with poor survival compared to those with low markers of NK cell exhaustion.

Figure 6.

Discussion

Immune cell dysfunction impacts the development, progression, and prognosis of diseases and serves as the basis for immunotherapy. The current understanding of immune cell dysfunction in inflammatory diseases is limited, which hinders the development of effective therapeutic strategies. Our retrospective clinical analysis indicates that WD contributes to the unfavorable prognosis of cholecystitis. We developed an atlas of immune cell dysfunction linked to specific metabolic disruption in WD patients, including 19 major cell types. The immune microenvironment change primarily encompasses increased proportion and strengthened function in innate immunocytes, including upregulated antigen processing and presentation and activation of immune responses. Furthermore, we conducted a comprehensive assessment and comparison of immune cell function in WD patients and control group. This includes an evaluation of pro-inflammatory and immune regulatory function in macrophages and KCs, phagocytosis and chemokine activity in neutrophils, and cytotoxicity in T cells, etc. We identified NK cell exhaustion in the WD patients, with upregulation of KLRC1 and TIGIT and downregulation of CD16, TNF, IFNG, GZMB, GNLY, and CCL4. The validation of our findings in patient liver tissues and blood samples was conducted using immunohistochemistry assays and flow cytometry. Our findings indicated that the NK cell exhaustion signature has the ability to predict and prognosticate inflammatory diseases such as cholecystitis, pancreatitis, lung adenocarcinoma, and hepatocellular carcinoma. These results unveil novel perspectives on the immune cell dysfunction role in inflammatory diseases.

WD represents a model of metabolic abnormality in hepatocytes caused by ATP7B mutation. Therefore, immune cell dysfunction was induced by the extracellular signals within the reshaped tissue microenvironment. Emerging research indicates that metabolic crosstalk in the tumor microenvironment disrupts metabolic and signaling pathways in immune cells, leading to dysfunction^97–101. Targeting metabolic pathways to enhance the efficacy of immunotherapy is a promising approach. Our study found multiple metabolic abnormalities in immune cells, including the widely upregulated ATP metabolic process in T and NK cells. Interestingly, glycolysis and oxidative phosphorylation contribute to maintain the cytotoxic ability and promote expansion of NK cells^102–104. However, high rate of oxidative phosphorylation, increased proportion, and exhausted status were simultaneously observed in NK cells of WD patients. Further studies are needed to elucidate the metabolic landscape, mechanisms underlying metabolic alterations, and their role in immune cell dysfunction.

Methods

Patients and clinical samples

All the clinical data were obtained from The First Affiliated Hospital of Anhui University of Chinese Medicine. The whole blood samples and liver tissue samples were collected from patients at The First Affiliated Hospital of Anhui University of Chinese Medicine. The study was approved by the Ethics Committee of the First Affiliated Hospital of Anhui University of Traditional Chinese Medicine and the Ethics Committee of the First Affiliated Hospital of USTC.

Sanger sequencing

The sanger sequencing on 20 WD patients were conducted by Sangon Biotech (Shanghai) Co., Ltd. The primers to detect mutations are: R778L-F, 5’-AGCCTTCACTGTCCTTGTCTTTC-3’; R778L-R, 5’-ATTAGCTGGGATTTCAGAAGTAGTG-3’; T935M-F, 5’-ATGCTTGTGGTGTTTTATTTCTTC-3’; T935M-R, 5’-AGAAGCAAGCAAATAAAATGTAATG-3’; P992L-F, 5’-TCTCAACCTGCCTCTGACTCTG-3’; P992L-R, 5’-TGGTGGCTACTCTGTTGCTACTG-3’.

Cell culture, gene knockout, and metabolic assays

WRL 68 cells were cultured using DMEM high-glucose media supplemented with 10% FBS at 37D°C with 5% CO₂. The lentivirus-mediated CRISPR-Cas9 KO vector targeting human ATP7B was obtained from OBiO Technology (Shanghai) Co., Ltd. Cells at 30% to 40% confluency were incubated in a medium containing optimal dilutions of lentivirus, then subjected to puromycin selection to obtain the stable transfected cells. OCRs and ECARs were measured using Seahorse XF Cell Mito Stress Test Kit (103015-100, Agilent) and Seahorse XF Glycolysis Stress Test Kit (103020-100, Agilent) according to manuals in Agilent XFe96 Analyzer, respectively.

Single-cell RNA sequencing

In the cell preparation section: For the quality check and counting of single cell suspension, the cell survival rate is generally above 80%. The cells that have passed the test are washed and resuspended to prepare a suitable cell concentration of 700∼1,200 cells/μL for 10x Genomics ChromiumTM. The system is operated on the machine.

Steps in GEM Creation & Thermal cycling: GEMs (Gel Bead in Emulsion) were constructed for single cell separation according to the number of cells to be harvested. After GEMs were normally formed, GEMs were collected for reverse transcription in a PCR machine for labeling.

Post Cycling Cleanup & cDNA Amplication: The GEMs were oil-treated, and the amplified cDNA was purified by magnetic beads, and then subjected to cDNA amplification and quality inspection.

Library Preparation & Quantification: The 3’Gene Expression Library was constructed with the quality-qualified cDNA. After fragmentation, adaptor ligation, sample index PCR, etc., the library is finally quantitatively examined.

Sequencing: The final library pool was sequenced on the Illumina Hiseq instrument using 150-base-pair paired-end reads.

Cell type annotation and differentially expressed genes (DEGs) analysis

The cell type identity of each cluster was determined with the expression of canonical markers found in the DEGs using SynEcoSys database. The markers and references are shown in Supplementary Table S4. To identify DEGs, we used the Seurat FindMarkers function based on Wilcox likelihood-ratio test with default parameters, and selected the genes expressed in more than 10% of the cells in a cluster and with an average log (Fold Change) value greater than 0.25 as DEGs. For the cell type annotation of each cluster, we combined the expression of canonical markers found in the DEGs with knowledge from literatures, and displayed the expression of markers of each cell type with heatmaps that were generated with Seurat Heatmap function. Doublet cells were identified as expressing markers for different cell types, and removed manually.

Pathway enrichment analysis

To investigate the potential functions of DEGs, the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis were used with the “clusterProfiler” R package 3.16.1¹⁰⁵. Pathways with p_adj value less than 0.05 were considered as significantly enriched. Gene Ontology gene sets including molecular function (MF), biological process (BP), and cellular component (CC) categories were used as reference.

Cell-cell Interaction Analysis

CellChat (version 0.0.2) was used to analyze the intercellular communication networks from scRNA-seq data¹⁰⁶. A CellChat object was created using the R package process. Cell information was added into the meta slot of the object. The ligand-receptor interaction database was set, and the matching receptor inference calculation was performed.

Trajectory analysis

Cell differentiation trajectory was reconstructed with Monocle2¹⁰⁷. Highly-variable genes (HVGs) were used to sort cells in order of spatialDtemporal differentiation. We used DDRTree to perform FindVairableFeatures and dimension-reduction. Finally, the trajectory was visualized by plot_cell_trajectory function. Next, CytoTRACE (a computational method that predicts the differentiation state of cells from single-cell RNA-sequencing data using gene Counts and Expression) was used to predict the differentiation potential¹⁰⁸.

Transcription factor regulatory network analysis

Transcription factor network was constructed by pyscenic v0.11.0 using scRNA expression matrix and transcription factors in AnimalTFDB¹⁰⁹. First, GRNBoost2 predicted a regulatory network based on the co-expression of regulators and targets. CisTarget was then applied to exclude indirect targets and to search transcription factor binding motifs. After that, AUCell was used for regulon activity quantification for every cell.

UCell gene set scoring

Gene set scoring was performed using the R package UCell v 1.1.0¹¹⁰. UCell scores are based on the Mann-Whitney U statistic by ranking query genes in order of their expression levels in individual cells.

Because UCell is a rank-based scoring method, it is suitable to be used in large datasets containing multiple samples and batches. The gene sets for signature scoring are listed in Supplementary Table S3.

Multiplex immunohistochemistry

The procedure of deparaffinization, dehydration, antigen retrieval, endogenous peroxidases quench, block, antibody incubation was conducted as standard immunohistochemistry. Solution of tyramide signal amplification (TSA) (WAS1101010, BioMed World) was applied on the slides and antibodies were removed by microwave heating. Repeat block, antibody incubation, TSA, and microwave heating to perform next cycle of staining. DAPI (WAS1301050, BioMed World) was used to stain nuclei and antifade mountant (WAS1302050, BioMed World) was applied. The primary antibody to CD56 (ab75813, Abcam) and KLRC1 (ab260035, Abcam) and second antibody (WAS1201100, BioMed World) were used. Image of immunofluorescence-stained slides were acquired with Nikon Eclipse Ci-L and 3DHISTECH Pannoramic MIDI D.

Flow cytometry

Lymphocytes were isolated from the whole blood using density gradient centrifugal method (P8610, Solarbio). Isolated cells were blocked (564219, BD Pharmingen) and stained with antibodies against CD3 (560835, BD Pharmingen), CD56 (555518, BD Pharmingen), KLRC1 (375114, Biolegend), and TIGIT (568672, BD Pharmingen). For intracellular cytokine staining, cells were stimulated for 4 h with Leukocyte Activation Cocktail and BD GolgiPlug (550583, BD Pharmingen) according to manuals. After stimulation, cells were blocked, stained for surface markers, fixed and permeabilized with Fixation/Permeablization Kit (554714, BD Pharmingen). Fixed cells were stained with antibodies against IFNγ (552887, BD Pharmingen).

Construction and verification of prognostic model

On the basis of the markers of NK cell exhaustion, univariate Cox analysis was conducted in cholecystitis, pancreatitis, lung adenocarcinoma, and hepatocellular carcinoma patients to screen survival-related markers of NK cell exhaustion, and markers with P values < 0.05 were retained. We conducted multivariate analysis to identify the optimal prognostic markers for the prognostic model. The risk scores of the patients were calculated based on the normalized gene expression levels and the Cox regression coefficients of selected markers. Time-dependent receiver operating characteristic (ROC) curves were utilized to verify the prognostic performance of the model for overall survival (OS). For Kaplan-Meier survival, we used the R software package maxstat (Maximally selected rank statistics with multiple p-value approximations version: 0.7-25) to calculate the optimal cutoff value for RiskScore. We set the minimum group sample size to be greater than 25% and the maximum group sample size to be less than 75%, ultimately obtaining the optimal cutoff value. Based on this, patients were divided into high and low groups, and further analyzed the prognostic differences between the two groups using the R software package’s survival function. The logrank test method was used to evaluate the significance of prognostic differences between different groups of samples, and significant prognostic differences were ultimately observed.

Statistical analysis

Statistical analyses were performed using appropriate tests as indicated in legends (unpaired two-tailed t-test, chi-square test, and Wilcox test).

Data availability

sc-RNA sequening data that support the findings of this study have been deposited in Gene Expression Omnibus (GEO) with the accession numbers GSE254082.

Additional information

Author contributions

Y.W., Q.F., H.W., and Q.Y. conceived the concept for this study. Y.J., J.X., L.J., and M.H. analyzed data and implemented experiments. C.D. conducted the bioinformatic analysis. L.J., Q.T., and Z.L. collected the clinical data and clinical samples. Y.J., J.X., and C.D. wrote the manuscript. W.Z. and Q. F. proofread the manuscript.

Fundings

This work was supported by grants from National Natural Science foundation of China (Project: 2021KY340; NSFC: 81773112 and 82073186 for Q.F.). This work was supported by Construction project of the State Administration of Traditional Chinese Medicine in 2021 (Teaching Letter [2021] No. 270 for Q.Y.). and Anhui Province Clinical Key Specialty Project in 2022 (Anhui Weichuan [2022] No. 297 for Q.Y.).

Competing interests

The authors declare no competing interests.