Introduction

Functional dyspepsia (FD) is a prevalent gastroduodenal disorder, with a global prevalence of 8.4%.¹ Clinically, FD is often manifested by symptoms such as epigastric pain or burning, postprandial fullness, and early satiety.² Although FD does not pose a survival threat to patients, it is characterized by a natural history of recurrent and remittent episodes. Research indicates that nearly two-thirds of dyspeptic patients continue to experience symptoms related to dyspepsia over a 10-year follow-up period.³ The chronic and fluctuating nature of FD symptoms significantly impacts the quality of life and mental health of patients, resulting in economic losses for both individuals and society, thus presenting a considerable public healthcare challenge.^4–7

Regarding pathogenesis, FD, as the most prevalent functional gastrointestinal disorder, lacks organic lesions that can explain the symptoms but is closely associated with central and gastrointestinal alterations. Consequently, the Rome IV consensus redefined “functional gastrointestinal disorders” as “disorders of gut-brain interaction” (DGBI).⁸ Over the past two decades, visceral hypersensitivity has supplanted motility disorders as the central theme in the pathophysiology of DGBI, with a particular focus on immune activation and neuro-immune interactions.⁹ The mucosal immune function in disorders of gut-brain interaction (DGBI) is believed to integrate and regulate bidirectional signaling between the brain, nervous system, and gut. Specifically, abnormal mucosal immunity can directly activate brain regions associated with stress and anxiety perception. Simultaneously, stress responses originating from higher cortical areas may profoundly influence immune system activity within the mucosal tissue, directly driving gastrointestinal symptoms.¹⁰ Recent studies have identified that pathological changes in the duodenum, leading to duodeno-gastric reflex activation, are associated with gastrointestinal sensory-motor dysfunction in patients with FD. Moreover, immune-inflammatory responses in the duodenal mucosa have been implicated as playing a pivotal role in the pathogenesis of FD through disruption of the brain-gut bidirectional interaction. Researchers are increasingly focusing on the critical contribution of enhanced mucosal permeability and micro-inflammation in the duodenum to the development of FD.^11,12 Several studies have observed that low-grade duodenal inflammation leads to increased mucosal permeability, changes in tight junction proteins, and other pathological alterations resulting from mild inflammation.¹³ These pathological changes are similar to the mechanisms underlying duodenal inflammation and the development of duodenal eosinophilia. Moreover, related studies have shown that activated eosinophil counts in the duodenum of FD patients are significantly higher than those in the control group. Additionally, these studies have identified changes in the function and structure of submucosal neurons in the duodenum of FD patients, which are associated with the accumulation of eosinophils and mast cells near neuronal structures.^14,15 The immune-inflammatory response in the duodenum is believed to potentially transmit harmful stimuli via afferent nerves, leading to gastric motility dysfunction and hypersensitivity, which in turn causes delayed gastric emptying and dyspeptic symptoms in FD patients.¹⁶ Consequently, low-grade inflammation and immune responses in the duodenum are emerging as novel targets in FD pathogenesis.¹⁷

While the role of duodenal immune-inflammatory responses in FD pathogenesis has gained attention, a systematic characterization of immune infiltration in the duodenal tissue of FD patients remains lacking. Recent advancements in microarray technologies and bioinformatics have significantly accelerated progress in the life sciences. By utilizing high-throughput data from public databases and employing integrative approaches combining bioinformatics and machine learning, researchers can identify high-dimensional features within gene expression data can be identified. This methodology is pivotal for biomarker discovery, the construction of diagnostic models, and further exploration of disease mechanisms.^18,19

In this study, we conducted a Mendelian randomization analysis, utilizing summary-level data from genome-wide association studies (GWAS) of functional dyspepsia (FD) and expression quantitative trait loci (eQTL) data to investigate genes associated with disease expression levels and complex traits. The results indicate that genes exhibiting a strong causal relationship with FD in terms of biological function are primarily enriched in aspects in immune response pathways. We retrieved the original microarray data for functional dyspepsia (FD) from the Gene Expression Omnibus (GEO) database to conduct differential gene expression identification and functional enrichment analysis. Subsequently, we employed the CIBERSORT and XCELL algorithms for immune infiltration analysis, followed by weighted gene co-expression network analysis (WGCNA) to identify co-expressed gene modules significantly associated with characteristic immune cells based on immune infiltration levels. From the key gene modules, we identified core immune-related genes using the cytoHubba application in Cytoscape software. Further, we applied the Least Absolute Shrinkage and Selection Operator (LASSO), an unsupervised machine learning method, to screen for nine immune-related hub genes in FD patients and constructed an immune-related diagnostic model capable of distinguishing high-risk from low-risk FD patients. Additionally, we employed AutoGluon, a deep learning approach, to train and test predictive models for classification. Finally, we validated the expression differences of hub genes and characteristic immune cells in FD animal models and control groups. We believe that immune inflammation activation is a key pathological change in the duodenum of FD patients, and the hub genes associated with the duodenal immune microenvironment could serve as potential new biomarkers for FD. Therefore, this study aims to elucidate the role of immune inflammatory responses in the pathological changes of the duodenum in FD, characterize the immune infiltration features in the duodenum of FD patients, and identify central genes related to duodenal immune responses. Moreover, we combine traditional machine learning and deep learning methods to construct a diagnostic model for FD patients, offering new insights for the accurate diagnosis and treatment of FD based on immune-related biomarkers.

Methods

Data Sources

Expression quantitative trait loci (eQTL) data were generated by the GTEx consortium and obtained from http://www.gtexportal.org/.²⁰ This dataset served as exposure data, with cis-eQTLs filtered at p < 5×10⁻⁸ for subsequent Mendelian randomization (MR) analysis. We established cis-regulatory regions within a 10,000 kb range on either side of coding sequences and performed linkage disequilibrium clustering using R² < 0.001. Only SNPs with valid allele frequencies and F-statistics ≥ 10 were retained, while corresponding palindromic SNPs were removed to mitigate confounding and weak instrument bias. The FD dataset (finn-b-K11_FUNCDYSP) was accessed through the IEU OpenGWAS database (https://gwas.mrcieu.ac.uk/). This dataset, comprising individuals of European ancestry, included 4,376 FD patients and 189,695 controls, with 16,380,380 SNPs.

Gene expression profiles were sourced from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) using the search terms “functional dyspepsia”, “human genome”, and “duodenum”. Datasets meeting the following inclusion criteria were considered eligible for this study: a) Homo sapiens with specific examination of duodenal mucosa; b) expression analysis by microarray; c) The study design includes FD patients who were clearly diagnosed based on symptom presentation and whose gastrointestinal organic lesions were excluded through endoscopic examination, along with healthy controls (HC); and d) a minimum of 20 individuals per dataset, with at least 10 subjects per group. The dataset GSE169304 (platform file GPL20301) used in this analysis comprised duodenal mucosal gene expression data from 24 healthy individuals and 40 FD patients.

Mendelian Randomization Analysis

Mendelian randomization (MR) analysis was conducted using the “TwoSampleMR” package to assess the causal relationship between biomarkers and functional dyspepsia (FD). We initially identified single nucleotide polymorphisms (SNPs) significantly associated with the target exposure (gene expression levels) from GWAS data. These SNPs (p < 5×10^–8) were selected as instrumental variables (IVs), and SNPs with an F-statistic < 10 were excluded to avoid weak instrument bias. This approach was used to explore the causal effects between the exposure data and outcome data.²¹ We employed the inverse variance weighted (IVW) method to investigate the relationship between specific genes and FD. Additional sensitivity analyses were performed using MR Egger regression, weighted median, fixed-effects IVW, and random-effects IVW methods.^22,23 The following criteria were applied to filter the results: (1) genes with p < 0.05 in the IVW method were selected; (2) genes were further refined based on the consistency of the direction of MR analysis results (Odds Ratio values) across three different methods, with all OR values being either >1 or <1; (3) genes showing signs of pleiotropy with p > 0.05 were excluded. In the sensitivity analysis, Cochran’s Q test was used to assess SNP heterogeneity, with p > 0.05 indicating no significant heterogeneity in the analysis.²⁴ MR-Egger regression was applied to test for horizontal pleiotropy, with p > 0.05 in MR-Egger indicating no significant intercept, suggesting the absence of horizontal pleiotropy.²⁵

Given the potential for strict p-value correction to overlook biologically meaningful exposure factors, and considering the sample size and specific research objectives, no p-value correction was applied in this section to allow for a more comprehensive exploration of potential biological links. Based on the MR results, we prioritized genes with strong evidence supporting a causal relationship with FD.

DEG Data Preprocessing

Gene expression profile data were downloaded from the GEO database and preprocessed and normalized using the R software (version 4.4.1, https://www.r-project.org/). To ensure data integrity, probes without matching gene symbols were excluded during preprocessing. When multiple probes corresponded to the same gene, the average value was taken as the final expression value for that gene to improve the accuracy and consistency of expression values. Sample preprocessing and outlier detection were performed using the goodSamplesGenes function from the WGCNA package. Outliers were identified based on distance and connectivity, filtering out samples with a z-score > 1. A sample connectivity plot was then generated through connectivity analysis. Further data refinement was conducted, with the counts matrix retained for differential analysis, while the standardized TPM matrix was used for other analyses. The TPM matrix was log-transformed using log2(exp + 1), followed by filtering out genes with zero variance and samples with more than 80% missing values. Genes with more than 50% missing values across samples were also excluded. Missing values were imputed using the impute R package, employing the nearest-neighbor imputation method. Principal component analysis (PCA) was then performed on the processed data to generate a PCA plot, highlighting the primary differences between experimental groups.

Identification of Differentially Expressed Genes (DEGs)

Differential expression analysis was performed on the corrected dataset using the DESeq2 package in R. The thresholds for differential expression were set as adj.p < 0.05 and |log2(fold change) (log2FC)| > 1, with the Benjamini-Hochberg (BH) method applied for correction. Volcano plots and heatmaps for visualizing DEGs were generated using the ggplot2 and pheatmap R packages, respectively.²⁶

Functional Enrichment Analysis

Gene Ontology (GO) analysis^27,28 is a widely used method for large-scale functional enrichment studies, including biological process (BP), cellular component (CC), and molecular function (MF) categories. The Kyoto Encyclopedia of Genes and Genomes (KEGG)²⁹ is a comprehensive database containing information about genomes, biological pathways, diseases, and drugs. To analyze the biological functions and pathways of related genes, GO and KEGG analyses were performed using the clusterProfiler R package, with p < 0.05 and FDR (q-value) < 0.1 considered as the significance thresholds for enrichment results. Gene Set Enrichment Analysis (GSEA) was used to evaluate the distribution of predefined gene sets within the gene expression data sorted by phenotypic correlation, thereby assessing their contribution to the phenotype. This powerful tool helps interpret gene expression data.³⁰ To identify the most significant functional terms between FD and control samples, we used the ClusterProfiler R package for GSEA to assess the distribution of DEGs in a phenotypic correlation sequencing gene table to determine their contributions to the FD phenotype.³¹ The gene set “c2.cp.all.v2022.1.Hs.symbols.gmt” was sourced from the MSigDB database. A significance level of adj.p < 0.05 was applied, with p-value correction performed using the Benjamini-Hochberg (BH) method.

Protein-Protein Interaction (PPI) Network Construction

The differentially expressed genes identified above were imported into the STRING database (v12.0, https://string-db.org)³² for constructing the protein-protein interaction (PPI) network. The minimum required interaction score was set to 0.4 in STRING to ensure moderate-confidence interactions.

Immune Cell Infiltration Analysis

In this study, we first used the CIBERSORT algorithm³³ and xCell algorithm³⁴ to assess the degree of immune cell infiltration in the gene expression profile based on the GSE169304 dataset. The IOBR package³⁵ in R was utilized to evaluate the proportions of 22 immune cell types in all samples using the CIBERSORT algorithm. Additionally, the xCell algorithm, which uses RNA sequencing (RNA-seq) data, was applied to calculate the abundance of 64 immune cell types in the samples. These methods, utilizing deconvolution algorithms, integrate the advantages of gene enrichment analysis to evaluate immune cell enrichment. Box plots generated using the ggplot2 package were used to present the abundance and proportion of infiltrating immune cells in each sample. The Wilcoxon test was applied to assess the differences in immune cell functional enrichment between the normal and FD groups, with p < 0.05 considered statistically significant.

Weighted Gene Co-Expression Network Analysis (WGCNA)

Weighted Gene Co-Expression Network Analysis (WGCNA) is a systems biology method used to describe gene association patterns across different samples. This approach helps identify gene sets with highly correlated changes and, by revealing the relationships between modules within these gene sets and their correlation with external sample traits, it uncovers potential candidate biomarker genes or therapeutic targets.³⁶ In this study, the WGCNA package in R was used to construct a co-expression network to identify modules most correlated with immune cells in patients with FD.³⁷

The pickSoftThreshold function was used to evaluate and select an appropriate soft threshold power (β), following the scale-free topology criterion to ensure the network structure adhered to standard scale-free features. The gene correlation matrix was then transformed into an adjacency matrix, which was further converted into a topological overlap matrix (TOM) to construct the weighted gene co-expression network. In the module detection stage, hierarchical clustering combined with dynamic tree-cutting algorithms was used to identify gene modules with similar expression patterns, with each module containing at least 30 genes. Genes that could not be classified were assigned to the gray module. Finally, the correlation between module eigengenes (MEs) and target traits was calculated to identify disease-related modules. Core genes within the modules were determined by selecting genes with a gene significance (GS) > 0.5 and module membership (MM) > 0.8.

Hub Gene Selection and Construction and Validation of Lasso Diagnostic Model

In this study, we uploaded the genes from the magenta module, which showed the most significant association with FD and key immune cells in WGCNA, to the STRING database (https://string-db.org/) to construct potential PPI relationships. A combined score > 0.4 was considered statistically significant. The PPI network was visualized using Cytoscape 3.7.2 (https://cytoscape.org). The cytoHubba plugin in Cytoscape was used to assign values to each gene based on the maximal clique centrality (MCC) algorithm, identifying the top 20 most important hub genes. Subsequently, machine learning methods were employed for data analysis. LASSO logistic regression was performed using the glmnet R package, and the optimal penalty parameter (λ) was adjusted via leave-one-out cross-validation. This procedure allowed the selection of feature genes for FD-related immune diagnostic biomarkers. To validate the model’s efficacy, the pROC R package was used to plot the receiver operating characteristic (ROC) curve of the hub genes and calculate the area under the curve (AUC) to assess their diagnostic performance for FD.³⁸ The AUC of an ROC curve typically ranges from 0.5 to 1. AUC values closer to 1 indicate better diagnostic performance. An AUC between 0.5 and 0.7 indicates low accuracy, between 0.7 and 0.9 indicates moderate accuracy, and above 0.9 indicates high accuracy. The glmnet R package was used to develop the optimal model. A scoring system based on 9 hub genes was constructed to calculate each patient’s risk score. The risk score formula is as follows: The risk score is as follows: (−2.83)* (AMBP expression) + (−7.21)* (CCK expression) + (−2.06)* (CHGA expression) + 0.26* (CLU expression) + (−2.28)* (GCG expression) + 5.23 *(RBP4 expression)+(−0.46)*(SST expression)+(−0.28) *(TTR expression). The R package “rms” was used to construct a nomogram, visualizing the interrelationships among the predictor variables in the predictive model. The calibration curve was employed to assess the discrepancy between the predicted probabilities and the actual probabilities, evaluating the model’s goodness of fit. The clinical utility of the model was assessed through Decision Curve Analysis (DCA). Precision-Recall (PR) analysis was conducted on the data using the pROC package, and the results were visualized with ggplot2 to construct the PR curve, thereby evaluating the model’s classification performance under class-imbalanced conditions.

Construction of a Diagnostic Model for Hub Genes Based on the AutoGluon Framework

This study employed the AutoGluon framework (version 1.2) for model training. AutoGluon is an open-source automated machine learning (AutoML) tool designed to streamline the machine learning process and enhance model performance. It integrates a variety of popular machine learning algorithms, which are used as base models for training. In this study, we used the presets=“best_quality” parameter in AutoGluon, which automatically selected several models, including XGBoost, LightGBM, CatBoost, Random Forest, Extra Trees, K-Nearest Neighbors, and Neural Networks. The dataset was split into 80% for training and 20% for testing. Given the small sample size of the training set (N = 38), we applied Leave-One-Out Cross-Validation (LOOCV) as the evaluation strategy for the training set. In each iteration, one sample from the training set was used as the validation set, while the remaining samples were used for training. The model was trained on the training subset and evaluated on the validation set. By setting num_bag_folds=None, we forced AutoGluon to use LOOCV for training and evaluating all models, ensuring consistency and fairness in model evaluation. After the LOOCV process, AutoGluon selected the model with the best cross-validation performance on the training set and retrained it using the entire training set (ie, 80% of the original data) to build the final predictive model. To obtain an unbiased estimate of the model’s generalization performance, we used the reserved testing set (20% of the original data) for the final performance evaluation. The model’s performance was evaluated using the following metrics: accuracy, ROC AUC, F1 score (the harmonic mean of precision and recall), recall, and precision. These metrics assess the model’s performance from multiple perspectives, providing a comprehensive evaluation of its effectiveness and robustness in practical applications.

Biological Function Analysis of Hub Genes

The GOSemSim R package was used for FRIENDS (Functional Relevance of Interconnected Gene Networks Derived from Semantic Similarity) analysis to investigate the importance of feature genes.³⁹ The chromosomal localization of feature genes was visualized using the RCircos R package to generate chromosome maps. To further explore the co-expression relationships between feature genes, other genes, and immune cells, psych R package was used for correlation analysis.

Based on the median model score, all samples were divided into high- and low-score groups for diagnostic prediction. To explore the biological functions of feature genes under different conditions, gene expression differential analysis and enrichment analysis were performed based on the model scores and high/low median score grouping. Gene expression differential analysis was carried out using the limma R package (Linear Models for Microarray and RNA-Seq Data), with empirical Bayes estimation to improve the robustness of the estimates. The p-values were adjusted for multiple hypothesis testing using the Benjamini-Hochberg method to control the false discovery rate (FDR), with genes having adj.p < 0.05 considered as differentially expressed. GO/KEGG enrichment analysis and GSEA were performed using the clusterProfiler package. The differentially expressed genes were then used for GO/KEGG analysis to identify associations with specific biological processes or signaling pathways. GSEA was conducted using all genes and their logFC values to assess the overall gene expression trend.

Construction of Regulatory Networks for Hub Genes

In this study, to explore post-transcriptional regulatory mechanisms, we constructed mRNA-transcription factor (TF), mRNA-miRNA, mRNA-RBP, and mRNA-drug regulatory networks. TFs regulate gene expression during post-transcriptional stages through interactions with hub genes. We retrieved relevant TFs from the ChIPBase v3.0 database⁴⁰ (http://rna.sysu.edu.cn/chipbase/)/hTFtarget database (http://bioinfo.life.hust.edu.cn/hTFtarget)⁴¹ to analyze their regulatory roles on hub genes. miRNAs play crucial roles in biological development and evolution by regulating multiple target genes, and a single target gene can be regulated by several miRNAs. To analyze the relationship between key genes and miRNAs, we retrieved relevant miRNAs from the ENCORI database (version 3.0)⁴² (https://starbase.sysu.edu.cn/). RNA-binding proteins (RBPs) control mRNA stability and translation efficiency during post-transcriptional stages by binding to mRNA. We also retrieved RBP data from the ENCORI database (version 3.0) to analyze the impact of RBPs on hub genes. The Comparative Toxicogenomics Database ⁴³ (http://ctdbase.org/) was used to obtain the targeting relationships between mRNA and drugs, constructing the mRNA-Drug network to analyze the potential regulatory effects of drugs on mRNA and their expression products, providing theoretical support for clinical drug selection. Finally, we used Cytoscape software to visualize the aforementioned regulatory network.⁴⁴

Establishment of the Rat FD Model

Healthy male SD rats (5 days old) were purchased from Beijing SPF Biotechnology Co., Ltd. Animal procedures were conducted in accordance with the ARRIVE guidelines and the Declaration of Basel. All animal experiments were approved by the Laboratory Animal Center Institute of Basic Theories of Chinese Medicine, China Academy of Chinese Medical Sciences (Approval No. IBTCMCACMS21-2404-02). The SD pups were randomly divided into two groups: the Normal Group and the Model Group (n=5). The FD animal model was induced using the widely recognized iodoacetamide gavage method, with the following protocol: Starting at postnatal day 10, rats in the Model Group were administered a mixture of 0.1% iodoacetamide and 2% sucrose solution (0.2 mL per gavage). Rats in the Control Group were given an equal volume of 2% sucrose solution once daily for 6 consecutive days.⁴⁵ The success of model induction was evaluated through gastric tissue pathological examination, rat body weight, gastric emptying, and intestinal propulsion rate. At the 8th week, rats were euthanized, and gastric antrum and duodenal tissues were collected for subsequent Hematoxylin and Eosin (H&E) staining, immunohistochemistry, and qPCR experiments.

Gastric Emptying Rate and Intestinal Propulsion Rate Assessment

Three rats from each group were randomly selected to assess gastric emptying rate and intestinal propulsion rate. The rats were fasted for 24 hours with free access to water. The methods for evaluating gastric emptying and intestinal propulsion are as follows: A mixture was prepared using 10 g sodium carboxymethyl cellulose, 250 mL distilled water, 2 g activated charcoal powder, 8 g sugar, 16 g milk powder, and 8 g starch, which was thoroughly mixed to form a 300 mL, approximately 300 g, black semi-solid paste. Each rat was gavaged with 3 mL of this mixture. Fifty minutes after gavage, the rats were anesthetized with 7% chloral hydrate and euthanized. The entire stomach was removed, weighed, and the gastric contents were washed out to determine the net weight of the stomach. Gastric emptying rate (%) was calculated as follows: Gastric Emptying Rate (%)= [1-(Stomach Total Weight-Net Weight)/ Gavage Volume] x 100%. The abdominal cavity was then opened, and the segment of the intestine from the pyloric section to the ileocecal junction was excised. This segment was placed on a clean glass plate without traction. The distance from the pylorus to the furthest point of charcoal progression was recorded, along with the total length of the intestinal segment. Intestinal propulsion rate (%) was calculated as follows: Intestinal propulsion rate (%) = Charcoal Propulsion Distance/Total Length of Small Intestine x 100%.

Hematoxylin and Eosin Staining (H&E Staining)

The gastric antrum tissue from the rats was fixed in 4% paraformaldehyde for 48 hours, followed by paraffin embedding. Paraffin-embedded tissues were sectioned into 4 µm slices and stained with Hematoxylin and Eosin (H&E) following standard procedures.

Quantitative Polymerase Chain Reaction (qPCR)

After collection, duodenal samples were immediately snap-frozen in liquid nitrogen and stored at −80°C. Total RNA was extracted from the duodenal tissue using the RNAExpress Total RNA Kit (G3013, Servicebio, China). cDNA was synthesized from RNA using the SweScript All-in-One RT SuperMix for qPCR (One-Step gDNA Remover) (Servicebio, China). Quantitative PCR (qPCR) reactions were performed with 2×Universal Blue SYBR Green qPCR Master Mix (Servicebio, China), and PCR amplification was conducted using the Bio-Rad CFX Connect Real-Time PCR System (Bio-Rad Laboratories, USA). The relative expression levels of RNA were quantified using the 2^−ΔΔCT method. GAPDH was used as the reference gene for quantitative analysis. The primer sequences for RT-qPCR are listed in Supplementary Table S1.

Immunohistochemistry (IHC) Staining

After collection, duodenal tissues were immediately fixed in 4% paraformaldehyde. Following paraffin embedding, tissue sections were subjected to IHC staining. We used CD138 antibody (Servicebio, GB115052, 1:150) and FOXP3 antibody (Servicebio, GB112325, 1:500) for IHC staining. The primary antibodies were incubated overnight at 4°C, and the secondary antibody (Servicebio, GB23302, 1:200) was incubated with the slides at room temperature (25°C) for 50 minutes. The slides were then stained using the DAB detection kit and counterstained with hematoxylin. Based on the scanning results of immunohistochemistry, the positive area-related parameters were analyzed using the HALO platform (Indica Labs, U.S.A). The positive staining levels were classified as follows: negative, no staining, 0 points; weak positive, light yellow, 1 point; moderate positive, brownish-yellow, 2 points; strong positive, brown, 3 points. We then calculated the average optical density and the percentage of the positive area. The staining intensity was evaluated using the histochemical scoring (H-score) method for each slide. The H-score is calculated as follows: H-Score (∑ (pi×i) = (percentage of weak intensity cells ×1)+(percentage of moderate intensity cells ×2)+(percentage of strong intensity cells ×3).⁴⁶

Results

Mendelian Randomization Study

We obtained a total of 19,942 genes and 41,157 strongly correlated SNPs (F-statistic ≥ 10, detailed in Supplementary Table S2) from the GTEx consortium. From the FD GWAS dataset, we retrieved summary results for 22,306 SNPs associated with FD and performed Mendelian Randomization (MR) analysis using instrument variables (SNPs) derived from exposure data. We applied the Inverse Variance Weighted (IVW) method, and then selected the appropriate MR method based on the heterogeneity test results (p-value of the Q statistic). The final results were retained for those with consistent effect directions and significant associations (p < 0.05) according to the IVW method. A total of 259 genes were identified as causally associated with the onset and development of FD, including 133 low-risk genes and 126 high-risk genes (Supplementary Table S3). High/low-risk genes refer to those whose increased expression is associated with either an increased or decreased likelihood of FD occurrence.

After inputting the genes causally associated with FD from the MR analysis into the PPI network, the results showed that these genes primarily function in the Biological Process category (Gene Ontology) related to defense responses. Reactome Pathways analysis further highlighted that these genes were associated with the immune system, particularly the Innate Immune System (Supplementary Figure 1A). Reactome Pathways enrichment analysis revealed that these causally significant genes in FD were mainly involved in immune-related functions, including neutrophil degranulation, immune modulation interactions between lymphocytes and non-lymphoid cells, and Toll-like receptor signaling (Supplementary Figure 1B).

Identification and Comprehensive Analysis of DEGs

Gene expression profiling data were first preprocessed and normalized. According to the results of PCA, there was a clear separation between the control and disease groups (Figure 1A). The sample connectivity plot revealed that most samples had Z-scores close to 0, indicating a high degree of connectivity among the samples (Figure 1B). Differential expression analysis was performed using the limma R package, and a total of 159 DEGs were identified. Among these, 90 genes were significantly upregulated, and 69 genes were significantly downregulated (Figure 1C). To further visualize the expression patterns of the top DEGs, a heatmap was constructed to display the top 20 upregulated genes and top 20 downregulated genes (Figure 1D). To comprehensively investigate the biological processes and pathways associated with these DEGs, GO and KEGG enrichment analyses were conducted (Supplementary Table S4). The most significant results were represented using a bubble chart (Figure 1E). The results indicated that the differentially expressed genes were primarily involved in various biological processes. In the BP assessment, DEGs were mostly engaged in detoxification of inorganic compound, stress response to copper ion, and other functions. DEGs changes associated with CC include immunoglobulin complex, circulating immunoglobulin complex, and multivesicular body. DEGs have been localized to the antigen binding, immunoglobulin receptor binding, aspartic-type endopeptidase activity in MF. In the KEGG pathway analysis, DEGs were enriched in pathways related to mineral absorption and ribosome biogenesis in eukaryotes. Additionally, GSEA was performed on the entire set of genes, and a total of 635 pathways were enriched using the Reactome gene set file (Supplementary Table S5). The six most significantly upregulated pathways included CD22-mediated BCR regulation, Fcgr activation, scavenging of heme from plasma, role of phospholipids in phagocytosis, role of LAT2 NTAL on calcium mobilization, and creation of C4 and C2 activators (Figure 1F). Using the KEGG gene set file, 109 pathways were enriched (Supplementary Table S5). Among these, the most significantly upregulated pathway was olfactory transduction, while the most significantly downregulated pathways included metabolism of xenobiotics by cytochrome P450, pentose and glucuronate interconversions, ascorbate and aldarate metabolism, drug metabolism by cytochrome P450, and renal cell carcinoma (Figure 1G).

Figure 1 Identification of DEGs between the FD and control groups. (A) PCA score plot of the training group dataset, showing the distinction between the control and FD groups. (B) Sample connectivity plot displaying the Z-scores of the samples. (C) Volcano plot of DEGs. Red dots represent upregulated genes, and blue dots represent downregulated genes. (D) Heatmap showing differential expression of the top 20 upregulated and top 20 downregulated genes. (E) Bubble plots for GO and KEGG pathway analyses of DEGs. (F) GSEA based on the Reactome database identified the top 6 upregulated and downregulated pathways in DEGs. (G) GSEA based on the KEGG database identified the top 6 upregulated and downregulated pathways in DEGs.

Immune Cell Infiltration Analysis

Using the GSE169304 dataset, we employed the XCELL and CIBERSORT algorithms to assess the immune cell infiltration in different groups within the dataset (Figures 2A and 3C). We observed significant differences in specific immune cell subtypes between FD and control samples. Specifically, based on the XCELL algorithm, the immune cells that were significantly increased in FD samples compared to the control group included astrocytes, B cells, cDC, class-switched memory B cells, CMP, endothelial cells, GMP, HSC, M1 macrophages, pDC, plasma cells, and Th2 cells. On the other hand, Tregs showed a significant decrease (Figure 2A and B). Using the CIBERSORT algorithm, we observed that plasma cells were significantly elevated and Tregs were significantly decreased in FD samples compared to the control group. Analysis of these immune cells using both CIBERSORT and XCELL algorithms indicated that the relationship between Tregs and plasma cells was consistent across both algorithms, suggesting that these immune cells play a key role in the differences observed between FD and control groups (Figure 2C and D). (Figure 3A and B) show the correlation of immune cells across different samples based on the XCELL and CIBERSORT algorithm. Based on the XCELL algorithm, it was observed that the cell types positively correlated with plasma cells were Preadipocytes (r = 0.68) and B cells (r = 0.79), while the cell types negatively correlated with plasma cells were NKT cells (r = −0.42) and Treg cells (r = −0.39). The cell types positively correlated with Treg cells were Basophils (r = 0.39) and NKT cells (r = 0.37), while those negatively correlated with Treg cells were GMP cells (r = −0.41) and plasma cells (r = −0.39). Furthermore, the CIBERSORT algorithm identified that the cell types positively correlated with plasma cells were Mast cells resting (r = 0.31) and B cells (r = 0.79), while those negatively correlated with plasma cells were CD4 memory resting T cells (r = −0.55) and Monocytes (r = −0.42). The cell types positively correlated with Treg cells included B cells naive (r = 0.34) and T follicular helper cells (r = 0.26), while plasma cells (r = −0.25) were negatively correlated with Treg cells.

Figure 2 Immune infiltration landscape between the control and FD groups. (A) Percentage of immune cell infiltration between the control and FD groups based on the XCELL algorithm. (B) Differential immune cell infiltration between the control and FD groups based on the XCELL algorithm. (C) Percentage of immune cell infiltration between the control and FD groups based on the CIBERSORT algorithm. (D) Differential immune cell infiltration between the control and FD groups based on the CIBERSORT algorithm. Green represents the control group, and red represents the FD group. *p <0.05; **p <0.01; ***p <0.001. Abbreviation: ns, not significant.

Figure 3 The correlated heatmaps representing the correlation between different immune cells based on XCELL (A) and CIBERORT (B) algorithms.

Identification of Immune-Related Genes

We performed WGCNA analysis and identified co-expression gene modules that were significantly associated with specific immune cells based on immune infiltration levels. A scale-free co-expression network was constructed, and the sample clustering results showed no outliers in the dataset (Figure 4A). The soft threshold power was set at 20 (scale-free topology fitting index R² = 0.90), and the cut-off height was set at 0.25 (Figure 4B). The hierarchical clustering dendrogram is shown in Figure 4C. WGCNA analysis identified a total of 6 modules (Figure 4D). Plasma cells showed significant positive correlations with the magenta module (cor = 0.38), black module (cor = 0.36), pink module (cor = 0.35), and tan module (cor = 0.32). In contrast, Tregs were significantly negatively correlated with the pink module (cor = −0.45), grey module (cor = −0.36), and tan module (cor = −0.28). Figure 4E presents the strong negative relationship between module membership in the pink module and gene significance for Tregs, while Figure 4F illustrates the strong positive relationship between module membership in the magenta module and gene significance for plasma cells.

Figure 4 Identification of co-expression modules based on WGCNA. (A) The samples were clustered to construct the sample clustering dendrogram to observe whether there are outliers. (B) The scale-free topology fit index (R 2) for soft threshold power (β) from 1 to 30. (C) Gene dendrogram and module colors; each color represents a module. (D) Correlation of six modules and two traits (plasma cells and Tregs). (E) The relationship of module membership and gene significance in MEpink module for Tregs. (F) The relationship of module membership and gene significance in MEmagenta module for plasma cells.

Hub Genes Selection and Lasso Diagnostic Model Construction

To identify the hub genes associated with immune processes in FD duodenum, we selected the magenta module, which exhibited the most significant positive correlation with the characteristic immune cell plasma cells (cor = 0.3, p = 1e-04). First, we used the Cytoscape software with the cytoHubba plugin to filter the top 20 most important DEGs based on the MCC (Maximal Clique Centrality) score and constructed a protein-protein interaction (PPI) network (Figure 5A). The extracted DEGs included GCG, SST, CCK, GAST, GHRL, GLP1R, CCKAR, CLU, VTN, AMBP, TTR, GC, CHGA, SERPINA3, RBP4, SOX9, MUC6, TFF2, PGC, and MUC1.

Figure 5 Screening for immune-related hub genes. (A) Analysis of the top 20 key genes based on the MCC score. (B) Coefficient profile plot derived using the LASSO algorithm. (C) Misclassification error plot derived using the LASSO algorithm. (D) Receiver operating characteristic (ROC) curve for the training set, with the area under the curve (AUC) indicating the feasibility of the model. (E and F) Performance of hub genes in the diagnosis of FD. (G) Nomogram of the model. (H) DCA curves of the model. (I) Calibration curves of the model. (J) PR curve of the model.

Next, we applied the LASSO algorithm for machine learning to further analyze the key DEGs. The coefficient profile plots (Figure 5B) and Misclassification Error (Figure 5C) derived from the LASSO algorithm led to the identification of 9 hub genes associated with FD, which were AMBP, CHGA, GCG, SOX9, TTR, CCK, CLU, RBP4, and SST. These genes were considered potential diagnostic biomarkers for FD immune-related processes.

Subsequently, we performed Receiver Operating Characteristic (ROC) analysis to assess the diagnostic potential of our model. The results showed that the Area Under the Curve (AUC) value of the nomogram was 0.94 (Figure 5D), indicating excellent diagnostic accuracy. Additionally, the ROC curve analysis evaluated the AUC values of the hub genes to determine their sensitivity and specificity in the immune-related diagnosis of FD (Figure 5E and F).

To validate the model, we first constructed a nomogram using hub genes to predict the risk of FD onset (Figure 5G). Both the decision curve (Figure 5H) and the calibration curve (Figure 5I) demonstrated the excellent performance of the nomogram in predicting FD. The PR curve (Figure 5J) confirmed the high predictive accuracy of the hub genes.

Construction of the Hub Genes Diagnostic Model Based on the AutoGluon Framework

The AutoGluon framework was used to train a diagnostic model based on the previously selected hub genes. The classification performance of nine different classifiers was evaluated (Supplementary Figure 2). Model performance was primarily assessed using the F1 score. The top 20 models ranked by the average F1 score are shown in Figure 6. The best-performing models were CatBoost, KNeighborsUnif, KNeighborsDist, ExtraTreesEntr, RandomForestGini, RandomForestEntr, and ExtraTreesGini, all of which achieved an F1 score of 0.91, demonstrating superior performance and strong classification ability. These models were effective in distinguishing functional dyspepsia (FD) samples.

Figure 6 Model Performance Comparison. Notes: Models sorted by average F1 score, showing top 20 performers.

Additionally, models such as LightGBMLarge and XGBoost achieved an F1 score of 0.83, also indicating strong classification ability. In terms of ROC AUC (Area Under the Curve), the LightGBMLarge model showed the highest ROC AUC of 0.925, indicating excellent ability to discriminate between positive and negative samples and capturing patterns and features in the data effectively. Furthermore, the XGBoost model achieved a ROC AUC of 0.8375, also demonstrating strong classification ability.

Regarding accuracy, most models had an accuracy between 0.71 and 0.86, with CatBoost and KNeighborsUnif achieving the highest accuracy of 0.86, indicating good overall performance in the classification task. The precision of most models ranged from 0.71 to 0.83, suggesting a high proportion of correctly identified samples among those predicted as positive. All models had a recall of 1.0, meaning that no samples of functional dyspepsia were missed in the identification process (Figure 7). Based on these results, the LightGBMLarge and XGBoost models performed the best in the overall diagnostic task and should be prioritized for use in diagnostic applications for clinical practice.

Figure 7 The classification performance of nine different classifiers.

Biological Function Analysis of Hub Genes

We further conducted a correlation analysis to explore the relationship between core genes identified through immune infiltration analysis and key immune cells (Figure 8A–C). It was found that, based on the XCELL algorithm, CHGA and TTR were significantly negatively correlated with Treg cells. However, using the CIBERSORT algorithm, no significant correlation was observed between hub genes and either Treg or plasma cells. Upon classification into disease and control groups, it was revealed that in the control group, Treg cells showed a significant negative correlation with RBP4. In FD patients, plasma cells exhibited a significant positive correlation with CHGA, SOX9, and SST. FRIENDS analysis revealed the significance of the core genes and indicated that the GCG gene played a central role among the hub genes (Figure 8D). These results suggest that CHGA may serve as an important biomarker contributing to the immune infiltration response in FD patients. To further elucidate the chromosomal distribution of the aforementioned genes, we visualized the characteristic genes (Figure 8E).

Figure 8 Biological significance of characteristic genes. (A) The correlated heatmap of immune cells and hub genes based on XCELL algorithm. (B) The correlated heatmap of immune cells and hub genes based on CIBERSORT algorithm. (C) The correlated heatmap of the correlation between immune cells and hub genes in control group and FD group based on CIBERSORT algorithm. (D) FRIENDS analysis of the hub genes. (E) Chromosomal distribution of hub genes. (F) GO and KEGG analysis between the high-risk and low-risk groups. (G) GO and KEGG analysis of the high-risk group; the color represents the log|FC| of the pathway. (H) GSEA analysis between the high-risk and low-risk groups.

We then stratified the samples in the GSE169304 dataset into high-risk and low-risk groups based on the median score of the hub genes, followed by differential gene expression analysis to further explore the potential biological mechanisms of the hub genes. The GO and KEGG enrichment analyses revealed that high-risk scores were associated with biological functions and pathways, including the hydrogen peroxide catabolic process, digestion, hydrogen peroxide metabolic process, endocytic vesicle lumen, haptoglobin binding, oxygen carrier activity, African trypanosomiasis (Figure 8F). The chord diagram illustrates the enrichment of various biological functions and pathways for the different genes (Figure 8G). The GSEA results demonstrated that the hub genes were highly enriched in several biological pathways related to immune inflammation and cell adhesion, such as the dectin 2 family, integrin cell surface interactions, and integrin3 pathway (Figure 8H).

Construction of the Hub Genes Regulatory Network

We retrieved transcription factors (TFs) associated with hub genes from the ChIPBase database to construct an mRNA-TF regulatory network, which was then visualized using Cytoscape software (Figure 9A). A total of 81 mRNA-TF interaction pairs were identified, including 9 hub genes and 38 TFs, with detailed information provided in Supplementary Table S6. Next, miRNAs related to hub genes were obtained from the StarBase database to build the mRNA-miRNA regulatory network, which was similarly visualized using Cytoscape (Figure 9B). We observed 25 mRNA-miRNA interaction pairs, involving 2 hub genes and 15 miRNAs, with detailed information in Supplementary Table S7. Additionally, using the ENCORI database, we retrieved RNA-binding proteins (RBPs) associated with hub genes to construct the mRNA-RBP regulatory network (Figure 9C). The results revealed that the network contained 2 hub genes and 7 RBPs, as detailed in Supplementary Table S8. Finally, based on the Comparative Toxicogenomics Database, we examined the mRNA-drug interactions and constructed an mRNA-drug interaction network (Figure 9D). Our findings showed that 5 genes formed 4 distinct drug-gene interaction pairs, involving 9 drug targets, providing a foundation for subsequent research. Specific mRNA-drug relationships are detailed in Supplementary Table S9.

Figure 9 Hub gene regulatory networks. (A) mRNA-TF regulatory network. (B) mRNA-miRNA regulatory network. (C) mRNA-RBP regulatory network. (D) mRNA-drug regulatory network.

Immune Cell Infiltration and Hub Gene Validation

To investigate immune cell infiltration in the duodenum of FD rats and quantify the expression of hub genes in this region, we induced the FD rat model through oral administration of iodine acetamide. The results revealed that, compared to the control group, FD rats exhibited no significant organic changes in the gastric antrum, but displayed notable weight loss, reduced weight gain, and delayed gastrointestinal transit, indicating the successful establishment of the FD rat model (Figure 10A–D).

Figure 10 Immune Infiltration Cells and Core Gene Expression in the Duodenum of FD Rats. (A) HE staining of gastric antrum tissue in rats. Scale bars, 20 μm. (B) Changes in body weight across different experimental groups. (C) Comparison of gastric emptying rates between groups. (D) Comparison of small bowel transit rates between groups. (E) Distribution and quantification of plasma cells in duodenal tissue, as observed by CD138 antibody immunohistochemical staining. Scale bars, 20 μm. (F) H-score assessment of plasma cell infiltration in duodenal tissue following CD138 antibody immunohistochemical staining. (G) Distribution and quantification of Treg cells in duodenal tissue, as observed by FOXP3 antibody immunohistochemical staining. Scale bars, 20 μm. (H) H-score assessment of Treg cell infiltration in duodenal tissue following FOXP3 antibody immunohistochemical staining. (I–Q) Differential mRNA expression of 9 hub genes between the FD and control groups. *p <0.05; **p <0.01; ****p <0.0001. Abbreviation: ns, not significant.

Immunohistochemistry (IHC) was employed to examine the infiltration of two key immune cell types—plasma cells and Treg cells. IHC analysis revealed that plasma cells were significantly more abundant in the duodenal tissues of FD rats compared to controls (Figure 10E and F). In contrast to the immune infiltration results, no significant difference was observed in the expression of Treg cells in the duodenal tissue of FD rats relative to the control group, although a slight downward trend was noted (Figure 10G and H).

Furthermore, we validated the expression of hub genes in the duodenum of FD rats. The results showed that the mRNA levels of CLU, RBP4, SOX9, AMBP, CHGA, GCG, and TTR were significantly different between the control and FD groups, while the mRNA expression of CCK and SST did not exhibit significant changes (Figure 10I–Q).

Discussion

Recent studies have increasingly recognized the duodenum as a central pathogenic site in the pathogenesis of FD.⁴⁷ Anatomically positioned between the stomach, jejunum, and ileum, the duodenum plays a crucial role in the digestion and absorption of nutrients. The duodeno-gastric reflex has been implicated in gastrointestinal motility dysfunctions, and some studies suggest that the duodenum, via afferent neurons, transmits harmful stimuli from the mucosa, which in turn may lead to gastric motility abnormalities and hypersensitivity.⁴⁸ Thus, the duodenum is thought to integrate the duodeno-gastric axis and central-peripheral organ interactions in the manifestation of FD symptoms.¹¹ As a key player in the pathological processes of FD, the duodenum contributes to upper gastrointestinal symptoms in FD patients, primarily due to the stimulation of duodenal contents, mild inflammation, and increased mucosal permeability. Notably, the immune-inflammatory response within the duodenal mucosa is considered to be closely associated with typical gastrointestinal symptoms and plays a critical role in the brain-gut axis dysfunction underlying FD. Therapeutic strategies aimed at reducing local immune activation have been suggested to offer potential benefits; however, currently, only a limited number of immune-targeted treatments are available.^9–11

In this study, we identified 159 DEGs between FD and healthy controls. GO and KEGG pathway enrichment analyses revealed that these DEGs were primarily involved in processes related to inorganic compounds, immunoglobulin complexes, antigen binding, immunoglobulin receptor binding, mineral absorption, and ribosome biogenesis in eukaryotes. Furthermore, Mendelian randomization analysis confirmed the significant role of immune responses in the pathogenesis of FD.

Current research generally acknowledges that the immune-inflammatory response in FD predominantly involves duodenal mucosal infiltration by eosinophils and mast cells.¹⁷ However, our study has identified plasma cells as a hallmark immune cell type in FD patients, a finding confirmed through animal model experiments. Plasma cells, the terminally differentiated form of B lymphocytes, are responsible for the secretion of large amounts of antibodies. Approximately 80% of the body’s antibody-secreting plasma cells reside within the gastrointestinal tract. The role of plasma cells in intestinal inflammation is complex, as they can either promote or suppress inflammation depending on their antibody secretion profile, microbiome composition, and the local immune environment.⁴⁹

Currently, there is limited research on the role of plasma cells in FD; however, studies have confirmed that patients with epigastric pain syndrome (EPS) associated with FD exhibit increased mast cells and plasma cells expressing IgE in the duodenal mucosa. In patients with EPS, duodenal mucosal susceptibility to acid leads to epithelial barrier dysfunction, with a response to weak acid that facilitates the penetration of luminal pathogenic antigens deeper into the mucosal layer, triggering plasma cell-driven immune responses and subsequent mast cell infiltration.⁵⁰ Previous studies have demonstrated significant plasma cell infiltration in the intestines of patients with conditions such as celiac disease, ulcerative colitis (UC), and Crohn’s disease.^51,52 Specifically, research has shown a marked increase in the proportion of IgG⁺ plasma cells in UC patients compared to healthy controls, which is associated with interferon (IFN)-γ-mediated type 1 immune responses, potentially promoting the progression of inflammation in UC.⁵³ Furthermore, the increased presence of intestinal-homing plasmablasts in the circulation of UC patients is linked to disease activity and an increased risk of long-term complications, suggesting that this may be a predictor of disease outcomes. In addition, studies have indicated that, in patients with mild to moderate duodenal inflammation, the number of duodenal immunoglobulin A (IgA) plasma cells is higher than in controls, while those with severe duodenal inflammation show microscopic ulceration, villous loss, and a reduction in all types of plasma cells, despite an overall increase in cell and inflammatory cell counts.⁵⁴ Plasma cells are also thought to have a close interaction with the gut microbiota.⁵⁵ The secretory antibody IgA produced by plasma cells plays a crucial role in mediating direct communication between host cells and the microbiota, maintaining intestinal homeostasis, and regulating immune response signaling.^56–58 Plasma cells that secrete antibodies are believed to be induced by various pathogenic and symbiotic microorganisms, and UC patients have been shown to exhibit microbial-specific colon plasma cell responses that contribute to intestinal inflammation.⁵⁹ Similarly, FD patients also present with dysbiosis in the duodenum and are found to have small intestinal bacterial overgrowth, which is thought to be a significant factor in duodenal immune activation and mucosal barrier disruption. This may serve as an important stimulus for the increased plasma cell count observed in FD patients’ duodenum.It is evident that plasma cells play a pivotal role in intestinal immune-inflammatory responses. Some researchers have proposed B-cell-targeted therapeutic strategies, such as the use of anti-CD20 antibodies (eg, rituximab), which may be effective in certain inflammatory bowel disease (IBD) patients. These treatments can indirectly affect plasma cell quantity and function through B-cell depletion, thereby reducing antibody production and alleviating inflammation.⁶⁰ However, such therapies may also interfere with the anti-inflammatory functions of B-cells, complicating treatment outcomes.

Interestingly, our study also found a reduction in Treg cells in the duodenum of FD patients. Treg cells are known to play a role in suppressing immune responses and maintaining peripheral tolerance and immune homeostasis.⁶¹ Characterized by the expression of Foxp3, Treg cells produce TGF-β, IL-10, and IL-35, which are specialized inhibitory factors that suppress various immune reactions and inflammation.⁶² The reduction in Tregs in the duodenum of FD patients may indicate a weakened regulatory capacity of the immune system, potentially contributing to chronic gastrointestinal inflammation and immune dysregulation, which could be a key factor in the pathogenesis of FD. In the gut, Treg cells interact with the microbiota, and their number and homeostasis are regulated by symbiotic microorganisms. Tregs play a crucial role in maintaining the peaceful coexistence of gut microbiota and food antigens.⁶³ The dysbiosis of the duodenal microbiota in patients with FD has increasingly become a focal point of research. The interplay between duodenal microbiota dysregulation and Treg cells in FD patients warrants further investigation. Studies have shown that the Treg cell subsets, Helios and Gata3, which dominate in the small intestine, are relatively less dependent on the microbiota but more responsive to tissue stress. These Tregs differentiate and function through the signaling of the alarm cytokine interleukin (IL)-33.^64,65 Related research has confirmed that transient depletion of Treg cells can lead to long-term autoimmune gastritis, which is associated with autoantibodies against H⁺K⁺ATPase and intrinsic factor.⁶⁶ A reduction in Treg cells leads to decreased Helicobacter pylori colonization in mice and negatively correlates with the severity of intestinal metaplasia in gastritis and peptic ulcer groups.^67,68 Furthermore, studies have indicated that Setd2 is critical for maintaining thymus-derived Treg cells, which influence the development and progression of colitis and colorectal cancer through various mechanisms. Loss of Setd2 exacerbates intestinal inflammation, likely due to its epigenetic and transcriptional regulation of Treg cell subsets.⁶⁹ However, animal studies have shown no significant changes in Treg cell populations between normal and model groups in the duodenum, and the precise role and infiltration characteristics of Treg cells in FD require further exploration.

Based on immune infiltration data, we performed WGCNA to identify key immune-related modules associated with FD. Subsequently, we employed machine learning techniques to pinpoint nine critical genes: SOX9, RBP4, CLU, AMBP, CHGA, GCG, CCK, TTR, and SST. Based on the LASSO logistic regression model, we confirmed the excellent diagnostic accuracy of these nine genes and further developed a diagnostic model using the AutoGluon automated deep learning framework.

SOX9, a member of the SOX (SRY-related HMG-box) gene family, is an important transcription factor crucial for the development of multiple organs, including the skeleton, testes, heart, lungs, pancreas, intestines, and nervous system.⁷⁰ SOX9 is essential for Paneth cell differentiation; in its absence, Paneth cells fail to form, leading to an enlargement of intestinal crypts and a significant increase in crypt cell proliferation, with proliferating epithelial cells replacing Paneth cells.⁷¹ SOX9 is thought to regulate Paneth cell function and maintain the intestinal barrier. An increase in SOX9 may indicate enhanced stem cell activity or function in the duodenum, potentially linked to the microinflammation and repair processes in FD patients’ duodenal regions.⁷² RBP4 is a lipophilic vitamin A transport protein primarily secreted by the liver and adipose tissue. Its main function is to transport retinol (vitamin A) to peripheral tissues.⁷³ Elevated levels of RBP4 are associated with both systemic and local inflammatory states. RBP4 can activate immune cells, such as macrophages, to increase the secretion of pro-inflammatory cytokines, potentially influencing the intestinal inflammatory response.⁷⁴ Moreover, vitamin A is crucial for maintaining epithelial cell health and promoting immune function. A deficiency in vitamin A can impair vitamin A-dependent immunity, including the homing of B and T cells to the intestine and the production of immunoglobulin A, leading to epithelial dysfunction and affecting intestinal health.⁷⁵ CLU is a molecular chaperone protein involved in various biological processes, playing essential roles in cell adhesion, immune response, and cell survival.^76,77 Elevated levels of CLU in the duodenum may reflect enhanced local inflammation or immune response. CLU plays a critical role in maintaining protein homeostasis under cellular stress, inhibiting cell death pathways, and regulating pro-survival signaling and transcriptional networks.⁷⁸ Increased expression of CLU in response to certain pathophysiological changes may provide a protective response in FD patients’ duodenal epithelial cells.⁷⁹

AMBP is an important biomarker involved in immune response, anti-inflammation, and cell protection. Its two main components—α1-microglobulin and bikunin—both participate in regulating immune responses and counteracting protease activity.^80,81 Currently, no studies have investigated the expression of AMBP in FD patients. However, based on the anti-inflammatory properties of AMBP, we hypothesize that a reduction in AMBP levels may impair both intestinal immune function and the integrity of the intestinal barrier. CHGA is a protein primarily expressed in endocrine cells, such as chromaffin and neuroendocrine cells in the gastrointestinal tract, involved in multiple biological processes within the endocrine, neuroendocrine, cardiovascular, and immune systems, and plays a regulatory role in the gut and nervous system.⁸² Reduced levels of CHGA may reflect a decline in neuroendocrine function in the duodenum and other gastrointestinal sites. Additionally, studies have shown that downregulated CHGA is a prominent central gene in colorectal cancer, suggesting its role in digestive system diseases.⁸³ Research also suggests that CHGA regulates changes in macrophages and intestinal epithelial cells and promotes the development of intestinal inflammation.⁸⁴ Therefore, the role of CHGA in FD pathogenesis warrants further investigation. GCG is expressed by specific intestinal endocrine cells (L cells), pancreatic α-cells, and a set of discrete neurons in the solitary tract. GCG encodes several peptides, including GLP-1.⁸⁵ After food intake, GCG rapidly increases, which is associated with blood glucose regulation, insulin secretion, and gastric emptying speed. Reduced GCG levels in the duodenum may affect GLP-1 secretion, thereby impacting gastric emptying and blood glucose regulation, potentially increasing the risk of functional dyspepsia.⁸⁶ CCK is a neuropeptide that regulates digestion, glucose levels, neurotransmission, and memory. Studies have found that both gastric and duodenal vagal afferent fibers are sensitive to CCK, and the endogenous release of CCK may activate the vagal afferent pathway in the stomach, which is related to the endocrine action of CCK originating from the intestine, associated with feeding signals.⁸⁷ CCK also has neurotrophic and anti-inflammatory properties. CCK analogs and liraglutide have been shown to restore tight junction integrity and reduce colon inflammation in Parkinson’s disease mouse models, inhibiting the reduction of dopaminergic neurons and accumulation of α-syn oligomers in the colon.⁸⁸ SST is a neuropeptide playing a key role in regulating cortical circuits and cognitive function and is widely distributed throughout the gastrointestinal system. SST analogs have a long-lasting inhibitory effect on gastric acid, pancreatic enzyme, bicarbonate secretion, and bile flow. SST also inhibits stimulated intestinal secretion, the release of neuropeptides from the gut and pancreas.⁸⁹ TTR is a protein primarily synthesized in the liver and choroid plexus of the brain, whose main function is to transport thyroid hormones and retinol-binding protein (vitamin A).⁹⁰ Currently, there are no studies on the role of SST and TTR in FD, which warrants further exploration. It is apparent that the key genes involved in FD pathogenesis predominantly play roles in immune inflammation, neuroregulation, and retinol transport, which may provide insights into the pathogenesis of FD and potential therapeutic targets.

This study provides valuable insights into the immune-related mechanisms of FD, particularly highlighting the role of immune cell infiltration in the duodenum. Our findings suggest that changes in plasma cells and regulatory Treg cells are crucial in the pathogenesis of FD, with a significant increase in plasma cells and a reduction in Tregs observed in FD patients. Furthermore, our research elucidates the immune-related hub genes in the duodenum of FD patients and constructs a diagnostic model based on these hub genes. However, there are several limitations to this study. First, our study is retrospective, with most of the analyses and conclusions based on publicly available data and a limited sample size. Further validation through prospective studies is required. Second, although our findings were validated using a machine learning model, external validation using independent datasets was not performed, which introduces the possibility of bias. Third, while our study identified immune cell infiltration characteristics and related hub genes in the duodenum of FD patients, the findings were validated only in a limited animal model, necessitating further experiments to clarify the expression patterns of these immune cells and characteristic genes in the disease and their underlying mechanisms. Furthermore, future studies could further investigate the manifestations of mild inflammation in the duodenum of FD patients based on endoscopic examinations and histopathological findings, as well as explore the association with changes in relevant immune cells.