eQTL Highlights the Potential Role of Negative Control of Innate Immunity in Kawasaki Disease

Introduction

Kawasaki disease (KD) as an acute, self-limited disease predominately affects children under 5 years.¹ KD is featured by the development of autoimmune-like polyangiitis and can cause damage to the coronary arteries (including coronary artery aneurysms and long-term cardiovascular sequences).² The incidence of reported KD is increasing over year, and the Asian population has a higher incidence. Moreover, the higher incidence of KD in siblings and twins suggests a genetic predisposition to KD.³ Some studies have investigated the KD-related genes from available databases and revealed some KD-related genes are mainly related to the immune response.⁴

Genome wide association study (GWAS) provides a tool for understanding the genetic changes in the KD. The single nucleotide polymorphisms (SNP) of inositol 1,4,5-trisphosphate 3-kinase C (ITPKC) have been found to be associated with the KD susceptibility and coronary risk in American and Japanese children. ITPKC can negatively regulate T-cell activation through a calcium signaling pathway, and the functional polymorphism of ITPKC may increase T-cell activation and cytokine expression.⁵ The coatomer protein complex beta-2 subunit (COPB2), a factor associated with T cell activity, is also related to the KD susceptibility.⁶ In humoral immunity, the polymorphism of IgG receptor gene FCGR2A has also been found to be associated with the KD susceptibility in European and Asian populations. Immune cells (such as dendritic cells [DCs], neutrophils, monocytes/macrophages) express Fc fragment of IgG receptor IIa (FCGR2A) on the surface, which can transmit activation signals into cells when the receptor binds to the immune complex, leading to the immunosuppression in some cases.⁷ The self-limitation of KD, the therapeutic effectiveness of intravenous immunoglobulin (IVIG), and the functions of genes identified by GWAS suggest the important role of immunoregulation in the pathogenesis of KD.

GWAS mainly focuses on genome-wide SNPs information. However, few studies have focused on the roles of immunoregulation-related genes in the pathogenesis of KD. In the target capture sequencing, a specific probe is designed to capture the genome region of interest and then used to hybridize with the whole genome DNA; subsequently, the DNA fragments of target genome region are enriched and processed for sequencing. It is cost-effective and capable of finding disease-related genes as well. Meanwhile, the expression quantitative trait loci (eQTLs) is a key to map SNPs that affect gene expression.⁸

The present study was to investigate whether eQTL highlights the potential role of negative control of innate immunity in the KD. First, the genes related to four immune regulation pathways were investigated, and their expressions were compared between KD patients and healthy children. Then, the inflammatory process was simulated to further confirm the differentially expressed genes between the above two groups. Thereafter, the eQTLs analysis was employed to obtain the SNPs affecting gene expression. After stimulation, the secretion of cytokines was also evaluated. Our results showed some of genes and cytokines associated with immunoregulation were involved in the pathogenesis of KD.

Patients and Methods

Patients and Healthy Controls

The patients in this study were enrolled from Shanghai Children’s Hospital between December 2015 and October 2017. All patients were recruited after informed consent was signed by their parents or guardians; this study was conducted in accordance with the Declaration of Helsinki and approved by the Ethics Committee of Shanghai Children’s Hospital. The KD was diagnosed according to the diagnostic criteria for KD proposed by the American Heart Association in 2017.⁹ Children who had not received IVIG therapy were excluded. In the control group, healthy children who received routine physical examination were recruited from the Department of Outpatient. The healthy controls had no prior immunological disease and heart disease.

The enrolled 309 Han Chinese children in this study were unrelated. All of the 190 patients with KD were included in the next-generation sequencing (NGS) phase and 20 in the peripheral blood mononuclear cells (PBMCs) stimulation phase. Of the 119 healthy controls, 99 were included in the NGS phase, and 20 in the PBMC stimulation phase.

Target Capture Sequencing and Data Analysis

Blood samples (2–4 mL) were obtained from KD patients before IVIG treatment. Genomic DNA was extracted from 300 μL of whole peripheral blood using Genomic DNA Extraction Kit (Cat. No. DP329, TIANGEN Bioscience Beijing) according to the manufacturer’s instructions. Two sources of information (one involved in 4 pathways and the other involved SNP previously published) were used to design hybridization probes in this study. First, 512 candidate genes among 4 pathways including Toll-like receptor signaling pathway, cytokine receptor interaction, TGF-beta signaling pathway and T cell receptor signaling pathway were selected. The capture region included all exons and the adjacent 15-bp intron of these genes. Then, 472 GWAS hotspots related to the KD according to GeneReviews (NCBI) were extracted. Sequences ±50 bp flanking the risk loci were included in this analysis. The final custom capture array was synthesized by NimbleGen, Roche. Target enrichment was performed using the custom capture array, and the enriched target fragments were subsequently sequenced on the Illumina HiSeq X10 platform.

The raw reads were mapped to the human reference genome (hg19) using the Burrows-Wheeler alignment (BWA v0.7.15).¹⁰ Then, Picard (v1.135) was used to process for PCR duplicates (http://broadinstitute.github.io/picard/). Variants in the target regions were detected using GATK (v3.7) HaplotypeCaller algorithm¹¹ and SNPs with “QD < 2.0 || FS > 60.0 || MQ < 40.0 || MQRankSum < −12.5 || ReadPosRankSum < −8.0” were excluded. Subsequently, functional annotation of filtered variants was performed using the ANNOVAR pipeline.¹²

Collection of Peripheral Blood Mononucleated Cells (PBMCs)

3 mL of blood was collected from 20 KD patients and 20 healthy children and then mixed with Ficoll-Paque PLUS (Cat. No.17544202, General Electric Co., USA), followed by centrifugation at 1000 g for 30 min at 4°C. The interphase PBMCs were collected and washed twice with RPMI 1640 medium (Cat. No.11875093, Life Technology Co., USA) including 10% fetal bovine serum (FBS) (Cat. No.10270106, Life Technology Co., USA). PBMCs were resuspended in plasma from the same person with 10% RPMI (No.11875093 Life Technology Co., USA) and then stored in the liquid nitrogen. The PBMCs were processed for further detection within a month.

LPS Stimulation and Immunoassays

Aliquots were thawed and centrifuged at 1000 rpm for 10 min at 24°C to precipitate cells. The cells were collected and incubated in RPMI 1640 supplemented with 10% FBS for 4 h at 37°C. The cells were collected after filtration and then re-suspended with RPMI at 1×10⁶/mL. The PBMCs were cultured at 1×10⁶ cells/mL in an environment with 5% CO₂ at 37°C for 24 h.

To induce inflammation, lipopolysaccharides (LPS) from Escherichia coli 0111:B4 (Cat. No. L4391, Sigma Aldrich, USA) was added to the cells at 1 μg/μL, followed by incubation for 22 h. Then, the supernatant and cells were collected independently. The expression of 13 cytokines (IL-8, IL-1β, IL-2, IL-6, IL-10, IL-4, IL-5, IL-12P70, IL-17A, IL-18, TNF-α, INF-γ and SCD25) was detected with the Luminex 200 ProcartaPlex Mix&Match Human-13-plex kit (Affymetrix eBioscience, UK). These 13 cytokines have been previously reported to be associated with KD pathogenesis and play important roles in four immunoregulation pathways.

RNA Sequencing and Data Analysis

Total RNA was extracted from the cells using Trizol reagent and then reversely transcribed into cDNA using the PrimeScript^TM Master mix (Cat. No. RR036a TaKaRa). The RNA library was constructed using NEBNEXT^® UltraTM II RNA Library Prep Kit for Illumina^® kit and submitted to Illumina HiSeq X10 platform.

The data analysis of RNA sequencing was conducted according to previously reported.¹³ Briefly, raw reads from each sample were quality controlled by Q30. After 3’ adaptor-trimming and removing low quality reads, the high-quality trimmed reads were mapped to the human reference genome hg19 (ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/grch37.tar.gz) with HISAT v2.1.0. Then, StringTie (version 1.3.3b) was used for transcript assembly and estimating the expression level of each gene with a gene annotation file (ftp://ftp.ensembl.org/pub/release-75//gtf/homo_sapiens/Homo_sapiens.GRCh37.75.gtf.gz). The Fragments Per Kilobase per Million mapped reads (FPKM) was employed to determine gene expression levels. Finally, differentially expressed genes (DEGs) and transcripts were detected by R packages: Ballgown and DEseq. A value of P<0.05 was considered statistically significant.

Variant Filtering and eQTL Analysis

In order to evaluate the impact of filtered SNPs on the KD, association analysis was performed by using the R package SNPassoc. SNPassoc estimates the numbers, percentages, odds ratios and P-value for the likelihood ratio test of association. A value of P < 0.05 was considered statistically significant.

To further confirm the SNPs associated with KD, an eQTL analysis was conducted to identify the cis eQTLs using FastQTL (v2.184) in 20 KD samples. If the SNP location was within 1 Mb upstream or downstream of a gene’ transcription starting site, this cis eQTL would be computed. FastQTL proposes an efficient permutation procedure to control for multiple testing and an adjusted P-value can be estimated. A significance threshold was set at 0.05, and therefore eQTLs with P-value <0.05 were considered statistically significant. The eQTL with the lowest FDR was selected for each gene to represent the best gene/eSNP pair.

Statistical Analysis

The Ingenuity Pathway Analysis (IPA) was carried out in order to analyze the molecular function, tox function and disease disorder of DEGs. The cytokine expression levels were analyzed using SPSS version 23.0 (IBM) and a value of two-sided P < 0.05 was considered statistically significant.

Results

Target Capture Sequencing of KD Patients

This study focused on the exon regions of 512 genes in the 4 KEGG pathways, including T cell receptor signaling pathway, Toll-like receptor signaling pathway, cytokine–cytokine interaction pathway and TGF-β signal pathway, and 472 risk loci of KD reported in previous GWAS studies (Supplementary Table S1). A total of 190 KD patients were recruited, including 34 patients with coronary artery lesions and 26 patients with IVIG resistance. The targeted regions were sequenced with an average depth of 440X (Supplementary Table S2).

A total of 12,470 SNPs (Supplementary File S1) were found in all samples. The annotation of SNPs showed that they consisted of 2247 nonsynonymous SNPs, 41 stopgain, 1 stoploss and others that do not affect the gene function (Supplementary Tables S3 and S4). To confirm the risk SNPs, SNPassco was employed to investigate the association of selected SNPs with KD.

A total of 2247 SNPs affected gene expression. Burden test showed there were 7 genes SNV sets that were closely associated with KD (P < 0.05), including 216 SNVs mutations (Table 1).

Table 1 The Number of SNVs Included in Differentially Expressed Genes Between KD Patients and Controls

RNA Sequencing of PBMCs After LPS Stimulation

In the present study, LPS was used to stimulate PBMCs to simulate the process of superantigen activation. After the stimulation, the PBMCs were separated from 20 KD patients (including 8 patients with coronary artery lesions and 1 patient with IVIG resistance) and 20 controls. Total RNA was extracted, transcriptome data were compared to the reference genome, and 10 KD patients and 10 controls were selected for RNAseq data with a comparison ratio of more than 90%. The ballgown was used to screen out the genes with significant differences in the transcriptome data between KD patients and controls (Table 2). IPA software was used to analyze the physiological functions of DEGs and their roles in clinical disease processes.

Table 2 Significant Differentially Expressed Genes Between KD Patients and Controls After RNAseq

IPA annotation showed these DEGs were related to the immune cell trafficking and cardiovascular system development and function (Figure 1). Then, the functions of genes with the most significant differences were further investigated. RUNX3 gene is related to the differentiation of regulatory T cells (Tregs), which is involved in the immune regulation.¹⁴ Results showed DEGs were mainly related to the connective tissue disorders and inflammatory diseases, which were consistent with the clinical manifestations of KD (Figure 2).

Figure 1 Enrichment analysis of DEGs related to the physiological development and function in KD patients and controls.

Figure 2 Enrichment analysis of DEGs in disease disorder of KD patients and controls.

eQTL Associated DNA and RNA Analysis

SNV sites in 20 KD patients that identified by the above-mentioned targeted sequencing showed correlation with the DEGs in the transcriptome sequencing data. A total of 356 SNVs were found to be associated with 313 DEGs (Table 3).

Table 3 eQTL Analysis for the Examination of DEGs

Then, 313 DEGs were subjected to IPA analysis. In terms of physiological functions, they were involved in the connective tissue development, cardiovascular system development and cell-mediated immune response (Figure 3). In terms of effects on the diseases, they were related to the connective tissue diseases, dermatological diseases and inflammatory diseases, which are consistent with the clinical phenotype of mucocutaneous damage in the KD at acute phase (Figure 4). The genes with the most significant changes in their expression (such as PIK3CD, ITPR3 and PRKCZ) also play an important role in the immune regulation.^15–17

Figure 3 DEGs associated with SNVs involved in physiological development and function.

Figure 4 DEGs associated with SNVs function in disease order.

Expression of Immunoregulation Related DEGs Among 40 Samples

After expression analysis and eQTL associated analysis, four immunoregulation-related DEGs (RUNX3, PIK3CD, ITPR3 and PRKCZ) were identified. Further examination showed their expression significantly decreased in the KD patients as compared to controls (Figure 5).

Figure 5 Expression of immunoregulation related DEGs decreased in the KD patients.

Cytokine Analysis After LPS Stimulation of PBMCs

The expression of 13 cytokines, including proinflammatory cytokines, anti-inflammatory cytokines, and cytokines previously reported to be associated with KD, was measured in the supernatant of cells after stimulation. Although the expression of immunoregulation-related cytokines IL-4 (Z=−2.39, P = 0.02), SCD25 (Z=−2.79, P = 0.02), and IL-10 (Z=−3.10, P < 0.01) significantly increased after LPS stimulation, their expression was still lower in the KD group than in the control group (Figure 6C,G and K). The expression of other proinflammatory cytokines (such as TNF-α, INF-γ and IL-6) was comparable between KD patients and controls (Figure 6E,L and M).

Figure 6 Expression of cytokines in the PBMCs from KD patients and controls before and after LPS stimulation (A–M).

Discussion

KD is a complex disease. Some studies have shown that the pathogenesis of KD is related to the immunological dysfunction. The genetic susceptibility is activated by antigens, causing immune activation and inflammatory reaction, which then produce corresponding clinical symptoms and vascular injury.¹⁸ However, the self-limitation in KD suggests the activation of anti-inflammatory pathways to regulate the immune system.

The mechanism underlying the IVIG-regulated inflammatory response in KD is still poorly understood. Studies have shown that, in vitro IVIG can negatively regulate T cell receptor signaling pathway, thereby inhibiting the proliferation and function of T lymphocytes.¹⁹ Thus, the present study aimed to investigate whether these genes associated with anti-inflammatory pathways were also involved in the pathogenesis of KD. Similarly, in Toll-like receptor pathway, stimulation of TLR2/TLR6 expression on mouse DCs induces their differentiation into resistant DCs which then secrete IL-10 to promote T-cell differentiation into Tregs.²⁰ Tregs production in mice is also induced by synergistic activation of T cells and smad3 nuclear factors in the TGF-β signaling pathway.^21,22 These may help reduce the inflammatory response. There is evidence showing that the cytokine storm is present in the acute phase of KD, and cytokine interactions also influence the inflammatory response.^23,24 In the present study, the chip covering genes from those four pathways was used. After sequencing, results showed the expression of some genes related to immune regulation was significantly different between KD patients and healthy controls. In the inflammatory response secondary to LPS stimulation, the DEGs between two groups mainly concentrated on the immune system and cardiovascular system. After eQTL association analysis, the genes that had mutation sites affecting their expression concentrated on those two systems. All of them were associated with vasculitis and the inflammatory response in the coronary arteries of KD.

RUNX3 was differentially expressed between KD group and control group. Although no evidence has indicated that its expression is affected by the mutation site, but, as a member of the RUNX gene, Rudra et al²⁵ found RUNX expression was up-regulated in the Treg cells, while that of Foxp3 (a gene regulating Treg) decreased when RUNX was silenced. Thus, RUNX may affect the number and function of T cell subsets. In tumor-induced CD8⁺ Treg cells, RUNX3 enhances the promoter activity of Foxp3.²⁶ In this study, our results showed that RUNX3 expression decreased in the KD group, which may affect the immunomodulatory function of KD and be related to the KD pathogenesis. The targeted capture sequencing showed the expression of PIK3CD, ITPR3 and PRKCZ in the KD group decreased and was affected by the mutant loci. Studies have shown that the PIK3CD family is involved in the regulation of cell proliferation, differentiation, apoptosis and glucose transport.^27–29 In the primary immunodeficiency disease, there is PIK3CD gene mutation, which leads to the T lymphocyte injury and aging, and compromises the immune regulation.¹⁷ ITPR3 receives the signals from second messenger ITP. ITP can be phosphorylated by ITPKC and negatively regulate the T cell signaling pathway by regulating Ca⁺/nuclear factors. Huang et al reported that SNP of ITPR3 gene (rs2229634) was related to KD in Taiwan children with concomitant coronary artery lesions (CAL).³⁰ This study showed that this SNV site was related to KD. In addition, ITPR3 gene is involved in the regulation of T lymphocyte apoptosis, can activate T lymphocytes and plays an important role in the progression of KD.¹⁵ It is also related to a variety of autoimmune diseases, such as systemic lupus erythematosus and rheumatoid arthritis.³¹ PRKCZ affects PKCζ expression. As an isozyme of PKC family, PKCζ can regulate T cell division, and affect subsequent T cell function.¹⁶ Pro-mitogen-activated protein kinases (MAPKs) also play a key role in the cytokine secretion in the Tregs via a PKC-dependent manner. In neonates with immature immunity or patients with immune deficiency, the defective PKCζ expression and IL-12 production may be two main events related to the T cell function.³² PKCζ may affect the secretion of cytokines in T cells, which also affects the susceptibility to autoimmune disease and is involved in the disease progression. The reduced expression of these genes may be associated with the hyperactive immune response and compromised immunoregulation in the KD.

The expression of these cytokines was further detected. Results showed, among 13 cytokines, the expression of IL-10 and SCD25 (IL-2R) related to immunoregulation was significantly lower in the KD group after LPS stimulation than in the control group. On the contrary, some pro-inflammatory cytokines have been previously reported to be related to the KD, but there was no significant difference in the expression of these cytokines. This may be ascribed to the difference between in vitro and in vivo experiments. Many studies have suggested that IL-10 is an important cytokine in the pathogenesis of KD.³³ The IL-10 gene not only showed different SNV sites between KD group and control group, but also had decreased protein expression. IL-10 plays an important role in the immunoregulation and immunotherapy through Treg cells.³⁴ SCD25 (IL-2R) is an important cytokine related to the Treg cells and effector T cells.³⁵ As a part of the IL-2 receptor, SCD25 (IL-2R) can recognize T cell activity and has become a biomarker of some immune diseases and tumors.³⁶ Treg cells can secrete IL-2 which together with SCD25 (IL-2R) plays an immunomodulatory role.³⁷ The decreased SCD25 (IL-2R) expression after LPS stimulation indicates that the decreased activity of this pathway compromises the regulatory effect on the KD. There were no significant differences in other proinflammatory cytokines (INF-γ, TNF-α and IL-6) after LPS stimulation. Although previous studies have reported that these cytokines are significantly elevated in a variety of inflammatory diseases including KD, no studies have shown the dynamic changes of these cytokines in the KD at acute stage, and these cytokines have not yet become specific biological markers for KD diagnosis. Therefore, we speculate that the immunomodulatory related pathways and cytokines may become a breakthrough for exploring KD pathogenesis and biological markers for the diagnosis of KD.

There were still limitations in the present study. Although some immunoregulatory molecules were found to be involved in the pathogenesis of KD, they mainly concentrated on the immune system and cardiovascular system. However, the DEGs affected by the SNVs have never been reported to be associated with KD, and the specific relationship between the hypothesized genes and cytokines has not been revealed. The genes related to anti-inflammatory pathways have no relationships with clinical phenotypes, disease severity, coronary artery lesions and alterations after IVIG therapy. Therefore, the mechanism by which these genes act in the pathogenesis of KD should be further investigated. In addition, the sample size was small, which limits the expansion of results.

Conclusions

Our study shows some genes and cytokines associated with immunoregulation are related to the pathogenesis of KD, which provides evidence for the development of therapies for KD.