Pan-genomic and phylogenetic characterization of the collected emm89 S. pyogenes isolates.

(A) emm genotyping of the 207 clinical isolates collected in Japan. (B) Pan-genome analysis of the Japanese and global cohorts. All genes detected in each cohort were classified into four groups, according to prevalence: core, soft-core, shell, and cloud genes. (C) Phylogenetic tree for the global cohort, based on the sequences of the core genes. From the inside, the color bars show clusters, phenotypes, MLSTs, clades, countries where strains were isolated, and global regions, in the order mentioned. The root of the tree was set as the mid-point. The scale located upper left indicates 0.001 times substitution of the bases on average. Arrowhead indicates strain belonging to the novel clade.

MLST, multilocus sequence typing.

Pan-genome-wide association study on SNPs/indels.

(A–B) Manhattan plots for the Japanese (A) and global (B) cohorts. The X-axis shows the location of each SNP/indel on the core gene alignment, while the Y-axis indicates the p-value. For each cohort, a permutation test was performed by iterating the calculations 1,000 times with randomly permuted phenotypes, with the significance level set at the 5th percentile of the 1,000 minimal p-values (p=5.75×10-4 and p=1.05×10-4 for the Japanese and global cohorts, respectively). Plots with lower p-values than genome-wide significant levels are colored magenta and blue, based on the direction of their effect size (positive and negative, respectively). (C) Distribution heatmap for the global cohort, with the strains possessing the significant SNPs/indels colored orange. Only the 20 SNPs with the lowest p-values are shown in this heatmap. Colored bars above indicate countries and phenotypes, and magenta bars represent invasive phenotypes. Using the Roary program, gene names starting with “Group_” were automatically assigned. Position indicates the location of each SNP/indel on the core gene alignment. The full results are shown in Table S6. SNP, single-nucleotide polymorphism; indel, insertion/deletion.

Pan-genome-wide association study on gene presence.

(A–B) Volcano plots for the Japanese (A) and global (B) cohorts. The X-axis shows the effect size, while the Y-axis indicates the p-value. Plots with a lower p-value than the genome-wide significant levels (p=1.09×10−4 and p=7.72×10−5 for the Japanese and global cohorts, respectively) have been colored magenta and blue, based on the direction of their effect size (positive and negative, respectively). (C) Distribution heatmap for the global cohort, with the strains possessing the significant genes colored orange. Only the 20 genes with the lowest p-values are shown in this heatmap. Colored bars above indicate countries and phenotypes, and magenta bars represent invasive phenotypes. Using the Roary program, gene names starting with “Group_” were automatically assigned. The full results are shown in Table S8.

K-mers related to pathology.

Detailed results of the GWAS on k-mers for the Japanese cohort in two genomic regions: covS (A, top left) and group_184 (B); and the global cohort in five regions: covS (A, top right), an intergenic region (C), group_141-143 (D), sagG (E), and fhuB (F). (A–F top) de-Bruijn graphs generated using DBGWAS. The respective nodes in the graphs indicate k-mers, with the significant nodes indicated using green arrowheads. The size of each node corresponds to the allele frequency. (A–F, bottom) Alignments of k-mers or maps of k-mers on the genome sequence of MGAS27061 around regions including significant mutations. Each base is colored according to the interpreted amino acids. (A) The arrowheads indicate significant k-mers that cause a frameshift mutation, truncating CovS protein to 35 amino acids. (B) The significant k-mers indicate that the presence of a sequence mapped on the first 26 bp of group_184 and its upstream 20 bp is related to the pathology. (C) Two significant k-mers indicate the same intergenic region of 270 bps. (D) The significant k-mers are indicated using arrowheads.

Predicted protein structure models.

(A) Snake-like plot of transmembrane regions of CovS estimated using SOSUI. Frameshift mutations, detected by the SNPs/indel- and k-mer-based GWASes, cause truncation of CovS at the indicated 35th and 45th residues, respectively. (B) Structural model of the CovS homodimer (ipTM + pTM=0.614). The putative transmembrane regions are colored orange. The upper part is the sensor domain, while the lower is the C- terminal kinase domain involved in the phosphorylation of the transcriptional regulator CovR. (C) Snake-like plot depicting the transmembrane regions of FhuB and FhuG. The 73rd residue of FhuB is indicated by an arrowhead and in magenta. (D) Structural model of the FhuBDCCG complex (ipTM + pTM=0.791). The putative transmembrane regions of FhuB and FhuG are colored green-yellow and peach, respectively. The upper part of the model is located in the extracellular region, while the lower part is in the cytoplasm. (E) The 73rd valine in FhuB, shown in magenta, was substituted with alanine. The molecular surface is illustrated with a wireframe and that of the predicted indentation is shown with an arrowhead. ipTM + pTM: Weighted combination of interface-predicted TM and predicted TM scores. ipTM is used to measure structural accuracy in the protein-protein interface, while ipTM is a metric for overall topological accuracy.

Transcriptome analysis of the fhuB T218C mutant strain in THY and human blood.

(A) Principal component analysis plot of RNA-seq data. (B) Differentially expressed genes in four comparisons. (C) Plot of gene expression in the WT strain vs. fhuB T218C mutant in human blood. Significantly upregulated and downregulated genes in the WT are colored red and blue, respectively. The shapes of the plots indicate relative transcriptional changes between THY and blood. Genes depicted with upward triangles are either significantly upregulated in the WT or downregulated in the mutant strain, in blood vs. THY. The downward triangle plots indicate genes that are either downregulated in the WT or upregulated in the mutant strain, in blood vs. THY. WT, wild-type; THY, Todd Hewitt broth supplemented with 0.2% yeast extract.

Effects of the fhuB T218C mutant strain on ferric ion uptake and bacterial survival in human blood.

(A) Intracellular ferric ion assay. The WT and fhuB T218C mutant strains were incubated in healthy human blood for 3 h, following which the intracellular ferric ion concentrations were measured. (B) Bactericidal assay in human blood. The WT and fhuB T218C mutant were mixed with healthy human blood to measure the ratio of bacterial counts at 1, 2, and 3 h to those at 0 h after infection. (C–F) Bactericidal assay in human erythrocyte-rich medium (C), polymorphonuclear cell-rich medium (D), plasma (E), and plasma inactivated by heating at 56°C for 30 min (F). (G) Bacterial growth in brain heart infusion broth. The data were pooled from three independent experiments, each performed in sextuplicate. Thick bars and error bars indicate means and quartiles, respectively. Statistical significance was determined using the Mann–Whitney test with Benjamini– Hochberg’s correction.

Comparison of clustering on the phylogenetic trees for the Japanese and global cohorts.

Complexes including k-mers significantly related to severe invasiveness detected in the Japanese cohort.

Complexes including k-mers significantly related to invasiveness detected in the global cohort.

Phylogenetic tree for the Japanese cohort.

From the inside, color bars show clusters, phenotypes, and MLSTs, in the order mentioned. The root of the tree was set as the mid-point. The scale located in the upper left indicates 0.001 times substitution of the bases on average.

Variations in the nga promoter region.

Sequences of the nga promoter are classified into four major genotypes, according to three key residues (shown as -27, -22, and -18 in the figure) affecting promotor activity. Turner et al. further determined subtypes using mutations other than the key residues. FS24 strain has the same allele A-27G- 22T-18 as the type 3 or clade 3 but has an SNP highlighted in yellow in the -10 box.

Pan-genome-wide association study on SNPs/indels.

(A–B) Scree plot of eigenvalues generated by decomposition of the genetic distance matrix for the Japanese (A) and global (B) cohort. To determine the number of principal components used as covariates for pan-GWAS, the plot in which the subsequent principal components seem to have relatively small eigenvalues was sought. Consequently, we used seven and three principal components as covariates for the Japanese and global cohorts, respectively. (C) Number of significant SNPs accumulated in single genes found in the pan-GWAS for the global cohort. Only the top 30 genes are shown in this graph.

Pairwise correlation of significant COGs.

Each cell indicates the correlation coefficient of each significant COG pair on its presence in the global cohort. The COGs are sorted based on hierarchical clustering. Annotations and phenotypes are represented by circles and cells located on the right and bottom of the heatmap. The dotted line on the tree divides the COGs into eight clusters.

Predicted protein structure model of LacE.

(A) Snake-like plot of transmembrane regions of LacE estimated using SOSUI. The 554th residue has been indicated with an arrowhead and in magenta. (B) Structure model of LacE generated with AlphaFold (pLDDT=0.902). The putative transmembrane regions are colored in orange. LacE has two domains: the transmembrane domain EIIC and intracellular domain EIIB. The 554th glycine substituted with valine is colored in magenta.

pLDDT: The best predicted local difference distance test indicating local structural accuracy.

Pipeline for bacterial pan-genome analysis, pan-GWAS, and GWAS constructed in this study.

We developed a pipeline for whole genome analyses working on the National Institute for Genetics (NIG) Supercomputer and SQUID at the Cybermedi Center, Osaka University.