Patterns of Global Genetic Diversity of Streptococcus pneumoniae Inferred Based on Archived Multilocus Sequence Typing Data

1. Introduction

Streptococcus pneumoniae is a Gram-positive, lancet-shaped bacterium that contributes significantly to microbial diseases worldwide. S. pneumoniae is commonly found in the nasopharynx of humans where 40-50% of healthy children and 20%-30% of healthy adults are carriers [1]. However, S. pneumoniae can cause a diversity of pneumococcal diseases (PD) such as Otitis Media (OM), pneumonia, sinusitis, septicemia, and meningitis [2,3,4]. For example, children ≤ 4 years old have a high incidence rate of OM [5], with certain local populations in Sub-Saharan Africa and Asia showing 100% incidence rates for children ages 1-4 [6]. S. pneumoniae can spread from the nasopharyngeal into the lower airways and cause pneumonia. The estimated worldwide incidence rate of community-acquired pneumonia ranges between 1.5-14 cases per 1000/year [7]. Pneumonia can be caused by multiple viral, bacterial, and fungal pathogens and is responsible for approximately 2.5 million deaths worldwide/year [8,9]. Among bacterial pneumonia, about 90% is caused by S. pneumoniae [10]. In addition, S. pneumoniae can pass the blood-brain barrier (BBB) to cause pneumococcal meningitis. Meningitis is estimated at over 1.2 million cases worldwide each year, and if untreated, fatality rate can reach 70% and leaving 1 in 5 survivors with permanent brain damage [11].

Among the >100 known serotypes in S. pneumoniae, 16 cause approximately 90% of PDs worldwide [12]. For example, Hausdorff et al. [13] found that serotypes 4, 6, 9, 14, 18, 19, and 23 cause approximately 70-88% of PDs in young children in North America, Europe, Africa, and several regions in other continents. While much is known about the association between serotypes and PD, the relatively low number of serotypes in this species offered limited discriminating power to distinguish strains and was insufficient for many types of epidemiological investigations about this pathogen. In 1998, a multilocus sequence typing (MLST) scheme was established for genotyping strains of S. pneumoniae. In this scheme, DNA sequences at the following seven housekeeping genes were recommended to identify the genotypes for individual isolates: d-alanine-d-alanine ligase (ddl), signal peptidase I (spi), transketolase (recP), shikimate dehydrogenase (aroE), glucose-6-phosphate dehydrogenase (gdh), glucose kinase (gki), and xanthine phosphoribosyltransferase (xpt) [12]. At each locus, a unique sequence was assigned a distinct allele and combinations of allele types (ATs) at the seven housekeeping loci were used to provide a sequence type (ST) for each isolate [12]. The MLST scheme showed much greater discriminatory power over serotyping and provided a consistent method for many researchers around the world to expand on the MLST database of S. pneumoniae over the following 25 years.

As of November 2022, there were 74346 S. pneumoniae isolates with MLST data at the seven loci deposited at PubMLST.org. These isolates came from diverse geographical populations and spanned over nine decades. Using MLST data, previous studies have described the associations between STs, serotypes, antibiotic resistance, and/or PD prevalence among various demographic groups [e.g., 14-22]. However, most of the studies have focused on local level, with a few included isolates and data from multiple countries. Interestingly, several studies identified the emergence of non- pneumococcal conjugate vaccine (PCV) serotypes following the introduction of PCVs. As of 2020, the WHO documented that 149 countries have introduced PCV [10]. However, the effects of PCV introduction on ST distributions at the global level have not been investigated.

In this study, we aimed to use the archived global MLST data on S. pneumoniae to address several population genetic and epidemiological questions. First, how are STs distributed at the continental and national/regional levels? Second, are geographic populations genetically differentiated at the continental and national/regional levels? Third, what is the effect of PCV implementation on the distributions of STs? Because PCVs were designed to target the dominant serotypes, we expect that certain STs will be more impacted by PCV implementation than other STs. And lastly, is there evidence for recombination within individual continental populations and the global population of S. pneumoniae?

2. Materials and Methods

2.1. Collection and Organization of Archived S. pneumoniae MLST data

Allelic profiles, geographical, temporal, and ecological information, and nucleotide sequences of all STs were retrieved from the S. pneumoniae PubMLST database (pubMLST.org). The selection of gene loci for MLST, the primers, and the criteria for identifying ATs and STs were described in the original MLST scheme by Enright and Spratt [12]. Meta data of the isolates associated with each ST, including their geographical location, ecological niche, and time of isolation were all organized into Microsoft Excel. In total, information for 74346 isolates were retrieved in November 2022 for this study [16].

2.2. Phylogenetic Analysis of STs and ATs

Allele sequences at each of the 7 housekeeping loci and allelic profiles for each ST were retrieved from PubMLST [16]. To determine the relationships among ATs at each housekeeping locus, DNA sequences of all ATs were imported into Molecular Evolutionary Genetics Analysis (MEGA) software and were aligned using the Multiple Sequence Comparison by Log-Expectation (MUSCLE) tool as described in the MEGA X-Help manual [23,24]. The aligned AT sequences for each locus were then utilized to construct a neighbour-joining (NJ) tree using a Kimura 2-parameter (K2P) model [23]. The newly constructed NJ trees for each locus were imported into the Interactive Tree of Life (iTOL) online tool for better visualization [25].

To infer the relationships among STs, their allelic sequences at the seven loci were concatenated using MobaXterm based on their AT numbers for each ST [26]. The newly concatenated sequences were then aligned through Multiple Alignment using Fast Fourier Transform (MAFFT) [27]. Similar to the AT tree construction, the aligned ST sequences were utilized to construct a NJ tree using a K2P model [23]. The newly constructed ST tree was imported into iTOL along with the continental, ecological, and temporal data for each ST [25].

2.3. Population Genetics Analysis

To determine genetic differences among populations of S. pneumoniae, three sample types were analyzed: (i) the entire S. pneumoniae population (as of November 2022), (ii) the STs associated with serotypes used in the PCV-13 vaccination program (which contained STs associated with serotypes 1, 3, 4, 5, 6A, 6B, 7F, 9V, 14, 18C, 19A, 19F, and 23F), and (iii) non-PCV-13 STs (which contained STs associated with serotypes 6C, 8, 9N, 11A, 12F, 15BC, 22F, 23B, 35B, and 38) that are part of the top 50 most frequent STs. Our focus on the top 50 most frequent STs was to identify how some of the STs might have changed due to PCV implementation.

To investigate the population genetic patterns, the retrieved S. pneumoniae allelic profiles were imported into Microsoft Excel for each of the three sample types and three datasets were created for each type: (i) ecological niche-based dataset, (ii) country and continent-based dataset, and (iii) temporal and continent-based dataset. For the country and continent-based datasets, countries with <10 isolates were combined with geographically the closest country within the same continent as one sample. For the temporal/continental datasets, isolates prior to the year 1999 were omitted due to small sample sizes. Furthermore, each dataset consisted of a non-clone corrected (NCC) type and a clone-corrected (CC) type. The NCC datasets contain all isolates pertaining to that set. The CC datasets contain only one isolate of each ST as a representative within each sub-population.

An analysis of molecular variance (AMOVA) test was performed to estimate genetic variation within each population and to quantify genetic differences among sub-populations within each dataset, as described in the GenALeX 6.5 manual [28]. Because GenALeX 6.5 as well as all other known population genetic analyses programs only allowed a maximum of 999 ATs for each locus, rare alleles at two loci with > 999 ATs each (xpt and ddl) were combined to fit the program requirement. Specifically, rare alleles that were closely related to each other based on their phylogenetic relationships were combined. In our analyses, only STs with represented isolates and required metadata as reported in the pubMLST website were included in our analyses. Furthermore, allelic profiles of all STs were inputted into Multilocus V1.3. Genotypic diversity and linkage disequilibrium analysis were determined using formulas as described in the Multilocus V1.3 manual [29]. Datasets were divided by continent and a PCV dataset was created for PCV-13 STs and non-PCV-13 (top 50) STs.

3. Results

3.1. Global Distribution of Isolates and Sequence Types

As of November 2022, 74346 isolates belonging to 17915 STs of S. pneumoniae were successfully extracted from PubMLST and imported into Microsoft Excel for analysis. The geographic distributions of these isolates and STs are shown in Table 1. The extracted information included isolates from 6 continents and 115 countries/regions. Based on continental origins, Europe contributed the highest proportion of isolates (32.403%), followed by Africa (25.68%), Asia (23.274%), North America (12.731%), South America (3.891%), and Oceania (1.246%). Among the 74346 isolates, 576 (0.775%) had no associated geographic information. At the country level, the highest number of isolates came from the USA (11.679%), followed closely by South Africa (11.55%) and the UK (8.569%). The USA, South Africa, and the UK also had the highest number of uniquely represented STs among the countries with 2005, 2013, and 2323 unique STs respectively. The lowest number of isolates were found in Sudan, Cuba, Greenland, Monaco, Azerbaijan, and Wallis and Futuna Islands, with one isolate from each country/region (Table 1).

Table 2 shows the breakdown of the top 50 globally most frequent STs. ST180 was the most frequently observed and made up 1.98% of the entire dataset. Overall, the top 50 most frequent STs made up 30.914% of isolates in the entire dataset. Among the top 50 STs, 39 were represented by those in PCV-13, including ST180, while 11 were non-PCV-13 STs. Similar to the skewed ST frequencies, the frequencies of ATs at each locus were also highly skewed (Table 2). For aroE, AT7 was the most frequent (21.80%, represented in 22% of STs). For gdh, AT5 was the most frequent (21.56%, in 24% of STs). For gki, AT2 was most frequent (21.39%, in 18% STs). For recP, AT1 was most frequent (27.57%, in 26% STs). For spi, AT6 was most frequent (37.40%, in 34% STs). For xpt, AT1 was most frequent (19.55%, in 14% STs). Finally, for ddl, AT14 was most frequent (33.90%, in 36% STs). We note that while for six of the seven loci, the most frequent AT at the whole population level was also found in the highest number of STs. However, for gki, AT4 was represented by more STs (n(ST)=11), but fewer isolates (n=4 568) than AT2.

Table 3 shows the geographical, ecological, and temporal distributions of 17915 STs of S. pneumoniae. At the geographical level, 36 STs (n=15263 isolates) were found in all six continents, 27 STs (n=4970 isolates) were found in five of six continents, 55 STs (n=5069 isolates) were found in four of six continents, 184 STs (n=6517 isolates) were found in three of six continents, 737 STs (n=7796 isolates) were found in two of six continents, and 15929 STs (n=33607 isolates) were found in one continent only. The remaining STs were represented by isolates with no associated geographical information. Europe had the highest number of unique STs (n=7986, n(ST)=4 813). At the ecological level, 21 STs (n=4665) were found in both ecological niches (clinical and veterinarian) and 16998 STs (n=69119) were found in one ecological niche only. The ecological niche information was unknown for the remaining 896 STs. Most isolates and the majority of STs were from clinical sources (n(ST)=16987, n=69098 isolates). Based on temporal information, the isolates were grouped into the following five time intervals: Int 1: <1949; Int 2: 1950-1969; Int 3: 1970-1989; Int 4: 1990-2009; Int 5: 2010-2022. At the temporal level, only 1 ST (n=717) was found in all five-time intervals, 12 STs (n=1390) were found in four of five intervals, 72 STs (n=10690) were found in three of five intervals, 1584 STs (n=31349) were found in two of five intervals, and 13627 STs (n=20660) were found in one interval only. Temporal information was unknown for the remaining 2619 STs. Majority of STs in the dataset were collected between 1990-2009 only (n(ST)= 6727). Majority of isolates were collected between 1990-2022 (n=31132).

Geographical and temporal distributions of top five STs in representative countries and time intervals were further highlighted (Table 4). Geographical distributions were shown at the country level where countries with the most collected isolate data on each continent were shown. Interestingly, we see diverse ST distributions among countries, with different countries often having different most frequent STs. However, there were some overlaps among countries. For example, ST199 was among the most frequent in the USA, UK, and Australia. Together, the top 5 most frequent STs in each country make up approximately 12-22% of their respective populations, with the highest proportion found in the USA (22.19%). We note that among the common STs in Table 4, ST4133 was associated with a nontypeable serotype and PCV-13 STs were represented in 17 of the 30 most frequent country-ST combinations. At the temporal level, ST199 and ST180 were among the most frequent, found in 4 of the 6 time intervals, and PCV-13 STs were found 18 of the 30 time interval-ST combinations.

Table 5 shows the global frequencies of ATs. The length of housekeeping loci ranges from 405-486 bp. ddl has the greatest number of ATs with 1157. The most frequent ATs for each housekeeping loci were found in both PCV-13 and non-PCV-13 (top 50) STs. We found AT6 of the spi gene to be represented the most among the dataset (35.31%).

3.2. Phylogenetic Analysis

In total, 17915 STs of S. pneumoniae were concatenated at seven housekeeping loci and aligned for phylogenetic analysis. Figure 1 shows a NJ tree of 17915 STs of S. pneumoniae using a K2P model. Geographical, ecological, and temporal data were respectively added through iTOL to demonstrate their distributions. The NJ tree shows both divergent and close relationships among STs. The most divergent STs were ST13466 and ST9066. Since this tree (Figure 1) contains too many STs for clear visualization of each ST, further phylogenetic analysis was performed for subsets of STs. Specifically, we analyzed the top 50 globally most frequent STs (Supplementary Figure 1). The most divergent STs among the top 50 most frequent were ST558 and ST416. Geographically, though highly abundant, ST4414 was found only from Asia. As the fifth most abundant, ST191 was the only ST found among all five time intervals. Among the top 50 STs, we found no clear phylogenetic clustering between those covered by PCV-13 and those belonging to non-PCV-13 STs. Phylogenetic analysis was also performed for PCV-13 STs and those not covered by PCV-13 (Supplementary Figure 2 and Figure 3). With only a few exceptions, in general, STs that were associated with the same serotypes were found to be phylogenetically closely related to each other among both PCV-13 ST and non-PCV-13 (top 50) ST trees. Phylogenetic trees showing the relationship between AT sequences among seven housekeeping loci can be seen in Supplementary Figures S4–S10.

3.3. AMOVA

Due to the large sample size of the non-clone-corrected data in the total sample, AMOVA test of the non-clone-corrected data among geographic populations were conducted separately for among countries within individual continents. For clone-corrected data, a two-level hierarchical AMOVA was conducted to assess population genetic differentiation at both the continental and country/region levels. Table 6 shows the results of AMOVA test among geographical populations. For clone-corrected data, genetic variation within countries contributed 98.515% to the total observed genetic variation. Genetic variation among countries and among continents contributed 1.074% and 0.411% respectively to the total genetic variations. For the non-clone-corrected data, genetic variation within countries contributed 97.914%, 96.975%, 98.086%, 97.322%, 97.731%, and 96.080% of the observed genetic variation in Africa, Asia, Europe, North America, South America, and Oceania, respectively. The remaining variances were attributed to among countries within individual continents. The among continents, among countries, and within countries contributions were all statistically significant with p<0.01.

Further, two-level hierarchical AMOVA testing was conducted for both PCV-13 and non-PCV-13 (top 50) STs at the geographical level (Supplementary Tables S1 and S2). Only non-clone-corrected data was used for both datasets due to lack of sample size for clone-corrected data. For PCV-13 STs, genetic variation within countries contributed 90.172% to the total observed genetic variation. Genetic variation among countries and among continents contributed 7.436% and 2.392% respectively (Supplementary Table 1). For non-PCV-13 (top 50) STs, genetic variation within countries contributed 90.568% of the total observed genetic variation. Genetic variation among countries and among continents contributed to 7.539% and 1.893% respectively (Supplementary Table 2). The among continents, among countries, and within countries’ contributions were all statistically significant with p<0.01.

Table 7 shows the results of AMOVA test at the temporal level. In this analysis, populations were defined at intervals beginning in 1999, then every three years until 2022. For clone-corrected data, genetic variation within time intervals contributed 99.180% to the total observed genetic variation. Genetic variation among time intervals and among continents contributed 0.189% and 0.631% respectively. For non-clone-corrected data, genetic variation within time intervals contributed 99.556%, 99.393%, 99.546%, 97.947%, 98.891%, and 98.519% to the observed genetic variation in Africa, Asia, Europe, North America, South America, and Oceania, respectively. The remaining amounts were attributed to among time intervals. The among continents, among time intervals, and within time intervals’ contributions were all statistically significant with p<0.01.

The two-level hierarchical AMOVA testing was also conducted for PCV-13 and non-PCV-13 (top 50) populations at the temporal level (Supplementary Table S3). Due to sample size issues, only six countries were used for analysis as described in Supplementary Table S3. In addition, only the non-clone-corrected data was used due to the small sample sizes for clone-corrected data. Here, genetic variation within time intervals contributed 94.784% to the total observed genetic variation for the six countries. Genetic variation among time intervals and among continents contributed 4.173% and 1.043% respectively. All values were statistically significant with p<0.01.

Table 8 shows the results of AMOVA test between the two ecological niches at the global level. For clone-corrected data, genetic variation within ecological niches contributed 98.478% to the total observed genetic variation, with between ecological niches contributing 1.522%. For non-clone-corrected data, only European samples were used due to small sample size for veterinary populations in other continents. Genetic variation within ecological niches and between ecological niches contributed 87.922% and 12.078% to the total observed genetic variation respectively in the European samples. All values were statistically significant with p<0.01.

3.4. Effects of PCV on Population Structure

Next, we analyzed the effects of PCV implementation on relative frequencies of the top 50 most common STs associated with PCV-13 and non-PCV-13 (Figure 2). The associations between STs and PCV-targeted serotypes are shown in Supplementary Table S4. We tested whether there was a rise of non-PCV-13 STs after PCV implementation. A line histogram was produced and as seen in Figure 2, A noticeable trend was observed following 1981-1985, where there was a decrease in relative frequency of PCV-13 STs, followed by a steady frequency following global PCV-13 implementation, and a notable decrease following global PCV-13 implementation. The opposite was observed with non-PCV-13 (top 50) STs where there was an increase after 1981-1985, followed by a steady frequency upon global PCV implementation, and a significant increase in relative frequency after global PCV-13 implementation. However, we note that the relative isolate frequencies taken from 1936-1975 and 2021-2022 did not follow the trends as expected based on PCV-13. This was likely due to small sample sizes collected during these time intervals (Figure 2, Supplementary Table S5).

Aside from the global analyses, changes in ST frequencies related to PCV-13 implementation was also conducted in countries where sample sizes were large (>600 isolates; Supplementary Figure S11). Since PCV-13 implementation differed among some of the countries and continents, we hypothesized that the changes in relative frequencies of STs over time may differ among geographic regions. Indeed, the frequencies of PCV-13 associated and not associated STs differed among geographic regions, with the implementation of PCV-13 showing different but notable effects among the US, South Africa, and UK (Supplementary Figure 11 and Supplementary Table 6 and Table 7). However, the effects of PCV-13 implementation on ST changes were not observed in Brazil and Australia (Supplementary Figure 11). At the continental level, the relative frequencies of PCV-13 associated STs increased from ~2.3% during the 2001-2005 interval to about 9% from 2006-2020 in Asia. However, in Europe, there was a large decrease from ~21% during 2001-2005 to ~10% during 2006-2020. In North America where PCV-13 was implemented the earliest, the decreases in PCV-13 associated STs happened earlier, from ~15% during 1996-2000 to 5.5% or less after 2000 (Supplementary Table 6).

3.5. Recombination & Linkage Disequilibrium

Evidence for potential recombination within individual geographic populations of S. pneumoniae was investigated using two indicators: phylogenetic incompatibility and index of association. In these two analyses, only clone-corrected data were used and due to the large sample size at the global level crushing the program Multilocus, only samples at individual continents were analyzed. An additional dataset was analyzed to test for recombination among STs associated with PCV-13 (Table 9). All datasets, with the exception of the PCV dataset, showed 0% phylogenetic compatibility, indicative of recombination between pairs of loci. However, all index of association values were statistically significant with p<0.01, rejecting the hypothesis of random recombination.

4. Discussion

In this study, we extracted and analyzed the genetic, geographical, ecological, and temporal data of the global S. pneumoniae isolates deposited within PubMLST up to November 2022. The extracted global population included 74346 isolates representing 17915 STs. The isolates were distributed among 6 continents and 115 countries. Among the six continents, Europe contributed the most isolates (32.403%). Among the 115 countries, USA had the most isolates (11.679%). The patterns of distribution of isolates among continents and countries likely reflected sampling and research efforts by scientists working on this pathogen as well as the differences in relative importance of S. pneumoniae to public health and in available resource for research among countries and continents.

Our analyses revealed a diversity of STs within most countries and continents. Indeed, most STs (15929 of 17915 STs) were found in only one continent. However, although only 36 STs were shared among all six continents, they represented over 20.53% of all the isolates in the database, consistent with recent gene flow and that a relatively small number of STs dominated the global S. pneumoniae population. Our results are similar to previous research results that showed a few serotypes causing most pneumococcal diseases [12]. A similar pattern was seen across temporal scales where a relatively small number of STs persisted across decades. Among these, ST191 was collected across all time intervals (Figure S1). ST191 is among the top 50 most frequent STs and is associated with serotype 7F, a PCV-13 target [30,31]. ST191 was a frequent ST prior to PCV-13 implementation. However, introduction of PCV-13 has resulted in a reduction of ST191 frequency over time. While incomplete sampling could have contributed to the low number of shared STs across the analyzed time frames, the persistence of ST191 suggests its significant adaptability in human populations. Interestingly, between 1990-2022, only 1532 STs were collected, representing about 10% of all STs in the database but these STs account for 41.87% of all the isolates in the dataset. The overrepresented STs indicate that a relatively small number of STs are highly transmissible and/or pathogenic to humans. Further genetic and genomic studies of isolates of these frequent STs (e.g., ST199, ST180, and ST81) could help reveal the potential mechanisms for their broad distribution and/or persistence in humans.

Among geographical regions, we observed differences in the prevalence of PCV-13 and non-PCV-13 (top 50) STs. In Thailand, 4/5 of the most frequently represented STs belonged to serotypes targeted by PCV-13, including ST4414 (of serotype 19F), a top 50 PCV-13 ST, found so far only in Asia, with 99.53% of the strains of this ST reported from Thailand (Figure S1). The high frequency of STs associated with PCV-13 is consistent with the lack of PCV immunization program in Thailand [32]. In South Korea, the introduction of PCV-13 resulted in the emergence of novel STs due to serotype replacement [33]. A similar pattern was found in other countries such as South Africa, the UK, and the USA where the implementation of PCV into their immunization programs [32] caused STs not associated with PCV-13 to increase in their relative prevalence in these countries.

Phylogenetic analyses of all STs revealed both divergent STs and closely related STs (Figure 1). Interestingly, our phylogenetic analyses revealed that two STs, ST14613 and ST17858, in the current MLST database for S. pneumoniae were very divergent from most other STs in the database. Further analyses through BLAST revealed that these two STs belonged to Streptococcus mitis, a close relative of S. pneumoniae. Significantly, S. mitis has been found to be a source for capsular polysaccharide variation for S. pneumoniae through horizontal gene transfer and can contribute to vaccine escape in S. pneumoniae [34].

Similar to the overall phylogenetic relationships among all STs, the top 50 most frequent STs in the database also showed both significant divergence and close relatedness among each other. Though small clusters of PCV-13 associated STs were found, STs associated with PCV-13 were overall inter-mixed with those not associated with PCV-13 (Supplementary Figures S2 and S3). Interestingly, 7 STs (ST156, ST81, ST63, ST193, ST320, ST695, and ST172) were each associated with more than one serotype covered by PCV-13. Three STs (ST199, ST162, and ST156) contained strains that were initially associated with PCV-13 prior to PCV implementation but changed their serotypes following PCV implementation. Previous literature has also identified ST199 switched from serotype 19A to 15B following PCV-7 implementation [35]. These findings suggest evidence of capsular switching, a process by which a new capsule operon is acquired through horizontal gene transfer (HGT) [36]. Due to flexibility and ease of DNA uptake in pneumococci, high rates of recombination via frequent HGT within the cps operon could have occurred, causing increases in non-PCV-13 serotypes [37]. Further genomic comparisons should help identify the relationship between STs and cps operon allele variation and how mutation and/or HGT might have impacted capsular switching and serotype changes in the global population.

Analysis of molecular variance revealed statistically significant genetic differentiations among geographical, temporal, and ecological populations of S. pneumoniae, rejecting the null hypothesis of no genetic differentiations among sub-populations in each the three types of analyses. However, in all cases, both the non-clone-corrected datasets and the clone-corrected datasets revealed that most genetic variations were found within subpopulations rather than between subpopulations. Among these, the highest inter-subpopulation differentiation (~12%) was observed for the non-clone corrected data between clinical and veterinary samples. Further investigation revealed that 46/112 veterinary isolates represented STs that were unique to veterinary niches. Among those 46 isolates, 26 belonged to ST6937 and 10 isolates belonged to ST6934. Localized selection and clonal expansion have likely contributed to the observed patterns. Differences between clinical and veterinary/environmental samples have also been reported in other microbial pathogens such as Campylobacter jejuni and Streptococcus agalactiae [38,39].

Population genetic analyses between pre- and post- PCV implementation revealed a decrease in relative frequency of PCV-13 associated STs and an increase in frequency of non-PCV-13 (top 50) STs as well as significant genetic differentiations between them based on AMOVA (Figure 2). These global patterns are consistent with previous literature in England and Wales which found that PCV implementation led to a rise in non-PCV serotype infections [40]. Interestingly in our analysis, we noticed a significant rise in some non-PCV-13 (top 50) STs that were not associated with new serotypes targeted by the upcoming PCV-20 (serotypes 35B, 23B, 9N, 15BC, and 6C). On the other hand, STs associated with serotypes 10A, 15B, and 33F, which are targeted by the new PCV-20, were not among the top 50 most represented STs. Thus, our results suggest that PCV-20 is not as optimized as it could be. We believe that further vaccine development and optimizations should take the global ST distributions into account to produce higher-valent PCVs. At the national level, we observed unique trends for countries that have implemented PCV and countries that have not (Figure S11). In Thailand, which has not implemented PCV into their national immunization program (NIP), the relative frequency of PCV-13 associated STs increased over time, reaching almost 100%. In countries that have implemented PCV into their NIP, a reduction in PCV-13 associated STs is typically seen after PCV-13/PCV-10 introduction together with an increase in non-PCV-13 (top 50) STs. This observation is consistent with the previous finding that ST4414, a unique PCV-13 ST only found in Asia, has contributed to increased PCV-13 ST frequency in Thailand (Figure S1).

S. pneumoniae has shown to be capable of importing DNA through transformation and homologous recombination, generating recombinant genotypes [41]. In this study, we found that strains of the same ST or closely related STs were often associated the same serotypes. However, differences were also found, consistent with recombination and horizontal gene transfer in natural populations of this species. In our analyses, all continent populations of S. pneumoniae showed evidence of recombination among the seven loci used for MLST. However, none of the analyzed populations showed evidence of random recombination, even in the clone-corrected samples.

In conclusion, global analyses of published MLST data for S. pneumoniae revealed great diversity and distribution of STs spatially, temporally, and phylogenetically. Analysis of molecular variance quantified the genetic variations within and among S. pneumoniae subpopulations. Implementation of PCV revealed an impact on both PCV-13 STs and non-PCV-13 (top 50) STs globally. We also demonstrated non-random association among alleles at the seven MLST loci. It is important to note that MLST data for S. pneumoniae only reflects data collected and deposited by researchers. The distributions observed in this study may underrepresent the actual distributions of S. pneumoniae. As new PCVs are being produced to reduce the incidence of pneumococcal diseases among global populations, the rise of non-PCV serotypes due to serotype replacement and capsular switching caused by recombination and horizontal gene transfer remains a major concern. At present, two new higher-valent PCVs, PCV-15 and PCV-20, are being implemented across the world to cover a wider range of serotypes causing pneumococcal diseases. However, information on the coverage and effectiveness of these new vaccines at the global level is scarce. At present, PCV-13 is still recommended by the US Centers for Disease Control and Prevention (CDC) for infants and younger children. However, as shown in our analyses, after PCV-13 implementation, new serotypes and new STs not associated with PCV-13 emerged and spread, reduced the effectiveness of PCV-13 [42]. New PCVs should be designed to target the most prevalent STs and serotypes that are not covered by previous PCVs in order to maximize the efficacy of the new vaccines. Indeed, genotype information from global isolates should be continuously deposited into the pubmlst database for monitoring the potential spatial and temporal patterns related to S. pneumoniae and pneumococcal diseases. Such data, in combination with those on host and environmental factors related to S. pneumoniae and pneumococcal diseases could help develop effective public health policies against this important human pathogen (Xu 2022).