We systematically searched for microarray studies from public GEO datasets using the following terms: “mucosal” and “ulcerative colitis”. The screening criteria were as follows: ① expression profiling by array; ② Homo sapiens; ③ dataset containing more than five active UC samples and healthy controls. Finally, five GEO datasets (GSE16879, GSE75214, GSE48958, GSE87473, and GSE92415) were included (27–30) ( Table 1 ). We downloaded series matrix file(s) from the GEO website and corresponding annotation documents from GEO or Bioconductor.

RRA analysis was used to integrate gene expression data from different datasets in an unbiased manner using a comprehensive ranking list algorithm (31). We combined the “limma” and “RobustRankAggreg” packages in R to obtain statistically changed genes (adj. p-value < 0.05) and differentially expressed genes (DEGs; adj. p-value < 0.05 and |logFC| > 1.5).

Gene Set Enrichment Analysis (GSEA), Gene Ontology (GO) analysis, and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were conducted using the “ClusterProfiler” R package (32). The reference gene list c2.cp.reactome.v7.0.symbols.gmt was downloaded from Molecular Signature Database (MsigDB) (33). GO analysis included the biological process (BP), cellular component (CC), and molecular function (MF) categories. FDR < 0.05 was considered significant.

The CIBERSORT deconvolution algorithm was employed to analyze the cellular composition of the tissues based on gene expression profiles, according to the known reference set LM22 (leukocyte signature matrix) (34). The permutation (perm) was set at 1000 to obtain more stable results.

WGCNA can be used as a data exploratory tool or screening method to identify key gene modules using an unsupervised clustering without priori defined gene sets (35). In our study, a total of 15162 genes (top 75% according to variation) were analyzed using the “WGCNA” R package. Scale-free network features were constructed when the power of β was equal to 12 (R² = 0.85). Dynamic tree cut algorithm was used to aggregate genes with similar expression profiles into the same module. Gene clusters in the module most related to the active UC traits were considered candidates for further validation.

STRING database (https://cn.string-db.org/) is a functional protein association network, assembling all known and predicted proteins (36). Multiple protein names were put into the list box, with one name per line. PPI network interactions file with medium confidence scores ≥ 0.4 was downloaded. Cytoscape software (version 3.7.2), a general-purpose, open-source software platform for network biology analysis and visualization (37). Molecular complex detection (MCODE) (K-core = 2, degree cutoff = 2, max depth = 100, and node score cutoff = 0.2) and cytohubba-MCC plugins with default parameters were used to explore important gene clusters and hub genes (38, 39). The iRegulon (Version: 1.3) plugin was used to screen key TF with the default cutoff criteria (40).

We used the R “corrplot” package to display a correlation matrix. We download the organized gene list of “Immunomodulator” from TISIDB database (http://cis.hku.hk/TISIDB/index.php), which is an online integrated repository portal for tumor-immune system interactions (41). p-value < 0.05 was considered to be statistically significant.

Animal experiments were approved by the Animal Ethics and Experimentation Committee of Soochow University and carried out according to the Guide for the Care and Use of Laboratory Animals. AGP or DMSO diluted in PBS was intraperitoneally injected into mice at a dose of 25 mg/kg on alternate days, with the first dose administered one day before DSS administration. The experimental colitis model was induced in 10-week-old C57BL/6 male mice by administering 3.5% DSS (MW: 36,000–50,000 Da; Yeasen Biotechnology Co., Ltd., Shanghai, China, 60316ES76) for seven days, followed by administration of normal water for three days. Mice were sacrificed on the tenth day, and colon tissues were obtained for hematoxylin and eosin (H&E) staining, immunofluorescence (IF) staining, and RT-PCR analysis.

Colon tissues were fixed with 4% paraformaldehyde (PFA) overnight and were then embedded in paraffin. Colonic sections of 5 μm were obtained and laid flat on a glass slide for H&E or IF staining. For IF staining, the tissue slices were permeabilized with 0.1% Triton X-100 for 15 min, blocked with 1% BSA for 1 h at 37°C, and then incubated with antibodies against NF-κB (Cell Signaling Technology, #8242, 1:500) and IκBα (Cell Signaling Technology, #4812, 1:500) at 4°C overnight. The next day, tissue slices were incubated with an Alexa Fluor^® 488-conjugated goat anti-rabbit IgG antibody (Life Technologies, A11088) for 1 h at room temperature. The nuclei were stained with Hoechst 33342 (1:10000 dilution, Beyotime) for 10 min prior to imaging. Representative cells were selected and photographed.

RNA was extracted from the intestinal tissues using TRIzol reagent (Invitrogen, 15596018) and reverse-transcribed to cDNA using the Hifair^® II 1st Strand cDNA Synthesis Kit (Yeasen Biotechnology, 11119ES60). Primer sequences were obtained from PrimerBank (https://pga.mgh.harvard.edu/primerbank/) and synthesized at GENEWIZ Biotech Co. Ltd. (Suzhou, China) ( Table 2 ). RT-PCR was performed using the Hieff^® qPCR SYBR Green Master Mix (High Rox Plus; Yeasen Biotechnology, 11203ES03) on an ABI Applied Biosystems. PCR amplification was conducted in triplicate for each sample, and the expression of target genes was normalized to GAPDH. Relative expression was determined using the 2^-ΔΔCt method.

TargetScan Human 7.2 (http://www.targetscan.org/vert_72/), miRDB (http://mirdb.org/), microT v5 (http://diana.imis.athena-innovation.gr/DianaTools/index.php?r=microT_CDS/index), and miRWalk 3.0 (http://mirwalk.umm.uni-heidelberg.de/) with default parameters were used to predict the miRNAs (42–45). StarBase (v2.0) (https://starbase.sysu.edu.cn/starbase2/index.php) was used to predict miRNA-lncRNA interactions with very high stringency (≥5) (46). The interaction networks were constructed and visualized using Cytoscape software.

The R “ggpubr” package was used to perform statistical analyses, and the R “ggplot2” package was used to generate images. Differences between two groups were assessed using the Student’s t-test. All data are shown as mean ± standard error of the mean (SEM). GraphPad Prism 7.0 (GraphPad Software, Inc., La Jolla, CA, USA) was used for statistical analysis and image construction. Adobe Illustrator (AI) CS6 software was used to edit the figures.

The study workflow is presented in Figure 1 . GSE16879, GSE48958, GSE75214, and GSE87473 expression data were processed and normalized. Boxplots show the normalized gene expression profiles, and principal component analysis (PCA) scatter plots show significant differences between active UC samples and normal controls ( Supplementary Figures 1A–D ). We performed RRA analysis to determine statistically changed genes in different datasets and then subjected them to GSEA. RRA analysis revealed 1838 statistically changed genes (adj. p-value < 0.05). GSEA indicated that the top annotated collection of genes was enriched in innate immune system, indicating that innate immune responses may play an important role during active UC ( Figures 2A, B ).

Several immune-related gene sets were enriched in the active UC group as compared with that in the normal group, especially with regards to the innate immune system. We explored the immune landscape in training datasets using the CIBERSORT algorithm. The abundance of 22 immune cell types is shown using stacked bar plots ( Figures 3A–D ). Cell types were filtered if not present in half of the samples, and their relative proportions are shown in the boxplots ( Figures 3A–D ). The results revealed that the active UC tissue was infiltrated by a higher fraction of activated DCs and macrophages (M0 and M1) and a lower fraction of M2 macrophages in all training datasets (FDR < 0.05), whereas other immune cells exhibited heterogeneity ( Figures 3A–D ). These results indicated that the inflammatory microenvironment reshapes the proportion and distribution of immune cells.

We identified 1072, 431, 490, and 303 DEGs between active UC patients and normal subjects in the GSE16879, GSE48958, GSE75214, and GSE87473 datasets, respectively ( Figure 4A ). Subsequently, RRA analysis integrated the four cohorts, revealing 64 up-regulated and 46 down-regulated DEGs. The top 15 upregulated and downregulated genes are shown using a heatmap ( Figure 4B ). Furthermore, GO and KEGG functional enrichment analyses were performed to determine the biological features of these 110 robust DEGs. GO functional enrichment analysis revealed 487 up-regulated and 221 down-regulated terms (FDR< 0.05) across BP, CC, and MF categories. DEGs were markedly enriched in granulocyte and neutrophil chemotaxis and migration in the BP category. In the CC category, genes were mainly enriched in vesicle and secretory granule terms. Enriched MF terms included cytokine activity and metallopeptidase activity ( Figure 4C ). KEGG pathway analysis revealed 16 up-regulated and 5 down-regulated items (FDR < 0.05), which included the IL-17 and NF-κB signaling pathways ( Figure 4D ). These results indicated that pro-inflammatory factors and pathways were enriched in active UC patients.

Correlation networks are used for identifying clusters of highly correlated genes across microarray samples. We employed WGCNA to construct and analyze active UC-associated networks. We analyzed GSE87473 microarray datasets, containing 106 active UC and 21 normal control samples, and the optimal β value to match approximate scale-free topology criterion compared to other datasets. The adjacency matrix was constructed based on the criterion of gene distribution conformed to a scale-free network when setting the soft-threshold power of β = 12 (R² = 0.85), retaining high connectivity information ( Figure 5A ). In this study, gene clusters were conducted using the hierarchical clustering method, and 11 consensus co-expression modules were identified ( Figure 5B ). The heatmap shows the association between each module and clinical traits ( Figure 5C ). The correlation analysis of module membership (MM) and gene significance (GS) indicated that the selected module genes exhibited good correlation with the red MM (R = 0.57, P = 1.2e-24, Figure 5D ), implying that genes in the red module were highly correlated with active UC. Genes within the red module were selected, and their biological function was inferred using GO and KEGG analyses ( Figures 5E, F ). Enrichment analyses yielded immune-related terms and pathways. A total of 55 overlapping DEGs were obtained via RRA analysis and WGCNA ( Figure 5G ).

An interaction network between proteins encoded by the 55 DEGs was constructed using the STRING database ( Figure 6A ). The interaction network comprised 49 nodes and 88 edges, visualized using the Cytoscape software. The MCODE plugin was used to identify gene cluster modules. According to the filter criteria, two cluster modules were identified. Cluster 1 had the higher score (score: 8.000, 8 nodes and 28 edges), followed by cluster 2 (score: 5.111, 10 nodes and 23 edges; Figure 6B ). Next, we used the cytoHubba-MCC plugin to identify hub genes and obtained eight hub genes when setting score >5000 ( Figures 6C, D ), consistent with cluster 1 in MCODE. We performed differential expression analysis of the eight hub genes in the validation cohort GSE92415. Consistent with our prediction, mRNA expression levels of the eight specifically expressed hub genes were significantly upregulated in active UC samples compared to those in the control samples ( Figure 6E ). Gene expression is spatiotemporally regulated via networks, which consist of interactions between TFs and their direct target genes, influencing development, homeostasis, and pathogenesis. Eight hub genes correlated with active UC traits were tested for TF binding motifs using the iRegulon plugin. The results indicated that six genes, excluding TIMP1 and MMP3, were regulated by NF-κB1 ( Figure 6F ).

In order to explore the expression patterns of hub genes, we performed correlation analysis of the expression data from the GSE92415 validation cohort. Positive correlation was observed between the hub genes ( Figure 7A ). We further explored the relationship between hub genes and immune cells, and the hub genes were significantly associated with pro-inflammatory cell types, such as neutrophils, M1 macrophages, and activated DCs. Hub genes showed a negative correlation with M2 macrophages and plasma cells ( Figure 7B ). We then analyzed the correlation between hub genes and different immune regulators, including immunosuppressive and immunostimulatory factors ( Figures 7C, D ). These analyses verified that hub genes were closely involved in the regulation of immune cell infiltration and inmmunoregulators, thus modulating the immune microenvironment.

To further validate the important role of hub genes in active UC, we employed the DSS-induced mouse colitis model. Studies have shown AGP and its derivates could ameliorate active UC symptoms and promote barrier integrity through the inhibition of NF-κB signaling (15, 47). NF-κB represents a family of TFs containing five members, with the NF-κB p65–p50 (NF-κB1) heterodimer being most widely studied. Mice were intraperitoneally injected with DMSO or AGP during the period of DSS administration ( Figure 8A ). AGP-treated mice exhibited milder colon blockage than that in the control mice ( Figure 8B ). Histologic analysis strengthened the protective effect of AGP against colitis, with less intestinal epithelial damage and more preserved goblet cells ( Figure 8C ). IF staining revealed that AGP attenuates inflammation by preventing NF-κB activation without affecting IκBα expression ( Figures 8D, E ). We also detected the expression levels of NF-κB and hub genes using RT-PCR. NF-κB activity was increased in the DSS group, whereas there was no significantly recused in the DSS-AGP treatment group. CXCL1, IL1B, MMP1, and MMP10 were significantly upregulated in the DSS group; however, the effect was rescued by AGP treatment. No significant changes were observed in MMP7 and PLAU expression ( Figure 8F ). Our results revealed four genes implicated in the progression of active UC, providing insight into the therapeutic mechanism of the traditional herb AGP.

MicroRNAs are widely studied refer to epigenetic regulators that have shown tremendous potential as therapeutic targets for various human diseases. Four online miRNA databases were used to predict the miRNAs regulating NF-κB activity. Interactions indicated by all three miRNA databases were selected. We obtained 55 regulatory miRNAs of NF-κB and constructed a network using the Cytoscape software ( Figure 9A ). The ceRNA network has attracted much attention in recent years and plays an important role in regulating gene expression. The ceRNA network was constructed with miRNA as a bridge to establish the relationship between target gene mRNA and lncRNA by combining miRNA response elements (MREs). We further explored the lncRNAs that interacted with regulator miRNAs using the StarBase 2.0 database and constructed three lncRNA–miRNA subnetworks ( Figure 9B ). Interestingly, we found that the lncRNA XIST mediated crosstalk between the three miRNAs. Thus, we speculated that the ceRNA network might play an important role in active UC by regulating NF-κB expression.