For the KD patient group, fecal samples were obtained at the time of admission (the acute phase), at the time of discharge (the convalescent phase), and at 4–6 months after the onset of KD (the non-acute phase). The study protocol was approved by the institutional medical ethics committee of the University of Tokyo, Tokyo Women's Medical University and the National Institute of Infectious Diseases in Japan (Approval No. 295), and it was conducted according to the Declaration of Helsinki Principles. Written informed consent was obtained from the parents of all children for publication of their individual details and accompanying images in this manuscript. The consent form is held by the authors' institution and is available for review.

Total DNA extraction was performed using a QIAamp® DNA Stool Mini Kit (QIAGEN, Tokyo, Japan) according to the manufacturer's instructions. To increase the recovery of bacterial DNA, particularly from Gram-positive bacteria, pretreatment with lytic enzymes was performed prior to extraction using the stool kit. Briefly, 100 mg of fecal sample was suspended in 10 mL of Tris-EDTA buffer (pH 7.5), and 50 μL of 100 mg/mL lysozyme type VI purified from chicken egg white (MPBIO, Derby, UK) and 50 μL of 1 mg/mL purified achromopeptidase (Wako, Osaka, Japan) were added. The solution was incubated at 37°C for 1 h with mixing, 0.12 g of sodium dodecyl sulfate (final conc. 1%) was added, and the suspension was mixed until it became clear. Next, 100 μL of 20 mg/mL proteinase K (Wako) was added, followed by incubation at 55°C for 1 h with mixing. The cell lysate was then subjected to ethanol precipitation. The precipitant was dissolved in 1.6 mL of ASL buffer from the stool kit and subsequently purified using a QIAamp® DNA Stool Mini Kit (QIAGEN).

A DNA library was prepared using a Nextera™ DNA Sample Prep Kit (Illumina-compatible, EPICENTRE Biotechnologies, Madison, WI, USA), and DNA clusters were generated on a slide using a Cluster Generation Kit (version 2) with an Illumina cluster station (Illumina, San Diego, CA, USA) according to the manufacturer's instructions. The general procedure described in the standard protocol (Illumina) was performed to obtain standard ~1.0 × 10⁷ short reads for 1 lane. All of the sequencing runs for generating 126-mers were performed with a Genome Analyzer IIx using an Illumina Sequencing Kit (http://www.illumina.com/systems/retired_gaiix/gaiix-kits.html). Fluorescence images were analyzed using Illumina base-calling pipeline (version 1.4.0) to obtain FASTQ-formatted sequence data. The short-read sequences have been deposited in DNA Data Bank of Japan (DDBJ; accession numbers: DRA000895 and DRA001171). All of the obtained DNA sequencing reads were aligned to a reference human genomic sequence using BWA-SW read-mapping software (Li and Durbin, 2010), with quality trimming to remove low-quality reads. The remaining sequence reads were subjected to a megaBLAST search against a nucleotide database. The results of this search were analyzed and visualized using MEGAN version 4.62.3 (Huson et al., 2011), with a minimum support of 1 hit and a minimum score of 150.

The sequenced reads were assigned to a taxonomic hierarchy using MEGAN software following a megaBLAST homology search. The raw read counts were normalized by the total number of reads, and then PCA was performed using the R “prcomp” and “plot” functions. Permutational multivariate analysis of variance (PERMANOVA) with “ADONIS” was performed using 10,000 × permutations and the “bray” method with R's vegan package (Anderson, 2001).

A metagenomic biomarker discovery approach, LEfSe, was used to identify the microbial components whose sequences were more abundant in the fecal samples of the KD patients during the acute phase than in those of the KD patients during the non-acute phase and the controls. For LEfSe, Kruskal–Wallis and pairwise Wilcoxon tests are performed, followed by LDA to assess the effect size of each differentially abundant taxon (Segata et al., 2011). In this study, a p-value of <0.05 was considered significant for both statistical methods. Bacteria with markedly increased numbers were defined as those with an LDA score (log₁₀) of over 2. Less than 0.01% of the total bacterial reads, corresponding with ≤10⁷ CFU/g feces, were omitted from further analysis because of low and unreliable read counts, although significant LDA scores were observed in LEfSe.

Cultivation of Streptococcus spp. was performed using phenylethyl alcohol agar with 5% sheep blood or chocolate agar under anaerobic conditions at 37°C for 48 h. The bacterial species present were determined by performing 16S-rRNA gene sequencing using the bacterial forward primer Bac27F (5′-AGAGTTTGGATCMTGGCTCAG-3′) and the universal reverse primer Univ1492R (5′-CGGTTACCTTGTTACGACTT-3′) (Eden et al., 1991). The obtained sequences were searched against SILVA ribosomal RNA gene database to identify the bacterial species (Quast et al., 2013).

A draft genome sequence was obtained by whole-genome sequencing using MiSeq with a NEXTERA XT library preparation kit (Illumina), followed by de novo assembly with A5-MiSeq pipeline (Tritt et al., 2012). The resulting scaffolds were annotated using RAST server (Aziz et al., 2008). Maximum likelihood phylogenetic analysis of Streptococcus 16S-rDNA was performed using MEGA 6.0 with 1000 bootstrap iterations (Tamura et al., 2013).

MIC testing was performed using an Etest (bioMerieux, France) on Muller-Hinton agar (Difco, Augsburg), according to CLSI guidelines (CLSI, 2013).

This study evaluated 28 KD patients (15 males and 13 females, aged 1–114 mo; median of 25 mo). All of these patients were enrolled within 4 days of the onset of illness, with day 1 defined as the first day of fever, and they all met the diagnostic criteria for KD established by the American Heart Association (Newburger et al., 2004). All of the KD patients in the study received intravenous gamma globulin (2 g/kg) and aspirin (30–50 mg/kg/day). One male patient (patient P2) had a persistent fever despite receiving these therapies and was administered additional intravenous gamma globulin (1 g/kg) and prednisolone sodium succinate (2 mg/kg/day). This patient had transient dilatation of the coronary artery, whereas the other 27 patients showed no evidence of cardiac abnormalities.

In this study, the time of admission was defined as the acute phase, while 4–6 months after the onset of KD was considered the non-acute phase. The profiles of the participants, including the age, sex, concomitant symptoms and empirical antimicrobial treatment received, are shown in Table 1.

A total of 56 samples (28 samples each for the acute and non-acute phases) were collected, including two samples from each KD patient (Figure 1 and Table 1). Extracted DNA was subjected to metagenomic sequencing using an Illumina GAIIx next-generation DNA sequencer, and more than 10 million short 126-mer reads were obtained for each specimen. The short reads were classified at the family level of bacteria, with a threshold megaBLAST homology score of = 150. Principal component analysis (PCA) was performed to elucidate the variations between the acute and non-acute phases of KD. The results suggested that the gut microbiota was more variable during the acute phase than during the non-acute phase based on family-level taxonomy (Figure 2A). PERMANOVA with 10,000 × permutations revealed significant dissimilarity of the bacterial communities at the family level between the acute and non-acute phases of KD (F-test = 3.7307, p = 0.0006) (Figure 2A). The components were highlighted based on the antimicrobial treatment, occurrence of diarrhea, and age group, indicating that such variations during the acute phase might be associated with patient-related factors. Further, the no antimicrobial treatment group during the acute phase had was clustered in the lower right area of the PCA plot (Figure 2B), whereas the diarrhea-positive group during the acute phase exhibited a relatively scattered cluster in the PCA plot (Figure 2C); however, both groups during the non-acute phase were clustered together at the relative center of the plot (Figures 2B,C). The patients in the antimicrobial treatment group did not always exhibit diarrhea symptoms (8 diarrhea/22 antimicrobial treatment), and PERMANOVA indicated no significant differences between the two subject groups, suggesting that the gut microbiota was not affected during the acute phase of KD, regardless of the presence of antimicrobial treatment or diarrhea. The age factor showed a possible association with gut microbiota composition because the p-value was close to 0.05 when comparing subjects who were less than 2 years of age with those who were over 2 years of age (Figure 2D).

To determine the variations in gut microbiota composition between the acute and non-acute phases, linear discriminant analysis (LDA) coupled with effect size measurements (LEfSe) was applied to determine which taxa were enriched in the different groups according to metagenomic analysis (see detailed parameters in Materials and Methods). LEfSe determines which features (organisms, clades, operational taxonomic units, genes, or functions) are most likely to explain differences between classes by coupling standard tests for statistical significance (between the acute and non-acute phases in this study) with additional tests of biological consistency and effect relevance (Segata et al., 2011). The obtained metagenomic reads were classified at the genus level (Figure 3). Although LEfSe revealed that the genera Rothia and Staphylococcus were the most abundant during the acute phase, this dominance was not observed in all of the patients (Figure 3). Regardless, relatively increased numbers of Ruminococcus, Blautia, Faecalibacterium, and Roseburia bacteria were observed during the non-acute phase of KD (Figure 3), indicating that these genera was possibly related to remission in KD patients.

The above genus classifications did not reveal which common features were most likely to explain the differences between the acute and non-acute phases in the tested KD patients (n = 28). Because both pathogenic and non-pathogenic species may be included within one bacterial genus, we speculated that genus classifications would not reveal certain potential pathogens involved in KD pathogenesis; thus, further taxonomic classifications at the species level may allow for effective detection of KD-related pathogens (Figure 4). Roseburia spp. were relatively abundant during the non-acute phase and could be identified at the genus level (Figure 3), whereas Streptococcus spp. were predominantly identified during the acute phase, suggesting that some Streptococcus spp., including S. pneumoniae, S. oralis and other strains, are candidate KD-related pathogens (Figure 4). Staphylococcus hyicus might have been misidentified due to an insufficient amount of reads (less than 0.01% of the population; <1000 of the assigned reads), although the LDA score was significant.

Indeed, Streptococcus spp. were highly abundant in the gut microbiotas of some of the KD patients; for example, 77% of the bacterial reads of one KD patient (the acute phase in P7) were from Streptococcus. However, the above-mentioned S. pneumoniae and oralis species were not cultured from feces during the acute phase of KD under conventional aerobic cultivation on a phenylethyl alcohol agar plate with 5% sheep blood in screening for Streptococcus spp., while anaerobic cultivation on chocolate agar resulted in positive Streptococcus colonies. In fact, fifty colonies of Streptococcus spp. were isolated from P7-feces on chocolate agar under anaerobic conditions with incubation at 37°C for 16 h, and then 16S-rDNA sequencing was performed to determine the bacterial species present. The results suggested that seven Streptococcus spp. were unique isolates (P7-Anaero4, P7-Anaero13, P7-Anaero24, P7-Anaero25, P7-Anaero36, P7-Anaero39, and P7-Anaero45) (Figure 5C). Using the draft genome sequences of seven P7-Streptococcus isolates, including publicly available Streptococcus spp. complete genomes (23 species), the metagenomic short reads of all 56 fecal samples were classified at the species level by a megaBLAST search and LEfSe. The results suggested that the amounts of the six P7-feces-related Streptococcus isolates (P7-Anaero4, 13, 24, 36, 39, and 45, but not P7-Anaero25) and the five detected Streptococcus spp. (S. pneumonia, pseudopneumoniae, oralis, gordonii, and sanguinis) were apparently increased during the acute phase in most of the KD patients, including P7, whereas S. pasteurianus was increased during the non-acute phase (Figures 5A,B). Intriguingly, the top 4 most abundant positive isolates were P7-feces-related Streptococcus spp. (P7-Anaero4, 45, 24, and 36) rather than defined pathogenic Streptococcus species, and all positively detected Streptococcus spp. were classified within a taxonomic lineage closely related to S. oralis or pneumonia (Figure 5C). All P7-feces-related Streptococcus isolates showed susceptibility to most antimicrobial agents, including cephem, indicating that the detection of abundant P7-feces-related isolates was most likely not correlated with antimicrobial selection.

Based on species-level identification, the first and second highest hits with identical BLAST scores constituted 17.9% of all of the short reads, which were mostly homologous to rRNA genes; thus, 82.1% of the short reads could be assigned to a unique top hit. Although the obtained 126-mer short reads might not have been sufficiently long for correct species assignments, the BLAST search results suggested that the above-mentioned P7-related Streptococcus groups were highly abundant during the acute phase of KD, in contrast with other pathogenic Streptococcus spp., such as S. pyogenes, dysgalactiae, mutans, and pasteurianus. The significant detection of unique isolates in KD patients implies a possible association of KD with uncharacterized Streptococcus spp.

The 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.