Identification of key genes and biological processes contributing to colitis associated dysplasia in ulcerative colitis

Data selection and processing

By searching the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) database, the microarray dataset GSE47908 based on the GPL570 platform was selected for analysis because it embraces the transcriptional profiles of colonic mucosa from both UC and CAD patients (Bjerrum et al., 2014). Specifically, it included data from 45 patients with UC (20 with left-sided colitis, 19 with pancolitis, and six with UC-associated dysplasia) and 15 healthy controls. After obtaining the data, principal component analysis (PCA) was performed to check and visualize the grouped data. For validation, another microarray dataset GSE87466 comprising the gene expression profiles of colonic mucosa from 87 patients with UC and 21 healthy controls was selected to check the gene expression trends in UC (Li et al., 2018). For comparison, we downloaded the transcriptome and clinical data of 480 patients with sporadic colorectal cancer (CRC) and 41 healthy controls from The Cancer Genome Atlas (TCGA) database (date limit: February 23, 2021), using “TCGAbiolinks” R package version 2.16.4. The gene expression profiles were normalized using the voom function in “limma” R package (Ritchie et al., 2015).

Differential expression analysis

Differential expression analysis was performed using the “limma” R package (Ritchie et al., 2015). The sum of the mean value of absolute log₂ fold change (FC) and the two standard deviations of absolute log₂ FC was used as the cut-off value of log₂ FC. Genes with absolute log₂FC >log₂ FC cut-off and adjusted p value <0.05, were considered as DEGs. The results are visualized as volcano plots.

WGCNA

The expression profiles of DEGs between CAD and control group were extracted to construct weighted gene co-expression network using R package “WGCNA” (Langfelder & Horvath, 2008; Zhang & Horvath, 2005). First, the candidate soft-thresholding powers (1 to 30) were used to calculate the scale independence and mean connectivity using the pickSoftThreshold function. The first candidate power, whose degree of independence was >0.8, was selected as the proper power. Then, a co-expression network was constructed and modules were identified using the blockwiseModules function, with the parameters mergeCutHeight set to 0.28 and minModuleSize set to 30. Each module was assigned a unique color. Finally, the Pearson correlation coefficient and corresponding p value between each module’s eigengene and phenotype were calculated. In addition to the UC and CAD phenotypes, another phenotype “Progress” was also included to quantify the dynamic process of UC-associated carcinogenesis. The phenotype “Progress” was assigned the values of 0,1, and 2, to represent the control, UC, and CAD groups, respectively; this reflected the progressive process from normal to UC and then to CAD. It means that if a module is significantly correlated to “Progress”, then this module is considered to be correlated to the progressive process from normal to CAD.

Identification of hub genes

For a specified module gene, the module membership (MM) calculated by signedKME function represents its degree of importance in the module, and gene significance (GS) calculated by the cor function represents the degree of correlation with the phenotype. By calculating MM and GS, the module genes with high MM and GS can be defined as hub genes. In the present study, thresholds of MM >0.9 and GS >0.7, were chosen to screen for the hub genes of each module (Zou et al., 2019).

Function enrichment analysis

Functional enrichment analysis was performed using the “clusterProfiler” R package (Yu et al., 2012). Significantly enriched Gene Ontology-Biological processes (GO-BP) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were identified using enrichGO and enrichKEGG functions, respectively. Gene set enrichment analysis (GSEA) was performed based on the gene list sorted by log₂FC obtained from differential expression analysis using gseGO and gseKEGG functions (Subramanian et al., 2005). Single-gene GSEA was conducted based on the gene list sorted by Spearman correlation coefficient between every gene and the specified hub gene to predict the significant biological processes and pathways associated with the hub gene. A p value <0.05 was considered significant.

GEO data overview

A flow chart of the present study is shown in Fig. 1A. As for GSE47908 dataset, 15 controls, 20 patients with left-sided colitis, and six patients with UC-associated dysplasia were enrolled in the present study, representing the control, UC, and CAD groups, respectively. Clinical data including sex, age, disease duration, Mayo score, Mayo endoscopic score, smoking status, and daily medication were extracted and compared using the chi-squared test (Bjerrum et al., 2014). Results showed that there were no significant differences in sex, disease duration, smoking status, or daily medication between the groups. As for age, Mayo score, and Mayo endoscopic score, only median and interquartile ranges (IQR) were known, and the corresponding p-values could not be calculated (Table S1). After processing, four outliers were removed and a total of 37 samples, including 13 controls, 18 patients with UC, and six patients with CAD, were enrolled. As shown in the PCA plot (Fig. 1B), the gene transcriptional profiles of these three groups were clearly distinct from each other. The PCA plot of gene expression profiles of 87 patients with UC and 21 healthy controls in GSE87466 dataset shown in Fig. 1C, also indicates a significant distinction between these two groups.

Differential expression analysis

To investigate the genetic changes during the progress from normal to UC and then to CAD, differential expression analyses between UC and control, CAD and UC, and CAD and control groups in GSE47908 were performed. The results are presented in Fig. 2. There were 679 upregulated and 288 downregulated genes in UC compared to those in the control group (Fig. 2A, log₂ FC cut-off of 1.092), 260 upregulated and 616 downregulated genes in CAD compared to those in the UC group (Fig. 2B, log₂ FC cut-off of 1.073), and 305 upregulated and 462 downregulated genes in CAD compared to those in the control group (Fig. 2C, log₂ FC cut-off of 0.804). Venn diagrams were then used to further detail the genetic changes during this dynamic process. There were 67 co-up DEGs (Fig. 2D) and 46 co-down DEGs (Fig. 2E) in UC vs. controls and CAD vs. controls. In addition, the Venn diagram shows that among 260 upregulated and 616 downregulated genes in CAD vs. UC, only 61 upregulated (Fig. 2D) and 180 downregulated genes (Fig. 2E) changed significantly compared to the control group. This indicates that some DEGs dysregulated during the progress from normal to UC may continue to change in CAD and even contribute to the initiation of CAD, while others might be only responsible for the development of UC and may change back to the normal level if the inflammation is controlled.

To gain further insights into the expression pattern of DEGs during the progress from normal to UC and then to CAD, the top 40 upregulated and top 40 downregulated DEGs between CAD and controls after sorting by adjusted p-value were visualized in a heat map (Fig. 2F). Row cluster analysis showed that these DEGs were grouped under two clusters: upregulated and downregulated genes in CAD compared to those in controls. However, when UC samples were taken into consideration, every cluster could be further divided into two sub-clusters, showing different changing trends of these DEGs during the progress from normal to CAD.

It should be noted that no co-up or co-down DEGs in UC vs. control groups and CAD vs. UC groups were identified in Venn diagrams, while the heatmap showed marked continuous changes in the expression patterns of some DEGs from control to UC and then to CAD. The traditional filter method for DEGs can potentially miss somekey information. To make our analysis more comprehensive, the cut-off criteria of log₂ FC was set to 0.3, leading to the identification of 4929 DEGs between the CAD and control groups, which is a suitable number of DEGs for WGCNA analysis.

WGCNA

WGCNA was performed using 4929 DEGs between the CAD and control groups. By setting the soft threshold to 12 (Figs. 3A–3B), a total of six functional modules, apart from the gray module, were identified (Fig. 3C). The turquoise, red, and blue modules were positively correlated with CAD, and the brown, yellow, and green modules were negatively correlated with CAD. Furthermore, the green and blue modules were correlated with UC and CAD in the same direction respectively, and were found to be correlated with “Progress” with the highest correlation coefficients (green module: correlation coefficient = −0.75, p-value = 9e−08; blue module: correlation coefficient = 0.71, p-value = 9e−07) (Fig. 3D). This indicates that genes belonging to blue and green modules and the relevant biological processes begin to change at an earlier stage and continue to change throughout the whole process of progression from normal to UC and then to CAD.

Function enrichment analysis

The GO-BP enrichment analysis of each module is shown in Table S2. Green module genes were significantly enriched in biological processes related to mitochondrial function and energy metabolism (Fig. 4A). Consistent with this observation, KEGG pathway enrichment analysis indicated that the green module genes were significantly enriched in energy metabolism-related pathways, such as the TCA cycle, carbon metabolism, and oxidative phosphorylation (Fig. 4B). Blue module genes were significantly enriched in biological processes related to cell–cell junctions and immune responses, especially neutrophil-mediated immunity (Fig. 4C). The KEGG pathways that blue module genes were significantly enriched in included infection-related pathways, focal adhesion, tight junction, Hippo signaling pathway, Notch signaling pathway, and ErbB signaling pathway (Fig. 4D).

To further explain the role/s of mitochondrial dysfunction during the process of UC-associated carcinogenesis, GSEA was performed to analyze the alterations in mitochondrial function-related GO-BP in the UC and CAD groups. The results indicated that mitochondrial function-related GO-BP, including cellular respiration, oxidative phosphorylation, respiratory electron transport chain, and mitochondrial gene expression, were all significantly downregulated in UC (Fig. S1A) and CAD (Fig. S1B) groups as compared to those observed for the control group.

To explore the differences between UC-CRC and sporadic CRC in terms of the alteration of mitochondrial functions, the gene expression profiles of 480 sporadic CRC patients and 41 healthy controls from the TCGA database were analyzed. After differential expression analysis, GSEA of GO-BP was conducted based on the gene list sorted by log₂ FC. The results showed that biological processes, including mitochondrial gene expression and mitochondrial translation, were upregulated in sporadic CRC patients, and no other mitochondrial function-related biological processes were identified (Table S3). This indicates that the carcinogenic mechanism underlying sporadic CRC is different from that of UC-CRC.

Identification and validation of hub genes

Based on MM >0.9 and GS >0.7, seven genes (SLC25A3, ACO2, AIFM1, ATP5A1, DLD, TFE3, and UQCRC1) were identified as hub genes in the green module (Fig. 5A) and six hub genes (ADIPOR2, SLC35D1, TOR1AIP1, PRR5L, ATOX1, and DTX3) were identified in the blue module (Fig. 5C). The differential expression patterns of these 13 hub genes during the progression from normal to UC and then to CAD are displayed in Figs. 5B and 5D. All of them were observed to be continuously up/downregulated significantly, supporting their potential persistent roles in UC-associated carcinogenesis.

Considering the small sample size of the GSE47908 dataset, validation using another dataset was considered indispensable. To validate the changes in the inflammatory phase, differential expression analysis of the 13 hub genes was performed in UC patients from GSE87466 dataset (87 patients with UC and 21 controls). The results showed that the expression changes of hub genes in UC patients of GSE87466 dataset were consistent with those observed for patients belonging to GSE47908 dataset (Fig. 6). To further explore the potential association of hub genes with UC-associated carcinogenesis, a brief literature review was performed to search for the proof of their carcinogenic effects (Bruggemann et al., 2017; Byeon et al., 2010; Chung et al., 2017; Ciccarone et al., 2020; Ding et al., 2020; Jana et al., 2020; Karginova et al., 2019; Khalil, 2007; Kim et al., 2019; Laiho et al., 2003; Li et al., 2019; Linehan et al., 2019; Liu et al., 2018; Millan & Huerta, 2009; Oehler et al., 2009; Perera et al., 2015; Rao et al., 2019; Seth et al., 2009; Thedieck et al., 2007; Viola et al., 2017; Wang et al., 2013; Wang et al., 2020). The results are presented in Table 1. With the exception of TOR1AIP1 and PRR5L, almost all other hub genes have been reported to be involved in cancers. Furthermore, five of the hub genes (AIFM1, ATP5A1, UQCRC1, ADIPOR2, and ATOX1) may play roles in the carcinogenesis of sporadic CRC. These known associations between hub genes and cancers suggest the possible carcinogenic effects of these hub genes in CAD.

Predication of the potential function of TFE3 in UC-associated carcinogenesis by GSEA

Among the proteins encoded by the hub genes, TFE3, a transcription factor, has been known to be involved in the onset and progress of cancers by regulating many biological processes, such as energy metabolism, lysosomal biogenesis, and immune response (Beckmann, Su & Kadesch, 1990; Perera et al., 2015; Willett et al., 2017). In this study, we found that TFE3 was continuously upregulated during UC-associated carcinogenesis. However, the plausible roles of TFE3 in this process remain unclear. Single-gene GSEA revealed that biological processes and KEGG pathways related to angiogenesis and immune response were positively correlated with the upregulation of TFE3 (Figs. 7A, 7C), whereas those related to mitochondrial function and energy metabolism were negatively correlated with the TFE3 upregulation (Figs. 7B, 7D).