Mitochondrial-related genes PDK2, CHDH, and ALDH5A1 served as a diagnostic signature and correlated with immune cell infiltration in ulcerative colitis

Weighted gene co-expression network construction

The GSE179285 and GSE107499 datasets were used as the training group, with 58 normal samples and 191 UC samples; and the GSE92415 and GSE48634 datasets were set as the testing group, with 47 normal samples and 111 UC samples. The 871 Mito-RGs isolated in the 302 UC samples were subjected to WGCNA analysis. In the following investigation, the value of the soft thresholding power was determined to be 7, as a result of the scale independence reaching 0.9 and having a reasonably high-average connectedness (Figures 1A and 2A). According to the merging of the highly linked modules at a clustering height limit of 0.25 (Figures 1B and 2B), 26 modules were chosen for in-depth analysis. Modules that had been primed and combined were finally shown in a clustering dendrogram (Figures 1C and 2C). To determine the connection between gene co-expression groups, eigen gene connectivity was evaluated. 14 modules were found to be divided into two clusters after we categorized the eigen genes (Figures 1D and 2D). Next, we established the connections between the recognized modules, and the analysis’s findings showed that gene expression was mostly independent between modules (Figures 1E and 2E). We explored significant correlations between modules and clinical features. In this analysis, the module darkgrey was most closely associated with UC (Figures 1F and 2F). A total of 1231 MRGs were found in the training group and 1140 MRGs in the testing group using WGCNA.

Figure 1. Construction of WGCNA co–expression network in training set. (A) The investigation involved the utilization of scale-free fit index and mean connectivity analysis across various soft-thresholding powers. (B) To identify and establish connections between related modules, the clustered dendrograms were truncated at 0.25 heights. (C) The WGCNA modules were assigned and represented by a colored horizontal bar in the cluster dendrograms. (D) The module feature genes were analyzed for collinearity and presented in a heatmap, where red indicates high correlation and blue indicates low correlation. (E) The topological overlap matrix (TOM) for each module was visualized in a WGCNA network heatmap. (F) Module-trait correlations. Positive correlations are represented by red and negative correlations by blue.

Figure 2. Construction of WGCNA co–expression network in validation set. (A) Scale-free fit index and mean connectivity analysis for various soft-thresholding powers. (B) Clustered dendrograms were cut at 0.25 heights to recognize and connect related modules. (C) Cluster dendrograms and module assignment for WGCNA modules. The colored horizontal bar represents the modules. (D) Collinear heat map of module feature genes. The color red represents a high correlation, whereas the color blue suggests the opposite results. (E) Visualization of the WGCNA network heatmap. The heatmap shows the topological overlap matrix (TOM) for each module. (F) Module-trait correlations. Positive correlations are represented by red and negative correlations by blue.

Screening genes for differential expression

2831 genes in the training group and 2883 genes in the testing group were identified as differential expression genes (DEGs) under the criterion of P-adjustment 0.05 and |log₂FC| >0.5. DEGs volcano plots are shown in Figure 3A, 3B, respectively, for the training group and the test group. After that, 610 intersecting genes were obtained in the training group (Figure 3C), and 437 intersecting genes were acquired in the testing group (Figure 3D) based on the intersections of the WGCNA analysis and the differential expression analysis.

Figure 3. Identification of DEGs. (A) and (B) show the Volcano plot for DEGs between healthy controls and UC in the training and test groups, respectively. (C) and (D) Venn diagrams for intersections of DEGs and the WGCNA module in the training and test groups.

Hub genes identification and functional enrichment analysis

Using Venn diagram, we detected 12 UC-Mito DEGs via overlapping MRGs and important DEGs (Figure 4A). To discover the biological roles played by the hub genes in the modules, we conducted a functional analysis. All of the 12 UC-Mito DEGs were shown to be functionally enriched, and the GO barplot showed 15 GO keywords according to P < 0.05 (Figure 4B). According to the data, the enrichment of the biological process (BP) was predominantly linked to the cellular amino acid catabolic process, small molecule catabolic process, and glutamate metabolic process. Oxidoreductase activity, acting on CH-OH group of donors, oxidoreductase activity, acting on paired donors, and sulfur transferase activity are all connected to enriched molecular function (MF). The enrichment of cellular components (CC) is connected to the mitochondrial matrix, oxidoreductase complex, and endoplasmic reticulum chaperone complex. In KEGG analysis, the signaling pathway of glycine metabolism, thiamine metabolism, and the steroid biosynthesis were related (Figure 4C).

Figure 4. (A) Functional enrichment analysis. (B) and (C) GO and KEGG analysis of DEGs.

Identification and evaluation of signature genes

We used two machine algorithms to identify the most significant feature genes with predictive values: LASSO regression analysis to select 5 predicted genes from statistically significant univariate variables (Figure 5A) and RandomForest combined with feature selection to determine the relationship between error rate, number of classification trees (Figure 5B), and 6 genes with relative importance (Figure 5C). Using the intersection of the two methods described above, we produced a Venn diagram to identify three overlapping genes, PDK2, CHDH, and ALDH5A1. (Figure 5D).

Figure 5. (A) Adjustment of feature selection in the minimum absolute shrinkage and selection operator model (lasso). (B) Random forest error rate versus the number of classification trees. (C) The top 12 relatively important genes. (D) Three algorithmic Venn diagram screening genes.

Using the Rms package, we constructed UC diagnostic column line graph models for the signature genes (PDK2, CHDH, and ALDH5A1) (Figure 6A, 6C), and we evaluated the predictive value of these models using Receiver operating characteristic curves analysis. ROC curves were constructed with an area under the curve that was 84.1305 in the training set and 72.8742 in the validation set respectively (Figure 6B, 6D).

Figure 6. Construction and validation of the UC diagnostic column line graph model. (A, B) Column line graphs are utilized to predict the occurrence of UC in the training and validation cohorts. (C, D) ROC curves were created, with an area under the curve of 84.1305 in the training set and 72.8742 in the validation set.

Analysis of immune infiltration

Given the importance of the immune response to the development of UC, we aimed to describe the immunological infiltration landscape in UC. Using the CIBERSORT method, we compared UC samples with healthy control samples to determine whether or not there was a difference in the infiltration of 22 different types of immune cells (Figure 7A). Meanwhile, Figure 7B depicts the relevant heatmap analysis of immune cells. Neutrophils were positively correlated with activated mast cells (r = 0.58), whereas mast cells at rest were negatively associated with activated mast cells (r = −0.78). When compared to the healthy control, UC patients had higher levels of T cells CD4 memory activated, M1 macrophages, mast cells activated, neutrophils plasma cells, according to violin plots (Figure 7C) of the immune landscape created from the CIBERSORT algorithm (P < 0.05). To gain additional insight into the function of PDK2, CHDH, and ALDH5A1 in UC immune response, the correlation between the expression of signature genes and various degrees of immune cell infiltration was analyzed using the corrplot software. The results showed that PDK2 linked closely with Mast cells at rest (correlation = 0.74, P < 0.001) and negatively with M1 macrophages (correlation = −0.54, P < 0.001), as shown in Figure 7D. Figure 7E demonstrates that CHDH showed a positive association with M2 macrophages (correlation = 0.63, P < 0.001) and a negative correlation with M1 macrophages (correlation = −0.51, P < 0.001). As shown in Figure 7F, ALDH5A1 was positively correlated with Mast cells at rest (correlation = 0.74, P < 0.001) and negatively correlated with M1 macrophages (correlation = −0.39, P < 0.001). Furthermore, we employed xCELL algorithm to validate the correlation between immune cells and PDK2, CHDH, and ALDH5A1. The result was demonstrated in Supplementary Figure 1. These data also revealed that signature genes were closely associated with immune cells in UC.

Figure 7. Analysis of immune infiltration. (A) Bar plot bar plot depicting the 22 subpopulations of immune cells in UC and healthy control (HC) samples. (B) Correlation heatmap of immune cells in UC samples. (C) Violin plot illustrating the varying proportions of immune cells between UC and normal control (NC). (D–F) Bubble plots demonstrate the relationship between immune cells and specific signature genes (PDK2, CHDH, and ALDH5A1).

Validation of signature gene expression

We evaluated the expression of these three genes in both the training and validation cohorts of UC patients and UC cells model (Figure 8). The ROC analysis was carried out so that we could further validate the diagnostic value of PDK2, CHDH, and ALDH5A1. It was discovered that PDK2 (AUC: 0.784), CHDH (AUC: 0.778), and ALDH5A1 (AUC: 0.813) all had better diagnostic values. The following findings were also confirmed by the validation datasets: PDK2 (AUC: 0.676), CHDH (AUC: 0.642), and ALDH5A1 (AUC: 0.694). (Figure 8A–8F).

Figure 8. (A–C) The mRNA expression levels of the feature genes in the training cohort, as well as their ROC curves. (D–F) expression levels and ROC curves for the feature genes present in the training cohort. (G–I) the mRNA expression of three feature genes were validated using the LPS-induced NCM460 cells.

Furthermore, the quantitative real-time PCR technique was employed to ascertain the downregulation of PDK2, CHDH, and ALDH5A1 expression levels in LPS-induced NCM 460 cells (Figure 8G–8I) and HIEC-6 cells (Figure 9A). The western-blotting method corroborated the decreased expression of ALDH5A1, PDK2, and CHDH in LPS-induced HIEC-6 cells (Figure 9B). To substantiate the involvement of ALDH5A1, PDK2, and CHDH in vivo, colitis was induced using a 3% DSS solution, resulting in a shorter colon length in the colitis group compared to the control (PBS) group (Figure 9C, 9D). Histological analysis revealed conspicuous infiltration of inflammatory cells in the intestinal mucosa and extensive epithelial erosion in DSS-induced mice (Figure 9E, 9F). Additionally, we assessed the protein expression levels of ALDH5A1, PDK2, and CHDH in the colon tissues of DSS-treated mice. Western blotting analysis demonstrated a reduction in the protein expression levels of ALDH5A1, PDK2, and CHDH in the colitis lesions of mice administered with DSS (Figure 9G).

Figure 9. Expression of PDK2, CHDH, and ALDH5A1 in LPS-induced HIEC-6 cells and DSS-induced colitis mice. The relative mRNA (A) and protein (B) expression levels of PDK2, CHDH, and ALDH5A1 were conducted using the LPS-induced HIEC-6 cells. (C) Photographs were utilized to assess the variations in colon lengths among the groups. (D) The length of the colon in different groups. (E, F) Hematoxylin and eosin staining images and colonic histological scores of colon cross-sections in mice (n = 5). Scale bar = 50 μm. (G) Expression of PDK2, CHDH, and ALDH5A1 in the colon tissues of mice was detected using Western blotting. (H) Double immunofluorescence staining of PDK2 (labeled in green), CHDH (labeled in green), and ALDH5A1 (labeled in green) with M1 macrophages (CD86, labeled in red) at DSS induced mice. The nuclei were stained by DAPI (blue). n = 5 per group. Scale bar = 50 μm. Values are expressed as the mean ± sd (n = 5). ^*P < 0.05, ^**P < 0.01, ^***P < 0.001. Abbreviation: DSS: dextran sodium sulfate.

Subsequently, immunofluorescence analysis was conducted to evaluate the association between immune cell infiltration and the expression of PDK2, CHDH, and ALDH5A1. The expression of CD86 (antigen-presenting surface molecules of M1 macrophage) was examined in relation to PDK2, CHDH, and ALDH5A1 expression. We found that CD86 was not detectable in control specimens. And the results also demonstrated that CD86 exhibited higher expression in mouse colitis lesions, whereas PDK2, CHDH, and ALDH5A1 displayed lower expression in these lesions (Figure 9H). Collectively, our findings indicate that the influence of PDK2, CHDH, and ALDH5A1 on the advancement of DSS-induced colitis in mice is potentially mediated through the regulation of macrophage infiltration. Nevertheless, the precise mechanism underlying the moderation of colitis remains to be elucidated.

Data processing and download

The gene expression array data and corresponding clinical datasets (GSE179285, GSE107499, GSE92415, and GSE48634) were retrieved from the Gene Expression Omnibus database. The training set for this study consisted of GSE179285 and GSE107499, while the validation set comprised GSE92415 and GSE48634. The former dataset included healthy controls (n = 58), ulcerative colitis patients with paired uninflamed colon (n = 48), and inflamed colon (n = 46), while the latter dataset included uninflamed colon (n = 22) and inflamed colon (n = 97). The GSE92415 dataset comprised healthy controls (n = 21) and ulcerative colitis patients (n = 87), and the GSE48634 dataset included healthy controls (n = 26) and ulcerative colitis patients (n = 24). Mitochondrial-related genes in the present study were downloaded from MSIGDB [3] (https://www.gsea-msigdb.org/gsea/msigdb/genesets.jsp?collection=H).

Weighted gene co-expression network construction

The present study employed the R package “WGCNA” to extract expression profiles of both UC and healthy controls, thereby establishing a weighted gene co-expression network [41]. The evaluation of scale independence and mean connectivity was performed using the “pick soft threshold function.” Subsequently, modules were identified through hierarchical clustering and the dynamic tree cut function. The modules were then linked to disease features, and gene significance and module membership calculations were conducted. Further analysis was carried out on the matching module gene details. Finally, we displayed the eigen genes network.

Identification of differentially expressed mitochondria-related genes

The MRGs expression profile was obtained from the transcriptome data that was downloaded, and a differential expression analysis was performed between healthy controls and UC patients using criteria of |Log fold change (FC)| >1 and false discovery rate (FDR) adjusted P < 0.05. The “ggplot” package was employed to generate a volcanic plot for visualization purposes.

The acquisition of key genes

Initially, the intersecting genes were chosen from the Weighted Gene Co-expression Network Analysis (WGCNA) and differential expression analysis (diff) in the training set. Subsequently, the intersecting genes were selected from the WGCNA and diff in the validation set. Following this, the common genes between the aforementioned steps and the mitochondrial genes were identified and deemed as pivotal genes for further analysis.

Functional enrichment analysis

Genes in interesting modules were taken out for further functional enrichment analysis. Essential genes’ biological properties were predicted and analyzed using a Gene Ontology analysis (GO) [42]. To find functional characteristics, KEGG pathway enrichment analysis was carried out [43]. The cut-off for significance was P < 0.05.

Identification and evaluation of signature genes

The top gene list was produced using the least absolute shrinkage and selection operator (LASSO) analysis, which was used to exclude confusing genes and screen important genes for relationships between them. Additionally, we employed random forest for feature selection, which has been commonly used and can accurately evaluate the significance of each feature in the dataset, to pick out more meaningful key genes. With the use of LASSO logistic regression and RandomForest, the intersection was then used to pick out the study’s most important feature genes. Training set and validation set further confirmed the expression of important genes. The “Rms” program in R was used to produce nomograms, which were then used to determine the independent prognostic factors. The “Rms” program in R was used to produce nomograms, which were then used to determine the independent prognostic factors. To calculate the area under the curve (AUC) and assess the diagnostic value of the signature genes, we utilized the receiver operating characteristic (ROC) curves in the R package.

Immune infiltration analysis

The CIBERSORT method was employed to investigate the proportion and composition of immune cells in UC and normal samples. Following that, immune cells correlation heatmap was constructed using the “corrplot” package. In order to show the dispersion of the Immune cells, the “vioplot” R package also was used. The Spearman rank correlation coefficient was calculated using the “corrplot” package to further investigate the link between immune cells and certain genes expression.

UC cells model construction

NCM460 and HIEC-6 cell line were obtained from the Key Laboratory of Hubei Province for Digestive System Disease (Wuhan, Hubei, China) and cultured in RPMI-1640 medium (SH30809; HyClone, USA) and Cytiva HyClone™ DMEM/F12 (SH30023.01; HyClone, USA) containing 10% fetal bovine serum (FBS, sijiqing, Hangzhou, China) at 37°C in a humidified atmosphere of 5% CO₂. NCM460 and HIEC-6 cells were seeded and then treated with LPS (1 μg/ml; Sigma) for 24 h to construct the UC cells model in vitro. After 24 h, LPS-induced NCM460 and HIEC-6 cells were harvested for further experiments.

Animals and induction of colitis

Male C57BL/6 mice, aged 6 to 8 weeks and weighing 15–20 g, were procured from Hubei Biont Biological Technology Co., Ltd. (Wuhan, China). The mice were housed in a controlled environment with a 12-hour light/dark cycle and maintained at a standard temperature of 20 ± 2°C. Prior to the commencement of the experiment, the mice were acclimatized to these conditions for a minimum of 7 days. All animal procedures were approved by the Animal Ethics Committee of Guizhou Provincial People’s Hospital (Animal, 2022-039). From day 8 to day 14, the mice were administered sterile drinking water supplemented with 3% (w/v) dextran sulfate sodium salt (DSS, MP Biomedicals, Irvine, CA, USA) to induce experimental murine colitis. The control group received sterile drinking water. Following a 1-week period of 3% DSS treatment, mice were anesthetized and euthanized. The length of the colon was measured after isolation. A portion of the distal colon was then fixed with 4% paraformaldehyde (Biosharp) at room temperature for 24 hours to facilitate hematoxylin and eosin staining. The remaining sections were promptly frozen in liquid nitrogen for future assays. The histological scores were determined based on established scoring criteria.

Quantitative real-time PCR

The total RNA of cells and tissues were extracted from Trizol (Cat#9109; TaKaRa, Japan), and the PrimeScript RT Master Mix kit was used to reverse it into cDNA (Cat#RR047A; Takara, Japan). All of the mRNA levels were determined by employing the SYBR Green PCR Mix (Cat#RR420A; Takara, Japan) in conjunction with the CFX Connect (Bio-Rad, USA) as previous study [44]. The 2^−ΔΔCT approach was utilized for data analysis, and each experiment was carried out with three separate sets of controls. The primer sequences used in this work were all synthesized by Sango Biotech (Shanghai, China), and they were all listed in the Table 1.

Histopathological and immunofluorescence

Tissue samples were fixed in a 4% paraformaldehyde fixation solution, embedded in paraffin, cut into 5-μm sections with a microtome, fully dewaxed, and hydrated. A hematoxylin-eosin staining solution was used to produce HE sections. Finally, we conducted the observation and evaluation using a microscope.

The immunofluorescence technique was conducted following the previously described protocol [45]. In summary, colon tissues were fixed in 4% paraformaldehyde at 4°C overnight, followed by immersion in 30% sucrose until they reached the bottom. Subsequently, the colons were embedded in optimal cutting temperature compound and frozen. Slices with a thickness of 10 μm were sequentially obtained at −20°C using a cryostat. These slices were then permeabilized in PBS containing 0.3% Triton X-100 for 30 minutes and blocked in PBS with 5% donkey serum at 37°C for 30 minutes.

Next, each section of the colon was subjected to overnight incubation at 4°C with the primary antibodies, namely ALDH5A1 polyclonal antibody (17319-1-AP), CHDH polyclonal antibody (17356-1-AP), PDK2 polyclonal antibody (15647-1-AP), and CD86 polyclonal antibody (13395-1-AP). Following this, the sections were washed and exposed to the corresponding secondary antibodies conjugated with fluorescence for a duration of 1 hour at 37°C. The nucleus was counterstained with 4,6-diamino-2-phenylindole (DAPI, Biosharp). The immunoreactivity density of ALDH5A1, CHDH, PDK2, and CD86 was quantified using ImageJ software (National Institutes of Health, MD, USA).

Protein extraction and western blotting

Cells and tissues were subjected to protein extraction. Subsequently, the extracted samples were subjected to electrophoresis using a 10% SDS-PAGE gel and transferred onto a polyvinylidene fluoride membrane. Following this, the membranes were blocked with protein-free rapid blocking buffer (EpiZyme, Shanghai, China) for 0.5 hours. The membranes were then incubated overnight at 4°C with the following antibodies: ALDH5A1 polyclonal antibody (17319-1-AP), CHDH polyclonal antibody (17356-1-AP), PDK2 polyclonal antibody (15647-1-AP), and GAPDH (60004-1-Ig) from ProteinTech group (Wuhan, China).

Data and statistical analysis

All statistical analyses were performed using R version 4.0.1 and GraphPad Prism 8.0 statistical software (La Jolla, CA, USA). The data in this study are presented as mean ± SD. Comparisons between two groups were conducted using Student’s t-test, Mann-Whitney U test, or Dunnett’s t-test, as appropriate. Differences between three or more groups were assessed using analysis of variance, and multiple comparisons between different groups were carried out using ordinary 1- or 2-way analysis of variance with a Sidak test. Statistical significance was defined as P < 0.05.

Availability of data and material (ADM)

The Supplementary Material for this article can be found online at this study, further inquiries can be directed to the corresponding authors. The datasets analyzed during the current study are available online from GSE179285, GSE107499, GSE92415 and GSE48634 dataset.