The study design is presented in Figure 1 .

The Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/info/datasets.html), curated by the National Center for Biotechnology Information (NCBI), functions as a comprehensive repository for gene expression data. From this resource, we acquired the single-cell transcriptome data corresponding to GSE189754, concentrating on 11 samples that provided complete single-cell expression profiles for single-cell analysis. Additionally, we downloaded the spatial transcriptome data for GSE189184, selecting two control groups (B10, C5) and two disease groups (B8, B4) for analysis. Furthermore, we procured the transcriptome data for GSE48958, encompassing data from 8 controls and 13 disease samples.

In this study, the processes of quality control and data standardization are essential to ensure the accuracy of subsequent analyses. We employed the Seurat package (11) for initial data processing. For cell quality control, we conducted screening based on the total number of unique molecular identifiers (UMIs) per cell, the number of expressed genes, and the proportion of mitochondrial gene expression. Typically, a high proportion of mitochondrial gene expression in a cell indicates low RNA expression levels and potential progression towards cell death, warranting the exclusion of such cells. Additionally, we utilized the median absolute deviation (MAD) for quality control, removing outliers that deviate from the median by more than three times the MAD to maintain data reliability. Subsequently, we applied DoubletFinder (version 2.0.4) (12) to individually filter doublet cells in each sample, thereby completing the comprehensive cell quality control process.

In the data standardization process, the LogNormalize method of global normalization is employed. This technique mitigates the impact of variations in total RNA content between cells on gene expression analysis by scaling the total expression level of each cell with a coefficient \(s_0\), adjusting it to 10,000, and subsequently normalizing it through logarithmic transformation. Cell cycle scores are computed using the CellCycleScoring function, and highly variable genes are identified via the FindVariableFeatures function. The ScaleData function is utilized to eliminate gene expression fluctuations attributable to mitochondrial gene expression, ribosomal gene expression ratios, and cell cycle differences. Linear dimensionality reduction is conducted on the expression matrix using the RunPCA function, with 20 principal components selected for further analysis. The Harmony algorithm is applied with default parameters to correct for batch effects, and finally, the RunUMAP function is employed with default parameters for nonlinear dimensionality reduction.

Cell types and corresponding marker genes were identified using CellMarker (13), PanglaoDB (14), and literature, supplemented by automated annotation with SingleR (15) software. The FindAllMarkers function was employed to filter marker genes within each category, with only positive markers expressed in at least 25% of the cells retained (only.pos = TRUE, min.pct = 0.25).

CellChat (16) is a sophisticated tool designed for the quantitative inference and analysis of intercellular communication networks derived from single-cell data. Employing network analysis and pattern recognition methodologies, CellChat facilitates the prediction of principal signaling inputs and outputs of cells, elucidating the mechanisms by which these cells and signals orchestrate their functions. In this study, we employed standardized single-cell expression profiles as input data, alongside cell subtype classifications obtained through single-cell analysis, to serve as cell-specific information. We conducted an in-depth examination of cell-related interactions, quantifying the strength and frequency of cell-to-cell interactions to observe the activity and impact of each cell type in the disease.

We utilized the Least Absolute Shrinkage and Selection Operator (LASSO) logistic regression and Support Vector Machine (SVM) algorithms to select features for diagnostic markers of diseases. The LASSO algorithm utilizes the “glmnet” package, while SVM-Recursive Feature Elimination (SVM-RFE) is a machine learning method based on support vector machines (17). SVM-RFE removes feature vectors generated by SVM to identify optimal variables, and establishes a support vector machine model through the “e1071” package to further assess the diagnostic value of these biomarkers in disease contexts.

The CIBERSORT method is a prevalent technique for assessing immune cell types within microenvironments (18). In this study, utilized the CIBERSORT algorithm was employed to analyze patient data, allowing for the inference of the relative proportions of 22 immune-infiltrating cell types. Furthermore, a correlation analysis was conducted to examine the relationship between gene expression and immune cell content.

Patients were categorized into high and low-expression groups based on the expression of key genes. Subsequently, Gene Set Enrichment Analysis (GSEA) was utilized to examine disparities in signaling pathways between these cohorts. The annotation gene set employed for the subtype pathway analysis was derived from version 7.0 of the Molecular Signatures Database (MsigDB). Differential expression analysis of pathways between the groups was conducted, and significantly enriched gene sets (adjusted p-value < 0.05) were ranked by consistency score. GSEA is frequently used to explore the correlation between disease classification and biological significance.

Gene Set Variation Analysis (GSVA) is a nonparametric, unsupervised method for assessing gene set enrichment in transcriptome data. GSVA assigns a comprehensive score to each gene set of interest, converting gene-level changes into pathway-level changes. This allows for the identification of potential biological function changes in different samples. In this study, gene sets were downloaded from MsigDB, and the GSVA algorithm was applied to comprehensively score each gene set, enabling the evaluation of potential biological function differences among the samples.

MicroRNAs (miRNAs) are small non-coding RNAs known to regulate gene expression by facilitating mRNA degradation or inhibiting mRNA translation. Consequently, we conducted an in-depth analysis to determine the presence of miRNAs associated with key genes involved in the transcriptional regulation or degradation of potentially deleterious genes. We identified miRNAs related to these key genes using the miRcode database and subsequently visualized the miRNA-gene interaction network utilizing Cytoscape software (19).

This study utilized the R package “RcisTarget” to predict transcription factors, with all computations conducted by RcisTarget being predicated on motif analysis. The normalized enrichment score (NES) of a motif depended on the total number of motifs in the database. In addition to the motifs annotated by the source data, we inferred further annotation files based on motif similarity and gene sequences. To estimate the overrepresentation of each motif in the gene set, we initially calculated the area under the curve (AUC) for each pair of motif-motif set. This was performed based on the recovery curve calculation of the gene set ranking of the motifs. The NES of each motif was calculated based on the AUC distribution of all motifs in the gene set.

To verify the expression of target genes in the diseased colon tissue of UC patients, tissue biopsy samples were collected from UC patients within the research cohort at the Digestive Endoscopy Center of Changshu Hospital Affiliated to Nanjing University of Chinese Medicine (Ethical Number: CZYLS-2024120). Patients with UC secondary to other diseases or with differing pathological results were excluded. Normal tissue samples for the control group were obtained from the periphery of pathological specimens diagnosed with colon cancer and subjected to Miles surgery in the General Surgery Department of Changshu Hospital Affiliated to Nanjing University of Chinese Medicine. The collection of all samples was approved by the hospital’s Ethics committee, and written informed consent was obtained from the patients.

Colon tissue sections fixed with paraformaldehyde were deparaffinized using xylene and incubated with primary antibodies (tissue inhibitor of metalloproteinase 1 (TIMP1):1:200, Absin, Shanghai, China; G protein subunit gamma 5 (GNG5):1:200, Abcam, Shanghai, China) at 37°C for 1.5 hours. After three washes with PBS, immunocomplex detection was performed using diaminobenzidine, and nuclei were counterstained with hematoxylin. The sections were examined under a microscope (Leica, Wetzlar, Germany) (20).

Male C57BL/6 mice (weighing 18–20 grams) were obtained from Beijing Vital River Laboratory Animal Technology Co., Ltd. (SCXK-2021-0011, Beijing, China). Before the experiments, the mice were provided with standard laboratory chow and water ad libitum under controlled conditions of 60 ± 5% humidity, 23 ± 1°C temperature, and a 12-hour light/dark cycle. The experimental protocol was approved by the Ethics Committee of the Experimental Animal Center at Nanjing University of Chinese Medicine (Ethical Number: NJUCCSHAE-2021-1123). The mice were randomly divided into two groups: a control group and a dextran sulfate sodium (DSS) group, each with 6 mice. The control group received drinking water, while the DSS group was given 3% DSS in drinking water for 7 days (21). On the 8th day, all mice were euthanized, and their colon samples were collected for further analysis.

Colon tissue was fixed in 4% paraformaldehyde, subsequently embedded in dehydrated paraffin, and sectioned at a thickness of 4.5μm. The sections were then stained with H&E. Pathological alterations in the tissue samples were examined using an optical microscope (Leica, Wetzlar, Germany).

Accurately weigh the colon tissue to achieve a weight (mg) to volume (µL) ratio of 1:9. Add nine times the volume of physiological saline and homogenize the mixture mechanically under ice water bath conditions to prepare a 10% homogenate. Centrifuge the homogenate at 2500–3000 rpm for 10 minutes and collect the supernatant for subsequent ELISA analysis. Following the manufacturer’s protocol, the concentrations of TNF-α (mIC50536-1, Mlbio, Shanghai, China) and interleukin-6 (IL-6) (ml098430, Mlbio, Shanghai, China), were quantified using a commercially available ELISA kit.

The sample slices were fixed in 10% formalin, embedded in paraffin, dewaxed, and subjected to antigen retrieval. After a one-hour blocking step at room temperature, the slices were incubated overnight at 4°C with primary antibodies TIMP1 (1:200, Absin, Shanghai, China) and F4/80 (1:50, Abcam, Shanghai, China). Following three 10-minute washes with PBS, the slices were incubated for one hour at room temperature with Alexa Fluor 488 and Alexa Fluor 594 secondary antibodies. After another three PBS washes, an anti-quenching medium was used to mount the cover glass onto the slide. The sections were then examined under a fluorescence microscope (Leica, Wetzlar, Germany) at a magnification of 80 for microscopic analysis and imaging (22).

Total protein was extracted from colon tissue samples of human or mouse origin using RIPA lysis buffer (Beyotime, Nanjing, China) and quantified using the Bicinchoninic Acid (BCA) protein assay kit (Beyotime, Nanjing, China). Subsequently, 20 micrograms of protein were separated on a 10% SDS-PAGE gel and transferred to a polyvinylidene fluoride (PVDF) membrane. The membrane was then blocked with 5% (w/v) bovine serum albumin (BSA) or skim milk at room temperature for 1 hour. Following the blocking step, the membrane was incubated overnight at 4°C with primary antibodies targeting GNG5 (1:1000, Absin, Shanghai, China), TIMP1 (1:1000, Absin, Shanghai, China), and GAPDH (1:5000, Proteintech, Wuhan, China). On the next day, the membrane was incubated with secondary antibodies (horseradish peroxidase-conjugated goat anti-rabbit or anti-mouse IgG, 1:5000, Cell Signaling Technology, Danvers, MA, USA) at room temperature for 1 hour. Visualization of the protein bands was performed using an ECL detection kit and a gel imaging system (Tanon, Shanghai, China). The intensity of the bands was then quantified using the densitometric analysis feature of Gel Pro 4.0 software (Tanon, Shanghai, China) (23).

All statistical analyses were performed using the R programming language (version 4.3.0), with a significance threshold set at p < 0.05.

During the initial processing of single-cell expression profile data, rigorous adherence to established quality control and standardized procedures was maintained. Following the screening process, cells with fewer than 200 captured genes were excluded, while those meeting the criteria were retained based on indicators such as nFeature-RNA, nCount-RNA, and percent.mt, resulting in a dataset of 22,345 high-quality cells. Concurrently, doublets were removed, and the top 2,000 highly variable genes were selected for subsequent analysis. The processed data demonstrated favorable distribution characteristics, as evidenced by violin plots and scatter plots, thereby establishing a robust foundation for precise cell subpopulation annotation and gene expression analysis in future studies ( Supplementary Figure S1A, C ). This approach effectively mitigates analysis bias associated with data quality issues.

The data underwent standardization, homogenization, and subsequent analysis using Principal Component Analysis (PCA), Harmony, and Uniform Manifold Approximation and Projection (UMAP) ( Supplementary Figures S1D-F , Figure 2A ). Each subtype was annotated to one of seven cell categories: CD4⁺ T cells, B cells, CD8⁺ T cells, Fibroblasts, Monocytes, Mast cells, and NK cells ( Figure 2A ). A bubble plot and histogram were generated to visualize the expression of classic markers and cell proportions for these categories ( Figures 2B, C ). The software package Cellchat was employed to examine ligand-receptor interactions within the single-cell expression profile, revealing intricate relationships between the cell subtypes ( Figure 2D ). Notably, Monocytes demonstrated a closer potential interaction with other cell types ( Figures 2E, F ).

To investigate the genetic underpinnings of UC, we retrieved expression profile data from the GEO database (GSE48958). This dataset comprised 21 patient samples. To identify key genes associated with this condition, we employed a two-step approach. First, we utilized a combination of LASSO regression and SVM algorithms to screen the 386 monocyte marker genes (p_adj < 0.05 & LogFC > 0.585) previously identified. LASSO regression yielded 18 characteristic genes ( Figures 3A, B ), while SVM identified 4 ( Figure 3C ). By intersecting these gene sets, we identified 2 genes, GNG5 and TIMP1, as the most promising candidates for further exploration in our research on UC ( Figure 3D ).

The microenvironment, a pivotal factor in disease progression, consists of a complex interplay between cellular and extracellular components. This intricate ecosystem includes fibroblasts, immune cells, extracellular matrix, growth factors, inflammatory factors, and unique physical and chemical properties. The microenvironment exerts a substantial influence on disease diagnosis, prognosis, and therapeutic response. Our investigation revealed distinct patterns of immune cell infiltration and correlation in various disease states ( Figures 4A, B ). Compared to the control group, the disease group exhibited significantly elevated levels of M1 Macrophages, resting NK cells, CD4⁺ memory activated T cells, and CD4⁺ memory resting T cells. Conversely, resting Mast cells and NK cells activated were significantly reduced in the disease group ( Figure 4C ). Further analysis of the relationship between key genes and immune cells demonstrated a strong positive correlation between TIMP1 and several immune cell types, including CD4⁺ memory resting T cells, CD4⁺ memory activated T cells, follicular helper T cells, resting NK cells, M0 Macrophages, M1 Macrophages, activated Dendritic cells, and Neutrophils. Conversely, TIMP1 was negatively correlated with Plasma cells, CD8⁺,T cells regulatory T cells (Tregs), activated NK cells, and resting Mast cells ( Figure 4D ). Moreover, our analysis of the correlation between key genes and different immune factors, including immunosuppressive factors, immunostimulatory factors, chemokines, and receptors, suggests that these genes are intimately involved in shaping the immune microenvironment ( Supplementary Figures S2A-E ).

To elucidate the specific signaling pathways involved in the key genes and explore their potential molecular mechanisms in disease progression, we conducted a comprehensive analysis. GSEA revealed that GNG5 was significantly enriched in signaling pathways such as propanoate metabolism, butanoate metabolism, and peroxisome proliferator-activated receptor (PPAR) signaling ( Figures 5A, B ). TIMP1, on the other hand, was enriched in pathways including B cell receptor signaling, interleukin-17 (IL-17) signaling, and NF-κB signaling ( Figures 5D, E ). Additionally, GSVA identified GNG5 as being enriched in pathways associated with protein secretion and adipogenesis ( Figure 5C ). TIMP1 was found to be enriched in pathways related to hedgehog signaling and epithelial-mesenchymal transition ( Figure 5F ). These findings collectively suggest that the key genes may influence disease progression through these identified signaling pathways.

Subsequently, we employed the miRcode database to conduct a reverse prediction of the key genes, resulting in the identification of 20 miRNAs and a total of 23 mRNA-miRNA regulatory relationships. These interactions were visualized using Cytoscape ( Supplementary Figure S3A ). By utilizing the key genes as a gene set for this analysis, we discovered that these genes were subject to regulation by common mechanisms, such as multiple transcription factors. To identify these transcription factors, we employed cumulative recovery curves and conducted motif-transcription factor annotation and selection analysis on the key genes. The motif with the highest standardized enrichment score (NES: 14) was determined to be cisbp:M6056. We have provided a comprehensive visualization of all the enriched motifs and their corresponding transcription factors associated with the key genes ( Supplementary Figure S3B, C ).

In the current study, the GeneCards database (https://www.genecards.org/) was utilized to identify genes potentially implicated in disease regulation. To assess inter-group expression differences amongst these genes, we analyzed the expression levels of 20 highly ranked genes (based on the Relevance*score) with confirmed expression within the transcriptome data. This analysis revealed significant expression differences between the two patient groups for genes including interleukin 23 receptor (IL23R), interferon gamma (IFNG), nucleotide-binding oligomerization domain 2 (NOD2), tumor protein 53 (TP53), transforming growth factor beta 1 (TGFB1), interleukin 1 receptor antagonist (IL1RN), interleukin-1 beta (IL1B), interleukin 8 (CXCL8), tumor necrosis factor (TNF) and ATP-binding cassette subfamily B member 1 (ABCB1) ( Figure 6A ). Furthermore, a correlation analysis was performed to investigate the relationship between key genes and disease regulation genes. The expression levels of these key genes demonstrated statistically significant correlations with the expression levels of disease regulation genes. Notably, TIMP1 exhibited a strong positive correlation (cor = 0.949) with IL1RN while displaying a significant negative correlation (cor = -0.807) with ABCB1 ( Figure 6B ).

We analyzed the spatial transcriptome data to assess the expression levels of two key genes. Compared with the control group, GNG5 expression was inhibited in the disease group, while TIMP1 expression was significantly upregulated in the disease group ( Figure 7A ). We assessed the differential expression levels of key genes across various groups utilizing bubble and violin plot visualizations. Our analysis revealed a downregulation of GNG5 and an upregulation of TIMP1 in UC ( Figure 7B ). IHC analysis was conducted on colon lesions from patients with UC to assess the expression levels of GNG5 and TIMP1. Results demonstrated a significant upregulation of TIMP1, and a significant downregulation of GNG5 in the disease group compared to the control group (p < 0.001), aligning with the previously presented spatial transcriptome data ( Figure 7C ).

A DSS-induced colitis mice model was generated successfully ( Figure 8A ). The distal colon tissue was collected and its length was measured and then imaged. It was found that compared with the control group, the colon of the DSS-induced colitis mice model was significantly shortened under inflammatory stimulation (p < 0.001) ( Figure 8B ). The colon tissue of mice was collected for H&E staining. As illustrated in Figure 8C , DSS induction resulted in the desquamation and necrosis of colonic epithelial cells, infiltration of inflammatory cells within the mucosal layer, and loss of crypt structures in DSS-induced colitis model mice. The levels of IL-6 and TNF-α in colon tissue were measured utilizing ELISA. The results revealed a significant increase in the levels of inflammatory cytokines IL-6 and TNF-α in the colon tissue of DSS-induced colitis model mice (p < 0.001) compared with the control group ( Figures 8D, E ). To further elucidate the expression of key genes in UC, we conducted Western blot to evaluate corresponding protein expression levels in colon samples obtained from a DSS-induced colitis mouse model. Results indicated that, compared to the control group, GNG5 protein levels were significantly downregulated and TIMP1 levels were significantly upregulated in the DSS group (p < 0.001, p < 0.01), which confirmed our spatial transcriptome predictions ( Figure 8G ).

Our single-cell sequencing data showed high monocyte expression in UC samples, with macrophages being crucial monocyte components. Lots of references indicates that macrophages are vital in UC inflammation and tissue repair (24). Immune infiltration analysis revealed a stronger correlation between TIMP1 and UC-related immune cells compared to GNG5. Thus, we used immunofluorescence co-localization to assess TIMP1 expression in macrophages. The results indicated that TIMP1 co-localizes with macrophage marker F4/80 in colon tissue, suggesting that TIMP1 may affect disease progression through functional expression in UC colon macrophages ( Figure 8G ).

To quantitatively assess the activity of immune metabolism genes in individual cells, we utilized AUCell. Bubble plots were employed to visualize the differential activity of key genes within these pathways. Our findings revealed that GNG5 and TIMP1 were significantly upregulated in oxidative phosphorylation, the unfolded protein response, and related pathways ( Figure 9A ). Furthermore, an analysis of classical exhaustion-related genes (LAG3, PDCD1, TIGIT, HAVCR2, CTLA4) in single cells indicated a pronounced T cell exhaustion phenotype ( Figures 9B, C ). We also investigated the correlation between exhaustion-related genes and immune infiltration as well as their differential expression in the transcriptome. CTLA4, LAG3 and TIGIT were found to be significantly up-regulated in UC ( Figure 9D ).Our analysis identified a significant positive correlation between five exhaustion-related genes and activated memory CD4⁺ T cells, T follicular helper cells, as well as other cell types. ( Figure 9E ). To delve deeper into the relationship between GNG5, TIMP1, and cellular depletion, a correlation analysis was conducted involving five depletion-related genes. This analysis identified a significant positive correlation between TIMP1 and TIGIT as well as CTLA4 ( Figure 9F ) (p = 3.5e-06, 1.1e-07).