Identification of potential biomarkers of tryptophan metabolism in Kawasaki disease and exploration of potential mechanisms

Introduction

Kawasaki disease (KD) is an acute, systemic vasculitis that is common in infants and children, with the onset generally occurring below 5 years of age. KD manifests as constant fever and a spectrum of mucocutaneous inflammation, including rash and conjunctival injection, alterations of the extremities, and cervical lymphadenopathy. Without timely and appropriate treatment, coronary artery dilatation and aneurysms can occur, which are significant causes of acquired heart disease in children (1).

Epidemiological studies have revealed a significant geographical and ethnic variation in KD incidence, with much higher frequencies in East Asia compared with Europe and North America (2). In pediatric populations in China, KD frequency shows differences between ethnic and regional subgroups, as well as a clear upward time trend (3). Currently, clinicians rely primarily on clinical findings and routine laboratory tests for diagnosis. Early disease detection is not currently possible, because there are no highly specific markers. In addition, children who do not present with obvious symptoms could be misdiagnosed due to the lack of these markers. The standard treatment for KD is high-dose intravenous immunoglobulin (IVIG) combined with aspirin; however, 10–20% of patients are IVIG resistant and are at higher risk for coronary artery lesions (CALs) (4). These issues highlight the need for more sensitive and specific biomarkers, which will facilitate early diagnosis, improve risk stratification, and guide therapeutic strategies for KD.

Many inflammatory diseases are regulated through immunometabolism. The kynurenine pathway of tryptophan (Trp) metabolism is central to immune modulation. L-kynurenine (Kyn) and kynurenic acid (Kyna) metabolites exert immunoregulatory effects, in part, by activating the aryl hydrocarbon receptor (AhR), which affects innate and adaptive immune signaling. The body primarily breaks down Trp through the Kyn pathway, which yields Kyn and other metabolites, many of which regulate inflammation and the immune system (5,6). Cardiovascular events and arterial inflammatory injury have been associated with these metabolites (7). Alterations in Trp-derived metabolites are associated with higher mortality in patients with cardiovascular diseases (8,9).

Metabolomic analyses on patients with KD have revealed marked Trp and derivative dysregulation. Plasma Trp and indole-3-acetic acid levels were significantly lower, whereas Kyn and Kyna levels were markedly higher than healthy controls (10). In a mouse model of coronary arteritis induced by a Lactobacillus casei cell wall extract, interruptions to metabolic processes were further confirmed. Analysis of the published single-cell RNA sequencing dataset revealed that peripheral blood mononuclear cells (PBMCs) of patients with KD showed increased expression of the AhR compared with healthy children. These results implicate the Trp-kynurenine-AhR axis in KD pathogenesis and highlight this axis as a source of novel diagnostic and therapeutic targets.

Based on these findings, we used bioinformatics methods to identify Trp metabolism-related genes (TMRGs) in KD. Using publicly available transcriptomic data, machine-learning algorithms [least absolute shrinkage and selection operator (LASSO) regression and random forest models] were used to identify candidate genes. Their functional relevance was determined by constructing protein-protein interaction (PPI) networks, conducting Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses, and assessing immune cell infiltration patterns. The results provide insight into the relationships between TMRGs and the inflammatory and immune processes in KD, which will provide novel mechanistic insights. We present this article in accordance with the TRIPOD reporting checklist (available at https://tp.amegroups.com/article/view/10.21037/tp-2025-1-905/rc).

Methods

Data collection

Transcriptomic data were obtained from the Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/). The training set, GSE68004 (platform GPL10558), consists of peripheral blood samples from 76 patients with KD and 37 healthy controls. To validate these findings, the GSE73461 dataset (also on GPL10558) was used, which consists of samples from 78 patients with KD and 55 controls. In addition, a curated list of 40 TMRGs was retrieved from the Molecular Signatures Database (https://www.gsea-msigdb.org/gsea/msigdb/) for subsequent analyses.

Identification of differentially expressed genes (DEGs) and key module genes

A differential expression analysis was conducted on GSE68004 using the “limma” package (version 3.54.2) to obtain DEGs in patients with KD compared with controls. The genes were significantly differentially expressed at an adjusted P value (P.adj) <0.05 and a |log₂ fold change| >1. Volcano plots were used to visualize the global expression profile of DEGs using the “ggplot2” package (version 3.4.4). Heat maps for the top 10 upregulated and downregulated genes were generated using the “ComplexHeatmap” package (version 2.14.0).

The “GSVA” package (version 1.46.2) was used to calculate Trp metabolism (TM) scores for the KD and control groups, which were then compared using the Wilcoxon test with a P<0.001. A weighted gene co-expression network analysis (WGCNA; version 1.72.1) was conducted on the GSE68004 samples with TM scores treated as an external trait. All samples were initially clustered together to identify and eliminate outliers. A soft-thresholding power (β) was selected, which yielded a scale-free topology fit index (R²) >0.9, with mean connectivity close to zero. Using dynamic tree cutting with a minimum module size of 30, we grouped the genes into modules, and similar modules were merged later at a height (MEDissThres) of 0.4. Pearson correlation analysis was used to identify modules associated with TM scores with a |cor| >0.3 and P value <0.05. A correlation heat map for module-trait relationships was generated. For the key module, gene-trait correlations were plotted individually, and the genes were filtered based on |gene significance (GS)| >0.2 and |module membership (MM)| >0.25 to obtain a key module gene set.

Functional analysis of key genes

We intersected the DEGs with the key module genes to obtain a set of genes using a Venn diagram generated using the “ggVennDiagram” package, version 1.2.2. Functional enrichment analyses were performed using the “clusterProfiler” (version 4.7.1) and “org.Hs.eg.db” (version 3.16.0) packages based on the GO (http://geneontology.org/docs/ontology) and KEGG (https://www.genome.jp/kegg/) databases. The criteria for determining pathway enrichment based on KEGG and GO terms were P<0.05. The top 10 terms or pathways were visualized using “GOplot” (version 1.0.2). To examine protein-level interactions among the key genes, a PPI network was constructed using the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database (https://string-db.org/) with a confidence score threshold of 0.4 and visualized in Cytoscape (version 3.7.1). The network modules were identified by applying the MCODE plug-in with the following parameters: K-core =2, degree cutoff =2, max depth =100, and node score cutoff =0.2. Within the CytoHubba plug-in, candidate target genes were further prioritized through four centrality algorithms: degree, closeness, maximal clique centrality (MCC), and maximum neighborhood component (MNC). The top 20 genes from each algorithm were intersected, and the candidate genes were visualized using a Venn diagram generated by “ggVennDiagram” (version 1.2.2).

Identification of biomarkers

A LASSO regression analysis was conducted using the glmnet package (version 4.1-7), and genes that corresponded to the optimal lambda value were selected as characteristic genes 1. A support vector machine (SVM) was used to build a model with a linear kernel using the caret package (version 6.0-86) and interpreted using the DALEX package (version 2.4.3). An assessment of feature importance was also performed. The top 10 genes were selected as characteristic genes 2 based on the importance score. In addition, a random forest model was applied using the “randomForest” package (version 4.7-1.1). Gene importance was assessed using a mean decrease in accuracy and also evaluated by the mean decrease in node purity (IncNodePurity). The top 10 genes were extracted and cross-validated, from which the top 5 were retained as characteristic genes 3. “ggVennDiagram” (version 1.2.2) was then used to visualize the intersection of characteristic genes 1, 2, and 3 using a Venn diagram, which resulted in the final set of characteristic genes. The diagnostic performance of the individual gene was assessed in the GSE68004 and GSE73461 datasets using the area under the curve (AUC). Genes exhibiting an AUC of >0.7 were considered potential biomarkers. A Wilcoxon test was performed to determine the expression difference between patients with KD and controls (P<0.05).

Gene set enrichment and immune infiltration analyses

To determine the putative biological function of the identified biomarkers, Spearman correlation analysis for each biomarker was done using the “psych” package (version 2.3.6) of the GSE68004 dataset, and the correlation coefficients were ranked. Gene set enrichment analysis (GSEA) was performed using the (version 4.7.1.003) “clusterProfiler” package with P.adj <0.05 and |normalized enrichment score| >1 as a threshold for significance. Enriched pathways and functions were visualized using the “GseaVis” package (version 0.0.5; https://CRAN.R-project.org/package=GseaVis). The CIBERSORT algorithm was used to estimate the relative proportions of 22 immune cell types in KD and control samples, and “ggplot2” (v 3.4.4) was used for visualization. The differences in immune cell composition between the patients with KD and healthy controls were determined using the Wilcoxon rank-sum test (P<0.05), and boxplots were generated. A Spearman correlation analysis was used to assess the relationship of biomarker expression to the differential immune cell populations (P<0.05).

Identification of N6-methyladenosine (m⁶A) modification sites and chromosome localization for biomarkers

The sequence-based RNA adenosine methylation site predictor (http://www.cuilab.cn/sramp) was used to predict m⁶A modification sites in the candidate genes, and the positions of these sites within the RNA secondary structure were visualized. The predicted m⁶A sites were annotated with confidence levels defined as very high (≥99%), high (≥95%), moderate (≥90%), and low (≥85%). The RCircos package was used to predict the chromosomal localization of the putative biomarkers.

TFs and drugs that modulate biomarkers

MicroRNAs (miRNAs) and transcription factors (TFs) targeting the candidate biomarkers were predicted using the miRDB (http://www.mirdb.org) and KnockTF (http://www.licpathway.net/) databases, respectively, to further examine the regulatory mechanisms of these genes in KD. miRNA-messenger RNA (mRNA) and TF-mRNA interaction networks were constructed and visualized using Cytoscape (version 3.7.1). Additionally, the DrugBank database (https://go.drugbank.com/) was queried to identify drugs that may target the biomarkers, and a drug-mRNA interaction network was generated using Cytoscape (version 3.7.1).

Validation of biomarkers by reverse transcription quantitative polymerase chain reaction (RT-qPCR)

Peripheral blood samples were obtained from 10 individuals, including 5 patients with KD and 5 age-matched healthy controls, at The First Affiliated Hospital of Hunan Normal University. The procedures were performed in accordance with the Declaration of Helsinki and its subsequent amendments, and written informed consent was obtained from all participants’ guardians. The Ethics Committee of The First Affiliated Hospital of Hunan Normal University (No. [2024]-135). RNA was extracted using TRIzol reagent, and the concentration was measured using a NanoDrop 2000 spectrophotometer. Complementary DNA (cDNA) was synthesized using a cDNA synthesis kit from SweScript. RT-qPCR was carried out with the primers listed in Table S1. The 2^−ΔΔCT method was used to calculate relative gene expression values, which were normalized to glyceraldehyde-3-phosphate dehydrogenase.

Statistical analysis

Version 4.2.2 of R software was used for the statistical analyses. Continuous variables that fit a normal distribution were used to compare two groups with a Student’s t-test, whereas non-normal continuous variables were analyzed using the Mann-Whitney U test (Wilcoxon rank-sum test). All P values were two-sided. A P value <0.05 was considered significant. GraphPad Prism 9 was used to generate a graphical representation of the RT-qPCR results.

Results

Identification of DEGs related to KD

The workflow for this study is presented in Figure 1. We identified 906 DEGs in the GSE68004 dataset, including 739 upregulated and 167 downregulated genes (Figure 2A,2B). The TM scores for the KD samples were significantly lower compared with the healthy individuals (P<0.001; Figure 2C). Hierarchical clustering of the samples was performed before WGCNA analysis, followed by confirmation of the absence of any outlier samples (Figure S1A). We determined the power of soft-thresholding to create a network with a scale-free topology with a scale-free fit index (R²) >0.9 and mean connectivity near zero. This yielded a power of 6 (Figure 2D). Using dynamic tree-cutting criteria, 17 modules were generated and subsequently merged into 10 modules with a merge threshold (MEDissThres) of 0.4 (Figure 2E, Figure S1B,S1C). Module-trait correlations were visualized with a heat map (Figure 2F). The MEbrown module displayed the greatest negative correlation with TM scores (cor =0.36, P<0.05) and therefore, we chose it as the key module. Using the thresholds |GS| >0.2 and |MM| >0.25, 1,194 key module genes were extracted (Figure 2G). The concurrent analysis of 906 DEGs and 1,194 key module genes identified 234 key genes which were used for subsequent analyses (Figure 2H).

Figure 1 Study flowchart. DEG, differentially expressed gene; GO, Gene Ontology; GSEA, gene set enrichment analysis; miRNAs, microRNAs; mRNA, messenger RNA; PPI, protein-protein interaction; qRT-PCR, quantitative reverse transcription quantitative polymerase chain reaction; ROC, receiver operating characteristic; TF, transcription factor; TM, tryptophan metabolism; WGCNA, weighted gene co-expression network analysis.

Figure 2 Identification of differentially expressed genes and TM-related module genes in KD. (A) Volcano plot showing the upregulated (red) and downregulated (blue) DEGs. (B) Heatmap of the top 10 DEGs; red and blue denote high and low expression levels, respectively. (C) Cloud plot showing significantly lower TM scores in KD samples compared with healthy controls (***, P<0.001). (D) Determination of the optimal soft-thresholding power for WGCNA. (E) Cluster dendrogram of all samples with corresponding phenotype heatmap. (F) Module-trait correlation heatmap, highlighting 10 modules and their association with TM scores. (G) Scatter plot of MM versus GS for the key module, showing correlation and statistical significance. (H) Identification of key genes. DEG, differentially expressed gene; GS, gene significance; KD, Kawasaki disease; MM, module membership; TM, tryptophan metabolism; WGCNA, weighted gene co-expression network analysis.

Identification of candidate genes and validation of biomarker

GO functional annotation of these 234 genes revealed 374 biological process (BP) terms, 29 cellular component (CC) terms, and 11 molecular function (MF) terms (Figure 3A). In particular, key genes were enriched for BPs related to the positive regulation of cytokine production and cytokine-mediated signaling pathways. The CCs included secretory granule membrane and specific granule, whereas the MFs included immune receptor activity and IgG binding. Moreover, osteoclast differentiation, tuberculosis, leishmaniasis, and other KEGG pathways were examined (Figure 3B). PPI network analysis revealed interactions between key genes, and the highest connected genes were IL1B and TLR4 (Figure 3C). The PPI network was divided into two subnetworks using the MCODE plug-in. IL1B, TLR4, SPI1, FCGR3A, and FCER1G were present in subnetwork 1, whereas MYD88 and FCGR3B were in subnetwork 2 (Figure S2). The top 20 scoring genes from four centrality algorithms were intersected simultaneously to yield 14 candidate genes (Figure 3D).

Figure 3 Selection of candidate characteristic genes. (A) Bubble plot showing the results of GO enrichment analysis of the core genes; bubble size indicates the gene count, and color represents the adjusted P value. (B) Bar plot showing the top 20 KEGG pathways significantly enriched among the core genes. (C) PPI network of the core genes constructed using STRING, with node size reflecting connectivity. (D) Identification of candidate genes. (E) Selection of candidate feature genes 1 by LASSO logistic regression analysis. (F) Selection of candidate feature genes 2 by SVM. (G,H) Selection of candidate feature genes 3 by Random Forest and Cross-validation curve. (I) Venn diagram of feature genes. FC, fold change; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; LASSO, least absolute shrinkage and selection operator; MCC, maximal clique centrality; MNC, maximum neighborhood component; PPI, protein-protein interaction; STRING, Search Tool for the Retrieval of Interacting Genes/Proteins; SVM, support vector machine.

Figure 3E shows the results of LASSO logistic regression analysis on 14 candidate feature genes. At a value of λ=0.0086, LASSO model selects 8 feature genes: FCGR3B, FCGR3A, SPI1, FCER1G, ITGAX, CCR1, CSF3R, FCGR1A, NCF2. At the same time, an SVM algorithm was executed to rank the importance of all the 14 genes, and those having importance score above grouping threshold and within the top 10 were defined as feature genes 2 (Figure 3F). A random forest model was also constructed, and gene importance values (lncMSE and IncNodePurity) were calculated; cross-validation identified the top five genes SPI1, FCER1G, ITGAX, NCF2, and FCGR2A as feature genes (Figure 3G,3H). Taking the intersection of the above three types of candidate feature genes, a total of 4 feature genes were obtained (Figure 3I).

The four feature genes (SPI1, FCER1G, ITGAX, NCF2) obtained from the intersection of the three machine learning algorithms were further verified for diagnostic efficacy. ROC curve analysis in the GSE68004 and GSE73461 datasets showed that the AUC values of all four genes were greater than 0.8, indicating favorable diagnostic performance for KD. These four genes were significantly upregulated in KD samples compared with healthy controls in both datasets, with statistically significant expression differences, which confirmed their potential as diagnostic biomarkers for KD.

Critical pathway analyses

GSEA revealed that all four biomarkers were positively correlated with neutrophil extracellular trap (NET) formation, ribosome formation, and ribosome biogenesis in eukaryotes. They were negatively correlated with lipid and atherosclerosis pathways and Fc gamma receptor (FcγR)-related signaling pathways (Figure 4A-4D). Moreover, FCER1G, ITGAX, and NCF2 were positively associated with the NOD-like receptor signaling pathway, and SPI1 was negatively associated.

Figure 4 GSEA enrichment analysis of the biomarkers. (A) GSEA of FCER1G, (B) ITGAX, (C) NCF2, and (D) SPI1. GSEA, gene set enrichment analysis.

Immune cell infiltration analysis

Analysis of immune cell infiltration showing the proportion of 22 immune cell types that were significantly different between patients with KD and the healthy controls. The abundance of CD8⁺ T cells, resting memory CD4⁺ T cells, naive CD4⁺ T cells, and 12 other immune cell subsets was altered (P<0.05) (Figure 5A,5B). The biomarkers showed robust interactions with various immune cell types. The neutrophils and monocytes showed the strongest positive correlation (r>0.7, P<0.05), whereas CD8⁺ T cells exhibited the strongest negative association (r<−0.8, P<0.05; Figure 5C). Negative correlations were observed between the biomarkers and other immune cells, including naive CD4⁺ T cells and natural killer (NK) cells.

Figure 5 Immune cell analysis of biomarkers. (A) Stacked histogram showing immune cell proportions across samples. (B) Box plot showing differences in immune cell proportions between patients with KD and healthy controls. (C) Heatmap depicting correlations between biomarkers and immune cell subsets. *, P<0.05; **, P<0.01; ***, P<0.001; ****, P<0.0001. KD, Kawasaki disease.

Identification of m⁶A modification sites and chromosome localization of the biomarkers

Table 1 shows the number of m⁶A modification sites predicted in the four biomarkers. FCER1G and ITGAX were predicted to have a total of four very high-confidence m⁶A modification sites. There was one site in FCER1G and three in ITGAX. The location of the sites in the RNA secondary structures of the two transcripts is depicted in Figure 6A,6B, and Figure S3. The positions of the m⁶A sites for NCF2 and SPI1 and their respective RNA secondary structures are shown in Figure S1. FCER1G and NCF2 were both located on chromosome 1, whereas SPI1 and ITGAX were located on chromosomes 11 and 16 (Figure 6C).

Figure 6 RNA methylation modifications of biomarkers. (A,B) Four highly confident m⁶A modification sites were predicted in the FCER1G and ITGAX genes. (C) Chromosomal localization circle for the biomarkers. m⁶A, N⁶-methyladenosine.

Computationally predicted TFs and candidate drugs targeting multiple biomarkers

The miRDB database was used to predict miRNAs targeting the four biomarkers (Figure 7A). To elucidate the regulatory mechanisms of these biomarkers in KD, TFs were predicted using the KnockTF database, and a TF-mRNA regulatory network was constructed, which consisted of 4 biomarkers, 20 TFs, and 28 regulatory edges (Figure 7B). Drug-gene interaction network analysis identified 76 compounds as computationally predicted candidates targeting the biomarkers. Among these, saxagliptin was predicted to potentially interact with NCF2 and SPI1; however, this prediction remains speculative and requires experimental validation (Figure 7C).

Figure 7 Interaction network between the biomarkers and miRNAs, TFs, and drugs. (A) miRNA-biomarker regulatory network. (B) TF gene regulatory network of the biomarkers. (C) Drug-biomarker interaction network. miRNAs, microRNAs; TF, transcription factor.

Validation of biomarkers by RT-qPCR

RT-qPCR was performed on peripheral blood samples from patients with KD and healthy controls to validate the expression of four biomarkers. SPI1, FCER1G, ITGAX, and NCF2 were significantly upregulated in the KD samples compared with the controls. The results were consistent with the bioinformatics analyses (Figure 8).

Figure 8 Validation of biomarkers in clinical samples. (A-D) Expression of SPI1, FCER1G, ITGAX, and NCF2 mRNA. *, P<0.05; **, P<0.01. KD, Kawasaki disease; mRNA, messenger RNA.

Discussion

KD is a febrile illness that affects children below 5 years of age. It is a systemic vasculitis with a predilection for the coronary arteries (11). Late intervention can result in coronary artery aneurysms in nearly 25% of patients with possible adverse consequences, such as myocardial infarction, heart failure, or death (12). Recent studies indicate that various metabolites of Trp are important for immune modulation, inflammation (13), and vascular endothelial repair (14), and the gut microbiota can affect their production. Thus, identifying biomarkers associated with TM will provide insight into the pathology of KD, which may lead to improved treatments. Through an integrated analysis, we identified DEGs and TMRGs that were significantly enriched in inflammatory pathways. SPI1, FCER1G, ITGAX, and NCF2 were identified as potential biomarkers for TM in KD. These genes may not only be used as diagnostic markers but also represent targets for therapy and reveal the underlying molecular mechanisms of KD.

SPI1, FCER1G, ITGAX, and NCF2 are known for their roles in immune cells. They perform a variety of functions and are involved in protein complexes (like NCF2). SPI1 activity plays an important role in myeloid and dendritic cell transcription. It activates myeloid populations, such as macrophages, monocytes, and dendritic cells, during the acute KD phase. This subsequently causes inflammation in the blood vessels. SPI1 attenuates phagocytosis and inflammatory gene expression and metabolic pathways that affect macrophage polarization and function, which are pathologically relevant to endothelial inflammation (15). High levels of SPI1 in the blood or skin suggest increased inflammation from immune cells. FCER1G is important for antibody-dependent cell functions, such as antibody-dependent cellular cytotoxicity, phagocytosis, and cell activation, because it associates with multiple activating Fc receptors. The pathophysiology of KD and IVIG therapy is involved in Fc receptor biology. Previous studies have linked KD susceptibility and IVIG responsiveness with Fc receptor gene expression and methylation (16,17). Altered FCER1G expression may modify the signaling of the antibody-mediated inflammation or the efficacy of IVIG. High ITGAX expression is a hallmark of activated dendritic cells and inflammatory myeloid subsets, such as CD11c⁺ monocytes/macrophages. Dendritic cells are important for initiating and sustaining adaptive immune responses and for producing proinflammatory mediators. Recruitment of CD11c⁺ cells to vascular tissue may enhance local inflammation and antigen presentation, thereby contributing to vascular injury in KD (18). Activated dendritic cells and inflammatory myeloid subsets, such as CD11c⁺ monocytes/macrophages, express high ITGAX levels. Dendritic cells are essential for the initiation and maintenance of adaptive immune responses. They also produce proinflammatory mediators. The accumulation of CD11c⁺ cells at vascular locations may enhance local inflammation and antigen presentation, which may lead to vascular injury in KD (18). NCF2 is involved in the assembly and function of NADPH oxidase and regulates reactive oxygen species (ROS) production. During the development of the disease, ROS may modulate the activation of neutrophils, the formation of NETs, and signaling in endothelial cells. The atypical expression or function of NCF2 may be the cause of endothelial damage in KD and increased inflammatory responses (19,20).

The GSEA results revealed that FCER1G, NCF2, ITGAX, and SPI1 were significantly enriched in ribosome, lipid and atherosclerosis, and NOD-like receptor pathways. Previous studies have indicated that the abnormal activation of the atherosclerosis pathway plays an important role in the progression of vascular injury in KD and exhibits synergistic effects with lipid metabolism disorders. Inflammation is a core pathological process in KD. Key molecules in the atherosclerosis pathway (e.g., TNF-α, IL-6, and VCAM-1) are significantly activated during the acute phase of KD. They induce vascular endothelial cell damage, monocyte infiltration, and foam cell formation (21). In addition, oxidized low-density lipoprotein generated by lipid metabolism abnormalities can activate the atherosclerosis pathway. They exacerbate vascular wall inflammation and fibrosis by upregulating inflammatory cytokine expression and promoting vascular smooth muscle cell proliferation and migration, which ultimately leads to the formation and progression of CALs (22,23). Moreover, long-term follow-up studies have found that lipid metabolism disorders and the abnormal activation of the atherosclerosis pathway may persist in patients with KD, even after the acute phase, which significantly increases the risk of atherosclerotic cardiovascular disease in adulthood. This suggests that lipid metabolism disorders and abnormal alterations in the atherosclerosis pathway serve as important predictors of long-term cardiovascular risk in KD (24). These results suggest that children with KD exhibit distinct lipid metabolism disorders and abnormal activation of atherosclerotic pathways, which synergistically exacerbate vascular damage.

This study revealed significant differences in 15 immune cell types between patients with KD and healthy control samples by immunohistochemical analysis. Neutrophils and monocytes showed the strongest positive correlation, whereas CD8⁺ T cells, CD4⁺ T cells, and NK cells exhibited negative correlations. As pioneer cells of innate immunity, neutrophils play an important role in initiating and amplifying the inflammatory response in KD. The study demonstrated that the expression of key enzymes in TM-indoleamine 2,3-dioxygenase (IDO) and Trp 2,3-dioxygenase was significantly increased in peripheral blood neutrophils of children with KD, and these levels were positively correlated with disease activity. This suggests that neutrophils may contribute to the inflammatory progression of KD by driving abnormal TM (25). The role of TM disorders in KD vascular injury has been verified in multiple studies, with neutrophil-mediated TM reprogramming likely a key regulatory mechanism. The study also revealed that the serum canthurianine/Trp ratio was significantly increased in children with KD, which serves as a marker of activated TM. Further mechanistic studies revealed that neutrophils can upregulate IDO activity by releasing ROS, thus facilitating the conversion of Trp to canthurianine. Canthurianine subsequently induces vascular endothelial cell damage and exacerbates the pathological changes of vasculitis (26,27). The association between neutrophil subpopulation functional heterogeneity and TM offers a unique perspective for understanding KD pathogenesis. Neutrophil-derived IDO can affect dendritic cell maturation by regulating TM, inhibiting the differentiation of anti-inflammatory regulatory T cells, and enhancing proinflammatory helper T cell 17 responses, thereby exacerbating immune dysregulation and promoting KD inflammatory progression (28). These results suggest that neutrophils play an important regulatory role in KD by driving Trp metabolic reprogramming, which contributes to inflammatory activation, immune dysregulation, and vascular damage. This provides a new theoretical framework for examining the pathogenesis of KD.

The strong positive correlations (r>0.7) of SPI1, FCER1G, ITGAX, and NCF2 with myeloid cell populations suggest that these genes act as plastic molecular biomarkers. They also indicate that altered immune modulation related to KD is characterized by hyperactivated innate immunity and suppressed adaptive immunity (29,30). Collectively, these four biomarkers are related to myeloid immune-related pathways and may drive the inflammatory phenotype of acute KD through various mechanisms. These include promoting myeloid cell development and activation (e.g., SPI1) (31), enhancing phagocytosis and FcγR-mediated signaling (e.g., FCER1G and ITGAX), and enhancing oxidative burst and NET generation (e.g., NCF2) (32).

In the present study, the regulatory network analysis revealed that these biomarkers predict 14 miRNAs (e.g., hsa-miR-653-3p and hsa-miR-532-5p) and 86 TFs (e.g., ZFX and IRF2). The expression of biomarkers may be temporally and spatially controlled by these TFs across acute versus convalescent or blood versus tissue. These may be considered candidate validation targets through chromatin immunoprecipitation sequencing (ChIP-seq) or single-cell transcriptomics. Our drug-mRNA network revealed 76 candidate drugs that target the biomarkers, with saxagliptin representing a notable hit that affects NCF2 and SPI1. Saxagliptin is a type 2 diabetes medication that blocks dipeptidyl peptidase 4 (DPP-4) and protects blood vessels. Saxagliptin ameliorates ox-LDL-induced endothelial dysfunction through the regulation of AP-1 and NF-kappaB and the reduction of the production of cytokines and vascular adhesion molecules induced by ox-LDL (33). Saxagliptin, which inhibits myeloid cell activation by SPI1 and NCF2, may avoid NCF2-induced ROS/NADPH oxidase activity or inhibit SPI1-induced inflammatory programming. Whether saxagliptin or other identified compounds can modulate SPI1 or NCF2 activity in the context of KD inflammation remains entirely speculative and requires rigorous experimental investigation, including binding assays, cellular functional studies, and validation in KD animal models. At present, these computational findings should be viewed solely as hypothesis-generating and not as evidence of therapeutic potential.

Several limitations to this study should be considered. The sample size was small, consisting of only 200 samples in the training and validation sets. Moreover, there were batch effects and differences between studies because all datasets were derived from GEO. Although RT-qPCR results showed that our bioinformatics analysis showed consistent expression, this small sample size was not sufficient to support strong biological conclusions or to extend the findings to a wider KD population. The execution of functional experiments, such as gene knockout or overexpression studies, is necessary to establish the mechanistic roles of these genes. Although immune infiltration analyses were estimated, they should be confirmed through flow cytometry and other immunophenotyping approaches. The information for predicting drug-biomarker interactions was obtained from a database. In vitro or in vivo studies are needed to evaluate its clinical efficacy and translational potential. Finally, the machine learning approach in this study did not fully nest the feature selection procedure within cross-validation, which may have introduced information leakage and thus led to an optimistic estimation of model performance. Therefore, external validation was needed in the future using completely independent multicenter cohorts that were not involved in any model construction process, to objectively evaluate its generalization ability.

Conclusions

This study identified SPI1, FCER1G, ITGAX, and NCF2 as potential biomarkers associated with TM in KD, which are implicated in immune regulatory pathways and transcriptional regulation. We present novel findings regarding the molecular mechanisms driving KD pathogenesis, which provide a valuable basis for the development of novel diagnostic and therapeutic approaches. It is necessary to apply larger multicenter cohorts to combine with multiomics approaches in future studies to validate these biomarkers. These molecular insights may lead to improved risk stratification, prompt diagnosis, and improved treatment for KD.

Acknowledgments

We thank Dr. Zheyan Chen for his kind help in graphical drawing and statistics.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tp.amegroups.com/article/view/10.21037/tp-2025-1-905/rc

Data Sharing Statement: Available at https://tp.amegroups.com/article/view/10.21037/tp-2025-1-905/dss

Peer Review File: Available at https://tp.amegroups.com/article/view/10.21037/tp-2025-1-905/prf

Funding: This work was supported by the Natural Science Foundation of Hunan Province of China (No. 2022J70016), Key R&D Project of Hunan Province (No. 2018SK21610), Hunan Provincial Key Laboratory of Pediatric Respiratory Medicine (No. 2019TP1043), Hunan Provincial Administration of Traditional Chinese Medicine (No. B2024101), National Key Specialty Construction Project (No. Z2023029) and Hunan Provincial Administration of Traditional Chinese Medicine (No. B2024101).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tp.amegroups.com/article/view/10.21037/tp-2025-1-905/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. This study was approved by the Ethics Committee of The First Affiliated Hospital of Hunan Normal University (No. [2024]-135). All participants’ guardians were informed about the study and provided informed consent.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Zeinali M, Frishman WH, Aronow WS. Kawasaki Disease. Cardiol Rev 2025; Epub ahead of print. [Crossref] [PubMed]
2. Lin MT, Wu MH. The global epidemiology of Kawasaki disease: Review and future perspectives. Glob Cardiol Sci Pract 2017;2017:e201720. [PubMed]
3. Mat Bah MN, Alias EY, Razak H, et al. Epidemiology, clinical characteristics, and immediate outcome of Kawasaki disease: a population-based study from a tropical country. Eur J Pediatr 2021;180:2599-606. [Crossref] [PubMed]
4. Rife E, Gedalia A. Kawasaki Disease: an Update. Curr Rheumatol Rep 2020;22:75. [Crossref] [PubMed]
5. Nishimura M, Yamashita A, Matsuura Y, et al. Upregulated Kynurenine Pathway Enzymes in Aortic Atherosclerotic Aneurysm: Macrophage Kynureninase Downregulates Inflammation. J Atheroscler Thromb 2021;28:1214-40. [Crossref] [PubMed]
6. Wang Q, Liu D, Song P, et al. Tryptophan-kynurenine pathway is dysregulated in inflammation, and immune activation. Front Biosci (Landmark Ed) 2015;20:1116-43. [Crossref] [PubMed]
7. Paz GA, Rangel MV, Farias CL, et al. Acute and chronic effects of physical exercise on atherosclerosis, kynurenine pathway, endothelial function and inflammation in patients with coronary artery disease: a clinical trial protocol. BMJ Open Sport Exerc Med 2025;11:e002432. [Crossref] [PubMed]
8. Zhang J, Jiang X, Pang B, et al. Association between tryptophan concentrations and the risk of developing cardiovascular diseases: a systematic review and meta-analysis. Nutr Metab (Lond) 2024;21:82. [Crossref] [PubMed]
9. Wang Y, Liu H, McKenzie G, et al. Kynurenine is an endothelium-derived relaxing factor produced during inflammation. Nat Med 2010;16:279-85. Erratum in: Nat Med 2010;16:607. [Crossref] [PubMed]
10. Fan X, Li K, Guo X, et al. Metabolic profiling reveals altered tryptophan metabolism in patients with kawasaki disease. Front Mol Biosci 2023;10:1180537. [Crossref] [PubMed]
11. Jone PN, Tremoulet A, Choueiter N, et al. Update on Diagnosis and Management of Kawasaki Disease: A Scientific Statement From the American Heart Association. Circulation 2024;150:e481-500. [Crossref] [PubMed]
12. van Stijn D, Schoenmaker NJ, Planken RN, et al. Myocardial infarction due to thrombotic occlusion despite anticoagulation in Kawasaki disease - a case report. BMC Pediatr 2022;22:85. [Crossref] [PubMed]
13. Seo SK, Kwon B. Immune regulation through tryptophan metabolism. Exp Mol Med 2023;55:1371-9. [Crossref] [PubMed]
14. Mannarino MR, Bianconi V, Scalisi G, et al. A tryptophan metabolite prevents depletion of circulating endothelial progenitor cells in systemic low-grade inflammation. Front Immunol 2023;14:964660. [Crossref] [PubMed]
15. Zhang G, Lu J, Zheng J, et al. Spi1 regulates the microglial/macrophage inflammatory response via the PI3K/AKT/mTOR signaling pathway after intracerebral hemorrhage. Neural Regen Res 2024;19:161-70. [Crossref] [PubMed]
16. Chang LS, Ming-Huey Guo M, Lo MH, et al. Identification of increased expression of activating Fc receptors and novel findings regarding distinct IgE and IgM receptors in Kawasaki disease. Pediatr Res 2021;89:191-7. [Crossref] [PubMed]
17. Uittenbogaard P, Netea SA, Tanck MWT, et al. FCGR2/3 polymorphisms are associated with susceptibility to Kawasaki disease but do not predict intravenous immunoglobulin resistance and coronary artery aneurysms. Front Immunol 2024;15:1323171. [Crossref] [PubMed]
18. Wu J, Wu H, An J, et al. Critical role of integrin CD11c in splenic dendritic cell capture of missing-self CD47 cells to induce adaptive immunity. Proc Natl Acad Sci U S A 2018;115:6786-91. [Crossref] [PubMed]
19. Beltran JVB, Lin FP, Chang CL, et al. Single-Cell Meta-Analysis of Neutrophil Activation in Kawasaki Disease and Multisystem Inflammatory Syndrome in Children Reveals Potential Shared Immunological Drivers. Circulation 2023;148:1778-96. [Crossref] [PubMed]
20. Paclet MH, Laurans S, Dupré-Crochet S. Regulation of Neutrophil NADPH Oxidase, NOX2: A Crucial Effector in Neutrophil Phenotype and Function. Front Cell Dev Biol 2022;10:945749. [Crossref] [PubMed]
21. Chen S, Savas A, Atici A, et al. Interleukin-1 Signaling on Vascular Smooth Muscle Cells Accelerates Atherosclerosis in a Murine Model of Kawasaki Disease. J Am Heart Assoc 2025;14:e040687. [Crossref] [PubMed]
22. Ferraro S, Benedetti S, Mannarino S, et al. Prediction of atherosclerotic cardiovascular risk in early childhood. Clin Chim Acta 2024;552:117684. [Crossref] [PubMed]
23. Motoji Y, Fukazawa R, Matsui R, et al. Kawasaki Disease-like Vasculitis Facilitates Atherosclerosis, and Statin Shows a Significant Antiatherosclerosis and Anti-Inflammatory Effect in a Kawasaki Disease Model Mouse. Biomedicines 2022;10:1794. [Crossref] [PubMed]
24. Yohannan TM. Heart Transplant, Kawasaki Disease, and Bone Marrow Transplant: Are There Consequences? Curr Atheroscler Rep 2022;24:243-51. [Crossref] [PubMed]
25. Fitz AR, Klimov DK, Lockhart C. Binding of Antimicrobial Peptide Indolicidin to DMPC Bilayer Using Replica-Exchange Molecular Dynamics. J Chem Inf Model 2025;65:9251-60. [Crossref] [PubMed]
26. Dickerhof N, Ashby LV, Ford D, et al. Dioxygenation of tryptophan residues by superoxide and myeloperoxidase. J Biol Chem 2025;301:108402. [Crossref] [PubMed]
27. Shahriari M, Jafari M, Doustdar F, et al. Comparative study of the interactions between C60 fullerene and SARS-CoV-2, HIV, eukaryotic, and bacterial model membranes: Insights into antimicrobial strategies with C60-peptide hybrids. Int J Biol Macromol 2024;271:132399. [Crossref] [PubMed]
28. Du Q, Zhang X, Sun H, et al. Value of the neutrophil percentage-to-albumin ratio in predicting intravenous Immunoglobulin resistance in Kawasaki disease: a retrospective cohort study. Ital J Pediatr 2025;51:330. [Crossref] [PubMed]
29. Wang N, Chen Z, Zhang F, et al. Intravenous Immunoglobulin Therapy Restores the Quantity and Phenotype of Circulating Dendritic Cells and CD4(+) T Cells in Children With Acute Kawasaki Disease. Front Immunol 2022;13:802690. [Crossref] [PubMed]
30. Wang S, Qian G, Liu Y, et al. Kawasaki disease: insights into the roles of T cells. Front Immunol 2025;16:1582638. [Crossref] [PubMed]
31. Huang J, Chen W, Jie Z, et al. Comprehensive Analysis of Immune Implications and Prognostic Value of SPI1 in Gastric Cancer. Front Oncol 2022;12:820568. [Crossref] [PubMed]
32. Schopf P, Mills MJ, Warshel A. The entropic contributions in vitamin B12 enzymes still reflect the electrostatic paradigm. Proc Natl Acad Sci U S A 2015;112:4328-33. [Crossref] [PubMed]
33. Ma S, Bai Z, Wu H, et al. The DPP-4 inhibitor saxagliptin ameliorates ox-LDL-induced endothelial dysfunction by regulating AP-1 and NF-κB. Eur J Pharmacol 2019;851:186-93. [Crossref] [PubMed]