Serum exosomal miR-328, miR-575, miR-134 and miR-671-5p as potential biomarkers for the diagnosis of Kawasaki disease and the prediction of therapeutic outcomes of intravenous immunoglobulin therapy.

The present study was conducted to screen serum exosomal microRNAs (miRNAs) for the early diagnosis of Kawasaki disease (KD) and to investigate their underlying mechanisms by analyzing microarray data under accession numbers GSE60965 [exosomal miRNA, including three pooled serum samples from 5 healthy children, 5 patients with KD and 5 patients with KD following intravenous immunoglobulin (IVIG) therapy] and GSE73577 (mRNA, including peripheral blood mononuclear cell samples from 19 patients with KD prior to and following IVIG treatment) from the Gene Expression Omnibus database. Differentially expressed miRNAs (DE-miRNAs) and genes (DEGs) were identified using the Linear Models for Microarray data method, and the mRNA targets of DE-miRNAs were predicted using the miRWalk 2.0 database. The functions of the target genes were analyzed using the Database for Annotation, Visualization and Integrated Discovery (DAVID). As a result, 65 DE-miRNAs were identified with different expression patterns between the healthy children and patients with KD and between patients with KD and patients with KD following IVIG therapy. The target genes of 15 common DE-miRNAs were predicted. Following overlapping the target genes of DE-miRNAs with 355 DEGs, 28 common genes were identified and further screened to construct a network containing 30 miRNA-mRNA regulatory associations. Of these associations, only miR-328-spectrin α, erythrocytic 1, miR-575-cyclic AMP-responsive element-binding protein 5/b-1,4-galactosyltransferase 5/WD repeat and FYVE domain-containing 3/cystatin-A/C-X-C motif chemokine receptor 1/protein phosphatase 1 regulatory subunit 3B, miR-134-acyl-CoA synthetase long chain family member 1/C-type lectin domain family 1 member A and miR-671-5p-tripartite motif containing 25/leucine rich repeat kinase 2/kinesin family member 1B/leucine rich repeat neuronal 1 were involved in the negative regulation of gene expression. Functional analysis indicated that the identified target genes may be associated with inflammation. Accordingly, serum exosomal miR-328, miR-575, miR-134 and miR-671-5p may act as potential biomarkers for the diagnosis of KD and the prediction of outcomes of the IVIG therapy by influencing the expression of inflammatory genes.

Introduction

Kawasaki disease (KD) is a common, acute, systemic vasculitis that occurs in children <5 years old in Asian populations, with estimated incidence rates of 264.8, 134.4, 66.24 and 71.9 per 100,000 children in Japanese (1), Korean (2), Taiwanese (3) and Chinese (4) populations, respectively. KD predominantly affects small- to medium-sized vessels, including coronary arteries (5). If diagnosis is delayed and KD is left untreated, coronary artery lesions may develop in 25% of patients, which, in turn, increases the formation risk of coronary artery aneurysms (CAAs) and, subsequently, can induce coronary artery thrombosis, myocardial infarction or even sudden death (6,7). Intravenous immunoglobulin (IVIG) infusion is an effective, first-line therapy for KD (8). However, ~30% of cases have been reported to be unresponsive to IVIG, and additional therapy is required (9). If patients who are unresponsive to IVIG therapy are not identified in a timely manner, CAAs may still develop (10). Therefore, the early diagnosis of KD and the prediction of the therapeutic outcomes of IVIG are important issues.

At present, the diagnosis and evaluation of therapeutic effects primarily depend on clinical symptoms (including fever for ≥ five days, conjunctivitis, erythema of the lips and oral mucosa, extremity swelling, rash and cervical lymphadenopathy) (11), and ultrasonic imaging results (including QT interval dispersion) (12). However, these presentations overlap with other febrile illnesses in childhood (including epistaxis, scarlet fever, bovillae or juvenile idiopathic arthritis) and, therefore, are not specific (13). In addition, accumulating evidence indicates that the development of KD may be associated with the abnormal activation of the immune system and inflammation (14,15). Several inflammatory cytokines in the serum [including tumor necrosis factor α (TNF-α), interleukin (IL)-6, IL-17, IL-22 and IL-23] have been suggested to be biomarkers for the diagnosis of KD and IVIG therapy (16–18). However, certain studies reported that inflammatory markers are also non-specific and usually decreased in patients who are unresponsive to initial therapy (19,20). Therefore, novel biomarkers for the diagnosis of KD and for the prediction of therapeutic effects are required.

MicroRNAs (miRNAs) are endogenous, small, non-coding RNAs (~19–22 nucleotides in length) that serve roles in the regulation of diverse physiological and pathological processes through the complementary binding of target genes in the 3′untranslated region, for cleavage or translational repression (21). Furthermore, miRNAs can be present in the serum due to resistance against ribonuclease digestion and serum miRNAs are stable with consistent levels among individuals of the same species (21). Therefore, serum miRNAs that regulate the expression of inflammatory genes may serve as potential biomarkers for the diagnosis of KD and the prediction of therapeutic outcomes. This hypothesis has been supported by recent studies. Using miRNA microarray assays, Yun et al (22) demonstrated that miR-200c and miR-371-5p were significantly upregulated in children with KD, compared with controls. Further, Zhang et al (23) observed that serum miR-200c and miR-371-5p levels were significantly increased in patients with KD who were unresponsive to IVIG therapy compared with patients with KD who exhibited a good response to IVIG therapy. The combination of serum miR-200c and miR-371-5p exhibited high predictive values for the diagnosis of patients with KD [area under the curve (AUC)=0.95] and for those with an excellent IVIG response (AUC=0.97) (23). Furthermore, Yun et al (22) predicted that miR-200c and miR-371 may be involved in KD by targeting the regulation of a series of inflammatory response genes. In addition to their intracellular presence, miRNAs can also be enveloped into nanoparticles, termed exosomes, which maintain the integrity of miRNAs and transfer the miRNAs to recipient cells, influencing their phenotypes (24). Therefore, serum exosomal miRNAs may represent important biomarkers for various diseases (24). To date, exosomal miRNAs have been demonstrated to be useful in the early diagnosis of several cancers (25,26). The studies on KD and the IVIG therapy remain rare, except for the study by Jia et al (27), which identified four exosomal miRNA biomarkers, including miR-1246, miR-4436b-5p, miR-197-3p and miR-671-5p. The roles and mechanisms through which these miRNAs may act were not investigated.

Several studies reported that the addition of serum induces the production of pro-inflammatory cytokines (IL-1Ra, IFNα, IL-6, TNF-α and C-C motif chemokine 20) by peripheral blood mononuclear cells (PBMCs) (28–31), indicating that serum-secreted factors (including exosomes) may be of importance in maintaining this function. This hypothesis was indirectly demonstrated by the following studies. Harshyne et al (32) demonstrated that exosome-enriched fractions from sera of patients with glioblastoma are capable of inducing scavenger receptor cysteine-rich type 1 protein M130 expression in normal monocytes. Zhou et al (33) demonstrated that serum exosomes primed macrophage polarization towards the M2 phenotype. Accordingly, the present study hypothesized that serum exosomal miRNAs may be involved in KD by influencing the genes expressed by PBMCs.

The aim of the present study was to further analyze the exosomal miRNA microarray data established by Jia et al (27) and to predict their functions by overlapping their predicted target genes with differentially expressed genes (DEGs) in PBMCs (34). Compared with the study by Jia et al (27), the cut-off value [log fold change (FC) ≥6 vs. FC >200] for screening differentially expressed miRNAs (DE-miRNAs) was broadened to identify an increased number of exosomal miRNAs.

Materials and methods

Microarray data

The miRNA microarray data were collected from the Gene Expression Omnibus (GEO) database (www.ncbi.nlm.nih.gov/geo), under the accession number GSE60965, which included the exosomal miRNA profiles of three pooled serum samples collected from 5 healthy children [the normal (N) group], 5 patients with Kawasaki disease (the KD group) and 5 patients with KD following IVIG therapy (the IVIG group) (27). The microarray platform was GPL16730 Agilent-039659 hs_miR_18_addvirus 038169.

The mRNA microarray dataset was also available from the GEO database, under the accession number GSE73577, in which PBMCs samples of 19 patients with KD were investigated prior to and following IVIG treatment (34). This dataset was dependent on a two-channel microarray platform (GPL4133, Agilent-014850 Whole Human Genome Microarray 4×44K G4112F), and two repeats were performed for each patient.

Data normalization and the identification of DE-miRNAs and DEGs

The raw data from the two datasets and the annotated symbols were downloaded from the corresponding platforms. All expression values were logarithmically transformed (base 2) and quantile normalized using the Bioconductor preprocessCore package (version 1.28.0; www.bioconductor.org/packages/release/bioc/html/preprocessCore.html) (35).

Since one pooled sample per group was used for the analysis of exosomal miRNA profiles, the DE-miRNAs between the N and KD groups and between the KD and IVIG groups were only screened by calculating the logFC value. |logFC|≥6 was set as the threshold value. Further, a Venn diagram was built using the Bioinformatics and Evolutionary Genomics online tool (bioinformatics.psb.ugent.be/webtools/Venn) to determine the overlap between DE-miRNAs that were upregulated or downregulated in the KD group compared with the N group and DE-miRNAs that were downregulated or upregulated in the IVIG group compared with the KD group, so that the DE-miRNA alterations that were observed in the KD group and were reversed by IVIG could be identified.

The DEGs of patients with KD prior to and following IVIG treatment were identified using the Linear Models for Microarray data (LIMMA) package (v.2.16.4; bioconductor.org/packages/release/bioc/html/limma.html) (36) in the Bioconductor R software environment (v.3.4.1; http://www.R-project.org/). The P-values of the DEGs were calculated using Student's t-test. P<0.05 and |logFC|>0.5 were considered the cut-off criteria.

A heat map of DE-miRNAs and DEGs was plotted using the pheatmap R package (v.1.0.8; cran.r-project.org/web/packages/pheatmap/index.html), based on a bidirectional hierarchical clustering analysis with Euclidean distance (37).

Target gene predictions for DE-miRNAs

mRNA targets of DE-miRNAs were predicted using the miRWalk database (v.2.0; zmf.umm.uni-heidelberg.de/apps/zmf/mirwalk2), which includes 12 prediction algorithms (DIANA-microTv4.0, DIANA-microT-CDS, miRanda-rel2010, mirBridge, miRDB4.0, miRmap, miRNAMap, PicTar2, PITA, RNA22v2, RNAhybrid2.1 and Targetscan6.2). Only the miRNA-target gene interactions predicted by ≥ four algorithms were collected to construct a regulatory network and visualized using Cytoscape software (v.2.8; www.cytoscape.org) (38).

In addition, the target genes of DE-miRNAs and DEGs were analyzed using a Venn diagram (bioinformatics.psb.ugent.be/webtools/Venn) to identify the overlapping genes. The shared DE-miRNA-target interactions were also visualized using Cytoscape software (v.2.8; www.cytoscape.org) (38).

Functional enrichment analysis

To analyze the potential functions of identified target genes, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed using the Database for Annotation, Visualization and Integrated Discovery (DAVID) 6.8 online tool (david.abcc.ncifcrf.gov) (39). P<0.05 was considered to indicate a statistically significant difference.

Results

Identification of exosomal DE-miRNAs

According to the threshold of |logFC|≥6, 86 exosomal DE-miRNAs were identified between the N and KD groups, including 49 that were upregulated and 37 that were downregulated, while 77 DE-miRNAs were identified between the KD and IVIG groups, including 30 that were upregulated and 47 that were downregulated.

Following overlapping the DE-miRNAs that were upregulated in the KD group compared with the N group with the DE-miRNAs that were downregulated in the IVIG group compared with the KD group, 45 common DE-miRNAs were obtained. A total of 20 common DE-miRNAs were identified when the DE-miRNAs that were downregulated in the KD group compared with the N group were overlapped with the DE-miRNAs that were upregulated in the IVIG group compared with the KD group (Fig. 1A; Table I). These results suggested that these 65 shared DE-miRNAs may be biomarkers for the development of KD and could be reversed by IVIG. Therefore, they may also be used as biomarkers for evaluating the effectiveness of IVIG. These DE-miRNAs could clearly distinguish between the three groups, and the N and IVIG groups were clustered identically (Fig. 1B).

Figure 1.

Differentially expressed serum exosomal miRNAs identified between healthy children and patients with KD, or between patients with KD and patients with KD following IVIG therapy. (A) Venn diagrams to analyze the common miRNAs between the two comparisons. (B) Heatmap indicating the differentially expressed miRNAs to distinguish between the three groups. KD, patients with Kawasaki disease; IVIG, intravenous immunoglobulin; miR, microRNA; N, negative control healthy children.

Table I.

DE-miRNAs in serum exosomes from healthy children, patients with KD and patients with KD following IVIG therapy.

miRNA	N vs. KD (logFC)	KD vs. IVIG (logFC)
hsa-miR-4695-3p	7.14	−6.98
hsa-miR-328	8.23	−8.07
hsa-miR-631	6.78	−6.62
hsa-miR-3190-5p	6.45	−6.29
hsa-miR-574-3p	6.30	−6.15
hsa-miR-1260a	7.85	−7.69
hsa-miR-3622a-3p	6.92	−6.76
hsa-miR-197-3p	8.40	−8.24
hsa-miR-664-3p	6.92	−6.76
hsa-miR-3622b-3p	6.75	−6.60
hsa-miR-3613-3p	6.46	−6.30
hsa-miR-33b-3p	7.05	−6.89
hsa-miR-149-5p	6.95	−6.80
hsa-let-7d-3p	7.41	−7.25
hsa-miR-4780	6.50	−6.34
hsa-miR-204-5p	6.95	−6.79
hsa-miR-3940-3p	7.28	−7.12
hsa-miR-1227	7.09	−6.94
hsa-miR-937	6.20	−6.04
hsa-miR-4763-5p	7.30	−7.15
hsa-miR-4713-5p	6.22	−6.06
hsa-let-7b-3p	7.47	−7.31
hsa-miR-211-5p	7.15	−6.99
hsa-miR-4730	7.12	−6.96
hsa-miR-885-5p	8.16	−8.00
hsa-miR-3184-3p	6.17	−6.01
hsa-let-7e-3p	6.25	−6.10
hsa-miR-644b-3p	7.33	−7.18
hsa-miR-1909-5p	6.39	−6.23
hsa-miR-483-3p	6.85	−6.69
hsa-miR-4436b-5p	7.77	−7.61
hsa-miR-5681b	6.72	−6.56
hsa-miR-4646-3p	6.27	−6.11
hsa-miR-1910	6.40	−6.24
hsa-miR-299-5p	6.69	−6.54
hsa-miR-4312	6.83	−6.68
hsa-miR-449b-3p	7.45	−7.29
hsa-miR-4701-5p	7.97	−7.82
hsa-miR-1296	6.91	−6.75
hsa-miR-4728-3p	7.45	−7.29
hsa-miR-4284	6.55	−6.39
hsa-miR-1976	6.65	−6.49
hsa-miR-4697-3p	7.00	−6.85
hsa-miR-3714	6.38	−6.22
hsa-miR-5010-3p	7.03	−6.88
hsa-miR-3141	−7.30	7.33
hsa-miR-134	−6.45	6.40
hsa-miR-4800-5p	−6.19	6.10
hsa-miR-19b-3p	−6.50	6.79
hsa-miR-4665-3p	−6.72	7.79
hsa-miR-483-5p	−6.70	6.71
hsa-miR-4463	−7.07	6.13
hsa-miR-1234	−6.38	7.60
hsa-miR-2392	−6.91	6.38
hsa-miR-22-3p	−6.71	6.44
hsa-miR-4689	−6.13	6.19
hsa-miR-575	−7.37	6.50
hsa-miR-3911	−6.30	6.75
hsa-miR-191-3p	−6.06	7.44
hsa-miR-1246	−7.85	7.04
hsa-miR-4271	−6.83	6.26
hsa-miR-671-5p	−8.16	6.14
hsa-miR-3610	−6.25	6.16
hsa-miR-3156-5p	−6.40	7.17
hsa-miR-4499	−6.39	7.07

These miRNAs were the shared differentially expressed miRNAs of the comparisons between N and KD, and between KD and IVIG. DE-miRNAs, differentially expressed miRNAs; N, normal; KD, Kawasaki disease; IVIG, intravenous immunoglobulin; FC, fold change; miRNA or miR, microRNA.

Target genes for DE-miRNAs

To understand how these 65 exosomal miRNAs may influence the pathogenesis of KD, their target genes were predicted using the miRWalk 2.0 database. As a result, only 1,192 potential targets for the 15 DE-miRNAs were identified and no targets were predicted for the remaining 50 DE-miRNAs.

To indirectly demonstrate that the potential target genes of the identified DE-miRNAs were indeed differentially expressed in KD, and, since the present study hypothesized that exosomal DE-miRNAs in serum may serve roles in KD by regulating immune cells, the mRNA expression profiles from PBMCs of patients with KD prior to and following IVIG treatment were also investigated. As a result, 355 DEGs were screened, including 35 that were upregulated and 320 that were downregulated (Fig. 2A). Following overlapping these DEGs with the predicted target genes of the DE-miRNAs, 28 common DEGs were obtained (Fig. 2B; Table II), which constituted a network containing 30 miRNA-mRNA regulatory associations. Of these, only miR-328-spectrin α, erythrocytic 1 (SPTA1), miR-575-cyclic AMP-responsive element-binding protein 5 (CREB5)/β-1,4-galactosyltransferase 5 (B4GALT5)/WD repeat and FYVE domain-containing 3/cystatin-A/C-X-C motif chemokine receptor 1 (IL8RA)/protein phosphatase 1 regulatory subunit 3B (PPP1R3B), miR-134-acyl-CoA synthetase long chain family member 1 (ACSL1)/C-type lectin domain family 1 member A and miR-671-5p-tripartite motif containing 25/leucine rich repeat kinase 2 (LRRK2)/kinesin family member 1B/leucine rich repeat neuronal 1 were involved the negative regulation of expression (Fig. 3).

Figure 2.

DEGs in peripheral blood mononuclear cells between patients with KD prior to and following IVIG therapy. (A) Heatmap indicating the differentially expressed genes to distinguish between the KD and IVIG groups. (B) Venn diagram to analyze the overlapping genes between the differentially expressed genes and the target genes of differentially expressed microRNAs. DEGs, differentially expressed genes; KD, Kawasaki disease; IVIG, intravenous immunoglobulin.

Table II.

DEGs in peripheral blood mononuclear cells from patients with KD prior to and following intravenous immunoglobulin therapy.

Gene	Log fold change	P-value
PIGL	1.73	1.13×10−18
CXCR3	1.16	5.66×10−18
MYL5	1.04	5.64×10−17
ADAT2	1.14	2.06×10−16
RBM3	1.91	7.68×10−16
HKDC1	1.03	1.99×10−15
TARBP1	1.17	3.04×10−15
LOC791120	1.15	4.04×10−15
IGH	1.19	2.94×10−14
TK1	1.04	5.89×10−14
SPTA1	1.09	3.88×10−6
RNF182	1.42	6.73×10−7
ANXA3	−2.60	4.80×10−24
UPP1	−1.75	2.85×10−23
RAB24	−1.10	3.35×10−23
CD177	−2.52	3.60×10−23
CARD17	−1.99	1.40×10−22
S100A9	−1.87	1.43×10−22
FCGR1B	−2.46	2.45×10−22
C19orf59	−2.70	4.38×10−22
S100A12	−2.25	7.53×10−22
LILRA5	−1.55	1.31×10−21
HSPA6	−1.50	4.27×10−20
CD59	−1.42	4.92×10−20
GALNT3	−1.38	9.15×10−19
TMEM49	−1.04	4.52×10−18
PPP1R3B	−1.43	7.36×10−18
TIFA	−1.91	5.18×10−17
TRIM25	−1.19	5.94×10−17
WDFY3	−1.33	2.45×10−16
KIF1B	−1.37	2.77×10−16
CREB5	−1.24	1.44×10−15
PSTPIP2	−1.48	2.53×10−15
B4GALT5	−1.24	3.15×10−15
GAS7	−1.29	4.45×10−15
ACSL1	−1.68	9.49×10−15
KREMEN1	−1.53	1.86×10−14
IL8RA	−1.07	1.41×10−13
CLEC2B	−1.49	1.71×10−13
LRRK2	−1.00	1.99×10−13
CSTA	−1.30	6.85×10−12
SCN9A	−1.04	1.36×10−11
CEBPD	−1.26	8.96×10−11
CLEC1A	−1.24	1.17×10−7
NECAB1	−1.05	6.29×10−7
LRRN1	−1.28	9.69×10−7
LYVE1	−1.31	1.02×10−5

These genes included top 10 differentially expressed genes and 28 differentially expressed genes that overlapped with target genes of differentially expressed microRNAs. DEGs, differentially expressed genes; KD, Kawasaki disease.

Figure 3.

The regulatory associations between DE-miRNAs and DEGs. miR, microRNA; DE-miRNAs, differentially expressed miRNAs; DEGs, differentially expressed genes; red, upregulated; green, downregulated; arrowhead, similar expression trend between miRNAs and target genes; vertical line, opposite expression trend between miRNAs and target genes.

The target genes of exosomal DE-miRNAs were subjected to analysis by DAVID for a functional enrichment analysis. The results indicated that 79 significant GO-biological process (BP) terms were enriched, including GO:0045893, positive regulation of transcription, DNA-templated (including CREB5); GO:0043406, positive regulation of MAP kinase activity (including LRRK2); and GO:0005975, carbohydrate metabolic process (including B4GALT5; Table III). A total of seven significant KEGG pathways were also enriched with target genes of exosomal DE-miRNAs, including hsa04390: Hippo signaling pathway (including WNT4, WNT8B); hsa05161: Hepatitis B (including CREB5); hsa04931: Insulin resistance (including CREB5 and PPP1R3B); and hsa04151: PI3K-Akt signaling pathway (including CREB5; Table IV).

Table III.

GO enrichment analysis for target genes of DE-miRNAs in exosomes of serum from healthy children, patients with KD and patients with KD patients following IVIG therapy.

A, Target genes of DE-miRNAs
Term	P-value	Genes
GO:0045893, positive regulation of transcription, DNA-templated	1.39×10−5	PPARD, GDF2, MITF, ZXDC, TGFB3, NFKB1, CTCF, GLI3, LGR4, ZIC3, WNT4, RRN3, ZNF281, TBL1XR1, EGR2, FOXJ2, SOX11, MED14, RB1, ESR2, MED13, SIX4, HIPK2, MAPK3, ZNF711, TFAP2B, ERBB4, HOXA11, SOX2, EHF, CDH1, NFYA, TFAM, NR1D2, NPAT, CREBL2, KLF6, IL5, TRIP4, TBX3, SMAD4, IGF1, MSTN, CREB5, ATMIN, FZD4, RLF, YWHAH, SFRP1, SP1, SETD7, NEUROD1, TP53INP2, F2R
GO:0045944, positive regulation of transcription from RNA polymerase II promoter	1.05×10−3	CCNT2, HLF, GDF2, RNASEL, E2F8, STAT5B, MITF, EDN1, ARID4B,TGFB3, CTCF, NFKB1, GLI3, ZIC3, ZBTB38, CRX, PGR, PCGF5, ZNF304, SERPINE1, OGT, FGF1, TBL1XR1, SATB2, EXOSC9, EGR2, FOXJ2, SOX11, MTA2, RXRA, MED14, RB1, SIX4, MED13, CD40, GRHL2, GTF2H1, ACVR2A, VEGFA, MAPK3, HIPK2, TFAP2B, TFAP2C, CRTC3, FGFR2, SOX2, ONECUT2, TAF9B, EHF, EGLN1, CDC73, NR2C2, ATF1, TFAM, CHD7, NIPBL, DDX3X, PKD2, MYF6, IKZF4, IKZF1, CEBPD, SMAD4, IGF1, EN2, CSRP3, TET1, PARK7, DDX58, RLF, ATRX, SP1, TRPS1, NEUROD1, IRF2, NHLH2, PBX3, PBX2, FOXI1
GO:0043200, response to amino acid	1.36×10−3	ICAM1, SLC1A2, GLRB, MTHFR, CDKN1B, GLRA3, EDN1, CDO1
GO:0043406, positive regulation of MAP kinase activity	1.49×10−3	RASGRP1, EDN1, VEGFA, PDE5A, ADRA2A, KITLG, PDGFC, KIT, CD40, FGF1, LRRK2
GO:0043268, positive regulation of potassium ion transport	1.57×10−3	DRD1, KIF5B, ADRA2A, STK39, DLG1
GO:0051091, positive regulation of sequence-specific DNA binding transcription factor activity	2.15×10−3	IL5, EDN1, TRIM27, TRIM14, TRIM25, ESR2, KIT, FZD4, TRIM21, PARK7, DDX58, TRIM32, HIPK2, NEUROD1, NHLH2
GO:0019228, neuronal action potential	3.97×10−3	DRD1, SCN1A, SCN3A, GRIK2, ANK3, KCNA1, SCN9A
GO:0014066, regulation of phosphatidylinositol 3-kinase signaling	4.01×10−3	FGFR2, C3ORF58, EREG, ERBB4, ERBB3, RASGRP1, MAPK3, KITLG, KIT, FGF1, PIP4K2A, PIP4K2C
GO:0050890, cognition	5.41×10−3	MAGT1, CHD7, NIPBL, PTCHD1, NF1, CHRNA4, CHRNB2, CHMP2B
GO:0007265, Ras protein signal transduction	5.42×10−3	ZNF304, PLD1, PLCE1, RASGRP1, NF1, ADRA2A, IGF1, RB1, SHC3, KSR1, PARK7
GO:0014070, response to organic cyclic compound	5.50×10−3	ICAM1, CD83, ACSL1, SFRP1, TRPA1, ABCC4, ABCD3, COMT, ATF1
GO:0047496, vesicle transport along microtubule	6.24×10−3	DYNC1I1, NDEL1, KIF5B, HTT, RASGRP1
GO:0010842, retina layer formation	6.48×10−3	PROM1, HIPK2, FJX1, TFAP2B, CALB1, DSCAM
GO:0061024, membrane organization	7.89×10−3	YWHAH, RAB14, YWHAB, TBC1D4, PMP2, RAB10, YWHAE
GO:0008585, female gonad development	8.14×10−3	WNT4, COL9A3, SFRP1, ZFP42, TIPARP
GO:0060021, palate development	9.66×10−3	WFIKKN2, ACVR2B, SATB2, CHD7, TBX3, TIPARP, TGFB3, SMAD4, COL2A1, C5ORF42, GLI3
GO:0090073, positive regulation of protein homodimerization activity	1.14×10−2	CRBN, TIRAP, PARK7, TRAF4
GO:0048565, digestive tract development	1.22×10−2	FGFR2, TRPS1, TGFB3, PDGFC, RB1, KIT, LGR4
GO:0007059, chromosome segregation	1.32×10−2	CIAO1, NDEL1, DDX3X, PPP1R7, SLC25A5, USP9X, NEK9, CTCF, SRPK1, MIS12
GO:0042552, myelination	1.34×10−2	EGR2, TSPAN2, MAL2, ATRN, CMTM8, XK, QKI, ACSBG1
GO:0001764, neuron migration	1.44×10−2	SATB2, TUBB2B, USP9X, CELSR1, PCM1, YWHAE, SEMA6A, NDEL1, NAV1, CCR4, NEUROD4, DCX, MYH10
GO:0043154, negative regulation of cysteine-type endopeptidase activity involved in apoptotic process	1.45×10−2	ARL6IP1, LAMP3, DDX3X, TNFAIP8, VEGFA, TFAP2B, USP47, RAG1, BIRC5, YWHAE
GO:0035136, forelimb morphogenesis	1.56×10−2	NIPBL, TBX3, RNF165, TFAP2B
GO:0007156, homophilic cell adhesion via plasma membrane adhesion molecules	1.59×10−2	PCDHA6, ME2, PCDHA2, CLSTN2, CADM2, PCDH9, CDH1, PTPRT, CELSR1, CDH2, IGSF9B, PCDHAC2, PCDHAC1, CDH9, PCDHA10, ROBO2, DSCAM
GO:0043372, positive regulation of CD4-positive, alpha-beta T cell differentiation	1.74×10−2	CD83, TNFSF4, SASH3
GO:0010951, negative regulation of endopeptidase activity	1.80×10−2	WFIKKN2, C5, CD109, PAPLN, FURIN, A2ML1, WFDC8, SERPINE2, SERPINE1, TFPI, PEBP1, ITIH5, CSTA, CRIM1
GO:0005975, carbohydrate metabolic process	1.83×10−2	GALNT3, GANAB, FUT9, GNPDA2, ST8SIA1, GPD1L, MAN2A2, PGM2, PGM3, GANC, ALDH1B1, SLC2A2, AKR1B1, ST8SIA5, FUT4, SPAM1, B4GALT5, PYGB
GO:0040007, growth	1.83×10−2	OPA3, BMP3, GDF2, VEGFA, BMP8B, FOXP2
GO:0006513, protein monoubiquitination	1.83×10−2	TSG101, DTL, KLHL12, RAD18, TRIM25, TRIM21
GO:0001894, tissue homeostasis	1.93×10−2	AKR1B1, TRIM32, COL2A1, RB1, TP53INP2
GO:0048745, smooth muscle tissue development	2.06×10−2	NF1, TIPARP, DLG1, FOXP2
GO:0046622, positive regulation of organ growth	2.06×10−2	ARX, IL7, RAG2, SASH3
GO:0061045, negative regulation of wound healing	2.06×10−2	WNT4, HMGCR, SERPINE1, CD109
GO:0007519, skeletal muscle tissue development	2.28×10−2	MYF6, CCNT2, CFL2, NF1, SIX4, FLNB, CSRP3, FOXP2
GO:0048839, inner ear development	2.56×10−2	CDKN1B, CXCL14, CEBPD, DUOX2, SOX2, TGFB3, NEUROD1
GO:0031954, positive regulation of protein autophosphorylation	2.73×10−2	RAP2B, VEGFA, PDGFC, RAD50, CALM2
GO:0097150, neuronal stem cell population maintenance	2.73×10−2	SOX2, CDH2, PCM1, HOOK3, MMP24
GO:0033157, regulation of intracellular protein transport	2.80×10−2	NDEL1, SH3TC2, LCP1
GO:0021631, optic nerve morphogenesis	2.80×10−2	CHRNB2, GLI3, EPHB1
GO:0048511, rhythmic process	3.03×10−2	HLF, SP1, NR1D2, SFPQ, PASD1, PRKAA2, NFYA, FBXL3
GO:0070911, global genome nucleotide-excision repair	3.13×10−2	SUMO3, UBE2N, DDB2, ERCC4, USP45, GTF2H1
GO:0098609, cell-cell adhesion	3.13×10−2	ZC3HAV1, KIF5B, CKAP5, RPL15, YWHAB, TRIM25, ARFIP1, FLNB, YWHAE, PARK7, MMP24, CHMP2B, EIF4G2, GAPVD1, DDX3X, FNBP1L, SERBP1, TMOD3, PCMT1, DNAJB1, MAPRE1, RAB10, UBAP2, AHNAK
GO:0039702, viral budding via host ESCRT complex	3.20×10−2	CHMP1A, TSG101, CHMP6, VPS37C, CHMP2B
GO:0060078, regulation of postsynaptic membrane potential	3.20×10−2	SCN1A, SCN3A, PKD2, SCN9A, SCN4B
GO:0006366, transcription from RNA polymerase II promoter	3.28×10−2	CCNT2, NCBP2, HLF, POLR2E, STAT5B, TAF9B, MITF, SOX2, ONECUT2, NFKB1, EHF, CTCF, NFYA, GLI3, ATF1, ZIC3, CRX, TFAM, MAX, DDX21, VEZF1, MYF6, ZNF831, EGR2, FOXJ2, CEBPD, SOX11, SNAPC3, SMAD4, CREB5, SIX4, GRHL2, GTF2H1, TRPS1, TFAP2B, IRF2, NEUROD1, TFAP2C, PBX3, FOXI1
GO:0006044, N-acetylglucosamine metabolic process	3.29×10−2	CHST7, GNPDA2, GNPNAT1, MGEA5
GO:0060134, prepulse inhibition	3.29×10−2	DRD1, SLC6A3, NRXN1, CTNNA2
GO:0032897, negative regulation of viral transcription	3.29×10−2	TRIM32, TRIM14, TRIM27, TRIM21
GO:0007399, nervous system development	3.30×10−2	PCDHA6, GLRB, FUT9, MOBP, PCDHA2, ERBB4, CAMK2G, ARID1B, IGSF9B, GAS7, NR2C2, PCDHAC2, PCDHAC1, SEMA6A, ATXN3, NDEL1, TPP1, VEGFA, MSI1, PCDHA10, DCX, CRIM1, DLG1, WNT8B, DSCAM
GO:0045892, negative regulation of transcription, DNA-templated	3.46×10−2	PPARD, GCLC, TSG101, CTCF, GLI3, LGR4, ZBTB38, WNT4, ZNF227, NIPBL, NR1D2, GATAD2A, ZNF425, PASD1, CRY1, BAHD1, MYF6, ZNF281, IKZF4, TNFSF4, TBX3, IKZF1, CEBPD, YWHAB, SMAD4, BIRC5, RB1, SIX4, FOXP2, CHMP1A, CDKN1B, SFRP1, TRIM33, EREG, SFPQ, RBAK, USP47, TFAP2B, XCL1
GO:0015758, glucose transport	3.52×10−2	PPARD, SLC2A10, SLC2A2, EDN1, SLC2A1, HK2
GO:0051402, neuron apoptotic process	3.52×10−2	MAX, USP53, GRIK2, ERBB3, RB1, NLRP1
GO:0006914, autophagy	3.58×10−2	TSG101, CHMP6, VPS41, VPS37C, PARK7, VTI1A, CHMP2B, TBC1D25, ATG5, FNBP1L, RB1CC1, ATG4A, LRRK2, VPS39
GO:0050680, negative regulation of epithelial cell proliferation	3.60×10−2	FGFR2, PPARD, EREG, SFRP1, SOX2, CDC73, RB1, DLG1
GO:0006479, protein methylation	3.70×10−2	PCMTD2, BHMT, PCMT1, ETF1, N6AMT1
GO:0045662, negative regulation of myoblast differentiation	3.70×10−2	PPARD, TBX3, CXCL14, MSTN, CSRP3
GO:0046854, phosphatidylinositol phosphorylation	3.75×10−2	FGFR2, EREG, ERBB4, ERBB3, PI4K2A, KITLG, PI4K2B, KIT, FGF1, PIP4K2A, PIP4K2C
GO:0045787, positive regulation of cell cycle	3.94×10−2	FGFR2, ANKRD17, CDKN1B, TBX3, TRIM32, TRIM21
GO:0007585, respiratory gaseous exchange	3.94×10−2	HNMT, TMPRSS11D, EDN1, CHRNA4, PBX3, TRAF4
GO:0006813, potassium ion transport	3.96×10−2	KCNS3, KCNMA1, KCNS1, CDKN1B, SLC12A2, ATP4B, SLC24A3, KCNA1, KCNA6, KCNJ12
GO:0051260, protein homooligomerization	4.03×10−2	CCDC88C, GLRA3, KCNA1, PRND, KCNA6, KCNA7, KCNS3, ANXA6, STOM, KCNS1, CLDN1, KCTD16, ZBTB1, SLC1A1, EHD3, SPAST, KCTD12
GO:0034454, microtubule anchoring at centrosome	4.05×10−2	KIF3A, PCM1, HOOK3
GO:0035020, regulation of Rac protein signal transduction	4.05×10−2	SSX2IP, OGT, CRK
GO:0010606, positive regulation of cytoplasmic mRNA processing body assembly	4.05×10−2	CNOT6L, CNOT2, CNOT6
GO:0010960, magnesium ion homeostasis	4.05×10−2	ANK3, KCNA1, TFAP2B
GO:0071910, determination of liver left/right asymmetry	4.05×10−2	PKD2, CCDC39, ZIC3
GO:0045165, cell fate commitment	4.22×10−2	FGFR2, WNT4, ERBB4, TRPS1, ONECUT2, NEUROD4, WNT8B
GO:0071456, cellular response to hypoxia	4.24 ×10−2	ICAM1, PPARD, PTGIS, TBL2, STC2, CPEB2, SFRP1, EDN1, VEGFA, BNIP3L, HIPK2
GO:0032456, endocytic recycling	4.25×10−2	STX6, VPS52, RAB14, ARL4C, EHD3
GO:0030307, positive regulation of cell growth	4.52×10−2	EIF4G2, EXTL3, EXOSC9, DDX3X, SFRP1, TRIM32, TAF9B, USP47, H3F3B, N6AMT1
GO:0042384, cilium assembly	4.54×10−2	KIF3A, DZIP1, ONECUT2, PCM1, C5ORF42, ACTR2, TTC30A, FNBP1L, C10ORF90, ABCC4, SSX2IP, EXOC5, EHD3
GO:0009636, response to toxic substance	4.81×10−2	GLYAT, MAPK3, SLC30A4, SLC6A14, SCN9A, CDH1, GUCY2C, HTR1D, NQO1, PON3
GO:2000679, positive regulation of transcription regulatory region DNA binding	4.82×10−2	NEUROD1, IGF1, RB1, PARK7
GO:0022408, negative regulation of cell-cell adhesion	4.82×10−2	NF2, TNR, EPB41L5, CDH1
GO:0007076, mitotic chromosome condensation	4.82×10−2	CHMP1A, NCAPH, NCAPG, CDCA5
GO:0006351, transcription, DNA-templated	4.96×10−2	IL16, ZNF451, ZXDC, CNOT2, ZNF250, MED22, CNOT6, ZNF254, PGR, ZNF304, EPC2, MIER3, ZNF445, CRY1, SAMD4B, ZNF449, IKBKAP, SATB2, RXRA, ARID1B, TRIM33, MAPK3, TGIF2, VGLL3, CRTC3, ERBB4, HOXA11, NR2C2, ARX, ZNF227, DDX3X, CNOT6L, ZNF697, ZNF425, ZNF124, CREBL2, IKZF4, KLF6, TRIP4, IKZF1, RFX5, SMAD4, ZNF521, ZNF320, ZNF585A, ZNF627, CSRP3, TET1, FOXP2, ZNF419, ZNF417, PNRC1, JAZF1, ZNF318, PHF6, CCNT2, PPARD, ZNF518B, ZNF81, ARID4B, ZFP42, E2F8, ZNF10, ZBTB38, PCGF5, HIF1AN, BRD9, ZNF281, NFKBIZ, TBL1XR1, ZNF33A, EGR2, ZNF354A, ZNF354C, ZFY, ZFX, SF1, RB1, ESR2, ZBTB26, PURB, GTF2H1, CHMP1A, BRWD1, HIPK2, ZNF711, ZNF480, LIN54, ZNF740, POLR2E, LIN9, SCML2, ZNF660, CHD7, ZSCAN22, NR1D2, RB1CC1, NPAT, GATAD2A, PRKAA2, ZNF470, BAHD1, ZNF267, TBX3, PPHLN1, CEBPD, NLK, ZNF770, ZFP1, ZNF667, TRIM27, BIRC5, ATMIN, ZNF665, RLF, ATRX, ATXN3, SFPQ, RBAK, ZBTB5, NHLH2, ZNF461, SETD7, NEUROD4, PBX2, ZNF766, TP53INP2
GO:0048661, positive regulation of smooth muscle cell proliferation	4.96×10−2	FGFR2, EREG, HMGCR, EDN1, AKR1B1, IGF1, ABCC4, CALCRL
GO:0045669, positive regulation of osteoblast differentiation	4.96×10−2	ACVR2A, WNT4, ACVR2B, GDF2, CEBPD, SOX11, IGF1, GLI3
GO:0000122, negative regulation of transcription from RNA polymerase II promoter	4.99×10−2	PPARD, IMPACT, E2F8, EDN1, MITF, CNOT2, NFKB1, CTCF, HSBP1, ZNF254, GLI3, CRY1, DLG1, ZNF281, TBL1XR1, SATB2, SOX11, MTA2, RXRA, HNRNPA2B1, RB1, ESR2, PURB, ACVR2B, TRIM33, VEGFA, HIPK2, TFAP2B, TGIF2, TFAP2C, FGFR2, USP9X, TAF9B, SOX2, CDC73, ARX, NIPBL, GATAD2A, IKZF1, TBX3, RFX5, PTPN2, SMAD4, TRIM27, FOXP2, DLX1, DKK1, SFPQ, TRPS1, JAZF1, IRF2, ZBTB1
B, DEGs regulated by DE-miRNAs
Term	P-value	Genes
GO:0060828, regulation of canonical Wnt signaling pathway	2.82×10−2	KREMEN1, LRRK2

DE-miRNAs, differentially expressed miRNAs; miRNA, microRNA; GO, Gene Ontology; KD, Kawasaki disease; IVIG, intravenous immunoglobulin; DEGs, differentially expressed genes.

Table IV.

KEGG pathway enrichment analysis for target genes of differentially expressed miRNAs in exosomes of serum from healthy children, patients with KD and patients with KD following IVIG therapy.

A, Target genes of differentially expressed miRNAs
Term	P-value	Genes
hsa04390: Hippo signaling pathway	6.54×10−4	NF2, SOX2, YWHAB, TGFB3, SMAD4, TEAD1, CDH1, BIRC5, YWHAE, FZD4, CTNNA2, WNT4, YWHAH, RASSF6, CCND3, SERPINE1, AXIN2, FGF1, BMP8B, WNT8B, DLG1
hsa05200: Pathways in cancer	2.19×10−3	FGFR2, ADCY1, PPARD, STAT5B, MITF, TGFB3, KITLG, NFKB1, EGLN1, CDH1, KIT, GLI3, CUL2, MAX, WNT4, RASGRP1, SLC2A1, GNG2, AXIN2, FGF1, TRAF4, WNT8B, COL4A4, PTGER3, RXRA, SMAD4, IGF1, BIRC5, RB1, FZD4, CTNNA2, CCDC6, LAMA4, CDKN1B, MAPK3, VEGFA, LAMC2, CRK, F2R
hsa05161: Hepatitis B	2.70×10−2	EGR2, STAT5B, YWHAB, TGFB3, TIRAP, SMAD4, BIRC5, CREB5, NFKB1, RB1, DDX58, CDKN1B, DDX3X, MAPK3, TICAM2, DDB2
hsa04931: Insulin resistance	2.79×10−2	PRKAB2, CREB5, NFKB1, PRKCQ, PPP1R3B, SLC2A2, SLC2A1, MGEA5, TBC1D4, GYS2, OGT, PRKAA2, PYGB
hsa04066: HIF-1 signaling pathway	3.23×10−2	CUL2, CDKN1B, CAMK2G, MAPK3, EDN1, VEGFA, SLC2A1, SERPINE1, HK2, IGF1, NFKB1, EGLN1
hsa04151: PI3K-Akt signaling pathway	4.10×10−2	FGFR2, PPP2R5A, KITLG, NFKB1, COL2A1, KIT, COL6A6, TNR, GYS2, PDGFC, GNG2, PRKAA2, FGF1, THBS2, COL4A4, SGK2, IL7, RXRA, YWHAB, IGF1, CREB5, YWHAE, LAMA4, CDKN1B, YWHAH, CCND3, VEGFA, MAPK3, LAMC2, F2R
hsa00500: Starch and sucrose metabolism	4.51×10−2	PGM2, GANC, HK2, GYS2, PYGB, PGM2L1
B, Differentially expressed genes regulated by differentially expressed miRNAs
Term	P-value	Genes
hsa00512: Mucin type O-Glycan biosynthesis	2.82×10−2	KREMEN1, LRRK2

KEGG, Kyoto Encyclopedia of Genes and Genomes; GO, Gene Ontology; KD, Kawasaki disease; IVIG, intravenous immunoglobulin. 1, for all target genes of differentially expressed miRNAs; 2, for differentially expressed genes regulated by differentially expressed miRNAs.

The target genes of DE-miRNAs that overlapped with DEGs were also subjected to analysis by DAVID for functional prediction. Only one significant GO-BP term (GO:0060828, regulation of canonical Wnt signaling pathway, including LRRK2; Table III) and KEGG pathway (hsa00512: Mucin type O-Glycan biosynthesis, including B4GALT5; Table IV) were enriched.

Discussion

The present study demonstrated that exosomal miR-328, miR-575, miR-134 and miR-671-5p in serum may be used as biomarkers for the diagnosis of KD and for the prediction of therapeutic outcomes of the IVIG therapy. miR-328 was upregulated, and miR-575, miR-134 and miR-671-5p were downregulated in patients with KD. These trends were reversed following IVIG treatment, leading to downregulation of miR-328, and upregulation of miR-575, miR-134 and miR-671-5p. The present study identified novel exosomal miRNAs (miR-575 and miR-134), in addition to those reported in the study by Jia et al (27) (miR-328 and miR-671-5p).

There have been a few studies that have investigated the miRNAs in KD prior to and following IVIG therapy (22,40,41). However, the roles of miR-328, miR-575, miR-134 and miR-671-5p remain unclear. Previous studies have suggested that some of these miRNAs are predictors for coronary artery diseases (42,43). He et al (42) demonstrated that elevated plasma miR-328 levels could distinguish between patients with acute myocardial infarction (AMI) and the control group, with an AUC of 0.887. Wang et al (43) also reported good diagnostic performance of miR-328 in plasma (AUC=0.810) and in whole blood (AUC=0.872) for patients with AMI. These results suggested that the miRNAs selected in the present study may identify the presence of cardiovascular lesions in KD. Although previous evidence indicated that miR-134 may serve a diagnostic role in AMI, this result was not consistent with the results of the present study (miR-134 was previously upregulated in AMI, but was downregulated in the present study) (42,44). This difference may be attributed to the difference in the miRNA spectrum between the plasma and serum (45), or between the whole serum and the exosomal fraction. Therefore, further clinical studies are needed to confirm the diagnostic value of the exosomal miRNAs identified in the present study.

The present study also predicted that the selected miRNAs may be involved in KD by regulating the inflammatory target genes expressed in PBMCs. CXCR1 was one of the targets and previous studies have reported that IL-8 and its receptors, CXCR1 and C-X-C chemokine receptor type 2, were upregulated in patients with KD (46) and in coronary artery diseases (47). By binding to CXCR1, IL-8 may promote the production of other inflammatory mediators through the activation of the p38-mitogen-activated protein kinase/extracellular signal-regulated kinase-nuclear factor (NF)-κB pathways (48), which may subsequently induce apoptosis in vascular endothelial cells, a potential mechanism for KD (49). In addition, transcription factor CREB was also reported to enhance inflammation by inducing IL-17A production and promoting coronary artery diseases, including atherosclerosis (50). B4GALT family genes encode enzymes for the biosynthesis of different glycoconjugates and saccharide structures and are involved in protein glycosylation (51). Lactosylceramide synthesized by B4GALT6 in astrocytes activated central nervous system (CNS)-infiltrating monocytes, in a non-cell-autonomous manner, by regulating granulocyte-macrophage colony-stimulating factor, resulting in chronic CNS inflammation (52). PPP1R3B is a gene encoding the hepatic glycogen-targeting subunit of protein phosphatase-1 (PP1), which targets PP1 to glycogen synthase, increasing the activity of this enzyme in the skeletal muscles and liver (53). It has been previously demonstrated that increased glycogen accumulation is associated with obesity-linked inflammation in humans (54). LRRK2 deficiency was reported to attenuate the lipopolysaccharide-induced expression of inducible nitric oxide synthase, TNF-α, IL-1β and IL-6 through inhibition of the p38 mitogen-activated protein kinase, and NF-κB pathways, and to alter neuroinflammation (55). ACSL1 enzyme catalyzes the thioesterification of fatty acids and is a marker of inflammatory activation. It has been reported that the inflammatory phenotype of diabetic mice is associated with the increased expression of ACSL1 (56). Myeloid-selective deletion of ACSL1 protects monocytes and macrophages from the inflammatory effects of diabetes and prevents accelerated atherosclerosis (56). Heterozygous loss of spectrin in mice lead to increased expression levels of IL-1α and IL-1β through the activation of signal transducer and activator of transcription 3 (57). Accordingly, the present study hypothesized that miR-328 may exhibit pro-inflammatory effects through the downregulation of SPTA1, and miR-575, miR-134 and miR-671-5p may exhibit anti-inflammatory effects, leading to the upregulation of CREB5, IL8RA, PPP1R3B, ACSL1 and LRRK2 in KD.

There are certain limitations that should be acknowledged when interpreting the results of the present study. Firstly, the sample size for screening exosomal miRNAs for KD was not large, and an increased number of samples should be analyzed in future studies. Secondly, although several exosomal miRNAs have been suggested as potential biomarkers for the diagnosis of KD and the prediction of therapeutic outcomes for IVIG therapy, further confirmation in clinical samples is necessary. Thirdly, negative associations between miRNAs and their target genes were revealed in the present study, however, in vitro and in vivo experiments are necessary to validate these results. Fourthly, the exosomal mechanisms of these four miRNAs remain to be elucidated.

In conclusion, the present study preliminarily revealed that exosomal miR-328, miR-575, miR-134 and miR-671-5p may serve as biomarkers for the diagnosis of KD and the prediction of therapeutic outcomes for IVIG therapy by influencing the expression levels of inflammatory genes.