Integrative Analysis of Clinical and Bioinformatics Databases to Reveal the Role of Peripheral Innate Immunity in Kawasaki Disease

,


Introduction
Kawasaki disease (KD), also known as cutaneous mucosal lymph node syndrome, was first reported from Japan in 1967 by Kawasaki [1].KD is a disorder characterized by abnormal inflammation and an atypical immune response.Genetic factors, particularly variations in genes associated with the immune system, contribute to the risk of developing this condition [2].Immune dysregulation and the activation of T cells and monocytes play a role in the disease [3,4].
Inflammation and the formation of aneurysms in the coronary arteries are caused by dysfunction in the endothelial lining of blood vessels [5].While the exact infectious agent responsible for KD remains unknown, some theories propose that specific viruses or bacteria may trigger the immune response [6].However, the exact cause of KD remains unknown.
Immunohistochemical analysis of postmortem tissue from patients with KD has revealed the presence of monocytes, macrophages, and neutrophils [3,4], as well as activated CD8 + T cells [7] and IgA + plasma cells [8,9], in the arterial wall.Infiltrating immune cells release proinflammatory cytokines (e.g., TNF and IL-1β) that then contribute to the development of endothelial lesions and CAL [10,11].In addition, endothelin autoantibodies can cause endothelial disease in KD [12].Recent research using single-cell RNA sequencing (scRNA-seq) has revealed alterations to immune cells in KD patients at the acute stage, including increases in immunomodulatory T cells, NK cells [13], plasma cells, and B cells [14].Immunomodulatory genes are involved in the pathogenesis of KD [15].Changes have also been observed in monocyte developmental trajectory [16], and CD14+CD16-monocytes were found to be expanded in KD.Multisystem Inflammatory Syndrome in Children (MIS-C), an inflammatory disorder associated with immune dysfunction, has clinical manifestations similar to KD which are sometimes difficult to distinguish.Immunological studies showed a decrease in the number of follicular B cells, an increase in the number of terminally differentiated CD4+T lymphocytes (LYM), and a decrease in IL-17A levels in MIS-C [17].These studies suggest that innate immunity plays an important role in KD pathogenesis.Furthermore, the most effective treatment for KD is currently high-dose intravenous gamma globulin (IVIG), which reduces CAL incidence.However, up to 15%-20% of patients do not respond to IVIG therapy, and CAL progression remains unaffected [17,18].An improved understanding of the mechanisms underlying the involvement of innate immunity could help to explain the nonresponsive patients.One promising area of research is immunity-related genes (IRGs), known to be critical during immune infiltration [19,20].No research has yet examined their regulation or expression characteristics in KD.
The objective of our study was to explore the role of innate immunity in the pathogenesis of KD and its potential in predicting treatment response.Here, we used scRNA-seq and RNA microarray data to investigate the expression characteristics of IRGs and their possible regulatory mechanisms.Furthermore, we validated the hub genes in Chinese KD patients.We demonstrated that CD14 + CD16 + monocytes are important effector cells in the acute stage of KD, and regulate cytokine levels through promoting the expression of inflammatory transcription factors.Moreover, CD14 + CD16 + and CD14 − CD16 + monocytes are closely related to the efficacy of IVIG treatment.Our findings provide insight into the role of innate immunity in KD pathogenesis and offer valuable biomarkers that can be useful for improving treatment efficacy.
2.2.Analysis of scRNA Sequencing Data.Processing of GSE168732 scRNA sequencing data followed methods from a previous study [14].The quality control standards are as follows: first, we assigned the CreateSeuratObject function various parameters, including min.Cell = 3 and min.features = 200.For most samples, the total UMI count is between 2,000 and 60,000, and the mitochondrial gene percentage was <5%.For P1 before therapy, a lower cutoff of total UMI count (1,000) was used owing to its lower median UMI count per cell.After quality control and filtering, 38,712 cells were selected for the analysis.Processing of the GSE152450 dataset also followed published methods [16], yielding 6,283 cells for further analysis after quality control and filtering.The "Findmarkers" function of the R package Seurat [23] was used to perform downstream analyses of differentially expressed genes (DEGs).Criteria for DEGs were adjusted P <0:05 and |log fold change| ≥ 0.25.R packages Monocle2 [24] and Dorothea [25] were used to determine the pseudo-time developmental trajectory and transcription factor activity of myeloid subpopulations, respectively.
2.5.Weighted Gene Coexpression Network Analysis.Gene coexpression networks were constructed with the R package weighted gene coexpression network analysis (WGCNA) [29].After analyzing the correlation modules and different sample periods, MEblue modules were selected for further analysis because they had the highest positive correlation with responders from acute samples.The code and tutorial were obtained online (https://rstudio-pubs-static.s3.amazona ws.com/687551_ed469310d8ea4652991a2e850b0018de.html).
2.6.Functional Analysis.Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed for DEGs using clusterProfiler [30].Differences were considered significant at adjusted P <0:05.The R package GSVA (gene set variation analysis) was applied on hallmark gene sets among cell clusters [31].All patients had acute KD.Individuals were excluded if they had other systemic diseases (e.g., kidney disease), or did not receive IVIG.The diagnostic criteria for KD were based on literature [33].Patients were divided into IVIGresponsive and IVIG-resistant groups [33].Inflammatory markers, including C-reactive protein (CRP), procalcitonin (PCT), white blood cells, neutrophils (NEUT), LYM, monocytes (MONO), and platelets (PLT) were collected; the neutrophil/lymphocyte ratio (NLR) was then calculated.These indicators were measured before IVIG administration.
2.9.Statistical Analysis.Wilcoxon tests or Kruskal-Wallis tests were used for between-group comparisons.All statistics and visualizations used R and the ggplot2 package [34].
Results from descriptive analyses are reported as percentages and medians (interquartile spacing), as appropriate.Significance was set at P <0:05.S1.Cells were assigned to known clusters based on marker genes [11].Using t-SNE analysis, we visualized nine clusters: myeloid cells, T cells, B cells, NK cells, erythrocytes, platelets, TABLE 2: Primers used for RT-qPCR with their sequence.

6
Mediators of Inflammation plasma cells, hematopoietic stem cells, and progenitor cells, and mixed clusters (Figure 1(a) and Table S1, for all marker gene expression, see Figure S2).The KD group exhibited an increase in the number of myeloid and B cells (Figure 1(a)).
We then identified DEGs (between KD and NC groups) in each cell cluster (the top five DEGs per cluster are shown in Figure 1(b)).
Compared with the NC group, T cells in the KD group decreased significantly, while myeloid cells increased, but not significantly (Figures 1(c) and 1(d) and Table 3).

Expression of IRG in PBMC Clusters from KD Patients.
We selected IRGs according to the import database, and then obtained them from DEGs of each cluster in the PBMCs, resulting in 381 IRGs for analysis (Table S2).
The uniform manifold approximation and projection analysis revealed that myeloid cells had the highest number of DEGs among the nine cell subgroups (Figure 2(a)).We confirmed this outcome with area under the curve (AUC) analyses of IRG activity, where cells expressing more genes had higher AUC values (Figures 2(b) and 2(c)).Next, GO and KEGG analyses of 381 IRGs revealed that they were mainly enriched in immune response-activating cell surface receptor signaling pathway and cytokine-cytokine receptor interaction (Figures 2(d) and 2(e)).These data suggest that myeloid cell clusters are heterogeneous and that the innate immune system plays an essential role in the development of KD.Therefore, we mainly focused on myeloid cells in subsequent analyses.

Characteristics of Monocyte Expression in KD.
Further reclustering analysis of myeloid cells based on cell-type marker genes showed that the main components were CD14 + CD16 + monocytes, leukocyte immunoglobulin-like receptor A4 + plasmacytoid dendritic cells (LILRA4 + pDC), CD14 + CD16 − monocytes, and CD14 − CD16 + monocytes   10 Mediators of Inflammation (Figure 3(a) and Figure S3).A previous study reported that peripheral blood monocytes play a central role in acute KD [35].We thus compared monocyte expression between the KD and NC groups, revealing that monocytes in the KD group expressed the calgranulin genes S100A8, S100A9, and S100A12 at significantly higher levels (Figure 3(b)).Dot plots showed that these inflammatory genes were mainly concentrated in CD14 + CD16 + monocytes (Figure 3(c)).Pseudo-time analysis of monocyte subsets indicated that the developmental trajectory started from CD14 − CD16 + monocytes and ended with CD14 + CD16 + monocytes; additionally, inflammatory gene expression increased during development (Figure 3(d)).Next, GSVA of significant hallmarks in myeloid cells indicated that immunomodulatory genes associated with the inflammatory response were significantly upregulated in CD14 + CD16 + monocytes from the KD group (Figure 4

Mediators of Inflammation
DEGs in KD (Table S3), including the inflammatory genes S100A8 and S100A12 (Figure 5(a)).Among these DEGs, we found 289 IRGs in the KD group (Figure 5(b)).According to GO analysis, these IRGs were enriched in pathways related to innate immunity, including neutrophil chemotaxis and myeloid leukocyte migration.(Figure 5

16
Mediators of Inflammation  Mediators of Inflammation genes were mainly enriched in innate immune pathways, such as neutrophil activation, neutrophil degranulation, and neutrophil-mediated immunity (Figure 7(c)).

Discussion
Increasing evidence suggests that the innate immune response may be involved in KD pathogenesis [13,14,36].Of particular note is the involvement of monocytes, innate immune cells from the bone marrow that have multiple functions, including tissue development and homeostasis, inflammation initiation and resolution, and tissue repair [37][38][39].Monocytes are involved in the pathogenesis of other cardiovascular diseases [40][41][42], and changes in monocyte subsets have been observed in KD patients [16].Other evidence to suggest that innate immunity plays a major role in KD pathogenesis includes research using mouse models of the disease.One such study identified a nucleotide-binding oligomerization domain-containing protein (NOD) 1 ligand as an important inducer of coronary arteritis [43].Another study using a mouse model found that inhibition of interleukin-1β attenuates vasculitis [44].
Here, we followed up on those previous reports in an effort to better understand the potential mechanisms of innate immunity in KD pathogenesis.Through an analysis of scRNA sequencing data and marker gene expression, we first found that T cells in patients with KD were significantly reduced (Figure 1(d)), and also confirmed that myeloid cells -the main source of monocytes-were the most heterogeneous cell group (Figure 2(a)).Our results corroborated earlier research showing changes to T cells in KD patients [13,14], but myeloid cells have rarely been studied in KD, although their involvement has been confirmed in various other inflammatory diseases [42,45].Here, functional analyses revealed that DEGs in myeloid cells were closely related to inflammatory response and cytokine regulation (Figures 2(d) and 2(e)), implicating them in the development of KD.Our subgroup analysis then revealed that myeloid cells in KD patients were mainly composed of monocytes and LILR4 + pDC.Monocytes can be broken down into three subtypes: CD14 + CD16 + , CD14 + CD16 − , and CD14 − CD16 + [46], or intermediate, classical, and nonclassical in humans.Our observed monocyte composition was consistent with previous studies [16].We also noted that the DEGs S100A12, S100A9, S100A8, and ITGAM were primarily expressed in CD14 + CD16 + monocytes (Figure 3(c)).Previous studies have found that patients with acute KD had higher circulating concentrations of S100A8/9 heterodimer and S100A12 than patients with fever caused by other diseases; IVIG treatment decreased these concentrations [47,48].Even after the acute stage, KD patients with large coronary aneurysms maintained higher S100A8/9 heterodimer levels [47].In addition, elevated S100A8, S100A9, and S100A12 levels have been found in inflammatory diseases associated with immune disorders, such as juvenile idiopathic arthritis [49].Hence, the increased levels of S100A8, S100A9, and S100A12 proteins cannot be exclusively attributed to KD.In addition, integrin ITGAM was upregulated in KD vasculopathy [50].Furthermore, we observed elevated expression of inflammatory regulation genes associated with CD14 + CD16 + cells (Figure 3(c)), and GSVA indicated that this monocyte subtype is important to the inflammatory response in KD (Figure 4(b)).Similarly, recent reports have shown that, in addition to having high antigen presentation capacity, CD14 + CD16 + cells highly express proinflammatory cytokines [51,52].Data from peripheral blood of KD patients also revealed a link between these cells and inflammation [53].Taken together, our study and previous research all indicate that CD14 + CD16 + monocytes are heavily involved in the inflammatory response of acute KD.
Changes in the trajectory of monocyte development are closely related to disease occurrence [54,55].Here, pseudotime analysis of monocyte subsets revealed a trajectory from CD14 − CD16 + monocytes to CD14 + CD16 + monocytes (Figure 3(d)), contradicting earlier research showing that CD14 + CD16 − monocytes are significantly elevated in acute KD [16].Moreover, during monocyte development, inflammation-related gene expression increased, but their expression in CD14 + CD16 + monocytes actually decreased at the end of the developmental trajectory.By contrast, the expression of inflammatory transcription factors was significantly higher in CD14 + CD16 + monocytes (Figure 4(b)).Therefore, we speculate that CD14 + CD16 + monocytes are the final effector cells in acute KD and that they regulate cytokine levels by promoting the expression of inflammatory transcription factors.
Although IVIG is an effective treatment for KD, drug resistance rates are high [56,57].Because IVIG-resistance is associated with an increased incidence of coronary artery aneurysms, IVIG-resistant patients should be identified before initiating treatment because they may benefit from additional anti-inflammatory therapy.In this study, we therefore analyzed multiple GEO scRNA datasets to identify potential genetic markers that could distinguish between IVIG-resistant and -responsive patients.Our results revealed that IVIG-responsive patients had highly synergistic differential genes (Figures 7(a) and 7(b)) that are involved in neutrophil activation and degranulation (Figures 7(c) and 8).Neutrophils are a critical part of innate immunity but can have harmful effects if excessively activated [58], causing immune diseases such as rheumatoid arthritis [59,60] and vasculitis [61,62].Elevated peripheral neutrophils in patients with KD are associated with coronary artery dilation and IVIG resistance [63,64].The NLR is a comprehensive indicator of neutrophil activation and immune disorders that also plays an important role in KD [64,65].Here, our clinical data revealed that IVIG-resistant patients had higher NLR than IVIG-responsive patients, validating the results of scRNA sequencing.
Further analysis showed that DEGs associated with neutrophils were mainly expressed in CD14 + CD16 + and CD14 − CD16 + monocytes (Figures 9(a) and 9(b)).In addition, the effect of IVIG on acute KD is closely related to the expression of immunomodulatory genes that activate neutrophils in these two monocyte subsets.The PPI network analysis then revealed that all hub genes in KD (Figure 9(c)) are inflammation-related genes that regulate the innate immune system.Broadly, these results corroborate a prior genome-wide transcriptome analysis demonstrating that increased CD177 transcript levels activate neutrophils and Mediators of Inflammation are closely related to KD [66].Autopsy results of patients who died during the acute phase of KD suggest neutrophilic involvement in damage to the coronary arteries [67].Consistent with our results, it has been found that the hub gene CCR1 is an important gene in the pathogenesis of KD [15].The hub genes TLR8 and IL1B are also closely related to neutrophil degranulation in patients with KD [68].All available data thus provide evidence of neutrophil activation being crucial to KD pathogenesis.Furthermore, it has been observed that IL1B, SOCS3, and IL1RN exhibit high expression levels in other autoimmune diseases that are closely linked to impaired innate immune function [69].
Although there are many clinical predictors of IVIG nonresponse in KD, the clinical application value is limited.In our study, we screened for hub genes by analyzing the expression of differential genes in monocytes of KD patients with IVIG response and IVIG nonresponse, and these results were further verified using RT-qPCR.The validated results indicated that the expression of IL1R1 (P <0:01), SOCS3 (P <0:01), TLR8 (P <0:01), CCR1 (P <0:01), IL1B (P <0:01), and IL10RB (P <0:001) was significantly up-regulated in the IVIG nonresponsive group compared with that in the control group and the IVIG responsive group, while IL1R2 (P <0:01), IL1RN (P <0:01), IL4R (P <0:01), and IFNGR1 (P <0:01) were significantly up-regulated in the IVIG nonresponsive group compared to that in the normal group, but there was no significant difference compared to the IVIG responsive group (Figure 10).The results suggested that IL1R1, SOCS3, TLR8, CCR1, IL1B, and IL10RB could be used as hub genes for screening IVIG responsive and IVIG nonresponsive patients of KD.
To our knowledge, this is the first study to elaborate on the role of innate immunity in the pathogenesis of KD and the mechanism of IVIG resistance.Combining both scRNAseq and microarray data analysis, our data strongly illustrates the involvement of innate immunity in the pathogenesis of KD and provides insights into the mechanism of IVIG resistance.Specifically, abnormal expression of immunomodulatory genes in CD14 + CD16 + and CD14 − CD16 + monocytes seems to trigger neutrophil activation, causing the worst symptoms of the disease and influencing the response to IVIG.The limitations of this study are as follows: first, as the study mainly relied on published RNA microarray datasets and RNA sequencing datasets for analysis, key clinical data could not be obtained.Second, the mechanism of action of these identified hub genes needs to be further studied.

FIGURE 1 :
FIGURE 1: Cluster distribution analysis.(a) Cell types identified in peripheral blood mononuclear cells (PBMC) using uniform manifold approximation and projection (UMAP).(b) Cluster heat map of top five differentially expressed genes (DEGs) per cluster.(c) Cluster distribution in each sample.(d) Cluster distribution in Kawasaki disease (KD) and normal control (NC) groups.* P <0:05.

FIGURE 2 :
FIGURE 2: Immunity-related gene (IRG) scores of PBMC clusters in KD.(a) DEGs in cell clusters from KD samples.Myeloid cells had the most DEGs.(b) Scores of 381 screened IRGs.The threshold was 0.11.(c) t-SNE plots of IRG scores in all clusters.Myeloid cells express more genes and exhibit higher AUC values.(d) Gene Ontology (GO) analysis of DEGs in myeloid cells, revealing enrichment in immune responseactivated cell surface receptor signaling pathway and immune response-activated signal transduction.(e) KEGG analysis of DEGs in myeloid cells, revealing enrichment in pathways associated with Th17 cell differentiation, cytokine-cytokine receptor interaction, and natural killer cell-mediated cytotoxicity.

3. 4 .FIGURE 3 :
FIGURE 3: Analysis of myeloid cell subsets.(a) Cell types identified in myeloid cells using UMAP.(b) Volcano plot of DEGs (|logFC| ≥ 0.25 and adjusted P <0:05).The top 20 DEGs are marked.(c) Expression of four major DEGs in myeloid cells.The dot size represents the percent expressed, and the color scale represents the average expression.The DEGs were mainly expressed in CD14 + CD16 + monocytes.(d) Pseudotime analysis of monocyte subsets in KD, using Monocle.The developmental trajectory began with CD14 − CD16 + monocytes and ended with CD14 + CD16 + monocytes; inflammatory gene expression increased during development.

FIGURE 4 :FIGURE 5 :
FIGURE 4: Functional enrichment analysis of DEGs in myeloid cells.(a) Gene set variation analysis (GSVA) in myeloid cell subsets.CD14 + CD16 + monocytes were the most abundant cells with DEGs related to inflammatory signaling pathways.(b) Heat map of transcription factors with functional enrichment in monocytes and LILR4 + pDC cells.CD14+CD16+ monocytes have the most immunomodulatory transcription factors.

FIGURE 7 :FIGURE 8 :
FIGURE 7: Analysis of DEGs across KD phenotypes (GSE63881): acute IVIG-responsive KD, acute IVIG-resistant KD, convalescent IVIGresponsive KD, and convalescent IVIG-resistant KD.(a) Cluster analysis of differentially expressed IRGs.Each color represents a module in the gene coexpression network constructed using weighted gene coexpression network analysis (WGCNA).(b) Membership in the blue module, showing DEGs that were positively correlated with IVIG-responsive KD phenotypes.(c) These DEGs are primarily associated with neutrophil activation, neutrophil degranulation, and neutrophil-mediated immunity.

TABLE 1 :
Enrolled datasets in the current study.
95°C for 10 min, followed by 95°C for 15 s and 60°C for 1 min for 40 cycles.GAPDH was selected as the internal control for mRNA, and the relative expression level of mRNA was calculated by the relative quantification (2 −ΔΔCt ) method.The primer sequences are listed in Table2.
[32]Protein-Protein Interaction Analysis and Validation ofHub Genes.Network analysis of protein-protein interactions (PPIs) was performed using STRING (https://string-db.org/)[32].Hub genes were identified in Cytoscape version 3.8.2.total RNA was extracted from each sample using TRIzol (Invitrogen, United States) according to the manufacturer's instructions.The cDNA was synthesized using the SuperScript III Reverse Transcriptase Kit (Invitrogen, USA).RT-qPCR was performed with Power SYBR Green PCR Master Mix (TransGen Biotech, China) on an ABI 7500 fast real-time PCR system.The amplification reaction procedure was as follows:

TABLE 3 :
Cell numbers of each cluster.
Cluster Myeloid cells B cells T cells NK cells Erythrocytes Platelet Plasma cells Mixed Hematopoietic stem and progenitor

TABLE 4 :
Comparison of the IVIG-responsive and -resistant patients with KD in terms of demographic and clinical variables.