The collection of S. pneumoniae clinical isolates analysed in this study was representative of the allelic (ST) diversity of CC156 strains. The following information on the 66 strains used in this study is listed in Table 1: name, serotype, ST, geographical origin, the number of MLST alleles in common with ST156 (the CC156 eBURST predicted founder) and with ST4945, data source and strain source. Overall, the strains analyzed belong to 41 CC156 STs and are associated with 14 capsular types. The isolate collection was composed of 8 strains for which either complete or draft genomic sequence was available and 58 additional isolates for which new sequence data were generated (see below). The 58 isolates already characterized for serotype and ST were kindly provided by many laboratories (complete list in the “Acknowledgments” section).

Bacteria were grown overnight at 37°C in 5% CO₂ on Tryptic Soy Agar plates (TSA) (Becton Dickinson) supplemented with 10 mg/l colistine, 5 mg/l oxolinic acid and 5% defibrinated sheep blood. Genomic DNA extractions for the 58 isolates were performed by using the Wizard Genomic DNA purification kit following the manufacturer’s instructions (Promega). PCR amplifications for 96-MLST were performed on genomic DNA as briefly described below, while the presence of PI-1 and the PI-1 clade were determined as reported elsewhere [27].

The 96-MLST loci and the amplification primers for the 96-MLST schema reported in Tables S1 and S2 were selected and designed as detailed elsewhere [34]. Briefly, the 96 loci were identified within the S. pneumoniae core genes by comparing the available complete and draft genomes of 39 S. pneumoniae invasive strains and selected for being short variable regions (ranging from 847 to 209 bp) flanked by conserved regions. For the 8 complete and draft genomes (Table 1) of S. pneumoniae the nucleotide sequences of the 96 loci were retrieved by in silico analysis; a locus was considered to be absent when a primer aligned against the genome sequence with more than three mismatches to guarantee consistency with the results of the PCR amplifications. For this reason, the two loci SP0180 and SP0181 were not extracted from the SP195 genome.

The 96-MLST loci of the 58 S. pneumoniae isolates were PCR amplified in a microtiter 96-well plate by using a unique amplification condition (55°C annealing temperature, 1 minute elongation) with the primer set reported in Table S2 [34]. The PCR products were then purified with paramagnetic beads (Agentcourt AMPure XP, Beckman Coulter Genomics, USA) and both strands sequenced with the same amplification primers by an ABI 3730xl DNA Analyzer (Applied Biosystems, USA). Chromatogram traces were edited and assembled with Vector NTI Advance 11 (Life Technologies). Consensus sequences were aligned using MUSCLE 3.8.31 [37]. The SP0278, SP1785, SP1909, SP2141 and SP2198 loci were not amplified on 2, 4, 1, 1 and 3 strains respectively. The complete nucleotide sequence dataset is provided in the File S1.

Sequences were converted into allelic profiles (Table S3), by assigning a unique ID number to each allele. When a primer pair did not amplify, the absent locus was assigned the ID number “0″. With this choice, strains sharing the same deletions or divergent sequences in the primer region were considered more similar than what would be by treating these loci as missing data, thus allowing the use of the information in the phylogenetic analysis.

Hierarchical clustering was performed using the package Cluster v1.13.1 [38] of the software package R v2.12.0 (www.r-project.org/). Distances between strains were computed using the function “Daisy” with Gower’s distance, counting the number of differences between allelic profiles. An agglomerative hierarchical clustering of the data was performed using the function “Agnes” with “average” (unweighted pair-group average method – UPGMA) method. Support of the clustering was assessed using bootstrap.

Minimum Spanning Tree analysis was performed using PHYLOVIZ [39].

For each of the 96 loci the sequences were aligned using MUSCLE [37]. Aligned sequenced where concatenated, and phylogenetic analysis was performed using Mega5 [40] with the Neighbor Joining method [41].

ClonalFrame V1.1 [42] was run on the aligned sequences of the 89 loci that where present in all strains. Seven independent runs of 10⁶ iterations were performed, and the first half of each run was discarded. When compared, the distributions of some of the parameters from the seven runs failed to fulfill the convergence requirement of a Gelman and Rubin statistics below 1.2, suggesting that the evolutionary model implemented by ClonalFrame gives a poor description of these data. The posterior samples of the tree topologies of the seven runs were combined, and a consensus phylogenetic network was generated using SplitsTree v4.10.

CC156 (predicted founder ST156) was identified as the largest clonal complex in the S. pneumoniae MLST database (accessed on 15th January 2012) by running the eBURST algorithm with the default settings on the complete dataset which, at the time of the analysis, comprised 7146 distinct STs (ST1-ST7160). CC156 presented a complex and heterogeneous structure (Figure 1A and [25]), which did not change by using other algorithms such as goeBURST [23], encompassing 13.8% (986 STs) of the total STs in the database (as a comparison the second largest CC, CC320, included only 235 STs, 3.3% of the total). Interestingly, based on analyses performed in 2008, CC156 comprises STs which formerly belonged to distinct complexes: CC124, CC146, CC162, CC176, and CC392. These CCs were associated with different serotypes and differed for the presence of PI-1 and for PI-1 clade [27]. As detailed elsewhere [25], [27], [43], CC162 strains were associated with serotypes 9V and 14 and were PI-1 positive (PI-1 clade I), CC124 and CC392 strains were associated with serotypes 14 and 17F and were PI-1 negative, CC146 strains were associated with serotype 6B and PI-1 positive (clade II) and CC176 strains were generally associated with serotype 6B and 23F and heterogeneous for the PI-1 presence (clade II).

To investigate whether the merger of multiple CCs into a single CC was due to the identification of one or multiple new STs, the eBURST algorithm was iteratively executed by progressively excluding the most recently identified STs from the analysis. This resulted in the observation that ST4945 (a new allelic combination), occupied a central position within the eBURST CC156 graphic representation (Figure 1 and Figure S1), and had been sufficient to induce the merger of three formerly distinct CCs: CC162, CC124 and CC176 into one larger CC (Figure 1A). At the time of this analysis, only two ST4945, both serotype 17F from two different countries were present in the MLST database. Noteworthy, the CC124 and CC176 identified by excluding ST4945 from the analysis comprised the formerly separated complexes CC124 and CC392, and CC176 and CC146, respectively [27].

To further investigate whether the strains comprising the newly formed CC156 had a common evolutionary descendent (ST156) or whether ST4945 strains contained a combination of genetic alleles from different lineages (as suggested by their MLST profile), we analyzed a panel of 66 representative strains of different STs belonging to CC156 (see materials and methods section, Table 1, Figure 1 and Figure S1). The strains isolated in different countries belonged to 41 STs and shared between zero and six MLST alleles with ST156 (Figure 1B) despite belonging to the same CC. As shown in Figure S1A, all of the strains displayed the PI-1 distribution expected based on the lineage partitioning present before the introduction of ST4945 (see Figure 1A and [26], [27], [43]). The two ST4945 strains were also part of the collection as well as single and double locus variant (SLV, DLV) strains of ST4945 (Figure S1B).

The 66 strains were typed by using the 96-MLST schema [34]. Following amplification and sequencing, the sequences were converted into allelic profiles (Table S3). The minimum number of distinct alleles was identified in the SP0841 locus (6 alleles), while the maximum number of alleles was identified in the SP334 and SP2194 loci (28 alleles). By performing hierarchical clustering of the 66 strains using the number of loci with different alleles as a measure of their genetic distance it was possible to identify ten genetically distinct lineages supported by high values of bootstrap (Table 1, Figure 2 and Figure S2A). As shown in Figure 2, strains with the same ST were always assigned to the same lineage, while being SLV was not predictive for the assignment to the same lineage. Indeed, ST392, ST789 and ST4404 were all SLVs (in different MLST alleles) of ST4945 strains but only ST392 was assigned to the ST4945 lineage; the SLVs ST273 and ST385 were assigned to different lineages as well as the SLVs ST2218 and ST176; finally, ST361 and ST2218 (in a DLV relationship) both SLV of ST171 were allocated into the same lineage, different from that of ST171 (Figure 2). In addition, similar lineages could be identified also by using ClonalFrame (Figure S2B). Small differences (e.g. the position of strains 08809744 and 08802447) with respect to the results obtained with hierarchical clustering could be attributed to the fact that, despite having run 7 independent simulations for a total of 7*10⁶ iterations, incomplete sampling of the tree topologies was achieved. Remarkably, a phylogenetic analysis performed with the same dataset by aligning the 96 loci concatenated sequences, identified the same lineage pattern (Figure S3).

As expected, by computing a hierarchical clustering analysis with the MLST allelic profiles for the same set of 41 STs, SLV STs were grouped in the same lineage (Figure S4). This was also true when the phylogenetic analysis was performed with the sequences of the seven MLST loci (Figure S5). However, for the 7-MLST analysis the groupings obtained with the two methods were not exactly the same (i.e. for the two SLV STs 4966 and 5420).

Interestingly, the most remarkable differences in the hierarchical clustering between 96-MLST and 7-MLST were that lineages “d” and “h” (the latter comprising ST4945), “f” and “e”, “a” “b” and “c”, and “i” and “j” clustered together in the 7-MLST, thus justified by their STs proximity within the CC156 eBURST representation.

The partitioning of CC156 into ten distinct lineages by hierarchical clustering was further evaluated by performing a Minimum Spanning Tree (MST) analysis of the 96-MLST alleles (Figure 3), and by creating a visual diagram of the allelic assortment within and across the identified lineages for both the 7 and the 96-MLST profiles (Figure 4). As shown in Figure 3, by applying a threshold of 75 loci (i.e. by cutting those links in the MST diagram that connected strains differing by more than 75/96 loci), seven distinct lineages could be identified. Notably, these lineages corresponded to those obtained with the hierarchical clustering analysis, although the strains belonging to lineages “a”, “b”, “c”, “e” and “f” were grouped by the MST analysis in two larger groups, “a+b+c” and “e+f”. Indeed, the lineages “a” and “b”, the lineages “b” and “c” and the lineages “e” and “f” were connected in the MST by strains differing in 56, 60 and 50 loci, respectively. In addition, the two ST4945 strains were unequivocally assigned to lineage “h” and the ST392 strain (SLV of ST4945 in the MLST classification) was the closest to the ST4945 strains according to 96-MLST, with 25 different alleles.

In order to visualize whether the allelic differences within and among lineages were concentrated in specific regions of the chromosome, and thus likely attributable to single recombination events, the strains were partitioned based on the lineages identified by hierarchical clustering analysis (Figure 4, columns). The alleles of the 7-MLST and 96-MLST loci (Figure 4, rows) were then colour coded by assigning the same colour to identical alleles within a locus. Since the rows in Figure 4 are ordered according to the position in the chromosome blocks of consecutive loci showing a colour pattern discordant with the rest of the chromosome or with other lineages are indicative of localized intra or inter-chromosomal differences, respectively. Overall, the distinction among genetic lineages was not due to differences present in specific hyper-variable regions of the chromosome. In fact, the differences were dispersed among the 96 loci and, apart from a few exceptions, each lineage contained a specific repertoire of alleles. Interestingly, lineages “d”, “e” and “h”, comprising ST4945 strains and SLV and DLV strains of ST4945 (black and blue arrows in Figure 4) shared several alleles in the loci probed by the 7-MLST typing schema, but were clearly distinct at the level of the 96-MLST loci. Besides, lineages “a”, “b” and “c” were the most heterogeneous, presented several unique loci and numerous alleles in common (thus justifying the fact they were classified into a unique cluster by MST and the lower support of the branching separating them in the hierarchical clustering and in the consensus network obtained from the posterior sampling of the tree topologies generated by ClonalFrame). A similar situation was shared by lineages “e” and “f”, although these two lineages presented a lower number of unique alleles.

With the exception of lineages “a” and “b”, all of the strains belonging to the same lineage were also closely related in the eBURST graphic visualization of CC156 (Figure 5). In detail, ST2218, spatially separated from the other lineage “a” strains, is both SLV of lineage “a” ST172 and of lineage “b” ST176, but had 70 and 31 96-MLST loci in common with the closest ST172 strain 08b09744 (lineage “a”) and ST176 strain PB011 (lineage “b”), respectively. Besides, ST5613, ST4966 and ST5420 were separated from the other lineage “b” strains. Noteworthy, based on MST analysis, strains belonging to these three STs were closer to one another than to other strains in the collection, being the ST171 strain BG02112 the closest to the ST4966 strain 1681 (46/96 alleles in common).

In addition, as shown in Figure 5, the majority of the lineages were homogeneous for the presence/absence of PI-1 and, when PI-1 was present, strains within the same lineage carried, (as a confirmation of their phylogenetic proximity) the same PI-1 clade.