Atlanta invasive S. pneumoniae isolates were obtained from the Georgia Emerging Infections Program and were isolated from patients in Atlanta, Georgia as part of an ongoing, prospective, population-based surveillance study (Farley et al., 2002). Macrolide resistant isolates were assayed for the presence of macrolide resistance genes mef(E), mel and erm(B) by PCR as previously described (Stephens et al., 2005). One hundred and twenty one isolates recovered from patients between 1994 and 2011 were non-randomly selected for whole genome sequencing (Table S1). Selection criteria included capsule serotype, antibiotic resistance, macrolide resistance genes, and date of isolation. In addition to macrolide resistance, clindamycin, tetracycline, and chloramphenicol were considered indicators of mobile elements. Sequence type was also considered if available (see below). Thirteen invasive pneumococcal isolates from U.S. states other than Georgia were provided by the Antimicrobial Bacterial Core Surveillance Program and the Centers for Disease Control and Prevention (CDC) (McCormick et al., 2003). These included four isolates from Maryland, two each from Connecticut, Oregon, and Minnesota, and one each from California, New York, and Tennessee. Three invasive isolates, two from Europe and one from the North Carolina, were provided by the Pneumococcal Molecular Epidemiology Network (PMEN) (McGee et al., 2001). Finally, the study was supplemented with 10 carriage isolates collected from the Atlanta population; five each isolated pre-PCV7 (1995) and post-PCV7 (2009) (Sharma et al., 2013). Carriage isolates were collected from nasopharyngeal swabs of healthy children between the ages of 5 months and 5 years old (Sharma et al., 2013). The total number of sequenced pneumococcal genomes was 147 (Table S1).

The capsule serotype of each isolate was determined by the institution from which it was acquired and verified by analyses of the capsule locus in the genome sequence. Minimal inhibitory concentrations (MIC) of erythromycin, clindamycin, tetracycline, and chloramphenicol were determined by microdilution or Etest (bioMerieux Inc., Durham, NC). The Multi Locus Sequence Type (MLST) of selected isolates was determined by Sanger nucleotide sequencing of PCR products with MLST allele-specific primers as previously described (Enright and Spratt, 1998). The sequence types of the remaining isolates were determined from allelic data extracted from their genomes. Isolates were assigned to a clonal complex based on eBURST analysis of the complete S. pneumoniae MLST database as of September, 2014 (http://pubmlst.org/spneumoniae/).

For each isolate sequenced, two shotgun libraries were constructed for 454 sequencing: one rapid non-paired-end library for coverage and one 3 kb paired-end library for assembly scaffolding. All isolates were sequenced using the 454 GS FLX Titanium to a depth of sequence coverage ranging from 16× to 57×. Shotgun reads were assembled into 4–32 contigs per isolate using Celera Assembler v6.1 (Myers et al., 2000). Genome sequences have been deposited in the GenBank Nucleotide Sequence Database (http://www.ncbi.nlm.nih.gov/) (Table S1) and the PathoSystems Resource Integration Center (PATRIC) (Wattam et al., 2014). Isolates are available through the American Type Culture Collection-managed Biodefense and Emerging Infections Resources repository (http://www.beiresources.org).

Two high throughput reference-based pipelines were used to perform single nucleotide polymorphism (SNP) discovery and validation in 166 genomes using the TIGR4 genome as a reference (Tettelin et al., 2001). The genomes included the 147 from the present project plus 15 publicly available, closed, pneumococcal genomes and four whole genome shotgun sequences from this project that were later excluded from this study due to the lack of associated meta-data (Table S1). The analysis pipelines made extensive use of Perl (https://www.perl.org) scripts and were implemented within the Ergatis framework (https://github.com/jorvis/ergatis) (Orvis et al., 2010). The alignment-based SNP discovery pipeline named Skirret mainly executes a series of utilities from the MUMmer 3.0 package (http://mummer.sourceforge.net/) (Kurtz et al., 2004) namely nucmer, delta-filter and show-SNPs to determine SNP positions. This pipeline was run on all 166 assembled genomes with a minimum alignment identity of 98%. A k-mer based algorithm called kSNP (http://sourceforge.net/projects/ksnp/) (Gardner and Hall, 2013) was also employed to identify SNPs in the assembled genomes using 41 bp as the length of the k-mer. Next, a BLAST-based validation pipeline was used to merge the SNP panels from the Skirret and kSNP pipelines and 41 bp sequence surrounding each SNP position was extracted from the reference TIGR4 genome. These short sequences were searched against all 166 genomes using all-vs.-all BLASTN searches (e-value cut-off 0.0001) to eliminate false positive SNP calls, only core SNPs in short sequences found in all 166 genomes were selected. Putative recombinogenic regions were not excluded. The set of 10,120 core SNP bases were concatenated from each of the input genomes to create a multi-FASTA file. This file was used for phylogenetic tree construction using the Geneious Basic software version 5.6.4 with neighbor-joining and the Hasegawa, Kishino and Yano (HKY) genetic distance model. The HKY model assumes that every base has a different equilibrium base frequency, and also assumes that transitions evolve at a different rate to transversions. A consensus tree was built based on 100 iterations to provide an estimate for the level of support for each clade in the final tree. We used a 100% consensus support threshold that resulted in a strict consensus tree where the clades included were those that were present in all the trees of the original set.

Mega and Mega insertion sites were located by comparing genomes with the 5,511 bp Mega element from strain GA17457 (Data File S1). ICE were identified by BLASTN (Altschul et al., 1997) search of the terminal direct repeats localized at the ends transposons ICESp23FST81 (81.0 kb; Accession no. FM211187), Tn5253 (66.2 kb; Accession no. EU351020) (Iannelli et al., 2014), Tn916 (18.0 kb; Accession no. U49939) and others. Elements were characterized by the int genes they carried. Integrase genes were identified by BLASTN search with query sequences from described ICE including those mentioned above. All query sequences used are supplied as supporting information (Data File S1). Comparisons of composite ICE were visualized with the Artemis Comparison Tool (Carver et al., 2005). Genomes were also searched for antibiotic resistance genes, including mef(E), mef(A), mel(msr(D)), erm(B), erm(A), tet(M) and cat (Data File S1). The nucleotide sequences of each Mega, Tn916-like and composite elements are provided as supporting information (Data Files S2, S3 and S4, respectively).

Isolates of S. pneumoniae were selected from the >13,000 strain collection of the Georgia Emerging Infections Program. Isolates were collected as part of an ongoing, prospective, population-based surveillance of invasive pneumococcal disease (IPD) in the Atlanta metropolitan area (Farley et al., 2002). One hundred and twenty one invasive S. pneumoniae isolates from 1994 to 2011 were non-randomly selected for genome sequencing. Isolates were selected to provide broad coverage of the genotypes existing in the population during the surveillance period, but with a focus on macrolide resistance, with 100 resistant to erythromycin (MIC >1 μg/ml). 74 displayed the M-phenotype as indicated by susceptibility to the lincosamide clindamycin and 24 were also resistant to clindamycin (MLS_B-phenotype) (Table S1). In addition to the Atlanta invasive isolates, 10 carriage isolates from Atlanta were included. Also included were 16 strains not from Atlanta. For more detail, please see the Material Methods and supporting information (Table S1). In total 147 pneumococcal genomes were sequenced including 23 different serotypes, 77 multi-locus sequence types (ST), 24 clonal complexes (CC) and two singletons (Table S1).

The 147 isolates from this study and 15 publically available genomes were compared using whole genome single nucleotide polymorphism (SNP) analyses. The resulting phylogenetic tree is shown in Figure 2. The genomes clustered into MLST clonal complexes with the exception of CC156 which was divided into three clades (Figure 2). This may be an example of a limitation of MLST typing as the three clusters appear to be unrelated. CC156 is a large complex consisting of multiple subgroups thus not all were as closely related as the label CC156 implied (Figure S1). For example ST90 has only one allele in common with ST156. However, it was interesting to observe that two of the CC156 subgroups appear to be linked by a common Mega insertion event (see below). Further, the eBURST-predicted subgroups (SG90, SG124 and SG156) are well-established independent clonal complexes (CC90, CC124 and CC156, respectively) indicating that their inclusion into a single clonal complex is an artifact of eBURST and not a true indication of common ancestry.

Of the Atlanta invasive isolates, 86 contained Mega integrated into the pneumococcal chromosome backbone or nested within larger mobile elements. Mega was found in 17 serotypes, 16 clonal complexes and two singletons (Table 1). Of these, 48 carried the 5,511 kb Mega type 1 (Mega-1) with a 119 bp mef(E)/mel intergenic region and 38 carried the 5,412 kb Mega type 2 (Mega-2) with a 99 bp deletion in the intergenic region (Table 1). The genome of GA07643 contained a previously unidentified 5,617 bp Mega variant, designated Mega-3, characterized by a 107 bp tandem duplication that included the first 90 bp of mel and 23 bp upstream containing its ribosomal binding site (RBS). Analyses of the 454 raw reads revealed that both copies of the repeated sequence were contained on single reads and the tandem duplication region had greater than 40× sequence coverage. This confirmed that the duplication was not an assembly error. The authentic duplication resulted in an intergenic region that contained the 119 bp that immediately follows mef(E) in Mega-1, followed by a 38 codon ORF that includes the first 30 amino acids of Mel. The ORF had the mel RBS and appeared suitable for protein expression. Theoretically this could have a polar effect of mel translation, however, GA07643 had an erythromycin MIC of >16 μg/ml so it does not appear to have had a negative impact on resistance (Table S1). The full-length copy of mel overlapped the 3′ end of the duplicated ORF by 10 bp and maintained the native mel RBS. There was no correlation between the size of the mef(E)/mel intergenic region and the level of efflux-mediated erythromycin expressed by that strain (data not shown).

Mega was found inserted into each of the four previously identified insertion sites within the pneumococcal chromosome backbone (classes I-IV) (Gay and Stephens, 2001), in Tn916-like elements (class V), and in one novel site (Table 1). Six genomes contained Mega integrated into the class I insertion site, five of which were clonal (6A, CC2090) indicating clonal dissemination (Figure 2, light blue dots; Table 1). The exception was GA17457, a serotype 19A, CC199 isolate that has been used as the genetic background for studying the regulation of mef(E) and mel (Zähner et al., 2010; Chancey et al., 2011). GA17457 was closely related to four isolates in CC199 that contained Mega-1 in class III sites (Mega-1.III) and was the only CC199 strain to have Mega-1 integrated into class I (Mega-1.I) (Figure 2). Class II Mega insertions were identified in 38 genomes (Table 1). The elements inserted into this site include Mega-1 (n = 2), Mega-2 (n = 35) and the novel Mega-3 in strain GA07643 (Table 1). Mega-1.III was found in ten genomes in diverse backgrounds: five serotypes and four clonal complexes (Table 1). Mega class III insertions clustered mainly in CC199 but single isolates containing Mega-1.III were identified in CC156, CC81, and CC62 (Figure 2). Five genomes contained Mega-1 or Mega-2 inserted into the class IV site which was divided into three subclasses, IVa, IVb, and IVc (see below).

Comparison of the nucleotide sequences flanking Mega inserted into each site revealed conserved sequences upstream of the Mega insertion site that may indicate a target sequence (Figure 3). The potential consensus was 5′ TTTCCNCAA 3′ and was located approximately six bp upstream of Mega in all sites (Figure 3). The six base pairs flanking Mega on either side resembled coupling sequences as described for conjugative transposons such as Tn916 (Figure 3). The coupling sequences and Mega flanking sequences illustrated in Figure 3 were representative of all Mega elements inserted into the respective sites. This implies that the Mega integrated into each site were descended from a single, or very few, ancestral transposition event(s). An example of this can be seen in the whole genome phylogeny (Figure 2). The class II Mega insertion was present in almost all CC156 isolates (Figure 2, dark blue dots). Genomes from isolates of this large complex clustered into three clades correlating to CC156 subgroups 90, 124, and 156 (Figure 2). Almost all genomes in subgroups 90 and 156 contained Mega (Figure 2). All 10 subgroup 156 genomes contained Mega; eight Mega-2, class II, one Mega-1 class III insertion and one integration into orf6 of Tn916 (Table 1; Figure 2). Further, CC2090, which branches from CC156, has two clades in one of which nine of nine isolates contain Mega class II. The other branch of CC2090 contains almost all isolates with Mega class I indicating a second conjugation event in this complex (Figure 2). Again this strongly implies that Mega inserted into the class II site early in the history of these lineages.

The Mega class IV insertion site, as previously described (Gay and Stephens, 2001), was near the 3′ end of rumA (SP_1029, TIGR4 annotation) that encoded a tRNA methyltransferase and was located at the left end of the Pneumococcal Pathogenicity Island-1 (PPI-1) (Brown et al., 2004) (Figure 4). Large deletions and multiple MGEs were associated with Mega integrated into this site, giving rise to the three class IV subclasses (Table 1). Starting with the simplest, the class IVb insertion contained Mega-1 integrated into rumA leaving Mega flanked by the fragments of the gene (Figure 4). Two class IVb isolates were identified, all serotype 33F, CC100 (Table 1). The class IVa insertion contained Mega-2 integrated into the same location in rumA and had an identical sequence to the left junction of Mega in class IVb (Figure 4). However, the downstream remnant of rumA was deleted along with the entire 30.7 kb PPI-1 (SP_1029 to SP_1067) (Figure 4). The island was replaced by Mega-2 and a 1,801 bp insertion sequence ISSmi2 found in S. mitis strain B6 (Accession no. FN568063.1) (Figure 4). The class IVc Mega-2 insertion was similar to class IVa in that PPI-1 was deleted and ISSmi2 was located at the right junction with the chromosome (Figure 4). However, a 41,995 bp transposon-like sequence most closely related to ICESz1 from S. equi subspecies zooepidemicus strain H70 (Accession no. FM204884) (Holden et al., 2009) was inserted between the rumA insertion site and the left end of Mega-2. Thus, the class IVc Mega site organization consists of, in tandem, a transposon previously undefined in pneumococci, Mega-2 and ISSmi2. The Atlanta invasive isolate GA17545 (6B, CC1536) contained the only class IVc Mega insertion identified (Tables S1, S2). Class IVa isolates are also typically 6B, CC100, however, a single Atlanta invasive isolate, GA04375, with serotype 19F, and in CC320 was also found to contain Mega-2.IVa (Table 1; Table S1). GA04375 was isolated in 1995 and was the earliest isolated CC320 in the genome collection and only one of two CC320 that contained Mega inserted into the pneumococcal chromosome backbone. The other was 3063-00, an invasive 19F, CC320 collected in 1999 in Tennessee (Table S1, S2). The serotype 14, CC156 isolate GA02254 from 1994 contained Mega-1 inserted into an unidentified transposon-like sequence and was most closely related to a sequence in S. suis 05HAS68 (Accession no. CP002007). Twenty-seven genomes contained Mega-1 inserted into orf6 of a Tn916-like element (class V) (Table 1). Tn916-like elements are discussed further in the next section. These data suggest that PPI-1 is a hotspot for integration by recombination and conjugation. This locus could be an entry point for non-pneumococcal DNA to integrate into the pneumococcal genome. Conversely, the conservation of the target site in rumA across non-pneumococcal streptococci makes it possible that pneumococcus has been the donor and a non-pneumococcal streptococcus the recipient. The high degree of recombination in the pneumococcus makes it difficult to determine directionality of gene flow, but it likely can flow in both directions. With just these few isolates, we showed evidence for HGT with three non-pneumococcal streptococci; S. equi, S. mitis, and S. suis.

The genomes were searched for Tn916-like elements by BLASTN using the conserved ends of Tn916 (Data File S1). Of the 86 macrolide resistant invasive isolates from Atlanta, 53 contained a full-length Tn916-like element (Table S3). A macrolide resistance determinant was nested within the Tn916 element in all but six, which instead contained Tn916 and an unlinked Mega-2.II (n = 5) or Mega-1.III (n = 1) (Table S3). Four distinct Tn916-like elements containing one or two macrolide resistance determinants were identified: Tn2009, Tn2010, Tn6002, and Tn3872 (Figure 1; Table 2). Tn2009 was found in 15 Atlanta invasive isolates spanning six CC and six serotypes (Table S3). Tn6002 was carried by nine Atlanta invasive isolates including five CC (and a singleton) and seven serotypes (Table S3). Tn2010 was found in 13 Atlanta invasive isolates, two Atlanta carriage isolates, and eight non-Atlanta invasive isolates (Table 2; Table S3).

To better understand the evolutionary relationships between these elements, the insertion site of each was determined. Tn916-like elements were integrated into four unique loci in the pneumococcal chromosome backbone and into several loci within Tn5253-like elements (Figure 5). Every Tn2010 element discovered was inserted into SP_1947 (TIGR4 annotation) encoding a hypothetical secreted protein (Table 2; Figure 5). This is the same location where Tn2010 has been previously reported to be integrated (Zhou et al., 2014). Six of 15 Tn2009 elements and four of nine Tn6002 elements found in invasive Atlanta isolates were inserted into the same locus as Tn2010 (Table 2). The 6-base coupling sequences were identical for all elements in this locus suggesting that the elements found integrated into this locus were derived from a single transposition event (Figure 5). This further implies that the macrolide resistance elements Mega-1, Tn917 and erm(B) were independently acquired by a Tn916 residing in the SP_1947 locus. In addition, elements in this locus were identified in two clonal complexes (CC320, CC156) indicating transformation and recombination of the element in this locus between clonal complexes.

Tn2009 was also identified integrated into SP_1638 (TIGR4), an iron-dependent transcriptional regulator, in a single serotype 9V, CC156 isolate from 2007. The fact that no element was identified in this locus from an earlier collected isolate may indicate that this integration event occurred more recently than the SP_1947 integration. Likewise, a single isolate (22F, CC433) from 2010 contained Tn6002 integrated into the ABC transporter gene SP_1438. Five serotype 6B, CC384 isolates contained Tn916 integrated into lytA (SP_1937) encoding an autolysin (Table 2). Each of these closely related isolates contained unlinked Mega-2.II (Table S2). The other Tn916-like elements identified were integrated into a composite ICE resembling Tn5253 and are discussed below.

Tn916-like elements with macrolide resistance determinants have been associated with large conjugative transposons also known as ICE (Henderson-Begg et al., 2009; Mingoia et al., 2014). To determine if ICE have been involved in macrolide resistance dissemination in Atlanta, the genomes of 100 macrolide resistant invasive pneumococcal isolates from Atlanta were searched by BLASTN for the conserved ends of known ICE (Data File S1). The genomes of 15 macrolide resistant isolates contained composite ICE-like elements that ranged in size from approximately 57 to 88 kb (Table 3). Comparison of the integrase genes encoded on the ICE revealed two distinct families of conjugative transposons. Nine belonged to the Tn5253 family and six belonged to the ICESp23FST81 family (Table 3). Each of the Tn5253-family transposons was inserted into the pneumococcal chromosome backbone in the gene for ribobiogenesis GTPase, rbgA (SP_1155, TIGR4 annotation) involved in assembly of 23S ribosomes. Downstream of Tn5253-like elements in this locus is the immunoglobulin A1 protease gene iga as has been described previously (Croucher et al., 2009). Each of the ICESp23FST81-like elements was inserted into the chromosome backbone gene rplL (SP_1354) which encodes ribosomal protein L7/L12 (Croucher et al., 2009). The correlation of ICE family and insertion sites reflects the specificity of the integrases Int_Tn5253 and Int_ICESp23FST81 (Mingoia et al., 2007; Wyres et al., 2013).

The five ICESp23FST81-like elements contained Tn916-like insertions, including a Tn916-like, a Tn3872-like and four Tn2009-like elements (Table 3). Each of the smaller elements was integrated between orf9 and orf8 of Tn5253 and contained identical flanking and coupling sequences, indicating that each was derived from a single conjugation event and subsequent transfer of Mega-1 between Tn916 and Tn2009 (Figure 6). GA44378 containing Tn916 integrated into the ICESp23FST81-like element was macrolide resistant due to an unlinked Mega-1.III insertion (Table 3; Table S2). ICESp23FST81 was originally described in a 23F, ST81 isolate (Croucher et al., 2009). Five of six ICESp23FST81 family elements were also in CC81, however only two were serotype 23F (Table 3). The others were serotype 19A (n = 1) and 19F (n = 2) and represent capsule switching from a PCV7 vaccine serotype to a non-vaccine serotype and a different vaccine serotype, respectively. This conclusion is supported by the close relationship between the isolates, all but one of which was also serotype 23F and CC81 (Figure 6). The lone exception was GA13494, a serotype 14, CC156 isolate that contained Tn2009, indicating horizontal transfer of the composite element between the clones (Figure 6). These were isolated over a 9-year period spanning the pre-PCV7 and post-PCV7 eras (1995–2004). This demonstrated clonal expansion and the long-term stability of the clone and the Tn5253-like element harboring Tn2009.

In the other composites identified, the Tn916-like element was inserted between orf43 and orf20 of Tn5253. Tn3872 was inserted in this locus in GA47597, a serotype 3, CC180 isolate from 2006 (Table 3). Tn3872 was also inserted downstream of orf20 and contained the same flanking DNA sequences, however, the right junction was flanked by orf18. This indicated that GA47597, a serotype 3, CC180 isolate from 2006, also contained Tn3872 which had a sequence identical to GA47502 flanking the left junction of Tn3872 (Figure 6). Tn6002 was also found flanked by orf20 on the left and orf18 on the right in two isolates (Table 3; Figure 5). Like GA47597, GA19101 was a serotype CC180 isolate suggesting a common integration event for a Tn916-like element (Figure 6). However, Tn6002 in both of these isolates was integrated in the opposite orientation from Tn2009 and Tn3872 (Figure 6). This suggests separate conjugation events (one for each orientation) followed by homologous recombination of one junction to create a chimeric insertion site (i.e., orf20/orf18). In both isolates with Tn3872, the element was integrated into orf20/orf18 in the same orientation as Tn2009 in this location (Table 3). Again, Tn2009 and Tn6002 in this location had identical coupling sequences indicating a single event integration event of a Tn916-like element followed by acquisition of nested macrolide resistance elements by homologous recombination. Finally, orf43/orf18 was identified as the site of insertion of Tn6002 in four isolates (Table 2). Tn5253 is also flanked by these genes and with conserved coupling sequences, indicating that each is descended from a single conjugation event followed by integration of the erm(B) element between orf43 and orf20. These data suggest that de novo integration Tn916-like elements into pneumococcal chromosome backbones or ICE was infrequent relative to the rate of horizontal transfer by transformation and recombination.

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.