Key genes of inflammation and Crohn’s disease severity identified via bioinformatics and clinical specimen analysis

Abstract

Crohn’s disease (CD) is a nonspecific inflammatory disease of digestive tract with limited known biomarkers. We aimed to identify potential molecules involved in CD severity. Weighted gene co-expression network and differentially expressed gene analyses were used to explore hub genes associated with inflammation and involved in pathological processes of CD in the GSE83448 cohort. The expression and relationship of genes with clinical characteristics in normal and CD intestinal mucosa were validated in GSE208303 and our research cohorts. Five genes related to disease and inflammation were identified in the clinically significant modules that were differentially expressed in healthy controls and CD patients. Among them, solute carrier family 39 member 5 (SLC39A5) expression was reduced in CD tissues compared with that in normal tissues; low SLC39A5 expression was associated with CD-simplifying endoscopic score and stenosis. Furthermore, we broadened the bioinformatic threshold with qRT-PCR verification and western blot, and found that the expressions of cytochrome P450 family 2 subfamily S member 1 (CYP2S1), fatty acid 2-hydroxylase (FA2H), fatty acid binding protein 1 (FABP1), and Rho family GTPase 3 (RND3) were reduced in CD tissues. Therefore, SLC39A5, CYP2S1, FA2H, FABP1, and RND3 may be novel biomarkers for CD and associated with severity.

Highlights

1. SLC39A5 expression was reduced in CD tissues.

2. Low SLC39A5 expression was associated with CD-simplifying endoscopic score and stenosis.

3. CYP2S1, FA2H, FABP1, and RND3 were reduced in CD tissues.

KEYWORDS:

• Crohn’s disease
• weighted gene co-expression network analysis
• differentially expressed gene
• biomarkers
• SLC39A5

Introduction

Crohn’s disease (CD) is a disease of the digestive system characterized by chronic segmental transmural inflammatory bowel disease that may affect the entire digestive tract, especially the distal ileum and adjacent colon (Ford 2021). The severity of CD varies among patients, and the main clinical symptoms include abdominal pain, diarrhea, fistula formation, and penetrating inflammation (Cushing and Higgins 2021). Microbial damage, genetic predisposition, changes in intestinal barrier permeability, stress, and immune abnormalities are believed to be involved in CD pathogenesis (Guo and Shen 2021; Hornschuh et al. 2021). Mesalamine, hormones, and immunosuppressants are the main therapeutic strategies for CD (Sulz et al. 2020). However, most patients experience recurrence, and the specific etiology of CD remains unclear. Therefore, exploring the molecular mechanisms involved in CD may contribute to its diagnosis and treatment.

In recent years, high-throughput sequencing and gene microarray, followed by bioinformatics analysis, have been employed to explore the key molecules involved in disease progression, contributing to the diagnosis and prognosis of disease (Pareek et al. 2011; McCombie et al. 2019). The Cancer Gene Expression Omnibus (GEO) is a major public database that stores various gene expression profiles of hundreds of diseases (McCombie et al. 2019; Albaradei et al. 2021). Differentially expressed gene (DEG) analysis and weighted gene co-expression network analysis (WGCNA) are two main bioinformatics analysis methods used to explore the key genes involved in disease progression. DEG analysis focuses on changes in gene expression levels among different samples (Costa-Silva et al. 2017). WGCNA is a systematic biological approach that focuses on the relationship between gene modules and clinical characteristics (Zhang and Horvath 2005). To date, a series of biomarkers for digestive system diseases have been identified using bioinformatics methods. For example, Xie et al. performed WGCNA and demonstrated that several genes associated with primary immunodeficiency, including CD4, CD3D, and CD40 ligands, are involved in the occurrence of inflammatory bowel diseases (Xie et al. 2018). By performing DEG analysis in a microarray cohort, Liu et al. demonstrated that four genes associated with m6A methylation regulation play key roles in colorectal cancer (Wei et al. 2021). However, effective biomarkers for CD remain unknown. Therefore, identifying the key genes involved in inflammation and CD severity has clinical significance in diagnosis and treatment.

In this study, bioinformatics analysis combined with qRT-PCR verification and clinical case review were used to identify key genes in CD. SLC39A5, cytochrome P450 family 2 subfamily S member 1 (CYP2S1), fatty acid 2-hydroxylase (FA2H), fatty acid binding protein 1 (FABP1), and Rho family GTPase 3 (RND3) were identified as core genes in CD, with lower expression in CD tissues. Among them, we found that low SLC39A5 expression was associated with CD, simplifying endoscopic score and stenosis. In conclusion, SLC39A5, CYP2S1, FA2H, FABP1, and RND3 may be key molecules involved in CD, which may act as biomarkers and targets for CD diagnosis and therapy.

Methods

Transcriptional profile preprocessing

The transcriptional profile (accession number: GSE83448) (Zabana et al. 2020) of intestinal mucosa biopsies obtained from healthy controls and CD patients and corresponding information of clinical characteristics were obtained from the GEO database (https://www.ncbi.nlm.nih.gov/gds) and set as the discovery cohort. The transcriptional profile was provided by Josep M, and was performed on the GPL18134 platform using GE Healthcare/Amersham Biosciences CodeLink Human Whole Genome Bioarray. In total, 14 intestinal mucosa biopsies obtained from healthy controls and 39 intestinal mucosa biopsies obtained from CD patients were included in the GSE83448 cohort. The gene expression profile of GSE208303 (Hong et al. 2017), containing nine intestinal mucosa biopsies provided by healthy controls and 26 intestinal mucosa biopsies obtained from CD patients, was used as the verification cohort. Before performing bioinformatics analysis, raw data were downloaded and pre-processed as follows: probe names were converted to gene names, the expression of the same genes was merged, genes with medium expression levels ≤1 were removed, and gene expression was normalized.

WGCNA construction and identification of significant modules

In this study, gene expression data of intestinal mucosa biopsies obtained from healthy controls and CD patients in the GSE83448 cohort were imported into R software (version: 4.0.5) (Chan 2018) and used to construct WGCNA using a standard procedure (Langfelder and Horvath 2008). First, hclust was used to cluster samples, and samples with heights >200 were set as outliers. After removing the outliers, a suitable soft threshold that could maintain scale-free topology >0.85 and mean connectivity close to 0 was set. A scale-free network was then constructed based on a suitable soft threshold. The weighted adjacency matrix of the gene pairs was then calculated and transformed into a topological overlap matrix (TOM). The hierarchical clustering method was used to construct the clustering tree structure of the TOM. Different branches of the cluster tree were represented by different modules. Modules with adjacency value <0.2 were merged, and the software automatically named each gene module with a color name. Genes with no co-expression relationships were placed in the gray module. The association between gene modules and clinical characteristics (disease state: with or without CD; inflammation: with or without inflammation) was calculated. Gene modules that were both associated with disease and inflammation with |correlation score| >0.3 and P<0.05, were set as potential clinically significant modules. Subsequently, the module membership (MM) and gene significance (GS) for each module were calculated; the MM of genes in the modules was significantly associated with their GS with a correlation score >0.3 and P<0.05, and the modules were set as clinically significant modules. Genes in clinically significant modules were used for further analysis.

DEG analysis and identification of potential hub genes

The DEGs between intestinal mucosa biopsies obtained from healthy controls and CD patients, as well as those with or without inflammation, were analyzed in R software using the edgeR package (Robinson et al. 2010). Genes with |LogFC| ≥ 1.0 or ≥0.5 and adjusted P-value < 0.05 were set as thresholds to consider as DEGs. DEGs and genes in clinically significant modules were used for intersection analysis, and intersecting genes were set as potential hub genes.

Protein–protein interaction (PPI) construction and identification of hub genes

Information on all potential hub genes was imported into the STRING online tool (https://string-db.org/) for PPI construction (Athanasios et al. 2017). The lines between the genes in the PPI network indicate connections.

Tissue collection and ethical statement

A total of 31 intestinal mucosa biopsies (from ileum and colon) of CD and adjacent tissues were collected during intestinal endoscopy in the gastroenterology department of Zhujiang Hospital Affiliated to Southern Medical University (Guangzhou, China). Tissue collection was approved by the Human Ethics Committee of Zhujiang Hospital Affiliated to Southern Medical University (approval number 2021-KY-150-01). Written informed consent was obtained from all patients enrolled in the experiments. The tissues were stored at −80°C prior to experimentation.

qRT-PCR

Total RNA was extracted from intestinal mucosa biopsies of CD and adjacent tissues using TRIzol reagent (Yeasen, Shanghai, China). The PrimeScript™ RT Reagent Kit (Takara, Japan) was used to reverse transcribe 800 ng of mRNA into cDNA, following concentration determination using a spectrophotometer. The final qRT-PCR was performed using 2 µL of cDNA, and SYBR Green (Thermo Fisher Scientific, USA) was used to detect the mRNA levels of the target genes. Gene expression was calculated using the 2^-△△CT method and normalized to GAPDH mRNA level. The primers used in this study were as follows:

GAPDH forward: 5′-GGAGCGAGATCCCTCCAAAAT-3′;

GAPDH reversed: 5′-GGCTGTTGTCATACTTCTCATGG-3′;

SLC39A5 forward: 5′-GGTCCCCAGTGAGTCATCTG-3′;

SLC39A5 reversed: 5′-GCCAGGTAATGGTTCTGCTCC-3′;

CYP2S1 forward: 5′-GCGCTGTATTCAGGGCTCAT-3′;

CYP2S1 reversed: 5′-CTTCCAGCATCGCTACGGTT-3′;

FA2H forward: 5′-CGCCTCTACGACCTCTCCA-3′;

FA2H reversed: 5′-TGGTTCCATAGCAGGATCTGTC-3′;

FABP1 forward: 5′-ATGAGTTTCTCCGGCAAGTACC-3′;

FABP1 reversed: 5′-CTCTTCCGGCAGACCGATTG-3′;

RND3 forward: 5′-GCTCCATGTCTTCGCCAAG-3′;

RND3 reversed: 5′-AAAACTGGCCGTGTAATTCTCA-3′;

ERBB3 forward: 5′-GGTGATGGGGAACCTTGAGAT-3′;

ERBB3 reversed: 5′-CTGTCACTTCTCGAATCCACTG-3′.

Western blotting

Total protein in intestinal mucosa biopsies of CD and adjacent tissues was extracted using RIPA reagent (Beyotime, Guangzhou, China), and the protein concentration of each sample was determined by the bicinchoninic acid method. Protein samples were separated by 12% Omni-PAGE (Epizyme, Shanghai, China) and transferred to PVDF membranes (Epizyme, Shanghai, China). The membranes were blocked with 8% skim milk powder (Epizyme, Shanghai, China) for 2 h and then incubated with anti-GAPDH (1:2000; 60004-1-Ig, Proteintech, Wuhan, China), anti-SLC39A5 (1:500; 17285-1-AP, Proteintech, Wuhan, China), anti-CYP2S1 (1:500, 19100-1-AP, Proteintech, Wuhan, China), anti-FA2H (1:1000, 15452-1-AP, Proteintech, Wuhan, China), anti-FABP1 (1:2000, 13626-1-AP, Proteintech, Wuhan, China), and anti-RND3 (1:1000, 66228-1-Ig, Proteintech, Wuhan, China) at 4 °C overnight. After washing twice with PBS and incubating with secondary antibodies for 2 h, the membranes were developed using the ECL reagent (Epizyme, Shanghai, China). The relative expression of SLC39A5, CYP2S1, FA2H, FABP1 and RND3 was determined using GAPDH.

Gene set enrichment analysis (GSEA)

According to the median values of SLC39A5, patients were categorized into high and low expression groups. GSEA was used to determine the enriched KEGG terms. Statistical significance was set at P < 0.05.

Statistical analyses

Gene expression in adjacent and CD tissues was analyzed using a paired t-test. The relationship between SLC39A5 expression in Crohn’s lesions and the SES-CD scores of patients was analyzed using Pearson’s correlation analysis. |R| >0.3 and P<0.05, were set as cut-offs to consider significance. The relationship between SLC39A5 expression and clinical characteristics of CD patients, including year, smoking history, upper gastrointestinal symptoms, stenosis (or penetration), ESR↑ (CRP↑), and perianal region symptoms, was analyzed using Fisher’s exact test. P<0.05 was set as the significance threshold.

Results

WGCNA construction

The transcriptional profile (accession number: GSE83448) of intestinal mucosa biopsies was obtained from healthy controls and CD patients, and corresponding information of clinical characteristics was obtained from the GEO database and used to construct WGCNA. The hclust function showed that there were no outliers (height >200) among these samples (Figure 1). Similarly, the clinical information (CD disease state: yes or no; inflammation: with or without) was thorough. Therefore, all the genes in these samples and their corresponding clinical information were used to perform WGCNA. We found that it ensured a scale-free topology >0.85 (Figure 2(A)) and mean connectivity close to 0 (Figure 2(B)); the soft threshold (β value) was set as 14. Furthermore, under β = 14, the topological overlap matrix met the scale-free topology criterion, with R² = 0.85 (Figure 2(C)). Therefore, we used β = 14 as a suitable soft threshold to construct WGCNA, and 15 co-expressed gene modules (darkturquoise, blue, brown, grey60, greenyellow, lightcyan, darkred, royalblue, turquoise, lightyellow, red, salmon, darkgreen, black, and magenta) were clustered (Figure 2(D)). In particular, when the genes had no co-expression relationship with each other they were clustered in the gray module.

Figure 1. Sample tree clustering of transcriptional profile of 53 intestinal mucosa biopsies and corresponding and clinical traits. Disease state: black indicates CD patients, and white indicates patients without CD. Inflammation: black indicates the tissues with inflammation, and white indicates the tissues without inflammation.

Figure 2. Construction of WGCNA. (A,B) Scale independence and mean connectivity of various soft threshold values (β). β with the scale independence >0.85 and mean connectivity closed to 0 was set as suitable. (C) Gene sets with the corresponding log10 and log10 P-values when the scale-free topology is set as β = 14. (D) Clustering dendrograms of all genes with dissimilarities based on topological overlap, together with their assigned module colors.

Identification of clinically significant modules

We then analyzed the relationship between gene modules and clinical characteristics. We found that the darkred (R = 0.68, P < 0.001; R = 0.42, P = 0.002), red (R = 0.36, P = 0.007; R = 0.35, P = 0.01), and royalblue (R = 0.62, P < 0.001; R = 0.44, P < 0.001) modules were positively and significantly associated with CD disease state and inflammation, whereas the darkgreen module (R = −0.7, P < 0.001; R = −0.51, P < 0.001) was negatively and significantly associated with disease state and inflammation (Figure 3(A)). We found that the MM in these gene modules was significantly associated with GS (Figure 3(B)). These four gene modules (darkred, red, royalblue, and darkgreen) were considered clinically significant gene modules. Therefore, 2752 genes (Supplementary Table 1) in these modules were used in further analysis.

Figure 3. Identification of clinical modules. (A) Heatmap showing the relationship between gene modules and clinical traits. (B) The relationship between module memberships of clinically significant modules and gene significance of inflammation and disease state are shown.

Identification of potential hub genes

DEG analysis was performed; a total of 489 DEGs with |LogFC| ≥ 1 were identified between the intestinal mucosa biopsies provided by healthy controls and those from CD patients (Figure 4(A); Supplementary Table 2); and 209 DEGs with |LogFC| ≥ 1 were identified between intestinal mucosa biopsies with or without inflammation from CD patients (Figure 4(B); Supplementary Table 3). After performing intersection analysis, 5 genes (CDHR5, SLC39A5, CYP2S1, TMEM61, and SCIN) which were in the clinically significant modules and differently expressed among intestinal mucosa biopsies from healthy controls, intestinal mucosa biopsies with inflammation from CD patients, and intestinal mucosa biopsies without inflammation from CD patients, were identified as hub genes (Figure 4(C)). Interestingly, these five genes were also downregulated in CD tissues, especially in the CD tissues with inflammation (Table 1). Furthermore, by performing the PPI network analysis, we found that there were relationships between CDHR5, CYP2S1, and SLC39A5 (Figure 4(D)).

Figure 4. Identification of differently expressed genes in intestinal mucosa biopsies obtained from healthy controls and CD patients with or without inflammation. (A) Volcano plot showing the differently expressed genes in the intestinal mucosa biopsies obtained from healthy controls and CD patients with the threshold as LogFC ≥ 1. (B) Volcano plot showing the differently expressed genes in the intestinal mucosa biopsies with or without inflammation obtained from CD patients with the threshold as LogFC ≥ 1. (C) Intersection analysis was performed to identify differently expressed genes which were also clustered in the clinical modules. List 1: differently expressed genes in the intestinal mucosa biopsies obtained from healthy controls and CD patients. List 2: the differently expressed genes in the intestinal mucosa biopsies with or without inflammation obtained from CD patients. List 3: Genes in the clinical modules. (D) PPI network of SCIN, TMEM61, SLC39A5, CYP2S1, and CDHR5.

Table 1. The change level of SCIN, TMEM61, SLC39A5, CYP2S1 and CDHR5 among intestinal mucosa biopsies from healthy control, intestinal mucosa biopsies with inflammation from CD patients and intestinal mucosa biopsies without inflammation based on the data from GSE83448.

Gene	LogFC (Crohn’s disease vs. normal)	LogFC (inflammation vs. non-inflammation)
SCIN	−2.388	−1.191
TMEM61	−1.352	−1.299
SLC39A5	−2.297	−7.839
CYP2S1	−2.763	−3.374
CDHR5	−1.433	−1.075

Verification of hub genes in the GSE208303 cohort

To verify the effectiveness of the five hub genes in CD, the GSE208303 cohort containing nine intestinal mucosa biopsies obtained from healthy controls and 26 intestinal mucosa biopsies obtained from CD patients was used as the verification cohort. Consistent with the results from the discovery cohort GSE83448, SLC39A5 expression was also reduced in intestinal mucosa biopsies obtained from patients in the GSE208303 cohort. However, CDHR5, CYP2S1, TMEM61, and SCIN expression levels did not change significantly in the GSE208303 cohort (Figure 5). Therefore, SLC39A5 was selected as the hub gene for CD.

Figure 5. SCIN, TMEM61, SLC39A5, CYP2S1, and CDHR5 expression in the intestinal mucosa biopsies obtained from healthy controls and CD patients based on the data from GSE208303.

Validation of clinical value of SLC39A5 in CD tissues from our cohort

We then determined the expression of the hub gene SLC39A5 in the 31 CD tissues from our cohort via qRT-PCR. The mRNA levels of SLC39A5 were significantly lower in CD tissues than those in adjacent mucosal tissues (Figure 6(A)). Matched-pair analysis revealed that the mRNA levels of SLC39A5 were significantly reduced in 48.4% of CD tissues compared to those in their corresponding adjacent tissues (Figure 6(B,C)). Furthermore, by reviewing the clinical traits of the 31 CD patients (Table 2), we found that mRNA levels of SLC39A5 were negatively associated with CD-simplifying endoscopic score (SES-CD; R = −0.51; Figure 6(D)). Moreover, we found that the protein levels of SLC39A5 were also lower in CD tissues compared with those in the corresponding adjacent tissues (Figure 6(E)). Furthermore, using Fisher’s exact test, we found that low SLC39A5 mRNA levels were associated with stenosis (Table 3). This evidence indicated that SLC39A5 may act as a biomarker for CD associated with the severity of CD.

Figure 6. Clinical value of SLC39A5. (A) SLC39A5 expression in the 31 pair CD tissues and adjacent mucosa biopsies. (B,C) Matched-pair analysis for SLC39A5 expression between CD tissues and corresponding adjacent mucosa biopsies. (D) Pearson analysis for the relationship between SLC39A5 mRNA levels in CD tissues and SES-CD score. (E) Representative protein levels of SLC39A5 in 5 pair CD tissues and adjacent tissues determined by western blot. **P < 0.05.

Table 2. The basic clinical traits of 31 CD patients who provided intestinal mucosa biopsies and expression of SLC39A5 in adjacent tissues and CD tissues.

ID	Year	Smoking history	Upper gastrointestinal symptoms	Stenosis (or penetrating)	CRP↑	ESR↑	Perianal region symptoms	SES-CD	SLC39A5 expression in adjacent tissues	SLC39A5 expression in CD tissues
1	80	No	No	Yes	Yes	Yes	No	18	1	0.582366793
2	28	No	No	Yes	Yes	Yes	No	15	1.109569472	0.450625231
3	24	No	Yes	Yes	Yes	Yes	Yes	21	0.456915725	0.574349177
4	25	No	No	No	Yes	Yes	No	8	1.028113827	1.071773463
5	18	No	No	No	Yes	Yes	No	10	0.432268616	1.283425898
6	17	Yes	No	No	Yes	Yes	No	8	1.602139755	0.864537231
7	17	No	Yes	Yes	Yes	Yes	No	20	1.310393404	0.489710149
8	44	No	No	No	No	No	No	3	1.424050196	0.882702996
9	49	No	No	Yes	No	No	No	12	1.494849249	0.574349177
10	14	No	No	Yes	No	No	No	11	1.2397077	0.514056913
11	19	No	No	No	Yes	Yes	No	14	1.453972517	0.602903914
12	23	No	No	Yes	Yes	Yes	No	18	1.474269217	0.473028823
13	30	No	No	No	No	No	No	10	1.301341855	0.426317446
14	19	No	No	No	Yes	Yes	No	11	0.812252396	0.574349177
15	44	Yes	No	No	Yes	Yes	No	16	1.347233577	0.386891248
16	29	No	No	Yes	Yes	Yes	No	18	0.993092495	0.469761375
17	41	No	No	No	No	No	Yes	18	0.489710149	0.50697974
18	15	No	No	No	No	No	Yes	8	0.476318999	0.870550563
19	35	No	No	No	Yes	Yes	No	6	0.687770909	0.888842681
20	24	No	No	No	No	No	No	8	1.197478705	0.846745312
21	75	No	No	Yes	Yes	Yes	No	17	1.310393404	0.607097442
22	33	No	No	Yes	Yes	Yes	No	12	1.580082624	0.737134609
23	38	No	No	No	Yes	Yes	No	15	0.840896415	0.628506687
24	19	No	No	No	Yes	Yes	No	21	0.742261785	0.316439148
25	31	No	No	No	No	No	No	3	0.907519155	1.064370182
26	47	No	No	No	Yes	Yes	No	15	0.655196702	0.939522749
27	34	No	No	Yes	Yes	Yes	No	19	0.716977624	0.343885455
28	13	No	No	Yes	Yes	Yes	Yes	22	1.474269217	0.721964598
29	22	No	No	No	No	No	No	6	1.443929196	0.61985385
30	34	No	No	Yes	Yes	Yes	No	24	1.140763716	0.888842681
31	26	No	No	Yes	Yes	Yes	Yes	22	0.309926925	0.650670928

Table 3. Fisher’s exact test was used to analyze the relationship between SLC39A5 expression in CD tissues and clinical traits of 31 CD patients.

Characteristics		Low expression	High expression	x^2	P
Year	≥30	7	7	0.03	0.87
	<30	9	8
Smoking history	Yes	1	1	0.01	0.96
	No	15	14
Upper gastrointestinal symptoms	Yes	2	0	2.00	0.16
	No	14	15
Stenosis (or penetrating)	Yes	10	4	4.02	0.04
	No	6	11
ESR↑ (CRP↑)	Yes	12	10	0.26	0.61
	No	4	5
Perianal region symptoms	Yes	2	3	0.32	0.57
	No	14	12

GSEA for the pathways affected by SLC39A5 in CD

GSEA was performed to explore the molecular mechanisms of SLC39A5 in CD. Results indicated that SLC39A5 expression deficiency induced the activation of pathways of ‘ribosome,’ ‘JAK-STAT,’ ‘protein export,’ ‘taste transduction,’ and ‘basal transcription factors’ (Figure 7(A)), and suppressed the pathways of ‘nitrogen metabolism,’ ‘alanine aspartate and glutamate metabolism,’ ‘aldosterone regulated sodium reabsorption,’ ‘glycerophospholipid metabolism,’ and ‘arginine and proline metabolism’ in CD (Figure 7(B)).

Figure 7. GSEA for SLC39A5 in CD tissues. (A) Pathways activated in the CD tissues with low SLC39A5 expression. (B) Pathways inhibited in the CD tissues with low SLC39A5 expression.

Broadened bioinformatic threshold identified CYP2S1, FA2H, FABP1, and RND3 as biomarkers for CD

Finally, to obtain more biomarkers for CD, we broadened the threshold (LogFC ≥ 0.5) in the DEG analysis. In total, 1858 DEGs were identified between the intestinal mucosa biopsies of healthy controls and those from CD patients (Figure 8(A)); 552 DEGs were identified between the intestinal mucosa biopsies with or without inflammation from CD patients (Figure 8(B)). After performing intersection analysis, 74 genes were identified that were in the clinically significant modules and differentially expressed among intestinal mucosa biopsies from healthy controls and CD patients, with or without inflammation (Figure 8(C)). The 74 genes were then imported into the STRING tool and a PPI network was constructed (Figure 8(D)). By calculating the degree score of all genes, it was found that FA2H (degree score = 10), ERBB3 (degree score = 9), FABP1 (degree score = 9), CYP2S1 (degree score = 9), and RND3 (degree score = 9) had the top 5 degree scores (Figure 8(E)). We also detected their expression in CD and adjacent tissues using qRT-PCR. We found that the expression of CYP2S1 (Figure 8(F)), FA2H (Figure 8(G)), FABP1 (Figure 8(H)), and RND3 (Figure 8(I)) was decreased in CD tissues, whereas the expression of ERBB3 (Figure 8(J)) was not significantly changed. Furthermore, western blotting demonstrated that the expression of CYP2S1, FA2H, FABP1 and RND3 was also reduced in CD tissues (Figure 8(K)). Therefore, CYP2S1, FA2H, FABP1, and RND3 may also act as biomarkers of CD. However, it is worth noting that the amplitude of variation of CYP2S1, FA2H, FABP1, and RND3 was lower than SLC39A5.

Figure 8. CYP2S1, FA2H, FABP1, and RND3 were identified as biomarkers for CD while broadening the threshold of bioinformatics analysis. (A) Volcano plot showing the differently expressed genes in the intestinal mucosa biopsies obtained from healthy controls and CD patients, with the threshold as LogFC ≥ 0.5. (B) Volcano plot showing the differently expressed genes in the intestinal mucosa biopsies, with or without inflammation, obtained from CD patients, with the threshold as LogFC ≥ 0.5. (C) Intersection analysis was performed to identify differently expressed genes, which were also clustered in the clinical modules. List 1: Genes in the clinical modules. List 2: the differently expressed genes between the intestinal mucosa biopsies with or without inflammation from CD patients. List 3: differently expressed genes between the intestinal mucosa biopsies provided by healthy controls and those from CD patients. (D,E) Construction of PPI network of the 74 genes, and selection of hub genes in PPI network based on degree score. (F–J) qRT-PCR analysis was performed to analyze the expression of CYP2S1, FA2H, FABP1, RND3, and ERBB3 in CD tissues and adjacent tissues. (K) Western blotting analysis was performed to analyze the expression of CYP2S1, FA2H, FABP1, and RND3 in CD tissues and adjacent tissues. **P < 0.01; *P < 0.05; ns, not significant.

Discussion

The incidence and clinical characteristics of CD differ among regions and races. Although the etiology of CD has not been fully ascertained, genetics, environment, and inflammation are considered key factors in the development of CD. Early identification of the etiology of CD is important for diagnosis and treatment. Some researchers have shown that CRP level and ESR are valuable in the diagnosis of CD (Daniluk et al. 2019; Lee et al. 2020). However, their sensitivities and specificities were low. Therefore, there is an urgent need to explore biomarkers with high sensitivity and specificity to accurately diagnose CD.

Bioinformatics is an important tool for identifying disease mechanisms (Wooller et al. 2017). Here, by performing WGCNA in GSE83448, four clinically significant modules (darkgreen, darkred, red, and royalblue) associated with both CD disease state and inflammation were identified, and a total of 2752 genes were identified in these modules. Via DEG and intersection analyses, five genes were identified that were in the clinically significant modules and differentially expressed among intestinal mucosa biopsies from healthy controls, intestinal mucosa biopsies with inflammation from CD patients, and intestinal mucosa biopsies without inflammation from CD patients. Among them, SLC39A5 was reduced in CD tissues from the GSE208303 cohort and 31 CD tissues from our research group. The SES-CD score is a simple, reproducible, and easy-to-use endoscopic scoring system for CD that can be used to evaluate the activity and prognosis of CD. Patients with high SES-CD scores have higher CD activity and poorer prognosis (Daperno et al. 2004; Zheng et al. 2020). We found that SLC39A5 expression was negatively associated with SES-CD scores. In particular, SLC39A5 expression was negatively associated with stenosis (or penetration) in CD patients.

The zinc transporter SLC39A5 is a member of the ZIP family. It transports zinc into the cell and plays a key role in managing zinc levels inside the cell (Hennigar et al. 2016). As zinc is a key element for cells, SLC39A5 plays various roles in diseases. Wang et al. demonstrated that SLC39A5 is involved in diabetes progression by regulating insulin secretion in pancreatic β cells (Wang et al. 2019). Dong et al. demonstrated that SLC39A5-deficiency contributes to the progression of high myopia (Dong et al. 2021). Under physiological conditions, SLC39A5 protects the pancreas from zinc toxicity by regulating intestinal zinc excretion (Geiser et al. 2013). However, the role of SLC39A5 in CD is unclear. Here, unlike the pathogenic effects of the cognate family member SLC3A8 (Melia et al. 2019), we provided the first evidence that SLC39A5 may act as a biomarker associated with the severity of CD.

Furthermore, GSEA revealed that SLC39A5 regulated a series of pathways, including the JAK-STAT pathway. In addition to its significant role in cellular signaling, the JAK/STAT pathway is associated with several immune-mediated diseases (Banerjee et al. 2017; Hu et al. 2021). Abnormal activation of the JAK/STAT pathway exists in most CD cases (Xu et al. 2021), and a series of JAK/STAT pathway inhibitors have exhibited significant therapeutic effects in CD (Ma et al. 2019; Moon et al. 2021). The present study provides the first evidence that SLC39A5 may be linked to the activation of JAK/STAT and CD state. The restoration of SLC39A5 expression in mucosal tissues may be a novel strategy for treating CD.

Interestingly, broadening the threshold (LogFC ≥ 0.5) in DEG analysis and performing PPI network analysis identified five genes, CYP2S1, FA2H, FABP1, RND3, and ERBB3, that may act as biomarkers for CD. Via qRT-PCR analysis and western blotting, we found that the expression of CYP2S1, FA2H, FABP1, and RND3 was reduced in CD tissues. CYP2S1 encodes a cytochrome P450 monooxygenase that is involved in the metabolism of retinoids and eicosanoids. It catalyzes the breakdown of the inflammatory factor PGH2 to hydroxyheptadecatrienoic acid and malondialdehyde (Bui et al. 2011). In physiological and pathological states, CYP2S1 acts as a protective factor. For example, under hypoxia, cells upregulate CYP2S1 to avoid DNA damage (Cabrera-Cano et al. 2021). Dysregulated CYP2S1 expression has been observed in cancer, atherosclerosis, and inflammatory diseases (Saarikoski et al. 2005). FA2H encodes a protein that catalyzes the hydroxylation of free fatty acids at the C-2 position to produce 2-hydroxy fatty acids and is stereospecific for the production of (R)-2-hydroxy fatty acids (Hirao-Suzuki et al. 2020). Previous studies have found that loss of FA2H expression leads to the accumulation of intracellular polyunsaturated fatty acids, which may inhibit cell proliferation and migration (Yao et al. 2019). As 2-hydroxy fatty acids are key mediators for inflammation, FA2H also has the potential to regulate the inflammation response under physiological and pathological conditions (Guo et al. 2019). FABP1 encodes a protein that plays a role in the lipoprotein-mediated cholesterol uptake. It regulates intracellular lipid transport and has the potential to bind free fatty acids and their coenzyme A derivatives (Tsai et al. 2020). FABP1 overexpression promotes intestinal inflammation and tumorigenesis (Abbasi et al. 2020). Inhibition of FABP1 can ameliorate hepatic steatosis and inflammation (Mukai et al. 2017). RND3 encodes a GTPase superfamily protein. The RND3 protein can bind to GTP but has no GTPase activity. It negatively regulates cytoskeletal organization, causing adhesion loss (Dankel et al. 2019). RND3 can bind to the P65 protein and activate the NF-κB pathway, thus leading to the production of inflammatory factors (Dong et al. 2021). Our study provides the first evidence that CYP2S1, FA2H, FABP1, and RND3 may act as biomarkers of CD.

Although we have provided some biomarkers for the diagnosis of CD, the present study has some limitations. First, bioinformatics has limitations. For the DEG analysis, the results were correlated with the threshold. High threshold settings will cause more false negative results, whereas low thresholds will cause more false positive results. Second, the threshold of DEG analysis is associated with the rangeability of the identified genes. Genes identified using a high threshold are often accompanied by a large change in their expression. Therefore, we only focused on the role of genes and ignored the effects of non-coding RNA on CD. Third, more cases should be examined to verify the effectiveness of SLC39A5, CYP2S1, FA2H, FABP1, and RND3 in CD. Moreover, the effects and molecular mechanisms of the biomarkers identified in our study on intestinal mucosal cells should be determined by additional cell or animal experiments.

In conclusion, we provided the first evidence to demonstrate that SLC39A5, CYP2S1, FA2H, FABP1, and RND3 may be novel biomarkers for CD, and that SLC39A5 has diagnostic value in estimating the activity of CD and the degree of inflammation. SLC39A5, CYP2S1, FA2H, FABP1, and RND3 may serve as novel biomarkers and targets for the diagnosis and treatment of CD.

Availability of data and materials

Transcriptional profiles (accession numbers: GSE83448 and GSE208303) of intestinal mucosa biopsies obtained from healthy controls and CD patients and corresponding information of clinical characteristics were obtained from the Gene Expression Omnibus database (https://www.ncbi.nlm.nih.gov/gds). All original data and R scripts are shown at doi:10.6084/m9.figshare.20632272.v6.

Acknowledgements

We appreciate the staff who participate in the daily maintenance and updating of the GEO database. Human Ethics Committee of Zhujiang Hospital Affiliated to Southern Medical University approved and administrated the tissue collection (approval number 2021-KY-150-01). LM, WQ, and ZZ designed the experiments and performed the analyses. ZQ, LY, MZ, and MH collected the data and constructed the figures. All authors reviewed the manuscript and approved its submission.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This study was supported by the Guangdong Medical Science and Technology Research Fund Project [grant number B2022084].