The mechanism underlying B-cell developmental dysfunction in Kawasaki disease based on single-cell transcriptomic sequencing

Background: Kawasaki disease (KD) is an acute systemic vasculitis that can lead to acquired heart disease in children mostly from in developed countries. The previous research showed that B cells in KD patients underwent a profound change in both the cell numbers and types after intravenous immunoglobulin (IVIG) therapy.

Methods: We performed the single-cell RNA-sequencing for the peripheral blood mononuclear cells (PBMCs) from three febrile patients and three KD patients to investigate the possible mechanism underlying B cell developmental dysfunction in KD. The pseudo-time analysis was employed to study the developmental trajectories of the PBMCs in febrile control and KD patients.

Results: Overall single-cell expression profiles show that the biological processes of immunity, B cell activation pathway and their related biological entities are repressed in KD patients before IVIG treatment compared to febrile patient and KD patients after IVIG treatment. The differentially expressed gene analyses further demonstrate that B cell signaling pathway is downregulated in B cells and plasma blast cells of KD patients before treatment while cell cycle genes and MYC gene are upregulated in dendritic cells (DCs) and hematopoietic stem and progenitor cells (HSPCs) of KD patients before treatment. The biological process of immune response is upregulated in the HSPCs of KD patients before treatment in our dataset while the biological process of inflammatory response is upregulated in the HSPCs of KD patients before treatment in GSE168732 dataset. Single-cell trajectory analyses demonstrate that KD patients before treatment have a shortened developmental path in which B cells and T cells are failed to differentiate into separate lineages. HSPD1 and HSPE1 genes show an elevated expression level in the early cell development stage of KD patients before treatment accompanied with the repression of MYC, SPI1, MT2A and UBE2C genes. Our analyses of all B cells from KD patients before treatment show most of B cells are arrested in a transitional state with an ill developmental path compared with febrile patients and KD patients after treatment.

Conclusion: Our results indicate that the immune premature HSPCs accompanied with the abnormal expression dynamics of cell cycle and SPI1 genes are the mechanism underlying B cell developmental dysfunction in KD patients.

1 Introduction

Kawasaki disease (KD) is an acute systemic vasculitis that occurs predominantly in children under 5 years of age, which is the most common cause of acquired heart disease in children in developed countries (1). Coronary artery lesion (CAL) is a major complication of KD. High-dose intravenous immunoglobulin (IVIG) combined with aspirin is the standard treatment for KD (1). The diagnosis of KD mainly depends on clinical features. Several febrile diseases have similar clinical manifestations to KD, such as scarlet fever, EB virus infection, juvenile idiopathic arthritis, measles, and adenovirus infection, leading to difficulties in the diagnosis and early treatment of KD (2). The etiology of KD so far remains unknown. The widely proposed theories have been in the categories of environmental fungal toxin exposure (3), autoimmune pathogenesis, and infectious diseases, but none of them have been proven (4). The most widely accepted etiologic hypothesis suggested that KD was caused by an infectious agent, which infected children with a genetic susceptibility (5).

In recent years, the hygiene hypothesis, initially proposed to explain the etiology of allergic diseases, had been introduced to explain the pathogenesis of KD (6). The hygiene hypothesis assumes that lack of early-life microbial exposure results in an impaired or delayed immune maturation, leading to an abnormal response to the microbes harmless to healthy population, thus triggering KD (7, 8).The KD incidence gradually increased after its emergence in 1967 (9). The epidemic data showed an important characteristic of KD that the disease is relatively rare in children less than 6 months or more than 5 years old, with a peak incidence in children between 6 and 24 months old in different ethnic groups (10–12). KD was followed by a prolonged period of antigen-specific “split” T-cell anergy, which reflected a maturational defect in immune responsiveness of KD patients (13). It suggests that KD may be caused by the immune maturation delay in early childhood, which is probably due to a nearly perfect sanitary environment. The B-cell development dysregulation may be the potential mechanism behind the immune maturation delay of KD (14). Genetic studies further showed that KD susceptibility genes are involved in B-cell development and function (15–17), underscoring the crucial role of B-cell immunity in KD. Although these studies related the hygiene hypothesis of KD to B-cell development, the possible pathogenic mechanism has not been investigated at the cellular level.

Single-cell transcriptomic sequencing (scRNA-seq) is a relatively new tool for quantifying the gene expression profile for an individual cell. Recent studies had utilized scRNA-seq data to study the immune disorders in KD (18–20). They revealed the dynamic changes of immune cells, differential gene expression, and even biological processes in KD compared with other conditions, providing new insights into the etiology of KD and the therapeutic mechanism of IVIG. In this study, we obtained the scRNA-seq of peripheral blood mononuclear cells (PBMCs) from febrile patients and KD patients before and after IVIG treatment to explore the possible mechanism of B-cell development dysfunction in KD. We found that the biological processes related to immunity and B-cell activation pathway were upregulated in the overall single-cell expression feature of febrile and KD patients after IVIG treatment, and the activation of cell cycle pathways was presented in dendritic cells (DCs) and hematopoietic stem and progenitor cells (HSPCs) pretreatment KD. The analysis of peripheral blood HSPCs confirmed the early immune activation propensity in pretreatment KD. Single-cell trajectory analyses further revealed that KD patients before treatment had a shortened cell development route and an ill cell differentiation outcome. The pseudo-time analyses indicated that a large proportion of B cells in KD patients before treatment are arrested in the transitional state. Our study indicates that the abnormal expression pattern of HSPD1 and HSPE1 in immune premature HSPCs would result in the dysregulated B-cell development and IVIG could rescue this process. Our study illustrates the probable etiology of KD and thus provides new insights into the diagnosis and treatment of KD.

2 Result

2.1 Clinical features and laboratory data of KD and febrile

In this study, we collected six peripheral blood samples from three KD patients (before/after IVIG treatment). All KD patients were sensitive to IVIG and did not develop CAL. We also collected peripheral blood samples from three febrile patients, which included two pneumonia patients and one pharyngitis patient. Table 1 summarizes the clinical and laboratory information of KD and febrile patients.

Table 1

Table 1. Clinical features and laboratory data for of KD and febrile patients.

2.2 Immunological parameters in KD and febrile

To validate the immune disorders in KD patients, we measured immune-related indicators in 28 KD and 28 febrile patients (Table 2). There were 15 men and 13 women in both KD group and febrile group. There was no difference in age between KD and febrile (28.36 ± 17.64 vs. 28.82 ± 10.89 months, p>0.05). There was no difference between fever days in KD and febrile (5.00 ± 0.00 vs. 4.82 ± 0.98 days, p>0.05). The percentage (17.41 ± 4.88% vs. 23.45 ± 6.51%) and absolute value (0.55 ± 0.35 10⁹/L vs 0.74 ± 0.36 10⁹/L) of CD8+T cells in KD were significantly lower than those in febrile (p<0.05). The ratio of CD4/CD8 (2.43 ± 1.28 vs. 1.59 ± 0.69) in KD was significantly higher than it in febrile (p<0.05). The percentages and absolute values of CD3+T cells, CD4+T cells, CD16 + 56+/CD3− cells, and CD19+ cells in the two groups showed no differences (p>0.05). We also measured the serum levels of IgG, IgM, IgA, and IgE in KD and febrile. The levels of IgG (6.29 ± 1.86 g/L vs. 8.59 ± 2.24 g/L) and IgM (0.91 ± 0.34 g/L vs. 1.14 ± 0.34 g/L) in KD were significantly lower than those in febrile (p<0.05), whereas the levels of IgA and IgE in KD and febrile showed no difference (p>0.05).

Table 2

Table 2. Immune-related indicators in KD and febrile patients.

2.3 Single-cell transcription profiling of PBMCs in febrile and KD patients before/after treatment

To describe the general features of single-cell expression features of febrile and KD patients before/after IVIG treatment, we collected the PBMCs from three febrile patients (used as control) and three KD patients before and after IVIG treatment. We also used a set of published KD data (GSE168732) to compare with our single-cell sequencing result (20). After quality control and filter, the total number of detected cells was 69,300, including 23,072 cells from febrile patients, 19,725 cells from KD patients before IVIG treatment (KD_BT), and 26,503 cells from KD patients after IVIG treatment (KD_AT); the total number of detected cells from the KD dataset of GSE168732 was 39,114, including 16,003 cells from KD patients before IVIG treatment (KD_BT_GSE168732) and 23,111 cells from KD patients after IVIG treatment (KD_AT_GSE168732). To compare the scRNA-seq profiles between control and KD patients, we first integrated nine of our samples together and clustered the cells across samples according to their expression features (Figure 1A, Supplementary Figure S1A). We also integrated our febrile samples with six GSE168732 KD samples according to their expression features (Figure 1B, Supplementary Figure S1C). The detected cells could be classified into 12 major cell types covering B cells, CD4 T cells, CD8 T cells, CD14 monocytes (CD14 mono), CD16 monocytes (CD16 mono), dendritic cells (DC), erythrocytes (Eryth), gamma-delta cells (gdT), hematopoietic stem and progenitor cells (HSPCs), natural killer cells (NK), plasma blast cells (plasmablast), and platelets. The canonical gene markers were used to validate major cell types (Supplementary Figures S1B, D). Figure 1C shows the percentage of major cell types for febrile patients and KD patients before and after IVIG treatment, as well as the KD samples from GSE168732 on the individual level. The scRNA-seq data demonstrated that KD patients after treatment had an increased percentage of CD8 T cells (p<0.05) and CD4 T cells (p<0.05) and a decreased percentage of monocytes (p<0.05) compared with those before treatment, which were well established by previous studies (20–22). Moreover, the number of B cells also shows a quite variance among febrile and KD patients before/after treatment. In both our dataset and the GSE168732 dataset, there is a trend of more plasma blast cells in KD patients after treatment than KD patients before treatment (Figures 1A, B).

Figure 1

Figure 1. Single-cell profiling of PBMCs in our dataset and the GSE168732 dataset. (A) The integration single-cell profiling analysis of nine samples in our dataset, which are three febrile patients (Febrile), three KD patients before treatment (KD_BT), and three KD patients after treatment (KD_AT). (B) The integration single-cell profiling analysis of three samples in our dataset and six samples in the GSE168732 dataset, which are three febrile patients (Febrile), three KD patients before treatment (KD_BT_ GSE168732), and three KD patients after treatment (KD_AT_ GSE168732). (C) The proportion of different cell types in each individual sample. The inferred cell types are marked with different colors.

2.4 Overall expression features of all single cells in febrile and KD patients before/after treatment

We first examined the expression features of all single cells in febrile and KD patients before/after treatment. The differentially expressed genes (DEGs) were identified in KD patients before treatment using febrile and KD patients after treatment as background. There are 303 upregulated genes and 516 downregulated genes (p<0.05) detected in KD patients before treatment from our dataset and the GSE168732 dataset, both of which are classified into four comparison groups (Figures 2A, B). A total of 112 out of 303 upregulated genes and 130 out of 516 downregulated genes are shared by both datasets (Supplementary Tables S1, S2). Next, we performed gene ontology (GO) analysis on these shared upregulated and downregulated genes (Figures 2C, D). There is a common biological process of immunity found in shared upregulated and downregulated genes, although the genes behind this biological process are totally different in KD patients before treatment. The shared upregulated genes in KD patients before treatment are mainly involved in innate immunity and the extracellular region, whereas the shared downregulate genes in KD patients before treatment are mainly involved in immunoglobulin, adaptive immunity, and cell membrane. It shows that KD patients before treatment and febrile patients/KD patients after treatment rely on different immune mechanisms. The latter is more inclined to utilize secreted immunoglobulin as the measure of adaptive immune response. Notably, we found that the B-cell receptor signaling pathway is downregulated in KD patients before treatment, which proposes B-cell-related pathology for KD patients before treatment and is consistent with the previous study (20).

Figure 2

Figure 2. Expression analyses of all single-cells in our dataset and GSE168732 dataset. (A) Venn diagram of upregulated genes in all single cells for KD patients before treatment in our dataset and the GSE168732 dataset. (B) Venn diagram of downregulated genes in all single cells for KD patients before treatment in our dataset and the GSE168732 dataset. (C) GO term enrichment analysis of 71 common upregulated genes in all single cells for KD patients before treatment. (D) GO term enrichment analysis of 29 common downregulated genes in all single cells for KD patients before treatment.

2.5 Distinct expression features of B cells, plasma blast cells, dendritic cells, and HSPCs in KD patients before treatment

We further compared the expression features of major immune cell types and HSPCs between febrile patients/KD patients after treatment and KD patients before treatment. Most expression features among major immune cell types are similar to the overall expression features (Supplementary Tables S3, S4). The genes related to the B-cell receptor signaling pathway (IGHC3, IGCHM, IGHG4, IGHG1, IGLC3, GO:0050853) are also downregulated in B cells and plasma blast cells in KD patients before treatment (Figures 3A, B). It is noticed that the genes related to the biological processes of acetylation, cell cycle, and cell division (UBE2C, HSPD1, HSPE1, MT2A, MYC, GO: KW-0007, KW-0131, and KW-0132) are selectively upregulated in DCs and HSPCs in KD patients before treatment (Figures 3C, D). We further compared the HSPC expression profiles between febrile control and KD patients before treatment, since HPSCs give the rise to the other types of cells. There are 122 upregulated genes in HSPCs in our dataset and 81 upregulated genes in HSPCs in the GSE168732 dataset (Figure 4A). A total of 14 HSPCs upregulated genes (FOS, BLVRB, MCM7, RASD1, HIST1H3H, SAT1, JCHAIN, NKG7, HSPA1A, SELENOS, S100A9, S100A8, FOSB, TYROBP) are shared by both datasets. They are mainly involved in the KEGG pathway of the IL-17 signaling pathway (GO: hsa04657, FDR value = 0.008). GO analyses show that there are common cellular components (GO:0070062~extracellular exosome, GO:0005576~extracellular region, GO: KW-0964~Secreted) shared by KD before treatment from our dataset and the GSE168732 dataset, although the genes underlying these cellular components are various in two datasets (Figures 4B, C). Interestingly, the biological process of immune response (GO:0006955) is upregulated in the HSPCs of our dataset whereas the biological process of inflammatory response (GO:0006954) is upregulated in the HSPCs of the GSE168732 dataset. Both immune response and inflammatory response are an innate function of immune cells, which proposes a premature immune activation propensity in the HSPCs of KD.

Figure 3

Figure 3. Dot plot for B-cell signaling pathway genes and cycle cell genes in our dataset and the GSE168732 dataset. (A) The expression level of B-cell signaling pathway genes in B cells and plasma blast cells of our dataset. (B) The expression level of B-cell signaling pathway genes in B cells and plasma blast cells of GSE168732 dataset. (C) The expression level of acetylation, cell cycle, and cell division genes in HSPCs and DC cells of our dataset. (D) The expression level of acetylation, cell cycle, and cell division genes in HSPCs and DC cells of the GSE168732 dataset.

Figure 4

Figure 4. Expression analyses of HSPCs in our dataset and the GSE168732 dataset. (A) Venn diagram of upregulated genes in HSPCs for KD patients before treatment in our dataset and the GSE168732 dataset. (B) GO term enrichment analysis of 122 upregulated genes in HSPCs for KD patients before treatment of our dataset. (C) GO term enrichment analysis of 81 upregulated genes in HSPCs for KD patients before treatment of the GSE168732 dataset.

2.6 Single-cell trajectory analyses of PBMCs in febrile and KD patients before/after treatment

UBE2C, HSPD1, HSPE1, MT2A, and MYC are the genes known to participate in the biological processes of cell cycle and cell division (GO: KW-0007, KW-0131, and KW-0132). Their abnormal upregulation in the HSPCs of KD before treatment could influence their descendant cells. We used pseudo-time analysis to reconstruct cell developmental trajectory of PBMCs in febrile and KD patients before/after treatment. Since HSPCs are the starting point of blood cell development, they are used to root the cell developmental trajectory of each dataset. The pseudo-time analyses show that the PBMCs in febrile and KD patients after treatment can be divided into five states whereas there are only three states of PBMCs in KD patients before treatment (Figures 5A–E). It proposes that KD patients before treatment have a shortened cell development route and an ill cell differentiation outcome. The number of DEGs in the cell developmental trajectory detected by the Monocle 2 package are 15,364, 14,819, 15,827, 12,819, and 11,867 in febrile patients, KD patients before treatment, KD patients after treatment, KD patients before treatment of GSE168732, and KD patients after treatment of GSE168732, respectively. Human PBMCs are known to form different cell lineages. We use a set of classic PBMC cell markers to identify the cell lineage for each state (Figure 5F) (23). For febrile patients, five states of their PBMCs are myeloid lineage (state 1), lymphoid lineage 1 (state 2), mixed lymphoid lineage 1 and 2 (state 3), erythro-megakaryocytic lineage (state 4), and lymphoid lineage 2 (state 5) (Supplementary Figures S2A–E). For KD patients before treatment, three states of their PBMCs are mixed myeloid and erythro-megakaryocytic lineages (state 1), myeloid lineage (state 2), and mixed lymphoid lineage 1 and 2 (state 3) (Supplementary Figures S3A–D). For KD patients after treatment, five states of their PBMCs are mixed myeloid and erythro-megakaryocytic lineages (state 1), lymphoid lineage 1 (state 2), mixed lymphoid lineage 1 and 2 (state 3), mixed myeloid lineage and lymphoid lineage 2 (state 4), and mixed lymphoid lineage 1 and 2 (state 5) (Supplementary Figures S4A–F). For KD patients before treatment of GSE168732, three states are lymphoid lineage 1 (state 1), mixed myeloid and erythro-megakaryocytic lineages (state 2), and mixed lymphoid lineage 1 and 2 (state 3) (Supplementary Figures S5A–D). For KD patients after treatment of GSE168732 are mixed lymphoid lineages 1 and 2 (state 1), mixed lymphoid lineages 1 and 2 (state 2), lymphoid lineage 1 (state 3), mixed lymphoid lineage 1 and erythro-megakaryocytic lineage (state 4), and mixed lymphoid lineage 1 and myeloid lineage (state 5) (Supplementary Figures S6A–F). The cell lineage analyses show that only febrile patients have a clear cell lineage division whereas KD patients before/after treatment have more than one mixed lineage. Compared with KD patients before treatment, KD patients after treatment have more B cell lineages (lymphoid lineage 1). These results propose that B-cell developmental dysfunction is the potential etiology of KD at the cellular level.

Figure 5

Figure 5. Pseudo-time analysis of cell developmental trajectory in our dataset and the GSE168732 dataset. (A) The differentiation trajectory of all cells in febrile patients by states in our dataset. (B) The differentiation trajectory of all cells in KD patients before treatment by states in our dataset. (C). The differentiation trajectory of all cells in KD patients after treatment by states in our dataset. (D) The differentiation trajectory of all cells in KD patients before treatment by states in GSE168732 dataset. (E) The differentiation trajectory of all cells in KD patients after treatment by states in the GSE168732 dataset. (F) The canonical markers of five cell lineages for state 1 in febrile patients in our dataset.

2.7 Pseudo-time expression dynamic analyses of cell cycle-related genes and SPI1 in febrile and KD patients before/after treatment

We have observed that the genes related to cell cycle and cell division are selectively upregulated in DCs and HSPCs in KD patients before treatment. Single-cell trajectory analyses further show that KD patients before treatment have a shortened cell development route and an ill cell differentiation outcome. Thus, we investigate the expression dynamics of these genes based on pseudo-time analysis. We also added the SPI1 gene in pseudo-time expression dynamic analyses, because it is a key gene in controlling myeloid and B-lymphoid cell lineage development (24). The expression dynamics of UBE2C, HSPD1, HSPE1, MT2A, MYC, and SPI1 were plotted along cell developmental trajectory based on pseudo-time analysis. We found that all six genes demonstrate a different expression dynamic among febrile patients, KD patients before treatment, and KD patients after treatment (Figures 6A–E). HSPD1 and HSPE1 show an elevated expression level in the early cell development stage of KD patients before treatment in our dataset and the GSE168732 dataset, which are both repressed in KD patients after treatment. In febrile patients, MYC has an expression peak in the middle-late stage of cell development and SPI1 is highly expressed in the early and late stages of cell development. In KD patients before/after treatment, the expression peak of MYC disappears in the middle-late stage of cell development and SPI1’s early expression is almost repressed. The expression of MT2A and UBE2C are repressed in KD patients before treatment. Interestingly, these six genes show a different expression dynamic between febrile patients and KD patients after treatment. Except MT2A, their expressions are repressed in KD patients after treatment compared with febrile patients. It indicates that IVIG treatment probably restores the B-cell development dysregulation in KD patients through suppressing the expression of UBE2C, HSPD1, HSPE1, MYC, and SPI1.

Figure 6

Figure 6. Pseudo-time analyses of expression dynamics of six cell-cycle-related genes in our dataset and the GSE168732 dataset. (A) Expression dynamics of six cell-cycle related genes in febrile patients of our dataset. (B) Expression dynamics of six cell-cycle-related genes in KD patients before treatment of our dataset. (C) Expression dynamics of six cell-cycle-related genes in KD patients after treatment of our dataset. (D) Expression dynamics of six cell-cycle-related genes in KD patients before treatment of the GSE168732 dataset. (E) Expression dynamics of six cell-cycle-related genes in KD patients after treatment of the GSE168732 dataset.

2.8 The developmental path of B cells in febrile patients and KD patients before/after treatment

The overall profiling of PBMCs in febrile patients and KD patients before/after treatment shows that febrile and KD patients after treatment have the larger number of plasma blast cells than KD patients before treatment in both our and GSE168732 datasets. The overall expression features also demonstrated that the B-cell activation pathway is repressed in KD patients before treatment. In order to elucidate the possible etiological mechanisms of KD, it is worth investigating the detailed developmental path of B cells in three sample datasets. Thus, we extracted all HSPCs and B cells from each sample dataset and used them to perform pseudo-time analysis. B cells could be divided into four subgroups including B naïve, B intermediate, B memory, and plasma blast cells. The B-cell developmental trajectories display seven states in febrile patients, three states in KD patients before treatment, seven states in KD patients after treatment, one state in KD patients before treatment of GSE168732, and three states in KD patients after treatment of GSE168732 (Figures 7A–E). This result indicates that the B cells are fully developed in febrile patients and KD patients after treatment than in KD patients before treatment. The different types of B cells were plotted along their developmental trajectories. The difference between the B cells in our dataset and the GSE168732 dataset is that KD patients in our dataset have more states than KD patients in the GSE168732 dataset. In both our and GSE168732 datasets, KD patients before treatment have a more shortened B-cell development path than in KD patients after treatment. It indicates that most B cells in KD patients before treatment might be in a transitional state compared with those in KD patients after treatment. CD24 is a marker gene usually expressed in transitional B cells (25, 26). We extracted B naïve, B intermediate, and B memory cells for five data samples and analyzed the CD24 expression level in all of them. Its expression range is wider in KD patients before treatment in both our and GSE168732 datasets. Our further analysis shows that CD24 is actually a cell maker gene for most B cells in KD patients before treatment, whose expression level is significantly higher in B naïve, B intermediate, and B memory cells compared with KD patients after treatment in our and GSE168732 datasets (Figure 7F). Furthermore, the CD24 expression level is higher in our dataset than in the GSE168732 dataset (Wilcox tests, p<0.05). It might explain the different number of states in these two datasets. In all five data samples, plasma blast cells are always at the end of the B development path. It is consistent with the current knowledge of the B-cell development process. Thus, the insufficient number of plasma blast cells in KD patients before treatment is due to B-cell developmental dysfunction whose main feature is the elevated expression of CD24. It proposes that B cells in KD patients before treatment are mostly arrested in a transitional state. We plotted the expression trajectory of CD24 in pseudo-time analysis as well (Supplementary Figure S7). It shows a conspicuous peak in KD before the treatment dataset and has an expression value in KD before treatment of the GSE168732 dataset, which is consistent with our violin plot results that a portion of cells are transitional B cells. The expression of CD24 on pro-B cells is also a herald signal of B-cell selection and development (27). It suggests that most B cells in KD patients before treatment could enter into their next developmental stage with proper stimulation.

Figure 7

Figure 7. Pseudo-time analysis of cell developmental trajectory of HSPCs and B cells in our dataset and the GSE168732 dataset. (A) The differentiation trajectory of HSPCs and B cells in febrile patients by cell types in our dataset. (B) The differentiation trajectory of HSPCs and B cells in KD patients before treatment by cell types in our dataset. (C) The differentiation trajectory of HSPCs and B cells in KD patients after treatment by cell types in our dataset. (D) The differentiation trajectory of HSPCs and B cells in KD patients before treatment by cell types in the GSE168732 dataset. (E) The differentiation trajectory of HSPCs and B cells in KD patients after treatment by cell types in the GSE168732 dataset. (F) The violin plot of CD24 expression in the B cells in our dataset and the GSE168732 dataset.

3 Discussion

KD was initially reported in 1967. Despite tremendous efforts, the etiology of KD remains unclear. Numerous etiology theories have been proposed for KD. The infection hypothesis suggests that an infectious agent leads to activation of the immune system in a genetically susceptible child who eventually developed KD (26). This theory is supported by the apparent seasonality of KD, the similar clinical features with other infectious diseases, and the increased levels of inflammatory markers (27–29). However, none of these potential KD infectious agents have been confirmed. A large number of studies have revealed the significant activation of innate immune response and adaptive immune response in KD, which support the autoimmunity hypothesis. However, its low recurrence rate and the absence of a family history of autoimmune diseases in KD are inconsistent with the typical presentation of autoimmune disorders (1, 30, 31). In short, none of these theories of KD have been fully validated and they only partially account for the characteristics of KD. Until now, there is no consensus on its etiology.

Then, hygiene hypothesis was introduced to explain the etiology of KD. Epidemiological and immunological data support it as a possible KD etiology, and the development of B cell is a crucial factor in this hypothesis (12). In previous studies, the increase of B-cell number in peripheral blood of acute KD patients has been documented. The infiltration of oligoclonal plasma cells producing IgA into the vascular wall of KD patients was also reported (32). A series of genome-wide studies identified that several KD susceptibility genes are involved in B-cell development and function (16–18). These findings have suggested that B cell-related immunity is crucial in the development of KD. The overall single-cell expression features show that B-cell receptor signaling pathway-related genes are downregulated in KD before treatment whereas innate immunity response-related genes are upregulated. These results confirm that B-cell development is dysregulated in KD and such dysregulation could lead to excessive innate immune response. B cells play a key role in adaptive immunity. Their undergrowth would shift KD patients’ immune system to rely on innate immune response for self-defense (28). Unfortunately, such shift would lead to overreaction of the immune system and further invoke the immune system to attack benign antigens. Thus, our findings are consistent with the concept of hygiene hypothesis. We further identified that HSPCs in our dataset and the GSE168732 dataset overexpress the genes related to the IL-17 signaling pathway which are mainly involved in inflammation. HSPCs in both datasets exhibit a premature immune activation propensity, which suggests that excessive immune response starts very early in KD patients and B-cell development dysregulation could begin in the bone marrow. It is known that HSPCs are the cellular origin of all immune cells, and thus we infer that any abnormality in HSPCs could affect their downstream cells’ development. Since the repression of B-cell development is a major cellular symptom of KD, the abnormality in KD patients’ HSPCs should play a part in this symptom. We also analyzed the expression of B-cell markers and B-cell ligands and receptors in KD patients’ HSPCs. The results are shown in Supplementary Figure S8. It shows that KD patients’ HSPCs not only have a relatively high expression level of B-cell markers (CD19, CD24, MS4A1, CD79A, CD79B, CD27) but also have a high expression level of B-cell ligands (TNFSF13, CD40LG, TNFSF13B) and their receptors (TNFRSF17, CD40, TNFRSF13B). B-cell ligands are information molecules binding to B-cell receptors, which promote lymphoid progenitor cells to differentiate into B cells (29). Thus, the proliferating HSPCs in KD patients might lead to the repression of B-cell development in two ways: 1) their premature B-cell propensity leads to the maldevelopment of their downstream cells, especially B cells; 2) their relatively high expression of B-cell ligands and receptors send the differentiation information to themselves and the other cells and keep them in the B-cell state.

The analyses of expression features of major immune cell types in KD before treatment further reveal the possible genes behind B-cell development dysregulation. First, the B-cell receptor signaling pathway-related genes are downregulated in B cells and plasma blast cells in KD before treatment. Their low expression level proposes that the upstream activation signal is repressed for the B cells in KD patients and thus the B cells in KD patients have difficulties to differentiate into downstream cells such as plasma blast cells. Second, the selective upregulation of five cell cycle- and cell division-related genes in dendritic cells and HSPCs lead to produce a large number of premature B cells in KD before treatment. The untimely expression of cell cycle- and cell division-related genes generates the premature B cells which have poor expression of the B-cell receptor signaling pathway-related genes. They mutually contribute to the symptom of a large number of B cells and a small number of plasma blast cells observed in KD patients before treatment. The B-cell developmental dysplasia in KD patients before treatment is also proven by clinical data, which show that the IgG/M level is low in KD patients compared with febrile patients.

The pseudo-time analyses of KD before treatment in our and GSE168732 datasets has an impaired cell development trajectory compared with febrile and KD after treatment. Canonical cell markers show that febrile patients have four clearly differentiated states, which are myeloid lineage (state 1), lymphoid lineage 1 (state 2), erythro-megakaryocytic lineage (state 4), and lymphoid lineage 2 (state 5), but KD patients before treatment have mixed myeloid and erythro-megakaryocytic lineage and mixed lymphoid lineage 1 and lineage 2 (30). It suggests that B-cell development dysregulation in KD actually has a deep root in the very early cell lineage differentiation, which affect not only lymphoid lineage but also myeloid lineage. SPI1 is also known as hematopoietic transcription factor PU.1. It is a very important regulator gene in controlling gene expression during myeloid and B-lymphoid cell development (24). The gene expression dynamic along the cell development trajectory shows that SPI1 is highly expressed in the early and late stages of cell development in febrile patients, but its early expression is repressed in KD before treatment and its expression is almost inhibited in KD after treatment in our dataset. Ikaros is one of known upstream regulators of SPI1, which represses SPI1’s expression (31). It is a gene family which encodes five IKZF zinc finger genes, named from IKZF1 to IKZF5. It is still unclear whether only IKZF1 functions as the repressor of SPI1 or all five IKZF genes might have the same function (32). If assuming that all of them can repress the expression of SPI1, the combination of different IKZF expressions would contribute to the control of SPI1’s expression. We plot the expression of all five IKZF genes (Supplementary Figure S9). In our dataset, the expression of IKZF1 and IKZF2 are always higher in KD after the treatment group, especially in HSPCs, but this phenomenon is not obvious in the GSE168732 dataset. It suggests that the regulation mechanism of SPI1 could be different in our dataset and the GSE168732 dataset, which is a part of KD heterogenicity that needs to be further studied. Its repression in KD before treatment is likely influenced by the early expression of HSPD1 and HSPE1, but how HSPD1 and HSPE1 repress the expression of SPI1 needs to be studied in future. The overall expression feature analyses show that the hematopoietic cell lineage-related genes are upregulated in KD before treatment. It further confirms that the SPI1 expression is repressed in KD before treatment, because it controls hematopoietic cells to differentiate into myeloid and lymphoid lineages. In hematopoiesis, HSPCs first give the birth to common myeloid progenitor cells and common lymphoid progenitor cells (30). The formers further split into myeloid lineage and erythro-megakaryocytic lineage, and the latter further split into T lymphoid lineage (lymphoid lineage 1) and B lymphoid lineage (lymphoid lineage 2) (30). Thus, the repression of SPI1 expression in KD patients before treatment render them to have only three developmental states in which myeloid lineage vs. erythro-megakaryocytic lineage and lymphoid lineage 1 vs. lymphoid lineage 2 are unable to separate from each other.

The examination of B-cell developmental path further confirms that KD before treatment has an ill B-cell developmental trajectory as well. The B cells in KD before treatment has a lower state than those in febrile patients and KD after treatment. In both our dataset and the GSE168732 dataset, IVIG treatment facilitates the B cells of KD to develop into more trajectory states and thus produce more plasma blast cells to perform normal adaptive immunity related functions. The expression dynamic analyses show that the expressions of UBE2C, HSPD1, HSPE1, MYC, and SPI1 are repressed in KD after treatment. UBE2C, HSPD1, and HSPE1 participate in the biological process of the cell cycle (GO: KW-0131), whereas MYC and SPI1 are two proto-oncogenes (33, 34). Their repression could stop cells from proliferating and thus render B cells unable to enter to the next developmental fate. This speculation is supported by the high expression of CD24 in the B cells of KD before treatment, which is a marker for transitional B cells. The expression level of CD24 in the B cells of KD before treatment clearly demonstrates that they are arrested in a transitional state, but the CD24 expression is also downregulated after IVIG treatment. These results indicate that IVIG performs a general suppressive signal for KD patients’ immune system, which would turn down its hyperactivity after receiving IVIG by facilitating the developmental fate of B cells. IVIG treatment is effective for 70%–80% of KD patients. Thus, it is the right signal for dendritic cells and HSPCs in most KD patients. We further examined the expression of stress response genes (HSP90A1, HSPH1, HSPA5, HSP90B1, HSPD1, HSPA1A, HSPE1) and hematopoietic stem cell differentiation genes (SMPD3, CITED2, FOXC1, CDK6, AP2A2, PRKDC) in dendritic cells and the other myeloid cells (CD8 T and CD4 T cells). Supplementary Figure S10 shows that stress response genes have a high expression level in the dendritic cells of KD before the treatment dataset and CITED2 has a low expression level in the dendritic cells of KD before treatment in our dataset and the GSE168732 dataset. It indicates that dendritic cells of KD patients before treatment are in a stress mode and insufficient to guide the differentiation of the other immune cells, since dendritic cells function as the messengers between the innate and adaptive immune systems (35). After IVIG treatment, dendritic cells are relieved from stress mode and their messenger function is rehabilitated. This rehabilitative message is documented by KD patients’ immune system and becomes a permanent memory of immune system. It can explain the low recurrence rate of KD, but the evidence of rehabilitation still lies in the patient’s bone marrow.

We observed some variances between our dataset and the GSE168732 dataset, especially in B-cell developmental path analyses. Before treatment, the total B cells from the KD patients of our dataset have three states whereas the total B cells from the KD patients of the GSE168732 dataset have only one state, which suggests that there exist some subtle pathological differences among different KD patients. There are approximately 10%–20% of KD patients insensitive to IVIG treatment (36). These differences among KD population need to be investigated by future studies. Although whether lacking of common pathogens in living environment sets off KD still needs to be investigated, our result unequivocally shows that KD has an immune system always ready for battle, i.e., a hyperactive immune state that could be easily activated by any antigen. Unfortunately, such hyperactive immune system is based on the disordered developmental fate of major immune cells, especially for B cells. There were several shortcomings in our study. Firstly, the sample size was limited, which affected the statistical power of this study. Secondly, there are very few public single-cell data for KD and we only found one available for this study. The future studies on KD’s pathological mechanisms requires more public data and endeavors from the scientific community. Thirdly, this study was performed on PBMCs and might not reflect the local inflammatory responses developing in the coronary artery and the maturation process of hematopoietic cells in the bone marrow. Lastly, we only performed scRNA-seq in this study, which only reflected the single-cell transcriptome level of KD. Single-cell multi-omics sequencing can be considered to be applied in the further studies of KD.

4 Materials and methods

4.1 Participants and statistical analysis

All participants were recruited from Shanghai Children’s Hospital. All manipulations were approved by the Ethics Committee of Shanghai Children’s Hospital (IRB number: 2022R121). Guardians of the participants had provided their informed consent in the study. KD was diagnosed according to the diagnosis criteria established by the American Heart Association (1). Two-dimensional echocardiography was used to examine whether there were coronary artery lesions during the acute and convalescent disease phases. Febrile patients were respiratory system infectious diseases admitted during the same period with KD patients.

Inclusion criteria for the febrile group are as follows: (1) febrile patients (respiratory system infectious diseases) admitted during the same period with KD patients; (2) age matching with KD patients; (3) Fever duration matching with KD patients; (4) children who had not received steroids or immunosuppressive drugs within 14 days.

Exclusion criteria of all participants are as follows: (1) combined with autoimmune diseases; (2) combined with congenital malformations; (3) combined with metabolic diseases, etc.; (4) incomplete clinical data or inability to cooperate with the study. Peripheral blood samples collected from KD patients were on the 5th day after the onset of fever before IVIG treatment and 24 h after IVIG treatment and subsidence of fever. Peripheral blood samples were also collected from three febrile patients immediately after admission.

Analysis was conducted using SPSS 27.0. Normally distributed continuous variables were expressed as mean ± standard deviation (s), with two independent samples t-test for between-group comparisons. Non-normally distributed measures were represented as median (M) (quartile 1 [Q1], quartile 3 [Q3]), and the Mann–Whitney U-test was used for the group comparisons. Count data were presented as the number of patients, using a four-compartment table χ² test for between-group comparisons. P<0.05 was considered statistically significant.

4.2 Single-cell RNA sequencing and data analyses

4.2.1 scRNA-seq library construction

Peripheral blood samples (2 mL each sample) were collected from the participants. PBMCs of participants were isolated according to standard density gradient centrifugation methods by using the Ficoll-Paque medium. The cell viability should exceed 90%. The single-cell library was constructed using the 5′ Library Kits. The cell suspension was loaded onto a chromium single-cell controller (10x Genomics) to generate single-cell gel beads in the emulsion (GEMs) according to the manufacturer’s protocol. Lysis and barcoded reverse transcription of polyadenylated mRNA from single cells was performed inside each GEM. Post-RT-GEMs were cleaned up, and cDNA was amplified. The barcoded sequencing libraries were generated using the Chromium Next GEM Single Cell V(D)J Reagent Kits v1.1 (10x Genomics) and were sequenced as 2×150-bp paired-end reads on an Illumina NovaSeq platform.

4.2.2 Published data acquisition and selection

A set of published KD patients’ single-cell data was also used in our study to expand KD sample size (20). The published KD data could be retrieved from Gene Expression Omnibus (GEO) repository with the accession number of GSE168732 through the NCBI website. Because we used three febrile patients as control in this study, three KD patients’ data were selected to get a paired dataset with our control from total six KD patients’ data of GSE168732. The selection procedure was as follows. First, we removed the low-quality samples from GSE168732. KD patient 1 before treatment of GSE168732 has a very low number of feature RNAs (average 700 feature RNAs) compared with the other five GSE168732 samples (average 1,200 feature RNAs). KD patient 2 after treatment has a relatively high percentage of expressed mitochondrion genes (approximately 10% of mitochondrion genes per cell) compared with the other five GSE168732 samples (approximately 7% of mitochondrion genes per cell). Thus, KD patient 1 and KD patient 2 were removed from this study for low quality. Second, we removed the sample with the excessive number of cells. KD patient 5 before treatment of GSE168732 has 8,783 cells (vs. average 5,585 cells in the other five before treatment datasets) and KD patient 5 after treatment of GSE168732 has 10,201 cells (vs. average 8,043 cells in the other five after treatment datasets). KD patients 5 of GSE168732 has approximately 25%–50% more cells than the other five datasets of GSE168732 and is an outlier to be removed. Finally, KD patients 3, 4, and 6 of GSE168732 passed our selection criteria and were used in this study.

4.2.3 scRNA-seq data processing

CellRanger (Version 6.0.0) software was used to process the raw FASTQ files, align the sequencing reads to the GRCh38 reference transcriptome, and generate a filtered UMI expression profile for each cell. Raw gene expression matrices were read into R (Version 4.2.1) and converted to Seurat objects. The number of genes, UMI counts, and percentage of mitochondrial genes were examined to identify outliers. The following criteria were applied for quality control: total UMI count between 2,000 and 60,000, and mitochondrial gene percentage <10%. After removal of low-quality cells, the count matrix was normalized by the SCTransform method, which is based on a negative binomial regression model with regularized parameters. Then, all datasets from the individual samples were integrated using the “FindIntegrationAnchors” and “IntegrateData” functions in Seurat (Version 4.1.1) (37). We identified “anchors” among all individual datasets with multiple canonical correlation analysis (CCA) and used these “anchors” to create a batch-corrected expression matrix of all cells, which allowed the cells from different datasets to be integrated and analyzed. The supervised principal component analysis (SPCA) was performed to reduction, and the weighted nearest neighbor (wnn) graph-based clustering was used to identify cell clusters. The cell identities were determined with multimodal reference mapping in Seurat (Version 4.1.1) (37).

4.2.4 Differential expression and functional enrichment analysis

DEG analysis for each cell type on the sample level followed the recommendation of Bioconductor (38). The differential expression analysis was conducted between conditions by using Libra (Version 1.0.0) (39), which implements a total of 22 unique differential expression methods that can all be accessed from one function. We used “run_de” functions with the pseudo bulk approach, implementing the DESeq2 (Version 1.36.0) (40) with a likelihood ratio test. GO analysis was performed using a DAVID online resource (41). P-values were adjusted to FDRs. FDRs <0.05 were chosen as the cutoff criterion indicating a statistically significant difference.

4.2.5 Pseudo-time analysis of cell differentiation trajectories

Pseudo-time analysis of cell differentiation trajectories for each sample dataset was performed with R package Monocle 2 (42). The expression feature and inferred cell type for each sample data from the Seurat result was used to construct the cell dataset for the Monocle analysis pipeline. We used the Monocle built-in approach named “dpFeature” to detect the variable genes that define the cell’s differentiation. Its advantages are needing no prior biological knowledge and discovering important ordering genes from data itself. Dimension reduction was performed with 2 max components and “DDRTree” method. HSPCs and B cells in each sample dataset for pseudo-time analysis were extracted according to the cell identities from the Seurat result.

4.2.6 Cell lineage verification

PBMCs are known to be differentiated into different lineages. The different cell lineage is verified with canonical cell markers in this study. The canonical cell markers were retrieved from the CellMarker 2.0 database (http://bio-bigdata.hrbmu.edu.cn/CellMarker/) (23). We defined five cell lineages in this study which are HSPC and progenitor lineage with CD34, CD38, SOX4, and ITGA6 as marker genes; erythro–megakaryocytic lineage with PPBP, PF4, GP9, HBB, HBA1, and HBA2 as marker genes; myeloid lineage with CD14, FUT4, CD16, CST3, CEBPA, LYZ, CD1C, and FCER1A as marker genes; lymphoid lineage 1 (B cell lineage) with MZB1, CD27, CD79A, CD40LG, JCHAIN, MS4A1, CD24, and IGHM as marker genes; and lymphoid lineage 2 (T cell lineage) with GNLY, CD8, CD4, GZMB, KLRD1, KLRF1, CD3D, and NKG7 as marker genes. A custom R script based on the vioplot library was developed to plot and verify all selected markers for each lineage.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.

Ethics statement

The studies involving humans were approved by Ethics Committee of Shanghai Children’s Hospital (IRB number: 2022R121). The studies were conducted in accordance with the local legislation and institutional requirements. Written informed consent for participation in this study was provided by the participants’ legal guardians/next of kin.

Author contributions

QL: Writing – original draft, Writing – review & editing. ZW: Data curation, Methodology, Writing – original draft, Writing – review & editing. GD: Conceptualization, Data curation, Writing – original draft, Writing – review & editing. GL: Data curation, Writing – original draft, Writing – review & editing. LC: Data curation, Formal analysis, Writing – original draft, Writing – review & editing. QQ: Data curation, Writing – original draft, Writing – review & editing. SS: Data curation, Writing – original draft, Writing – review & editing. WL: Data curation, Writing – original draft, Writing – review & editing. XJ: Data curation, Writing – original draft, Writing – review & editing. MH: Resources, Writing – original draft, Writing – review & editing. LS: Conceptualization, Data curation, Formal analysis, Software, Writing – original draft, Writing – review & editing. TX: Conceptualization, Resources, Writing – original draft, Writing – review & editing. LX: Writing – original draft, Writing – review & editing, Funding acquisition, Resources, Supervision, Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Project administration, Software, Validation, Visualization.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This work is jointly supported by the National Natural Science Foundation of China (82170518), the Shanghai Science and Technology Committee research Funding (22Y11909700), and the Seventh Period Jinshan Medical Key Specialty A-level-Pediatrics (JSZK2023A04).

Acknowledgments

We would like to thank all the donors who contributed samples.

Conflict of interest

Author GL was employed by the company Daozhi Precision Medicine Technology Shanghai Co., LTD.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2024.1438640/full#supplementary-material

Supplementary Figure S1 | Integrated single-cell profiling of PBMCs for our dataset and GSE168732 dataset. (A) Integrated single-cell profiling of PBMCs for our dataset. (B) Expression of canonical gene markers for each cell type in our dataset based on integration analysis. (C) Integrated single-cell profiling of PBMCs for GSE168732 dataset. (D) Expression of canonical gene markers for each cell type in GSE168732 dataset based on integration analysis. The inferred cell types are marked with different colors.

Supplementary Figure S2 | Pseudo-time analysis of cell developmental trajectory in febrile patients of our dataset. (A) The differentiation trajectory of all cells in febrile patients by cell types in our dataset. (B) The canonical markers of five cell lineages for state 1 in febrile patients in our dataset. (C) The canonical markers of five cell lineages for state 2 in febrile patients in our dataset. (D) The canonical markers of five cell lineages for state 3 in febrile patients in our dataset. (E) The canonical markers of five cell lineages for state 4 in febrile patients in our dataset. (F) The canonical markers of five cell lineages for state 5 in febrile patients in our dataset.

Supplementary Figure S3 | Pseudo-time analysis of cell developmental trajectory in KD patients before treatment of our dataset. (A) The differentiation trajectory of all cells in KD patients before treatment by cell types in our dataset. (B) The canonical markers of five cell lineages for state 1 in KD patients before treatment in our dataset. (C) The canonical markers of five cell lineages for state 2 in KD patients before treatment in our dataset. (D) The canonical markers of five cell lineages for state 3 in KD patients before treatment in our dataset.

Supplementary Figure S4 | Pseudo-time analysis of cell developmental trajectory in KD patients after treatment of our dataset. (A) The differentiation trajectory of all cells in KD patients after treatment by cell types in our dataset. (B) The canonical markers of five cell lineages for state 1 in KD patients after treatment in our dataset. (C) The canonical markers of five cell lineages for state 2 in KD patients after treatment in our dataset. (D) The canonical markers of five cell lineages for state 3 in KD patients after treatment in our dataset. (E) The canonical markers of five cell lineages for state 4 in KD patients after treatment in our dataset. (F) The canonical markers of five cell lineages for state 5 in KD patients after treatment in our dataset.

Supplementary Figure S5 | Pseudo-time analysis of cell developmental trajectory in KD patients before treatment of GSE168732 dataset. (A) The differentiation trajectory of all cells in KD patients before treatment by cell types in GSE168732 dataset. (B) The canonical markers of five cell lineages for state 1 in KD patients before treatment in GSE168732 dataset. (C) The canonical markers of five cell lineages for state 2 in KD patients before treatment in GSE168732 dataset. (D) The canonical markers of five cell lineages for state 3 in KD patients before treatment in GSE168732 dataset.

Supplementary Figure S6 | Pseudo-time analysis of cell developmental trajectory in KD patients after treatment of GSE168732 dataset. (A) The differentiation trajectory of all cells in KD patients after treatment by cell types in GSE168732 dataset. (B) The canonical markers of five cell lineages for state 1 in KD patients after treatment in GSE168732 dataset. (C) The canonical markers of five cell lineages for state 2 in KD patients after treatment in GSE168732 dataset. (D) The canonical markers of five cell lineages for state 3 in KD patients after treatment in GSE168732 dataset. (E) The canonical markers of five cell lineages for state 4 in KD patients after treatment in GSE168732 dataset. (F) The canonical markers of five cell lineages for state 5 in KD patients after treatment in GSE168732 dataset.

Supplementary Figure S7 | Pseudo-time analysis of expression dynamics of CD24 in cell developmental trajectory in our dataset and GSE168732 dataset. (A) Expression dynamics of CD24 in febrile patients of our dataset. (B) Expression dynamics of CD24 in KD patients before treatment of our dataset. (C) Expression dynamics of CD24 in KD patients after treatment of our dataset. (D) Expression dynamics of CD24 in KD patients before treatment of GSE168732 dataset. (F) Expression dynamics of CD24 in KD patients after treatment of GSE168732 dataset.

Supplementary Figure S8 | Dot plot for HSPC marker, B cell markers, B-cell ligands and receptors in our dataset and GSE168732 dataset. (A) The expression of HSPC marker, B cell markers, B-cell ligands and receptors in HSPCs, B cells and plasma blast cells in our dataset. (B) The expression of HSPC marker, B cell markers, B-cell ligands and receptors in HSPCs, B cells and plasma blast cells in GSE168732 dataset.

Supplementary Figure S9 | Violin plot for five IKZF genes in our dataset and GSE168732 dataset. (A) The expression of five IKZF genes in HSPC, DC, plasma blast and B cells in our dataset. (B) The expression of five IKZF genes in HSPC, DC, plasma blast and B cells in GSE168732 dataset.

Supplementary Figure S10 | Dot plot for stress response and hematopoietic stem cell differentiation genes in our dataset and GSE168732 dataset. (A) The expression of stress response and hematopoietic stem cell differentiation genes in DC, CD8 T and CD4 T cells in our dataset. (B) The expression of stress response and hematopoietic stem cell differentiation genes in DC, CD8 T and CD4 T cells in in GSE168732 dataset.

References

1. McCrindle BW, Rowley AH, Newburger JW, Burns JC, Bolger AF, Gewitz M, et al. Diagnosis, treatment, and long-term management of kawasaki disease: A scientific statement for health professionals from the American heart association. Circulation. (2017) 135:e927–e99. doi: 10.1161/CIR.0000000000000484

PubMed Abstract | Crossref Full Text | Google Scholar

2. Wright VJ, Herberg JA, Kaforou M, Shimizu C, Eleftherohorinou H, Shailes H, et al. Diagnosis of Kawasaki disease using a minimal whole-blood gene expression signature. JAMA Pediatr. (2018) 172:e182293. doi: 10.1001/jamapediatrics.2018.2293

PubMed Abstract | Crossref Full Text | Google Scholar

3. Manlhiot C, Mueller B, O’Shea S, Majeed H, Bernknopf B, Labelle M, et al. Environmental epidemiology of Kawasaki disease: Linking disease etiology, pathogenesis and global distribution. PloS One. (2018) 13:e0191087. doi: 10.1371/journal.pone.0191087

PubMed Abstract | Crossref Full Text | Google Scholar

4. Nagata S. Causes of Kawasaki disease-from past to present. Front Pediatr. (2019) 7:18. doi: 10.3389/fped.2019.00018

PubMed Abstract | Crossref Full Text | Google Scholar

5. Rowley AH. Kawasaki disease: novel insights into etiology and genetic susceptibility. Annu Rev Med. (2011) 62:69–77. doi: 10.1146/annurev-med-042409-151944

PubMed Abstract | Crossref Full Text | Google Scholar

6. Strachan DP. Hay fever, hygiene, and household size. BMJ. (1989) 299:1259–60. doi: 10.1136/bmj.299.6710.1259

PubMed Abstract | Crossref Full Text | Google Scholar

7. Burgner D, Carter K, Webster R, Kuijpers TW. Kawasaki disease, childhood allergy and the hygiene hypothesis. Pediatr Allergy Immunol. (2011) 22:751. doi: 10.1111/j.1399-3038.2011.01184.x

PubMed Abstract | Crossref Full Text | Google Scholar

8. Lee KY, Han JW, Lee JS. Kawasaki disease may be a hyperimmune reaction of genetically susceptible children to variants of normal environmental flora. Med Hypotheses. (2007) 69:642–51. doi: 10.1016/j.mehy.2006.12.051

PubMed Abstract | Crossref Full Text | Google Scholar

9. Singh S, Vignesh P, Burgner D. The epidemiology of Kawasaki disease: a global update. Arch Dis Child. (2015) 100:1084–8. doi: 10.1136/archdischild-2014-307536

PubMed Abstract | Crossref Full Text | Google Scholar

10. Holman RC, Belay ED, Christensen KY, Folkema AM, Steiner CA, Schonberger LB. Hospitalizations for Kawasaki syndrome among children in the United States, 1997-2007. Pediatr Infect Dis J. (2010) 29:483–8. doi: 10.1097/INF.0b013e3181cf8705

PubMed Abstract | Crossref Full Text | Google Scholar

11. Nakamura Y. Kawasaki disease epidemiology and the lessons from it. Int J Rheum Dis. (2018) 21:16–9. doi: 10.1111/1756-185X.13211

PubMed Abstract | Crossref Full Text | Google Scholar

12. Holman RC, Christensen KY, Belay ED, Steiner CA, Effler PV, Miyamura J, et al. Racial ethnic differences in the incidence of Kawasaki syndrome among children in Hawaii. Hawaii Med J. (2010) 69:194–7.

PubMed Abstract | Google Scholar

13. Kuijpers TW, Wiegman A, van Lier RA, Roos MT, Wertheim-van Dillen PM, Pinedo S, et al. Kawasaki disease: a maturational defect in immune responsiveness. J Infect Dis. (1999) 180:1869–77. doi: 10.1086/315111

PubMed Abstract | Crossref Full Text | Google Scholar

14. Lee JK. Hygiene hypothesis as the etiology of Kawasaki disease: dysregulation of early B cell development. Int J Mol Sci. (2021) 22:1–20. doi: 10.3390/ijms222212334

Crossref Full Text | Google Scholar

15. Johnson TA, Mashimo Y, Wu JY, Yoon D, Hata A, Kubo M, et al. Association of an IGHV3-66 gene variant with Kawasaki disease. J Hum Genet. (2021) 66:475–89. doi: 10.1038/s10038-020-00864-z

PubMed Abstract | Crossref Full Text | Google Scholar

16. Lee YC, Kuo HC, Chang JS, Chang LY, Huang LM, Chen MR, et al. Two new susceptibility loci for Kawasaki disease identified through genome-wide association analysis. Nat Genet. (2012) 44:522–5. doi: 10.1038/ng.2227

PubMed Abstract | Crossref Full Text | Google Scholar

17. Onouchi Y, Ozaki K, Burns JC, Shimizu C, Terai M, Hamada H, et al. A genome-wide association study identifies three new risk loci for Kawasaki disease. Nat Genet. (2012) 44:517–21. doi: 10.1038/ng.2220

PubMed Abstract | Crossref Full Text | Google Scholar

18. Fan X, Zhou Y, Guo X, Xu M. Utilizing single-cell RNA sequencing for analyzing the characteristics of PBMC in patients with Kawasaki disease. BMC Pediatr. (2021) 21:277. doi: 10.1186/s12887-021-02754-5

PubMed Abstract | Crossref Full Text | Google Scholar

19. Geng Z, Tao Y, Zheng F, Wu L, Wang Y, Wang Y, et al. Altered monocyte subsets in Kawasaki disease revealed by single-cell RNA-sequencing. J Inflammation Res. (2021) 14:885–96. doi: 10.2147/JIR.S293993

Crossref Full Text | Google Scholar

20. Wang Z, Xie L, Ding G, Song S, Chen L, Li G, et al. Single-cell RNA sequencing of peripheral blood mononuclear cells from acute Kawasaki disease patients. Nat Commun. (2021) 12:5444. doi: 10.1038/s41467-021-25771-5

PubMed Abstract | Crossref Full Text | Google Scholar

21. Matsubara T, Ichiyama T, Furukawa S. Immunological profile of peripheral blood lymphocytes and monocytes/macrophages in Kawasaki disease. Clin Exp Immunol. (2005) 141:381–7. doi: 10.1111/j.1365-2249.2005.02821.x

PubMed Abstract | Crossref Full Text | Google Scholar

22. Furukawa S, Matsubara T, Yabuta K. Mononuclear cell subsets and coronary artery lesions in Kawasaki disease. Arch Dis Child. (1992) 67:706–8. doi: 10.1136/adc.67.6.706

PubMed Abstract | Crossref Full Text | Google Scholar

23. Hu C, Li T, Xu Y, Zhang X, Li F, Bai J, et al. CellMarker 2.0: an updated database of manually curated cell markers in human/mouse and web tools based on scRNA-seq data. Nucleic Acids Res. (2023) 51:D870–D6. doi: 10.1093/nar/gkac947

PubMed Abstract | Crossref Full Text | Google Scholar

24. Oikawa T, Yamada T, Kihara-Negishi F, Yamamoto H, Kondoh N, Hitomi Y, et al. The role of Ets family transcription factor PU.1 in hematopoietic cell differentiation, proliferation and apoptosis. Cell Death different. (1999) 6:599–608. doi: 10.1038/sj.cdd.4400534

Crossref Full Text | Google Scholar

25. Sims GP, Ettinger R, Shirota Y, Yarboro CH, Illei GG, Lipsky PE. Identification and characterization of circulating human transitional B cells. Blood. (2005) 105:4390–8. doi: 10.1182/blood-2004-11-4284

PubMed Abstract | Crossref Full Text | Google Scholar

26. Zhou Y, Zhang Y, Han J, Yang M, Zhu J, Jin T. Transitional B cells involved in autoimmunity and their impact on neuroimmunological diseases. J Trans Med. (2020) 18:1–12. doi: 10.1186/s12967-020-02289-w

Crossref Full Text | Google Scholar

27. Mensah FFK, Armstrong CW, Reddy V, Bansal AS, Berkovitz S, Leandro MJ, et al. CD24 expression and B cell maturation shows a novel link with energy metabolism: potential implications for patients with myalgic encephalomyelitis/chronic fatigue syndrome. Front Immunol. (2018) 9:2421. doi: 10.3389/fimmu.2018.02421

PubMed Abstract | Crossref Full Text | Google Scholar

28. Hara T, Nakashima Y, Sakai Y, Nishio H, Motomura Y, Yamasaki S. Kawasaki disease: a matter of innate immunity. Clin Exp Immunol. (2016) 186:134–43. doi: 10.1111/cei.12832

PubMed Abstract | Crossref Full Text | Google Scholar

29. Ketchum C, Miller H, Song W, Upadhyaya A. Ligand mobility regulates B cell receptor clustering and signaling activation. Biophys J. (2014) 106:26–36. doi: 10.1016/j.bpj.2013.10.043

PubMed Abstract | Crossref Full Text | Google Scholar

30. Monga I, Kaur K, Dhanda SK. Revisiting hematopoiesis: applications of the bulk and single-cell transcriptomics dissecting transcriptional heterogeneity in hematopoietic stem cells. Briefings Funct Genomics. (2022) 21:159–76. doi: 10.1093/bfgp/elac002

Crossref Full Text | Google Scholar

31. Carotta S, Wu L, Nutt SL. Surprising new roles for PU.1 in the adaptive immune response. Immunol Rev. (2010) 238:63–75. doi: 10.1111/j.1600-065X.2010.00955.x

PubMed Abstract | Crossref Full Text | Google Scholar

32. Spooner CJ, Cheng JX, Pujadas E, Laslo P, Singh H. A recurrent network involving the transcription factors PU.1 and Gfi1 orchestrates innate and adaptive immune cell fates. Immunity. (2009) 31:576–86. doi: 10.1016/j.immuni.2009.07.011

PubMed Abstract | Crossref Full Text | Google Scholar

33. Finver SN, Nishikura K, Finger LR, Haluska FG, Finan J, Nowell PC, et al. Sequence analysis of the MYC oncogene involved in the t(8;14)(q24;q11) chromosome translocation in a human leukemia T-cell line indicates that putative regulatory regions are not altered. Proc Natl Acad Sci United States America. (1988) 85:3052–6. doi: 10.1073/pnas.85.9.3052

Crossref Full Text | Google Scholar

34. Ray D, Culine S, Tavitain A, Moreau-Gachelin F. The human homologue of the putative proto-oncogene Spi-1: characterization and expression in tumors. Oncogene. (1990) 5:663–8.

PubMed Abstract | Google Scholar

35. Banchereau J, Steinman RM. Dendritic cells and the control of immunity. Nature. (1998) 392:245–52. doi: 10.1038/32588

PubMed Abstract | Crossref Full Text | Google Scholar

36. Piram M, Darce Bello M, Tellier S, Di Filippo S, Boralevi F, Madhi F, et al. Defining the risk of first intravenous immunoglobulin unresponsiveness in non-Asian patients with Kawasaki disease. Sci Rep. (2020) 10:3125. doi: 10.1038/s41598-020-59972-7

PubMed Abstract | Crossref Full Text | Google Scholar

37. Hao Y, Hao S, Andersen-Nissen E, Mauck WM 3rd, Zheng S, Butler A, et al. Integrated analysis of multimodal single-cell data. Cell. (2021) 184:3573–87 e29. doi: 10.1016/j.cell.2021.04.048

PubMed Abstract | Crossref Full Text | Google Scholar

38. Amezquita RA, Lun ATL, Becht E, Carey VJ, Carpp LN, Geistlinger L, et al. Orchestrating single-cell analysis with Bioconductor. Nat Methods. (2020) 17:137–45. doi: 10.1038/s41592-019-0654-x

PubMed Abstract | Crossref Full Text | Google Scholar

39. Squair JW, Gautier M, Kathe C, Anderson MA, James ND, Hutson TH, et al. Confronting false discoveries in single-cell differential expression. Nat Commun. (2021) 12:5692. doi: 10.1038/s41467-021-25960-2

PubMed Abstract | Crossref Full Text | Google Scholar

40. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. (2014) 15:550. doi: 10.1186/s13059-014-0550-8

PubMed Abstract | Crossref Full Text | Google Scholar

41. Sherman BT, Hao M, Qiu J, Jiao X, Baseler MW, Lane HC, et al. DAVID: a web server for functional enrichment analysis and functional annotation of gene lists (2021 update). Nucleic Acids Res. (2022) 50:W216–W21. doi: 10.1093/nar/gkac194

PubMed Abstract | Crossref Full Text | Google Scholar

42. Qiu X, Mao Q, Tang Y, Wang L, Chawla R, Pliner HA, et al. Reversed graph embedding resolves complex single-cell trajectories. Nat Methods. (2017) 14:979–82. doi: 10.1038/nmeth.4402

PubMed Abstract | Crossref Full Text | Google Scholar