All strains used in this study were collected between 1997 and 2017 and stored in the repository at Gyeongsang National University Hospital (GNUH) in Gyeongnam Province, Korea. Forty-five S. pyogenes strains were selected according to the common emm types in three years: 1997, 2006, and 2017. Forty-eight strains were randomly selected based on their invasiveness between 1997 and 2017. An “invasive strain” was defined as one isolated from a normally sterile body fluid, such as blood, cerebrospinal fluid, pleural fluid, pericardial fluid, joint fluid, bone aspirate, or a deep-tissue abscess [6]. Fig. 1 shows a flow chart of the strain selection procedure. In total, 87 strains were included in this study. For the first analysis to evaluate the genetic evolution over a 20-year time span, non-invasive strains (N=45) were selected every 10 years: 1997, 2006, and 2017. In the second analysis to evaluate assocations between invasiveness and genetic profiles, non-invasive and invasive strains were compared. Sixty-three non-invasive strains were isolated from the throats of carriers who did not have any symptoms or signs of tonsillitis. The other 24 invasive strains were isolated from blood (N=21) or joint fluid (N=3) of patients (Fig. 1).

Bacteria were identified using a Vitek-2 automated identification system (BioMérieux Inc., Marcy l’Étoile, France). All strains were inoculated in 30% glycerol in Todd-Hewitt broth and stored at –70°C. They were recovered on blood agar plates for genetic analysis.

The study protocol was approved by the Institutional Review Board of GNUH (approval number: GNUCH 2018-01-008). Informed consent was waived because of the retrospective nature of the study.

Genomic DNA was extracted using a Wizard Genomic DNA Isolation Kit (Promega, Madison, WI, USA). The DNA and potential culture contamination were checked by 16S rRNA gene sequencing using an ABI 3730 DNA sequencer (Applied Biosystems, Foster City, CA, USA). A draft genome sequence of each strain was generated by MiSeq sequencing (300-bp, paired-end) using a MiSeq Reagent Kit v3 (Illumina, San Diego, CA, USA). Sequencing libraries were prepared using the TruSeq DNA LT sample Prep Kit (Illumina). The Illumina sequencing data were assembled with SPAdes v3.13.0 (Algorithmic Biology Lab, St. Petersburg Academic University of the Russian Academy of Sciences, St. Petersburg, Russia). The EzBioCloud genome database was used for gene finding and functional annotation of the whole-genome assemblies (https://www.ezbiocloud.net). Protein-coding DNA sequences (CDSs) were predicted using Prodigal 2.6.2 [7]. The CDSs were classified based on their roles, with reference to orthologous groups (EggNOG v4.5; http://eggnogdb.embl.de). For more detailed functional annotation, the predicted CDSs were compared with those from the Swiss-Prot (https://www.uniprot.org), KEGG (http://www.genome.jp/kegg/), and SEED (http://pubseed.theseed.org) databases using the UBLAST (https://www.drive5.com/) program.

CG analyses comprised two steps (Fig. 1). The first analysis was conducted according to the isolation year (strains of 1997 vs. those of 2006 vs. those of 2017). The second analysis was conducted according to invasiveness (strains from the throat vs. those from blood/joint fluid). CG analysis was conducted by comparing functional genes based on the clustering of orthologous genes. The genome sequences of all strains were obtained from the EzBioCloud database (http://www.ezbiocloud.net/), and average nucleotide identity (ANI) values were calculated. For ANI calculation, the query genomes were cut into small fragments (1,020 bp), and high-scoring pairs between two genome sequences were selected using the USEARCH program (http://www.drive5.com/usearch). Using the calculated ANI values, a dendrogram was constructed using the unweighted pair group method. Homologous regions in a target genome to query open reading frames were determined using the USEARCH program and were aligned using pair-wise global alignment. The matched regions in the subject contig were extracted and saved as homologs [8].

We used DDBJ Fast Annotation and Submission Tool v1.2.4 (DFAST; https://dfast.nig.ac.jp) for annotation and searched the sequences around the mga annotation encoding multiple virulence gene regulators based on the annotation data [9]. If the sequences around mga were not found, the sequences around the emm1 primer (forward: 5΄-TATT(C/G)GCTTAGAAAATTAA-3΄) were searched throughout the contig sequences using the FASTA format. We extracted the emm sequences between emm1 and emm2 primers (reverse: 5΄-GCAAGTTCTTCAGCTTGTTT-3΄) using the corresponding sequences recovered from the contig data. By inserting the extracted sequences into the Centers for Disease Control and Prevention (CDC) database (https://www2.cdc.gov/vaccines/biotech/strepblast.asp), emm genotypes (including subtypes) were assigned to the extracted sequences.

When sequences were incompletely matched with an emm genotype in the CDC database, we directly PCR-amplified emm using bacterial DNA templates and the emm1/emm2 primer set and sequenced the amplicons after purification using an AccuPrep Purification Kit (Bioneer Corp., Daejeon, Korea). The emm genotypes were assigned directly to the amplified sequences based on the CDC database [10, 11].

Phylogenetic analysis was accomplished using ~1.4 million bps of orthologous protein-coding regions for 45 strains according to the isolation year and 48 strains according to invasiveness, respectively (data not shown).

To determine five target genes (speA, speB, speC, ssa, and smeZ) encoding the superantigens (also known as exotoxins), we conducted PCR simulation analysis using the Serial Cloner (http://serialbasics.free.fr/Serial_Cloner.html) application with the contig sequences, as previously reported [12-14]. The primer sets used to amplify the speA, speB, speC, ssa, and smeZ are listed in Table 1. The speB product was included as an internal control in the PCR simulation analysis because all S. pyogenes strains possess the speB sequence (955 bp). Superantigen gene profiles were determined for each strain.

We determined the STs using allelic profiles consisting of seven housekeeping genes (gki, gtr, murI, mutS, recP, xpt, and yqiL) by inserting the contig sequences obtained into the online application MLST v2.0 (https://cge.cbs.dtu.dk/services/MLST/), which is managed by the Center for Genomic Epidemiology at the Technical University of Denmark [15]. The STs were grouped into clonal complexes (CC), whereby related STs were classified as single locus variants, differing in only one housekeeping gene. An expansion of the goeBURST program implemented in PHYLOViZ was used to produce a minimum-spanning tree representing possible relationships among the STs [16].

For novel allelic numbers/STs, we submitted the data (i.e., bacterial genotypic/phenotypic data and patient backgrounds) to the S. pyogenes PubMLST (http://pubmlst.org/organisms/streptococcus-pyogenes) database. The PubMLST curator assigned novel allelic numbers/STs to our strains.

We used Fisher’s exact test (two-sided) to determine significant differences in categorical variables, and the chi-square test to compare the proportions in each emm genotype/cluster between invasive and non-invasive strains. SPSS Statistics v22.0 (IBM Corp., Armonk, NY, USA) was used for the analysis. P<0.05 was considered significant.

The emm genotypes and cluster types are presented in Supplemental Data Table 1. In total, 21 emm genotypes were identified, with emm1 (emm1.00, emm1.18, emm1.30, and emm1.76), emm4 (emm4.00), and emm12 (emm12.00, emm12.19, and emm12.49) accounting for 19.5%, 13.8%, and 20.7%, respectively. In total, eight emm cluster types were identified, among which the A-C3, A-C4, and E1 types were the most common. Fig. 2 shows the differences among the emm clusters according to invasiveness. There were significantly more invasive strains in cluster A-C3 than in the others (P<0.05). The A-C3 type included the genotypes emm1.0, emm1.18, emm1.3, and emm1.76 (Supplemental Data Table 1, Fig. 2).

The STs are presented in Supplemental Data Table 1. The 87 strains comprised 21 STs with exact loci matched against the PubMLST database. ST36, ST28, and ST39, accounting for 20.7%, 18.4%, and 11.5%, respectively, were the most frequent. There were strong associations of genetic characteristics within the MLST complex. The goeBURST diagram is shown in Fig. 3. There were 17 singletons in the CG analysis, and ST28 showed a clonal distribution of invasive strains in the second analysis. The predominant emm1 lineage belonged to ST28 (Fig. 3).

The phylogenetic tree based on the periodic comparison showed a sporadic distribution (data not shown). The second analysis revealed the genetic relationships among emm genotypes or STs. Superantigen profiling revealed that speB was present in all strains. speZ-speB and speZ-speB-speC profiles were present in 37.9% and 28.7% of the total strains, respectively. We found no significant association between the coexistence of different superantigen genes and invasiveness.

When comparing gene origins by pan-genome orthologous group (POG) analysis, we found csn1, ispE, nisK, and citC were more significantly present in invasive strains than in non-invasive strains (all P<0.05) (Table 2).

We looked for common virulence-associated CDSs among all 87 strains by searching for annotated CDSs based on functional annotation of the whole-genome assemblies. We found 25 CDSs associated with bacterial virulence. Among them, 12 (lactocepin, oleate hydratase, putative glycoslytransferases, capsule biosynthesis protein [CapA], regulatory protein MsrR, internalin-I, deoxyribonuclease, biofilm-regulatory protein, listeriolysin-regulatory protein, streptokinase, C5a peptidase, and M protein) were identified in all strains. CDSs encoding exotoxin A and procollagen-proline 3-dioxygenase were frequently detected in invasive strains (all P<0.05) (Table 3).