Transcriptome Remodeling Contributes to Epidemic Disease Caused by the Human Pathogen Streptococcus pyogenes

Abstract

For over a century, a fundamental objective in infection biology research has been to understand the molecular processes contributing to the origin and perpetuation of epidemics. Divergent hypotheses have emerged concerning the extent to which environmental events or pathogen evolution dominates in these processes. Remarkably few studies bear on this important issue. Based on population pathogenomic analysis of 1,200 Streptococcus pyogenes type emm89 infection isolates, we report that a series of horizontal gene transfer events produced a new pathogenic genotype with increased ability to cause infection, leading to an epidemic wave of disease on at least two continents. In the aggregate, these and other genetic changes substantially remodeled the transcriptomes of the evolved progeny, causing extensive differential expression of virulence genes and altered pathogen-host interaction, including enhanced immune evasion. Our findings delineate the precise molecular genetic changes that occurred and enhance our understanding of the evolutionary processes that contribute to the emergence and persistence of epidemically successful pathogen clones. The data have significant implications for understanding bacterial epidemics and for translational research efforts to blunt their detrimental effects.

Importance: The confluence of studies of molecular events underlying pathogen strain emergence, evolutionary genetic processes mediating altered virulence, and epidemics is in its infancy. Although understanding these events is necessary to develop new or improved strategies to protect health, surprisingly few studies have addressed this issue, in particular, at the comprehensive population genomic level. Herein we establish that substantial remodeling of the transcriptome of the human-specific pathogen Streptococcus pyogenes by horizontal gene flow and other evolutionary genetic changes is a central factor in precipitating and perpetuating epidemic disease. The data unambiguously show that the key outcome of these molecular events is evolution of a new, more virulent pathogenic genotype. Our findings provide new understanding of epidemic disease.

FIG 1

Temporal and geographic distribution of the emm89 strain cohort. Shown is the temporal distribution of the emm89 strains by clade. The inset shows the geographic distribution of the isolates by country. The colored horizontal bars at the bottom of the figure show the temporal distribution of the strains by country. A single isolate from Italy is not illustrated. The reduced numbers of cases in 2014 are due to United States isolates not being available for study rather than to a decline in the frequency of infections. Clade 3 strains emerged in 2003 and expanded greatly in number, displacing the predecessor clade 1 and 2 strain samples studied in all 3 populations (United States, Finland, and Iceland).

FIG 2

Genetic relationships among emm89 strains. Genetic relationships were inferred by the neighbor-joining method based on concatenated core chromosomal SNP data using SplitsTree. (A) Genetic relationships based on 28,425 SNPs identified among the members of the entire population of 1,200 emm89 strains. (B) Genetic relationships based on 11,846 SNPs identified among the 1,193 major population strains. Isolates are colored by cluster as determined using Bayesian analysis of population structure (BAPS) as indicated in the hierarchy below the figure. Three major clades (C1, C2, and C3) are defined at the first level of clustering. Subclade 3D (SC3D), a recently emerged and expanding population of strains in Finland, is defined at the second level of clustering. Indicated for the inferred phylogenies are the mean genetic distances (MGDs), both inter- and intraclade, measured as differences in core chromosomal SNPs. The mean genetic distance among strains within clades is less than the MGD to strains of the nearest neighboring clade(s). Bootstrap analysis with 100 iterations gives 100% confidence for all of the clade-to-clade branches (i.e., C1-C2, C2-C3, and C3-SC3D). (C) Genetic relationships based on 8,989 SNPs identified among the major population of 1,193 strains, filtered to exclude horizontally acquired sites as inferred using Gubbins. Exclusion of sites attributed to horizontal gene transfer events collapses the MGD strain-to-strain both within and between the clades. The MGD within the clades remains less than the MGD to the nearest neighboring clade(s). Trees in panels B and C are illustrated at the same scale.

FIG 3

Genetic relationships between strains of various Emm/M protein serotypes. Genetic relationships were inferred among 49 GAS strains of 20 M types based on 75,184 concatenated core chromosomal SNPs by the neighbor network method. The analysis is based on 42 closed genomes and 7 whole-genome-sequenced emm89 genetic outlier strains (indicated in italics). The MGD interserotype consists of 16,340 SNPs. emm89 strains are the only emm-type strains with two distinct lineages (L1 and L2) in the interserotype network. The MGD of 14,247 SNPs between the emm89 L1 and L2 genomes is greater than the MGD of 11,548 SNPs among the serotype M5, M6, M18, and M23 genomes. Of note, the emm89 L1 to L2 MGD is greater than the emm89 L1 to M53 genome MGD of 14,194 SNPs.

FIG 4

Distribution of SNPs and regions of horizontal gene transfer. Illustrated in the genome atlas of clade 2 strain MGAS23530 from the 1st (outermost) ring to the 7th (innermost) ring are the following. (Ring 1) Genome size in megabase pairs (black). (Ring 2) Landmarks: rRNA, 23S 16S-5S rRNA; FCT, fibronectin/collagen/T-antigen; SLS, streptolysin S; SRT, streptin; SAL, salvaricin; MGA, mga operon; HAS, hasABC capsule synthesis operon. (Rings 3 and 4) Coding sequences on the forward (light green) and reverse (dark green) strands. (Ring 5) Clade 1 strain MGAS11027 SNPs (n = 1,915, light blue) relative to clade 2 strain MGAS23530. (Ring 6) Clade 3 strain MGAS27061 SNPs (n = 415, red) relative to clade 2 strain MGAS23530. (Ring 7) Predicted regions of horizontal gene transfer separating clade 1 and 2 strains (light blue), clade 2 and 3 strains (red), and clade 3 and subclade 3D strains (dark blue) as listed in Table 1. SNPs are nonrandomly distributed. Regions of elevated SNP density correspond to predicted horizontal gene transfer/recombination blocks.

FIG 5

Prophage content of the emm89 strains. Shown is the phylogeny inferred by neighbor-joining for the 1,193 clade 1, 2, and 3 isolates based on 8,989 core SNPs filtered to exclude SNPs acquired by horizontal gene transfer events. Isolates are colored by phage genotype (PG) as indicated in the index. PGs were assigned in order of prevalence of occurrence in the strain sample. With the exception of PG02 (absence of phage), most of the PGs are exclusive to a single clade. PG01 was first present in the strain sample in 2003 in two isolates, one each of clades 2 and 3. The year 2003 is also when epidemic clade 3 strains were first present in the strain sample.

FIG 6

Transcriptome analysis of genetically representative preepidemic and epidemic emm89 strains. RNAseq analysis was done in triplicate for six genetically representative strains. The strain index provided in panel C applies to all of the panels. (A) Growth curves. The graph shows the averages of growth curves analyzed in triplicate. The growth curves were closely similar for all strains. Cells were harvested for RNA isolation at mid-exponential growth (ME = optical density at 600 nm [OD₆₀₀] of 0.5) and early stationary growth (ES = 2 h post-exponential phase). (B and C) Principal component analyses. Illustrated are transcriptional variances among the strains expressed as the primary and secondary principal components, the two largest unrelated variances in the data. Strain replicates cluster, illustrating good reproducibility.

FIG 7

RNAseq and qPCR expression analyses. (A) Genes significantly altered in transcript level at 1.5-fold change or greater in RNAseq. The numbers shown refer to the total number of differentially expressed genes for each comparison. The representative strains of each clade analyzed are C2/C1 (MGAS23530/MGAS11027), C3/C2 (MGAS26844/MGAS23530), and SC-3D/C3 (MGAS27520/MGAS26844). (B) Transcript levels for the nga-ifs-slo operon. The transcript levels of nga, ifs, and slo were significantly (4- to 8-fold) (P < 0.05) greater in the epidemic strains than in the preepidemic strains. The index in panel B applies to panels B, C, and D. (C) Transcript levels of speC and spd1. Transcript levels of speC and spd1 were significantly greater in the preepidemic strain at the exponential growth phase (P < 0.01). (D) Transcript levels for the hasABC operon. Transcription of hasABC was very weak for the clade 2 strain at both growth phases and was significantly less than for the clade 1 strain at the mid-exponential growth phase (P < 0.002). (E, F, and G) Relative transcript levels as measured by qPCR for nga, slo, and hasA, respectively. Differences in strain-to-strain expression were assessed by one-way ANOVA. Expression of nga and slo was significantly greater for all 5 clade 3 strains than for all 6 clade 1 and 2 strains (P < 0.001). All 3 clade 1 and 2 strains with hasABC weak/repressed promoter pattern A expressed significantly less hasA than all 3 clade 1 strains with strong/derepressed promoter pattern B (P < 0.001). Levels of expression of hasA were not significantly different among all 3 strains with weak/repressed promoter pattern A and all 5 genetically acapsular clade 3 strains. RB, recombination block; RPKM, reads per kilobase per million reads mapped.

FIG 8

Virulence assays. (A) Kaplan-Meier survival curve for mice (n = 25/strain) inoculated intramuscularly in the right hind limb with 2.5 × 10⁸ CFU. The genetically representative epidemic strain (MGAS26844) was significantly more lethal than the preepidemic strains throughout the period of observation. The index of the strains compared in panel A applies to panels A to G. (B) Histopathology scores for muscle tissue sections as determined by pathologists blind to the infecting strain. Data represent means (n = 5 assessments/strain) ± standard errors of the means (SEM). (C) Cynomolgus macaques were inoculated intramuscularly in the anterior thigh with 1.0 × 10⁹ CFU/kg of body mass. Shown at the same magnification are micrographs of muscle tissue sections from the site of inoculation. (D and E) Epidemic strain MGAS26844 caused significantly larger lesions (D) with greater tissue destruction (E) than preepidemic strain MGAS11027. (F and G) Although the bacterial burdens were similar at the site of inoculation (F), they were significantly greater for the epidemic strain than for the preepidemic strain at the distal margin (G) showing greater dissemination. P values for panels B, D, E, F, and G were determined with the Mann-Whitney test. (H) Viability of naturally occurring variant strains MGAS28980 CovR (S130N) and MGAS27552 LiaS (K214R) in human saliva persisted for 2 and 4 weeks longer, respectively, than that of wild-type strain MGAS27520. No growth, <10 CFU/ml for a 1:10 dilution. IM, intramuscularly; NHP, nonhuman primate.