Aberrant peripheral immune responses in acute Kawasaki disease with single-cell sequencing

Kawasaki disease (KD) is the most common cause of acquired heart disease in children in developed countries. Although diverse immune aberrance was reported, a global understanding of immune responses underlying acute KD was lacking. Based on single-cell sequencing, we profiled peripheral blood mononuclear cells from patients with acute KD before and after intravenous immunoglobulin therapy and from healthy controls. Most differentially expressed genes were derived from monocytes, with upregulation of immunoglobulin receptors, complement and receptors and downregulation of MHC class II receptors before therapy. The percentage of B cells was significantly increased before therapy and rapidly returned to normal after therapy. There was also an increased abundance of B-cell receptors with IGHA and IGHG after therapy, accompanied by massive oligoclonal expansion. The percentage of CD8 T cells was remarkably decreased during acute KD, especially the subset of effector memory CD8 T cells. All lymphocyte compartments were characterized by underexpressed interferon response pathways before therapy. The identification of unique innate and adaptive immune responses suggests potential mechanisms underlying pathogenesis and progression of KD.


Introduction
Kawasaki disease (KD) is an acute, systemic febrile illness and vasculitis of childhood that can lead to coronary artery lesions (CALs) 1 . This disease predominantly affects children who are younger than 5 years and has become the most common cause of acquired heart disease among children in many developed countries 2 . The etiology of KD remains elusive, and diagnosis continues to depend on principal clinical features including fever, rash, conjunctivitis, changes in the oral mucosa and in the extremities, and cervical lymphadenopathy 3 . Because the signs and symptoms of KD resemble those of other childhood febrile illnesses, prompt diagnosis of KD is still challenging. High-dose intravenous immunoglobulin (IVIG) within the first 10 days after fever onset is the standard therapy for KD, which remarkably reduces the rate of CALs 4 . However, the mechanism of IVIG in the treatment of KD is unknown. Approximately 10% to 20% of KDs are IVIG resistant and are at increased risk for CALs 3 .
The most widely favored theory of KD etiology is that it is caused by a ubiquitous infectious agent in a small subset of genetically predisposed children 5 . Genome-wide wide association studies (GWAS) have identified a number of susceptibility genes for KD 6 . Many of the susceptibility genes including ITPKC, CASP3, FCGR2A, CD40 and MHC class II are related to immune functions. Nonetheless, the causative agents initiating the disease have been still widely debated over 50 years. Various pathogens have been proposed as the trigger including superantigen toxin, Epstein-Barr virus, coronavirus and retrovirus, but none has been confirmed by subsequent studies 7 .
Despite the uncertainty of the cause, activation of the immune system provides important evidence for its pathogenesis. Plasma levels of pro-inflammatory cytokines such as TNF-α, IL-1β, and IFN-γ are elevated during the acute phase of KD 8 . Innate immune cells including neutrophils and monocytes are elevated in the peripheral blood and have been identified in the arterial wall early in the disease 9 . Activation of peripheral blood lymphocytes seems also responsible for the development of KD, though their roles remain controversial 10 . In addition, antigen-specific IgA plasma cells and CD8 T cells infiltrate inflamed tissues, implying an immune response to an intracellular pathogen entering through the respiratory tract 11 . Recently, a high incidence of Kawasaki-like disease has been reported during the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) epidemic 12,13 , suggesting a possible link between coronavirus infection and KD.
The immune system is composed of numerous cell types with varying states. Most of previous studies performed immunophenotyping based on flow cytometry, which is constrained to a few selected cell types and markers. Several studies performed transcriptome analyses of KD on bulk cell populations 14,15 , but the heterogeneity of the immune system cannot be resolved. Over the past few years, the revolution in single-cell sequencing has enabled an unbiased quantification of gene expression in thousands of individual cells, which provides an more efficient tool to decipher the immune system in human diseases 16 . In this study, we profiled peripheral blood mononuclear cells (PBMCs) isolated from acute KDs at single-cell resolution. Our results revealed global and dynamic immune responses unique to each cell compartment before and after IVIG therapy.

Study design and single-cell profiling of PBMCs
We collected eight fresh peripheral blood samples derived from four patients with acute KD (P1-P4, Supplementary Table S1). The patients were diagnosed according to the criteria proposed by the American Heart Association 3 . Three patients met criteria for complete KD and one for incomplete KD. For each patient, the first blood sample was taken on the 5th days after the onset of fever before IVIG therapy. The second sample was obtained at 24 hours after completion of IVIG therapy and subsidence of fever (Supplementary Table S2). All patients were IVIG sensitive responders and had not developed CALs. We also collected fresh peripheral blood samples from three age-matched healthy donors as controls (H1-H3, Supplementary Table S1).
We used the 10× Genomics platform for single-cell RNA sequencing (scRNA-seq) of PBMCs isolated from the samples. The total number of detected cells passing quality control were 57,881, including 19,319 cells for patients before IVIG therapy, 23,874 cells after therapy, and 14,688 cells for healthy controls (Supplementary Table S3).
We also performed single-cell B cell receptor sequencing (scBCR-seq) and single-cell  Table S4 and Table S5).
Based on the scRNA-seq profiles, we clustered the cells across samples and visualized them in two-dimensional space ( Figure 1A). Major cell compartments comprising PBMCs could be characterized by canonical marker genes ( Figure 1B CD68, FCGR3A or CD16, 7.55%) and plasmacytoid dendritic cells (pDCs) (LILRA4, 0.23%). There were also some residual erythrocytes (HBB, 2.49%) and megakaryocytes (PPBP, 0.83%) mixed in the PBMCs, which were removed in the following analyses.

Overall dynamics of cell abundance
We compared the percentage of each cell compartment between KDs before and after IVIG therapy, as well as healthy controls (Figure 2A). To verify the scRNA-seq results, we also retrieved a clinical dataset of large sample size, which comprised flow-cytometric estimation of peripheral blood lymphocytes from 125 acute KDs before IVIG therapy ( Figure 2B). The scRNA-seq data demonstrated that the B cell abundance increased significantly in KDs before therapy, compared with those after therapy (P = 0.01) and in healthy controls (P = 0.04, two-sided t-test, Figure 2A). The increase of B cells in acute KD was observed in several previous studies 10, [17][18][19] , and our flow-cytometric dataset also confirmed a significantly higher levels of B cells than the reference range in both percentages (P = 4.74×10 -20 , Figure 2B) and absolute counts (P = 1.59×10 -7 , two-sided t-test, Supplementary Figure S1). The scRNA-seq data also showed a significantly lower abundance of CD8 T cells in KDs than healthy controls (P = 0.02, two-sided t-test), though no significant difference were found before and after IVIG therapy (Figure 2A). The decrease of CD8 T cells was consistent with previous reports 10, [17][18][19] , which could be confirmed by our flow-cytometric dataset in both percentages (P = 2.84×10 -34 , Figure 2B) and absolute counts (P = 5.00×10 -13 , two-sided t-test, Supplementary Figure S1).
A reduction of NK cells was reported in acute KD 18,20,21 and can be verified by our flow-cytometric dataset (P = 6.70×10 -19 , two-sided t-test, Figure 2B and Supplementary Figure S1). The tendency of NK cell alternations was consistent with the scRNA-seq data, though statistical significance was lacking ( Figure 2A). There remained a controversy in literature regarding whether the level of CD4 T cells was increased or decreased in acute KD 10,[17][18][19] . While our flow-cytometric dataset showed a significant increase of CD4 T cells in percentages of lymphocytes (P = 2.57×10 -6 , two-sided t-test, Figure 2B), they did not differed significantly in absolute counts compared with the reference range (Supplementary Figure S1). We did not find clear tendency of CD4 T cell alternations in the scRNA-seq data ( Figure 2A).
It was reported that the abundance of monocytes was decreased following IVIG therapy 10,22,23 . However, one patient (P1) in the scRNA-seq data showed remarkable increase of monocytes after therapy (Figure 2A). To exclude any technical bias, we retrieved the routine blood test dataset of the patients, as well as other 80 acute KDs before and after IVIG therapy ( Figure 2C and Supplementary Figure S2). The dataset confirmed a significant decrease of monocytes after therapy in general (P = 2.35×10 -7 , two-sided t-test), whereas P1 displayed a reverse tendency as in the scRNA-seq data.
It was reported that CD16 monocytes were more closely associated with many inflammatory conditions 24 , including acute KD 10,23 . The abundance of CD16 monocytes remained significantly lower in patients after therapy than healthy controls (P = 0.04, two-sided t-test), but we did not find significant difference in patients before therapy and healthy controls ( Figure 3B).
To identify more genes involved in the immune responses, we performed differential expression analysis for each cell type. The most differentially expressed genes (DEGs, false discovery rate [FDR] < 0.1) resulted from monocytes among all PBMC compartments (Supplementary Figure S4). Gene ontology (GO) and KEGG pathway 25 analyses for the DEGs suggested that they were significantly enriched in immunoglobulin receptors, cytokine receptors, complement and receptors and MHC class II receptors (FDR < 0.1, Supplementary Figure S5). Most DEGs encoding immunoglobulin receptors were upregulated in patients before IVIG therapy compared with those after therapy and healthy controls, including FCGR3A, FCGR3B, FCGR2B and FCEG1G ( Figure 3C). The reduced expression of FCGR3A and FCGR2B after IVIG therapy in monocytes were in agreement with previous studies 26,27 . Most DEGs of complement and receptors were also upregulated before therapy, including C1QA, C1QB, CFD, C3AR1 and C5AR1 ( Figure 3C). In contrast, most DEGs encoding MHC class II receptors were downregulated in acute KD than healthy controls regardless of therapy, including HLA-DQA1, HLA-DQB1, HLA-DRA, HLA-DPB1 and HLA-DMB ( Figure 3C). The expression patterns of DEGs encoding cytokine and chemokine receptors were more diverse. While some were upregulated in patients before therapy including IFNAR2, IL3RA, IL10RB and CX3CR1, others were overexpressed after therapy or in healthy controls including IL13RA1, IL6ST, CXCR3 and CXCR4 ( Figure 3C).
We then performed the gene set enrichment analysis (GSEA) 28 to detect changes of pathway activity, which did not depend on predefined DEGs but coordinated changes of functionally related genes. In comparison with patients after therapy or healthy controls, a large number of pathways related to cytokine signaling and inflammatory responses were upregulated before therapy (FDR < 0.1, Figure 3D). Especially, TNF signaling via NF-kB, complement and coagulation pathway were activated in both comparisons. The activation of NF-kB in monocytes of acute KD was previously reported 29 .

B cell subsets and BCR repertoires
We re-clustered B cells and identified three distinct B cell subsets ( Figure 30 . After IVIG therapy, the proportion of plasma cells were increased compared with that before therapy, although the result just reach the statistical significance level (P = 0.05, two-sided t-test, Figure 4B). As few DEGs could be detected in B cells, we performed GSEA to assess the pathway activity (FDR < 0.1, Figure 4C). We found that interferon response pathways were underexpressed before IVIG therapy compared with those after IVIG therapy and in healthy controls. In addition, many pathways related to inflammatory responses were upregulated after therapy, which was in contrast to the cases in monocytes. Pathways associated with cell cycles were upregulated before therapy including G2M checkpoint and E2F targets ( Figure 4C), which might be responsible for the high abundance of B cells in acute KD.
We analyzed the scBCR-seq data to explore dynamics of BCR repertoires during acute KD. The abundance of IGHA and IGHG were significantly elevated after IVIG therapy (P = 0.03, two-sided t-test, Figure 4D and Supplementary Figure S7), which implicated isotype switching of BCRs from IgM/IgD to IgG/IgA resulted from B cell activation. We then examined the clonality of BCRs based on unique VDJ sequences for each patients. In comparison with patients before therapy, the percentage of clonal BCRs (clonal size ≥ 3) was greatly increased after therapy ( Figure 4E). Particularly, in every patient, we could observe remarkable oligoclonal expansions after therapy. The significantly expanded clones were dominated by the isotype IgA and IgG (FDR < 0.05, Fisher's exact test, Figure 4E), suggesting that antibody-mediated responses might be involved in the recovery of acute KD. The IGHV gene usage of the expanded clones was diverse among the patients ( Figure 4E), indicating specific B cell reactions rather than responses to a superantigen. We also did not find obvious bias in IGHV usage across patients and healthy controls (Supplementary Figure S8).
Oligoclonal IgA plasma cells infiltrate vascular tissue in acute KD were reported in previous studies 31,32 , and synthetic antibodies were designed to detect pathogens of KD 33,34 . Our analyses suggested that the antibodies might also be identified from peripheral blood of patients recovered from acute KD.

T and NK cell subsets and TCR repertoires
We subdivided CD4 T, CD8 T and NK cells based on expression of canonical genes, respectively ( Figure 5A). Both CD4 T and CD8 T cells comprised three subsets were performed for peripheral blood, but the signals were overwhelmed by the activation of innate immune system due to the dominate abundance of neutrophils 14,15 .
In this study, we used single-cell sequencing to dissect the complex immune responses of PBMCs in acute KD. We observed remarkable temporal changes in the cell abundance before and after IVIG therapy, which were largely consistent with previous studies and our flow-cytometric dataset. In particular, the abundance of B There were several shortcomings in our study. The sample size for single-cell analysis was limited, which affected the statistical power in differential abundance and differential expression analysis. Because all patients diagnosed with acute KD would be treated with IVIG, the intrinsic immune responses to potential pathogens were mixed with the immunomodulatory effects of IVIG. For example, further efforts should be made to distinguish clonally expanded BCRs specific to unknown pathogens and exogenous immunoglobulins after therapy. Also, more patients with various clinical presentation are needed to determine the relationship between immune responses of different cell types and disease outcomes.

Patients
All patients were recruited with informed consent from Shanghai Children's Hospital between December 2019 and June 2020. The study was approved by the ethics committee of the Shanghai Children's Hospital. The diagnosis of KD was made by using the criteria proposed by the American Heart Association 3 , including fever of unknown etiology lasting for ≥ 5 days, and the presence of five associated symptoms (conjunctivitis, oral changes, extremity changes, rash, cervical lymphadenopathy).
Complete KD was diagnosed if at least four of the symptoms were fulfilled, otherwise incomplete KD. The patients received high-dose IVIG (1 g/kg per day) for two consecutive days combined with oral aspirin after diagnosed with KD. CALs were regularly monitored by two-dimensional echocardiography.

Single-cell preparation and sequencing
For each donor, 2 ml venous blood was collected in EDTA tubes and transferred to the

Differential expression and functional enrichment analysis
We performed differential expression analysis for each cell compartment on the sample level following the tutorial of Bioconductor 50 . A pseudo-bulk expression profile was generated by summing counts together for all cells with the same combination of compartment and sample. Then the differential expression analysis was conducted between conditions by using edgeR (v3.30.3) 51 , which modelled the overdispersed count data based on a negative binomial generalized linear model. We

BCR and TCR data analysis
The scBCR-seq and scTCR-seq data were assembled by the Cell Ranger VDJ pipeline (v3.0.1, 10x Genomics), and reference sequences in IMGT 54 were fetched for annotation. Only cells with productive and paired chains (IGH and IGL/IGK for BCRs, TRA and TRB for TCRs) were preserved. For cells with more than one consensus sequence detected for the same chain type, the one with higher UMIs was chosen. Clonotypes represented by CDR3 sequences and gene usages were explored by immunarch (v0.6.5) 55 . The abundance of a clonotype was compared between preand post-therapy samples by the Fisher's exact test. A clonotype must occur ≥ 3 times in at least one sample to be considered. Clonotypes and expression profiles were matched based on the same cell barcode.

Data Availability
The raw sequence data reported in this paper have been deposited in the National