Integrated Bioinformatics Analysis and Experimental Verification of Immune Cell Infiltration and the Related Core Genes in Ulcerative Colitis

Ulcerative colitis (UC) is an autoimmune disease characterized by recurrent diarrhea, abdominal pain, and mucous, pus, and blood in the stool. Although various drugs are available to relieve UC, up to 20% of the patients require colectomy.¹ Therefore, it is important to find new core genes related to UC and further investigate the pathogenesis of the disease as a clear understanding is still lacking.

Research has shown that UC is closely associated with the activation of the immune system, especially the Th2 immune response. Furthermore, immune cell infiltration contributes to the occurrence and progression of UC. For instance, macrophages are polarized to the M1 proinflammatory phenotype,² hundreds of neutrophils accumulate and damage the epithelium of the intestinal mucosa.³ In patients with IBD, intestinal inflammation affects the development of tolerogenic macrophages and results in the accumulation of proinflammatory monocytes.⁴

Various immune pathways and inflammatory genes are known to be involved in the pathogenesis of UC. It is unclear whether other novel genes play a role in the pathogenesis of UC apart from the known genes and whether the new genes are associated with immune cells. Recently, the advancements in bioinformatics, especially the continuous updates in the Gene Expression Omnibus (GEO) database, R software, machine learning, CIBERSORT algorithm, etc, have provided a new means for the study of UC related novel genes and immune cell infiltration. Hence, in this study, the novel genes of UC, infiltration of UC immune cells, and the relationship between the new genes and immune cells were analyzed using bioinformatics tools, which are of immense importance for elucidating the molecular mechanism of UC and identifying new targets for immunotherapy.

In this study, the UC microarray datasets were obtained from the GEO website. After screening the differentially expressed genes (DEGs), the UC core genes were identified by combining machine learning methods. Subsequently, CIBERSORT was used to study the differences in immune cells in 22 immune cell subsets between the UC and the normal group. The association between core genes and immune cells was also studied to understand the molecular immune mechanism involved in the pathogenesis of the disease. Finally, the infiltration of immune cells and the expression of core genes in UC were verified with in vivo experiments.

Materials and Methods

Screening Gene Chip Datasets

GEO is an international database that includes microarray gene expression profiles and high-throughput sequencing data. The term “ulcerative colitis” was used as the search term, while “Series” “Expression profiling by array “Homo sapiens” were used as a qualifier. The screening standards included the following: the microarray dataset involved whole-gene colonic gene expression profiles; the microarray dataset contained UC and healthy samples; drugs were not administered to any of the samples included, and the total number of samples was >20. After data screening, GSE87466, GSE179285, and GSE75214 were selected. Our analysis was limited to samples that were active and healthy and met the inclusion criteria. It was not necessary to obtain ethics approval for the present study, as it was a secondary bioinformatics analysis that was based on a public database (GEO database). All the data of this study are from public databases, an ethical review was exempted from the ethics committee by the Ethics Committee of The First Affiliated Hospital of Zhejiang Chinese Medical University (2023-YBK-001-01).

Differentially Expressed Genes (DEGs) Analysis

The study utilized two different datasets from varying platforms and personnel, resulting in heterogeneity and latent variables. To address this, we batch normalized the combined data sets using the “sva” package in R software. We then analyzed the genetic differences between normal people and UC patients using “Limma”, and identified DEGs using P<0.05 and |log FC|≥1 as screening conditions. Then, the “ggplot2” “ggrepel” “dplyr” and “pheatmap” packages were used to visualize the volcano plot and the heat map of DEGs.

Functional Enrichment of DEGs

DEGs’ biological functions were identified using GO and KEGG pathway enrichment analyses, and the enrichment analysis was performed using clusterProfiler, org.Hs.eg.db, enrichplot, and ggplot2 packages with a cut-off criterion of P < 0.05. DEG-disease associations were identified by the DO analysis.⁵ The Gene Set Enrichment Analysis (GSEA) program was applied to identify significantly enriched pathways and functional hallmarks, and the top 5 pathways were selected for graphical plotting.

Key Gene Screening and ROC Analysis

Several articles have utilized a method that combines the Least Absolute Shrinkage and Selection Algorithm (LASSO) as well as Support Vector Machine Recursive Feature Elimination (SVM-RFE) algorithm to screen key genes.^6–9The key gene screened based on GSE87466 and GSE179285 were validated in GSE75214. In this study, LASSO logistic regression analysis was conducted using the R package “glmne”, with the smallest lambda value considered as optimal. Feature gene selection was performed using the RFE function in the caret package, based on five-fold cross-validation. SVM classification was then constructed using the “e1071”, “kernlab”, and “caret” packages in R. The ROC curves were generated using the R package pROC, and the area under the curve (AUC) was calculated to evaluate the diagnostic accuracy. AUC > 0.9 indicated that the model fit well.

Immune Infiltration Analysis

On the CIBERSORT website, the expression matrix of individual immune cell subtypes can be deconvolved using the linear support vector regression method, which provides a set of gene expression features of 22 immune cell subtypes. The immune cell infiltration matrix was obtained with a cut-off criterion of P < 0.05. Analysis of the correlation between 22 infiltrating immune cells was performed using the “corrplot” package. The results were visualized using “limma”, “reshape2”, “ggpubr” ggExtra, and “vioplot” packages. The characteristic diagnostic marker was compared with the immune-infiltrating cells using these analytical tools.

Animal Experiments

An approval was obtained by the Animal Experimentation Ethics Committee at Zhejiang Chinese Medical University for this study (IACUC-20210517-24). All animal experiments were conducted in accordance with the guidelines of Zhejiang Chinese Medical University. In this study, dextran sulfate sodium (DSS) was used to construct a UC mouse model and explore the expression of immune cells and core genes in UC. Meanwhile, due to the immune regulation of Azathioprine (AZA), the correlation between immune cells and the core genes as well as the pathogenesis of UC was further analyzed by observing the changes in immune cells and core genes.

Male wild-type mice (C57BL/6; 9–12 week-old) were purchased (Slaughter Shanghai). Three groups of mice were used in the experiment: control/WT, DSS/WT, and AZA/WT (n = 6/group). Then, 3.0% DSS was administered in drinking water for 7 days, and saline or 0.01 mL/g AZA was applied via gavage for 10 days, starting from the first day of modeling. The mice were then sacrificed after the experiments, and the colonic tissues were retained for subsequent analyses.

Quantitative Real-Time PCR Assay

According to the manufacturer’s instructions, total RNA was extracted from samples with the Trizol reagent (TAKARA, Japan). Reverse transcription was performed by using the Reverse Transcription Kit (ThermoFisher Scientific, USA). Quantitative PCR was performed with the PowerUp SYBR Green Reagent (ThermoFisher Scientific) on the Roche LightCycler^® 480. The primer sequences used in the study are shown in Table S1. The 2^−ΔΔCt method was used and the relative levels of gene expression (mRNA) were calculated.

Immunohistochemical Analyses

The distal colon tissues were fixed with 4% paraformaldehyde before embedding in paraffin, sectioned at a thickness of 4 um, dewaxed, and subjected to antigen repair with citric acid (VECTOR Laboratories, USA) at a high temperature, and then blocked with 5% bovine serum albumin (Absin, China) for 30 min. The dilution solution of the primary antibody was incubated at 4°C overnight, and the dilution ratios of the anti-Ly6C antibody, anti-F4/80 antibody, and anti-Ly6G antibody were 1:300. The antibodies were purchased from HUABIO (Zhejiang, China). Immunohistochemical secondary antibodies were incubated with enhanced HRP-labeled sheep anti-mouse/rabbit IgG polymer (ZSBB-Bio, CHINA) for 30 min at room temperature, followed by treatment for 10 min with DAB, counterstaining with hematoxylin, and finally sealing with neutral resin. ImageJ analysis was used to calculate the positive results of the images captured with the VS120 microscope (Olympus, Tokyo).

Results

Screening of DEGs in Different Datasets

The sample distribution of the three genes in the microarray datasets is shown in Table S2. A total of 36 DEGs were identified after considering the intersection of GSE87466 and GSE179285, including 20 upregulated genes and 16 downregulated genes (Appendix 1), and the DEGs were plotted using volcano and heat maps (Figure 1).

Figure 1 (A) Volcano map of DEGs in GSE87466 and GSE179285, red dots represent the upregulated genes and the green dots represent the downregulated genes; (B) Heat map of DEGs.

Functional Correlation Analysis

In the present study, 36 DEGs were predicted to serve specific functions. GO enrichment analysis showed that BP was mainly enriched in antibacterial humoral response and humoral immunity, MF in secretory granules and vesicles, and CC in chemokine receptors, transmembrane transport, and other molecular functions (Figure 2A). The pathways involved in the KEGG enrichment pathway were chiefly the IL-17 signaling pathway, bile secretion, NOD-like receptor signaling pathway, and PPAR signaling pathway (Figure 2B). Moreover, the enrichment analysis of DO revealed a significant connection between intestinal disease and inflammatory bowel disease (Figure 2C). The GSEA results indicated that the UC active group was predominantly enriched in the pathways of chemokine signaling, cytokine–cytokine receptor interaction, and toll-like receptor signaling (Figure 2D).

Figure 2 (A) Results of GO enrichment analysis; (B) results of KEGG enrichment analysis; (C) results of DO analysis; (D) GSEA analysis results displaying the top 5 signaling pathways.

Screening Core Genes and Data Set Validation

In this study, 13 key genes were identified from the DEGs using the LASSO logistic regression algorithm (Figure 3A). Eight genes were recognized as the key genes from the DEGs using the SVM-RFE algorithm (Figure 3B). AQP8, VNN1, MME, HMGCS2, and LOC389023 were the overlapping genes in both algorithms (Figure 3C). ROC analysis was conducted and the AUCs for AQP8 (Figures 4A), HMGCS2 (Figures 4B), VNN1 (Figures 4C), and MME (Figures 4D) were found to be 0.963, 0.955, 0.963, and 0.399, respectively. However, LOC389023 was not present in the dataset GSE75214. In the validation set, AQP8 (Figures 5A), HMGCS2 (Figures 5B), and VNN1 (Figures 5C) showed significant differences, while MME did not (Figures 5D). These findings suggest that only AQP8, VNN1, and HMGCS2 may be considered as the core genes of UC.

Figure 3 (A) LASSO logistic regression algorithm used to analyze the feature genes; (B) SVM-RFE algorithm used to analyze the feature genes; (C) Venn diagram of the intersection of two sets of feature genes.

Figure 4 ROC analysis of AQP8 (A), HMGCS2 (B), VNN1 (C), and MME (D) in GSE75214.

Figure 5 The validation results of AQP8 (A), HMGCS2 (B), VNN1 (C), and MME (D) in GSE75214.

Infiltration of Immune Cells

The results of GSE87466 and GSE179285 (44 normal controls and 110 patients with UC) were summed up using the CIBERSORT algorithm. The proportions of the immune cell composition showed clustering and individual differences (Figure 6A). Based on the gene expression matrix comprising 22 immune cells, their abundance in GSE87466 and GSE179285 in the relevant samples was calculated using the R software. See Appendix 2A–D for original files and codes. The strongest positive correlation was observed between neutrophils and monocytes, with a correlation coefficient of 0.52, and the strongest negative correlation was noted between mast cells and CD8 T cells with a correlation coefficient of 0.64 (Figure 6B). In comparison with normal samples, UC samples contained a higher proportion of memory B cells and M1 macrophages (P < 0.05) (Figure 6C).

Figure 6 (A) Results of immune cell infiltration in the colon of UC patients and normal subjects; (B) heat map of the correlation of 22 immune cell species in UC; (C) the violin plot of the ratio of 22 immune cell species.

Correlation Analysis of Core Genes and Infiltration-Related Immune Cells

According to the results of correlation analysis, AQP8 (Figure 7A) was positively correlated with M2 macrophages (R=0.52, p=7.7E-09), eosinophils (R=0.33, p= 0.00045), and resting mast cells (R=0.29, p=0.0017). On the contrary, the gene was negatively correlated with M0 macrophages (R=−0.51, p=8.8e−09), neutrophils (R=−0.48, p= 7e−08), and M1 macrophages (R=−0.34, p=3e−04). HMGCS2 (Figure 7B) was positively correlated with M2 macrophages (R=0.62, p=2.2e−16), eosinophils (R=0.36, p=0.00011), and resting mast cells (R=0.33, p=4e−04). The gene was negatively correlated with M0 macrophages (R=−0.65, p=5.6e−15), neutrophils (R=−0.44, p=1.5e−06), and M1 macrophages (R=−0.38, p=4e−05). VNN1 (Figure 7C) was positively correlated with neutrophils (R=0.27, p=0.0035), plasma cells (R=0.26, p=0.0063), and M0 macrophages (R=0.21, p=0.026) and negatively correlated with CD8⁺ T cells, B cells, and M2 macrophages, but without statistical difference.

Figure 7 Correlation among AQP8 (A), HMGCS2 (B), VNN1 (C), and infiltrating immune cells.

Verification of Core Genes

To further verify the relationship between AQP8, HMGCS2, VNN1 and UC, the mouse model of UC was constructed with the DSS-free drinking water method and treated with AZA gavage. The results showed that the expressions of AQP8 (Figures 8A) and HMGCS2 (Figures 8B) mRNA decreased significantly in the UC mouse model group and increased after AZA treatment. However, VNN1 (Figures 8C) and TNF-α (Figures 8D) mRNA increased in the UC mouse model group and decreased after AZA treatment. Furthermore, the infiltration of macrophages, monocytes, and neutrophils improved (Figure 9).

Figure 8 Expression of AQP8 (A), HMGCS2 (B), and VNN1 (C) and TNF-α (D) mRNA in UC mice. *P<0.05, **P<0.01, ***P<0.001.

Figure 9 Expression of Ly6G, Ly6C and F4/80 in mouse colon. Scale bar: 50μm. Each data represents at least three independent experiments. *P<0.05, **P<0.01, ***P<0.001.

Discussion

UC is a disease of the colon that causes inflammation, and its exact cause is not yet fully understood. The inflammation affects the mucosa and can lead to chronic intestinal inflammation, colectomy, and even cancer. Research suggests that the infiltration of immune inflammatory cells plays a key role in the development and progression of UC. Identifying novel core genes associated with UC and investigating their correlation with immune cells is crucial for understanding the underlying immune molecular mechanisms of UC. The effects of TWP and AZA interventions can shed light on the differential regulation of core genes and immune cell infiltration.

For the integration of two groups of microarray datasets, this study selected 36 DEGs including 20 upregulated and 16 downregulated genes. The results revealed that, after enrichment, the DEGs were closely related to immune responses, inflammation, and signaling pathways. The presence of a strong immune activation process in UC was observed, and the participation of several immune cells eventually caused intestinal inflammation and damage to the intestinal mucosal barrier. After machine learning screening and verification, AQP8, HMGCS2 and VNN1 were found to be the likely core genes of UC, and they exhibited different degrees of correlation with neutrophils, macrophages, and monocytes.

Aquaporins (AQP) belong to a family of integrated membrane proteins that reside in the apical membrane of colonic epithelial cells and mediate passive transmembrane transport of water.^10–12 Hence, AQP8 plays an important role in the absorption of intestinal water. When the expression of AQP8 in the UC intestine is downregulated,¹³ its distribution is altered,^14,15 and the polarity of the epithelium is destroyed, which ultimately leads to intestinal barrier damage and water absorption impairment.¹⁶ However, when AQP8 expression is upregulated, intestinal barrier repair is promoted.¹⁷ AQP8 is positively correlated with M2 macrophages, eosinophils, and resting mast cells and negatively correlated with M0 macrophages, neutrophils, and M1 macrophages. A study has noted that AQPs are correlated with macrophage infiltration,¹⁸ but their relationship with other immune cells is less studied and needs to be further explored.

Mitochondrial 3-hydroxy-3-methylglutaryl CoA synthase 2 (HMGCS2) belongs to the HMG-CoA synthase family, and this protein is a key mitochondrial enzyme that is involved in ketogenesis.¹⁹ HMGCS2 is expressed in the most differentiated areas of the intestinal mucosa and promotes the differentiation of intestinal cells into goblet cells, Paneth cell, and epithelial cells. The decrease in the level of HMGCS2 in the intestine induces the development of UC, thereby resulting in an imbalance in the differentiation and proliferation of enterocytes.^20–23 More importantly, HMGCS2-mediated fatty acid β-oxidation-ketogenesis acts as a protective mechanism, inhibiting M1 macrophage activation and attenuating TNF-α-induced intestinal cell apoptosis, CXCL3 expression, and ROS production.^24,25 A negative correlation was perceived between HMGCS2 and M0 macrophages, neutrophils, and M1 macrophages in this study. M2 macrophages, eosinophils, and resting mast cells were positively correlated with HMGCS2 but M0 macrophages were negatively correlated with HMGCS2, which further confirms that HMGCS2 has a protective effect on UC. Therefore, HMGCS2 is expected to become one of the potential targets for the treatment of UC.

Vascular noninflammatory molecule 1 (VNN1) is a glycosylphosphatidylinositol anchored pan-enzyme, which is expressed by intestinal epithelial cells and induces macrophage polarization. When catalyzed to cysteamine, VNN1 in turn promotes chemokine secretion by neutrophils, which plays pertinent role in oxidative stress, inflammation, and cell migration.^26–29 VNN1 is highly expressed in UC colon, and three single nucleotide polymorphisms in VNN1 were found to be associated with genetic susceptibility to IBD.³⁰ VNN1-KO mice can resist oxidative stress-induced intestinal injury,^31,32 with significantly decreases the expression of inflammatory factors and reduces intestinal inflammation, which may be related to the inhibition of the PPAR-γ pathway.³³ This study showed that VNN1 was positively correlated with neutrophils, plasma cells, and M0 macrophages but not negatively correlated with other immune cells. Therefore, the proinflammatory effect of VNN1 may be one of the potential mechanisms in the pathogenesis of UC. Hence, treatment targeting VNN1 may alleviate the intestinal immune inflammatory response in UC.

Currently, there is no conclusive evidence indicating a direct functional relationship between AQP8, HMGCS2, and VNN1 in relation to ulcerative colitis (UC). However, each of these proteins represents distinct functions within the context of UC. AQP8 is involved in intestinal water absorption, HMGCS2 serves as a protective factor against UC, and VNN1 acts as a pro-inflammatory factor. Notably, KEGG analysis has identified the IL-17 signaling pathway as being significantly enriched in UC. IL-17 is a crucial effector factor produced by Th17 cells that can stimulate the release of pro-inflammatory factors and neutrophil chemokines, ultimately regulating tissue inflammatory response. Its role in the immune regulation network and its association with inflammation and autoimmune diseases is well-established. Therefore, it is plausible to assume that a correlation exists between the three genes and the IL-17 signaling pathway, although further confirmation is required.

UC is an autoimmune disease caused by the over activation of the immune system. To better illustrate the impact of immune cell infiltration in UC, immune infiltration was assessed using CIBERSORT. Increased infiltration of neutrophils, activated mast cells, monocytes, and M0/M1 type macrophages and decreased infiltration of M2 type macrophages, resting mast cells, and CD8⁺ T cells were noted in this study. The above immune cells may be associated with the development of UC. In the early stage of UC inflammation, neutrophils are activated and get accumulated in the damaged area, thus forming neutrophil extracellular traps to destroy the intestinal pathogens, relieve the intestinal mucosal inflammation, repair the intestinal mucosa, and maintain the intestinal homeostasis.³ If inflammation persists, excessive accumulation and activation of neutrophils can impede the recovery of the intestinal mucosa. Mast cells (MC) are usually associated with type I hypersensitivity, but MCs also have a key role in UC. MCs induce IL-33 secretion, thereby creating a proinflammatory microenvironment in the intestine by producing histamine, aggregating neutrophils, and exacerbating inflammation.^34,35 Blood monocytes play an important role in innate immune defense by phagocytosing pathogens and digesting, processing, and presenting antigens. Moreover, monocytes are the sole source of intestinal macrophages.³⁶ Macrophages are crucial for innate immunity and microenvironmental homeostasis in the intestine.³⁷ They are polarized into different phenotypes upon stimulation, with M1-type macrophages playing a significant proinflammatory role in UC.

In this study, macrophages, monocytes, and neutrophils were selected for verification in vivo. Ly6G is a marker of neutrophils, Ly6C is a marker of proinflammatory monocytes, and F4/80 is a surface marker of mouse macrophages. Several studies have confirmed that the expressions of the three markers are increased in the UC colon.^38–41 This finding was demonstrated in the present study too, where the expressions of neutrophils, monocytes, and macrophages were increased in the UC intestine. This increase was accompanied by an increase in the inflammatory factor TNF-α. Simultaneously, the expressions of AQP8 and HMGCS2 decreased in the UC model group, whereas that of VNN1 increased. As a commonly used immunosuppressive agent in the treatment of UC, AZA has a good immunosuppressive effect and can effectively reduce the infiltration of various immune cells.^42,43 This study showed that after AZA treatment, the infiltration of neutrophils, monocytes, and macrophages was significantly lowered, TNF-α was decreased, AQP8 and HMGCS2 were increased compared with the model group, and VNN1 was decreased compared with the model group. This finding indicates that immune cells and AQP8, HMGCS2, and VNN1 are involved in the pathways of AZA. This observation further suggests that immune cells, AQP8, HMGCS2, and VNN1 are involved in the pathogenesis of UC.

Intestinal inflammation in UC is commonly associated with immune inflammatory cell infiltration. This study highlights AQP8, HMGCS2, and VNN1 as core genes in UC that are linked to immune cells. TWP has been found to have a significant therapeutic effect on UC and also exhibits immunoregulatory properties. AZA, a commonly used immunosuppressant for UC treatment, has been observed to reduce the infiltration of various immune inflammatory cells (33–35). In vivo studies have shown that TWP and AZA have an impact on immune cells, but they have different mechanisms of immune regulation. TWP targets pro-inflammatory genes like VNN1 and has limited regulatory effects on protective genes such as AQP8 and HMGCS2. On the other hand, AZA affects all three genes and regulates them. The study highlights the differences in immune regulatory mechanisms between TWP and AZA in reducing immune cell infiltration. It also suggests that AQP8, HMGCS2, and VNN1 may be associated with immune cell infiltration, indicating a potential pathogenesis of UC immune inflammation.

Hence, further experimental validation is required. Although lacking direct experimental verification, the significance of the three genes in immune inflammation warrants further investigation. Our future research will focus on exploring their relationship with immune cells and underlying mechanisms.

Conclusion

This study identified AQP8, HMGCS2, and VNN1 as the central genes associated with UC. These genes are linked to the infiltration of neutrophils, monocytes, and macrophages, and potentially play a role in the IL-17 signaling pathway, contributing to the development of UC. Identifying the immunotherapy targets and optimizing immunomodulatory therapy for patients with UC may be achieved with an in-depth exploration of immune cells.