Identifying genetic variations in emm 89 Streptococcus pyogenes associated with severe invasive infections

Abstract

Streptococcus pyogenes causes mild human infections as well as life-threatening invasive diseases. Since the mutations known to enhance virulence to date account for only half of the severe invasive infections, additional mechanisms/mutations need to be identified. Here, we conducted a genome-wide association study of emm89 S. pyogenes strains to comprehensively identify pathology-related bacterial genetic factors (single-nucleotide polymorphisms [SNPs], indels, genes, or k-mers). Japanese (n = 311) and global (n = 666) cohort studies of strains isolated from invasive or non-invasive infections revealed 17 and 1075 SNPs/indels and 2 and 169 genes, respectively, that displayed associations with invasiveness. We validated one of them, a non-invasiveness-related point mutation, fhuB T218C, by structure predictions and introducing it into a severe invasive strain and confirmed that the mutant showed slower growth in human blood. Thus, we report novel mechanisms that convert emm89 S. pyogenes to an invasive phenotype and a platform for establishing novel treatments and prevention strategies.

Keywords: S. pyogenes; Streptococcus pyogenes; ferric iron transporter; infectious disease; invasive infections; microbiology; pangenome analysis; pangenome-wide association study.

Figure 1.. 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, and 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.

Figure 1—figure supplement 1.. Phylogenetic tree for the Japanese cohort.

From the inside, color bars show clusters, phenotypes, and multilocus sequence typing (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.

Figure 1—figure supplement 2.. 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 promoter activity. Turner et al. further determined subtypes using mutations other than the key residues. FS24 strain has the same allele A₋₂₇G₋₂₂T₋₁₈ as the type 3 or clade 3 but has a single-nucleotide polymorphism (SNP) highlighted in yellow in the –10 box.

Figure 1—figure supplement 3.. Pipeline for bacterial pan-genome analysis, pan-genome-wide association study (pan-GWAS), and GWAS constructed in this study.

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

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

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 1000 times with randomly permuted phenotypes, with the significance level set at the 5th percentile of the 1000 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 Supplementary file 1, table S6. SNP, single-nucleotide polymorphism; indel, insertion/deletion.

Figure 2—figure supplement 1.. Pan-genome-wide association study (pan-GWAS) on single-nucleotide polymorphisms (SNPs)/indels.

Scree plot of eigenvalues generated by decomposition of the genetic distance matrix for the Japanese (A) and global (B) cohorts. 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.

Figure 3.. Pan-genome-wide association study on gene presence.

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 Supplementary file 1, table S8.

Figure 3—figure supplement 1.. Pairwise correlation of significant clusters of orthologous genes (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.

Figure 4.. 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 bp. (D) The significant k-mers are indicated using arrowheads.

Figure 5.. Predicted protein structure models.

(A) Snake-like plot of transmembrane regions of CovS estimated using SOSUI. Frameshift mutations, detected by the single-nucleotide polymorphisms (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 pTM is a metric for overall topological accuracy.

Figure 5—figure supplement 1.. 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.

Figure 6.. 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 versus fhuB T218C mutant in human blood. Significantly up- 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 versus THY. The downward triangle plots indicate genes that are either downregulated in the WT or upregulated in the mutant strain, in blood versus THY. WT, wild-type; THY, Todd Hewitt broth supplemented with 0.2% yeast extract.

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

(A) Intracellular ferric ion assay. The wild-type (WT) and fhuB T218C mutant strains were incubated in healthy human blood for 3 hr, 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 hr to those at 0 h after infection. 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.