The UC-related flow cytometry standard (FCS) files were preprocessed and then downloaded from Cytobank (https://www.cytobank.org/) (Chen and Kotecha 2014). The channel of each FCS file was scaled according to the minimum and maximum proportion to avoid the data deviating from the recognizable area. Then, the cell population was identified by manual gating. The positive or negative cell subsets of two visual markers were selected and further stratified until the target population in a series of marker combinations was captured. The gating subsets were captured in Cytobank, including T cells (CD45⁺CD3⁺ cells), Tregs (CD45⁺CD8a⁻CD4⁺CD25⁺CD127^- cells), B cells (CD45⁺CD3⁻CD19⁺ cells), innate immune cells (CD45⁺CD3⁻CD19⁻ cells) and NK cells (CD3⁻CD45⁺CD56⁺CD16⁻ cells). To ensure channel consistency, the CyTOF FCS files were uploaded to FlowJo software (FlowJo, Ashland, OR).

All statistical calculations were conducted in R (version 4.0). The R packages employed in the analysis process included cytoWorkflow (Nowicka et al., 2017), CATALYST, flowCore (Hahne et al., 2009). FlowSOM (Van Gassen et al., 2015) was used for clustering, while t-stochastic neighbor embedding (t-SNE) (Kobak and Berens 2019; Zhou and Jin 2020) was used to reduce the dimensionality. Heatmaps were used to visualize differences in clusters, and boxplots were used to visualize differences in selected immune cell subsets (nodes). The formation of heatmaps and boxplots was accomplished using ggplot2 (https://cran.r-project.org/web/packages/gplots/index.html). Kruskal-Wallis test was used for comparison among multiple groups, and Wilconxon rank sum test was used for comparison between two groups, and p ≤ 0.05 was considered a statistically significant difference.

The detailed methods of visualization and calculation were as follows. 1) To visualize the differences in clusters, the proportion of 39 nodes was determined by the abundance of cluster cells. Analysis of variance of cell population abundance was used to compare the proportions of cell types to highlight populations with different proportions. 2) To visualize the differences in the selected nodes, the p-value of selected cell subsets was calculated, and the significant differences among them in individuals were compared using box line diagrams. 3) t-SNE dimension reduction of selected nodes and their important cytokines was used for revalidation.

To verify the CyTOF results, two datasets were collected from the GEO database (GSE125527 and GSE116222) (Boland et al., 2020; Mitsialis et al., 2020). Moreover, we also collected samples (colon tissues) for single-cell RNA sequencing to further verify the CyTOF results.

The Ethics Committee of The First Affiliated Hospital of Guangzhou University of Chinese Medicine approved the study (No: ZYYECK [2019]160). All the participants signed the informed consent form. The inclusion criteria were as follows: 1) participants were diagnosed with Uca or UCin or were healthy control subjects based on “The Asia-Pacific consensus on ulcerative colitis” published by APAGE on IBD in 2010 (Ooi et al., 2010); and 2) participants were aged between 18 and 75 years old. Then, all collected tissue samples will be sent to Huada Gene Technology Co., Ltd. (Shenzhen, China) for processing.

Colon tissues were placed in RPMI 1640 medium supplemented with 100 μg ml⁻¹ streptomycin, 100 U ml⁻¹ penicillin and 10 mm HEPES (Gibco, Gaithersburg, MD) on ice. Biopsies were thawed, washed twice in PBS and digested with collagenase type II (Worthington Biochem) under stirring (250 rpm) at 37°C for 20 min. The cell suspension was washed in PBS and passed through a 70-μM filter. Cells were counted with an automated cell counter, and their viability was confirmed with Trypan blue. Single cells were resuspended in freezing buffer (90% FBS and 10% DMSO freezing media), placed into a freezing container (Thermo Scientific Nalgene Mr Frosty) at −80°C, and then transferred to liquid nitrogen for long-term storage.

Cellular suspensions were encapsulated in droplets using GemCode technology and processed by Huada Gene Technology Co., Ltd. (Shenzhen, China). In brief, gel beads with barcodes, UMIs (unique molecular identifiers), primers and enzymes were mixed with single cells by using a GemCode Single-Cell Instrument (10X Genomics, Pleasanton, CA, United States). In droplets, the cells broke, and the released mRNA was linked to the cell tag sequence on the gel beads to form a single-cell GEM structure (gel bead in emulsions). cDNA for sequencing on Illumina platforms was generated by reverse transcription using a GemCode Single-Cell 3’ Library Kit. Cell mRNA was reverse transcribed in droplets to form cDNA and then demulsified to construct indexed sequencing libraries. Sequencing libraries were then subjected to scRNA-seq using the Illumina platform.

The Seurat package (3.2.2) was employed to merge and cluster datasets and to analyze differential gene expression (Satija et al., 2015). The entire work flow was as follows:

Quality control (QC) and merging Seurat objects: The QC process was performed using Seurat (version 3.2.2) (Butler et al., 2018). Single cells with the following characteristics were considered as low quality cells. The single cells are removed if they met the following criteria: 1) single cells containing less than 400 genes; 2) single cells with more than 4,000 genes; 3) single cells with more than 50% mitochondria-specific genes. For the GSE116222 database, the IntegrateData function was used to eliminate batch effects among the patients.

Dimensionality reduction: Scaled and log-normalized data were subsequently used to enable dimensionality reduction, and the main cell clusters were identified using the FindClusters function and visualized using t-SNE.

Cluster cell type annotation: To identify the markers of each cluster, the FindAllMarkers function was used. The main cell types were then recognized based on the markers obtained from the CellMarker database (Zhang et al., 2019). Eventually, T cells, B cells and NK cells were selected for further verification.

Verification of differential expression of specific factors: The variation in the expression of specific subsets in different clusters obtained by CyTOF analysis was visualized using a box graph. It was verified that these specific genes were veritably expressed in specific cell clusters. Differentially expressed gene-GO/KEGG analysis: To further investigate UC disease-associated cell functional states and potential molecular regulators, GO/KEGG analysis of differentially expressed genes (DEGs) in UC mucosa was conducted. Functional enrichment analysis was performed using the “ClusterProfiler” package in R software with a significant level of q-value < 0.05, for GO biological processes and KEGG pathways.

In addition, we also performed a differential abundance test using miloR (Dann et al., 2022), which will test for significant differences in cell abundance between the two different conditions, and the associated gene signatures.

To discover the differences among T cell populations across the three groups, CD3⁺CD45⁺ T cells were clustered in a dedicated analysis (Figure 1A,B). Multiple T cell subsets mainly contained effector memory [EM] (CD45RA⁻CD45RO⁺CCR7⁻CD27^+/-), central memory [CM] (CD45RA⁻CD45RO⁺CCR7⁺CD27⁺), and naïve (CD45RA⁺CD45RO⁻CCR7⁺CD27⁺) cells.

A population of EM T cells (node 11, containing clusters 35, 36 and 40) was found to express CCR6, TNF, CD161, IFNG and IL-17A. There was a statistically significant difference between the abundance levels of these cells among the UCa, UCin and HC mucosa groups (p = 0.011). These EM T cells were enriched in both UCa and UCin mucosa compared to HC mucosa (p < 0.05, Figure 1C). However, none of these differences were statistically significant among the three groups in peripheral blood.

Tregs were manually identified based on the expression of CD45, CD4, CD25, and CD127 and were then used for dedicated automated analysis (Figure 2A,B). Treg subsets were allocated on the basis of the same guidelines as for T cell subsets.

As shown in the figure, the abundance level of EM Tregs across node 15 (clusters 10, 11 and 17) and node 27 (clusters 25, 20, 24, 21 and 26) was statistically significant, with p = 0.00079 and p = 0.042, respectively. Node 15, which coexpressed CD38 and TNF, was comparatively enriched in Uca mucosa but decreased in HC and UCin mucosa (p < 0.05, Figure 2C). Node 27, which coexpressed CXCR3, CD27, CTLA4 and IL1B, was specifically increased in UCin peripheral blood compared with UCa and HC peripheral blood (p < 0.05, Figure 2D).

CD3⁻CD45⁺CD19⁺ B cells were clustered in an analysis dedicated to exploring the differences in B cell subset populations among the UCa, UCin and HC groups (Figure 3A,B). Multiple B cell subsets were categorized, including CD38⁺CD27⁺ plasmablasts, CD27⁻ naïve cells, and CD27⁺ memory cells.

The difference in the number of plasmablasts (node 6 in Figure 3B, containing clusters 37 and 36) among the three groups was statistically significant (p = 0.024). A decreasing trend of these plasmablasts in UCa mucosa samples was discovered (Figure 3C). In addition, for node 15 (clusters 8 and 9) and node 22 (clusters 28 and 32), the differences in naïve B cells were statistically significant (p = 0.039 and p = 0.036, respectively). CXCR3⁺CCR4⁺ naïve B cell clusters (node 15), expressing HLA-DR, TNF and CCR7, were found to be expanded more in UCa mucosa than in HC mucosa (p < 0.05, Figure 3D). Node 22 (expressing TBET, CD38, CD161 and IFNG) was decreased in UCin peripheral blood compared with HC and UCa peripheral blood (p < 0.05, Figure 3E).

Innate immune cells (CD3⁻CD45⁺CD19⁻ cells) were manually gated with CyTOF (Figure 4A). Several innate immune cell subsets were identified according to their high expression factors, mainly macrophages/monocytes (CD11C⁺CD14⁺), dendritic cells (DCs; HLA-DR⁺CD11C⁺CD14⁻), and innate lymphocytes (ILCs; HLA-DR^-IFNG⁺CD161⁺CD127^+/−) (Figure 4B).

There were statistically significant differences in the abundance levels of HLA-DR⁺IL21⁺ macrophages/monocytes among the three groups in both the mucosa (p = 0.027) and peripheral blood samples (p = 0.043). HLA-DR⁺IL21⁺ macrophages/monocytes (node 5, covering clusters 5, 6, 7, 9, 11, and 14) were enhanced in UCa mucosa (p < 0.05, Figure 4C). Nevertheless, in peripheral blood samples, these cells were enriched the most in the UCin group in contrast to the UCa and HC groups (p < 0.05, Figure 4C).

For node 14 (a cluster of ILCs), there was a distinction in the cell frequency among the three groups in the mucosa sample (p = 0.032). Node 14, which expressed AHR, IL1B, CXCR3 and IL23A, was specifically reduced in UCin mucosa (p < 0.05, Figure 5A). Although there was a higher abundance of AHR⁺IL1B⁺ ILCs (node 14) in UCin peripheral blood, the difference was not statistically significant (p = 0.12, Figure 5A). There was a decrease in CD45⁺IFNG⁺ ILCs (node 20) in UC peripheral blood compared to the levels in HC peripheral blood (p < 0.05, Figure 5B), but none of these differences were statistically significant in mucosa.

For node 25 (clusters 13, 15, 23, and 31), the differences in the number of HLA-DR⁺CCR7⁺ DCs among the three groups were statistically significant (p = 0.0032). In detail, these cells were increased more in UCa mucosa than in UCin and HC mucosa (p < 0.05, Figure 5C).

Natural killer (NK) cells (CD45⁺CD3⁻CD56⁺CD16⁻ cells) were identified by manual gating (Figure 6A). Several NK cell subsets were identified according to their high expression factors (Figure 6B), mainly including regulatory NK cells (rNK, CD11⁺CD27⁺), tolerant NK cells (tNK, CD11⁻CD27⁻) and cytotoxic NK cells (cNK, CD11⁺CD27⁻) (Fu et al., 2014).

1) The differences in the level of rNK cells (node 3, containing clusters 32, 37 and 38) were statistically significant (p = 0.00099). Specifically, rNK cells expressing CD38, CXCR3, CTLA4, IL22 and TNF were diminished in UCa mucosa compared to HC and UCin mucosa (p < 0.05, Figure 6C). 2) For node 17 (clusters 36 and 39), CXCR3⁺CCR4⁺ tNK cells among three groups were significantly different (P = 0.0017), and there was a lower level in UCa mucosa than in HC and UCin mucosa (p < 0.05, Figure 7E). (3) For AHR⁺CD14⁺ cNK cells (node 11, covering clusters 13, 15 and 18), there was a significant difference among the three groups (p = 0.00095). AHR⁺CD14⁺ cNK clusters were prominently increased in UCa mucosa (p < 0.05, Figure 7D). In addition, there were statistically significant differences in the abundance level of node 18 among the three groups in either mucosa (p = 0.018) or peripheral blood samples (p = 0.0041). IFNG⁺CD8A⁺CTLA4⁺ cNK cells (node 18, including clusters 1, 2, 9, 12, and 14) were significantly expanded in UCa mucosa and conversely declined in UCa peripheral blood compared with levels in HC samples (p < 0.05, Figure 7F).

Finally, intestinal mucosal biopsy samples, including 2 UCa patients, 2 UCin patients and 2 HC patients, were obtained (our dataset). scRNA-seq on mucosal samples from 6 subjects was performed to validate the previous findings. We identified 5 immune cell lineages (including NK cells, T cells, B cells, DCs, and mast cells) and 4 nonimmune cell lineages (including epithelial cells, endothelial cells, neuroepithelial cells, smooth muscle cells and tissue muscle cells, Supplementary Figure S1).

GSE125527 and GSE116222 were also used to validate the results of CyTOF. A total of 30 samples, including 16 HC samples and 14 UC samples, were included in GSE125527. For GSE125527, we identified 5 immune cell lineages, including NK cells, T cells, B cells, DCs, and macrophages (Supplementary Figure S2). A total of 9 samples, including 3 HC samples, 3 UCa samples and 3 UCin samples, were included in GSE116222. For GSE11622, we identified 5 immune cell lineages (including NK cells, T cells, B cells, DCs, and mast cells) and 2 nonimmune cell lineages (including epithelial cells and neuroepithelial cells, Supplementary Figure S3).

The results of GSE125527 and our dataset showed that CCR6⁺TNF⁺CD161⁺ EM T cells (node 11) were enriched in UCa mucosa, and the abundance levels between the UCa and HC groups were significantly different (p < 0.05, Figures 8A, 9A). In addition, node 11 had a comparatively increased trend in UCa mucosa compared to UCin and HC mucosa (Figure 8E). For CD38⁺TNF⁺ Tregs (node 15), the abundance in UCa was significantly higher than that in HC (p < 0.05, Figure 8B). Similarly, there was a trend in which node 15 was highest abundant in UCa mucosa and lowest abundant in HC mucosa (Figure 8F).

Single-cell analysis of three datasets verified that IFNG⁺CD8A⁺CTLA4⁺ cNK cells (node 18) were obviously increased in UCa mucosa and diminished in UCa peripheral blood (p < 0.05, Figures 8C,D, 9F). Moreover, an enrichment in cNK cells (node 18) in UCa in contrast with both HC and UCin mucosa was confirmed (Figure 8G). Increased abundance of CXCR3⁺CD38⁺ rNK cells (node 3) in UCin was also observed, and the difference between UCin and the other two groups was statistically significant (p < 0.05, Figure 9E).

The results for B cells demonstrated that the abundance of node 6 (CD38⁺CD27⁺ plasmablasts) was the highest in HC mucosa (Figures 8H, 9D), and the abundance levels between HC and UCin mucosa were significantly different (p < 0.05). In addition, scRNA-seq found that there was significantly higher abundance of node 15 (CXCR3⁺CCR4⁺ naïve B cells) in UCa mucosa than in HC and UCin mucosa (p < 0.05, Figure 9C).

HLA-DR⁺CCR7⁺ DCs (innate node 25) measured by scRNA-seq among the three groups showed the same trend as the results with CyTOF, and the abundance of DCs in UCa mucosa was greater than that of the HC and UCin groups (p < 0.05, Figure 9B).

CCR6⁺TNF⁺CD161⁺ T cells (T node 11), HLA-DR⁺CCR7⁺ DCs (innate node 25), CXCR3⁺CCR4⁺ B cells (B node 15), and CD8A⁺IFNG⁺ NK cells (NK node 18), which were increased in the UCa group and decreased in the HC group, were selected for gene-GO/KEGG analysis (Figures 9G,H).

GO analysis indicated that the primary functions of these four subsets were similar. It revealed the activation and differentiation of many immune cells in UC, including T cell activation, B cell activation, regulation of T cell activation, and lymphocyte differentiation. These subsets also revealed a higher proportion of biological processes, consisting of positive regulation, response to stress, cellular response, and cell adhesion.

KEGG pathway analysis also demonstrated disease enrichment in each subset, including coronavirus disease (COVID-19) and rheumatoid arthritis.

In Supplementary Figure S4, at a 10% FDR, we identified 6635 DA neighbourhoods. As shown in the figure (Supplementary Figure S4C), immune cells were differentially expressed in UC and HC conditions, which is consistent with the previous results.