Introduction

Once termed “the captain of the men of death” by Sir William Osler in the 20^th century, Streptococcus pneumoniae, or the pneumococcus, remains a significant cause of life-threatening invasive diseases, such as sepsis, pneumonia, and meningitis¹, despite the widespread implementation of life-saving serotype-specific polysaccharide pneumococcal conjugate vaccines (PCVs)². The high capsular diversity, amounting to ∼100 serotypes with varying potential for colonising the nasopharynx and causing invasive diseases^3,4 and mortality^5–7, demonstrates the remarkable genetic diversity of the pneumococcus. Such high genetic diversity reflects the plasticity of the pneumococcal genome, which readily swaps genetic material with other strains through recombination and horizontal gene transfer⁸. However, the impact of pneumococcal genomic variation on its phenotypes, especially those associated with disease pathogenesis and clinical manifestations, remains less understood.

The pneumococcus possesses an arsenal of virulence factors. The extracellular polysaccharide capsule is the main virulence factor critical for transmissibility, survival, and pathogenicity in the human host^9–11. Recent studies in the pneumococcus and other bacterial pathogens suggest that variability in growth kinetics impacts disease pathogenesis, clinical manifestation, and antibiotic therapy. For example, differential pneumococcal growth features appear to correlate with fitness for colonisation and invasive disease,¹² clinical manifestation^12–14 and response to environmental conditions^14–16. Studies of other bacterial pathogens suggest that growth rate influences gene and protein expression^17,18, adhesion¹⁹, bacterial competition²⁰, plasmid replication²¹, bacteriophage dynamics²², and antibiotic accumulation and efficacy^23–25. Similarly, experimental and mathematical modelling studies have demonstrated that bacterial growth correlates with bacterial load and host-mediated cell lysis and impacts dispersal to different host cells and tissues during systemic infection²⁶. Previous studies of pneumococcus and other bacteria have suggested the influence of environmental factors^15,16 and bacterial-specific factors, including capsular serotype^14,27 and the presence of plasmids²⁸ on bacterial growth phenotypes. These growth phenotypes may consequently impact invasive pneumococcal disease pathogenesis^12–14. However, the specific pneumococcal genomic variations affecting growth kinetics remain less understood.

Studying bacterial growth kinetics sheds light on the influence of genes, genomic variation, and growth media conditions on bacterial fitness^29,30. However, studying bacterial growth kinetics in vivo remains intrinsically challenging due to the complexity of the conditions in the hosts, although innovative approaches are being developed^31,32. On the other hand, in vitro methods provide a valuable alternative for assessing the impact of specific bacterial and environmental factors on bacterial fitness. Although these in vitro methods have revealed remarkable insights on bacterial fitness, these experiments have mostly yielded strain-specific rather than population-level insights due to the analysis of a small number of isolates from few genetic backgrounds.

Here, we studied the impact of pneumococcal genetics on in vitro growth kinetics to understand the influence of bacterial genetic variation on pneumococcal pathogenesis. Specifically, we performed high-throughput quantification of the in vitro growth kinetics and whole genome sequencing of 348 invasive pneumococcal isolates from the Netherlands to determine the impact of genomic variation on growth kinetics. We conducted growth curve experiments on 348 isolates in blood-like media conditions, and absorbance was monitored over 15 hours. Growth features were determined from these curves using statistical models, and we assessed the impact of capsular serotype and genetic background on the estimated growth features and tested their association with the phylogeny of the isolates. In addition, we used a genome-wide association study (GWAS) approach to identify specific genomic loci associated with the growth features independent of the genetic background. Finally, we calculated the variability in the growth features explained by pneumococcal genetics. Our study suggests that the capsular serotype, genetic background or lineage, and a combination of genetic loci may influence pneumococcal growth kinetics, a fundamental trait potentially influencing nasopharyngeal colonisation, invasive disease pathogenesis, and clinical manifestations.

Results

Population structure of a diverse collection of Dutch invasive pneumococcal isolates and derivation of their growth kinetics

We aimed to understand the impact of the capsular serotype, genetic background, and specific genomic loci, independent of the clonal population structure, on the pneumococcal growth features (Fig. 1 and Supplementary Data 1). Therefore, we investigated the in vitro growth kinetics of 348 pneumococcal isolates collected from patients with invasive pneumococcal diseases via the Pneumococcal Bacteraemia Collection Nijmegen (PBCN) cohort study at two hospitals in Nijmegen, the Netherlands before (2000-2006) and after (2007-2011) the introduction of the seven-valent pneumococcal conjugate vaccine (PCV) [PCV7] into the Dutch paediatric immunisation programme³³. Our isolates represent 34 capsular serotypes, 40 pneumococcal lineages, defined based on the global pneumococcal sequence cluster (GPSC) nomenclature³⁴, and 99 sequence types (STs) defined by the pneumococcal multilocus sequence typing (MLST) scheme³⁵. An interactive and annotated maximum likelihood phylogeny of the isolates used in this study is available online at the Microreact platform for open data visualisation and sharing for genomic epidemiology (https://microreact.org/project/wxmdblfgprepwugvbbtaq7). We estimated a Simpson diversity index of 0.932 (95% confidence interval [CI]: 0.924 to 0.940), 0.937 (95% CI: 0.928 to 0.945), and 0.958 (95% CI: 0.949 to 0.966) for the serotypes, lineages, and STs, respectively. The high richness and diversity of pneumococcal serotypes and lineages highlighted a substantial population-level diversity of the study isolates in the Netherlands.

Fig. 1.

Capsular serotype and genetic background significantly impact pneumococcal growth kinetics

To explore the impact of capsular serotype and genetic background or lineage on the pneumococcal growth features estimated from a fitted growth curve, namely lag-phase duration, maximum growth density (H_max), maximum growth change (ΔH), and average growth rate (r) (Fig. 1, 2, and Supplementary Data 1), we assessed the relationship between these growth features and serotypes and lineages associated with the study isolates (Fig. 3a). In general, we found significant variability of the growth features across serotypes: average growth rate, maximum growth density, maximum growth change, and lag phase duration (P<0.0001; Kruskal-Wallis test) consistent with previous studies^14,15 (Fig. 3b-e and Supplementary Data 2). We observed similar patterns across different pneumococcal lineages (Fig. 3f-i and Supplementary Data 2).

Fig. 2.

Fig. 3.

Recombination and horizontal gene transfer may lead to the replacement of part or the entire capsule biosynthesis locus resulting in the production of a polysaccharide capsule distinct from the original one, a phenomenon called capsule-switching³⁶. Consequently, some lineages are associated with multiple serotypes, and similarly, a single serotype can be expressed in distinct genetic backgrounds³⁷. Examples of lineages containing multiple serotypes included GPSC3, GPSC4, and GPSC7, and likewise, serotypes, including 3, 19A, and 19F, were expressed in multiple genetic backgrounds, which also express other serotypes (Fig. 3a). Overall, the average growth rate for serotypes expressed in different genetic backgrounds broadly showed similar patterns. For example, lineages expressing serotype 14 had consistently higher average growth rates. In contrast, lineages expressing serotype 3 showed lower rates, consistent with the notion that capsule production interferes with pneumococcal growth¹⁴ (Fig. 3b). Therefore, we investigated the impact of the capsular serotype and genetic background on the estimated pneumococcal growth features.

First, we identified lineages that expressed multiple serotypes and serotypes found across multiple lineages. To ensure robust statistical comparison, we only focused on serotype-lineage combinations with at least three sequenced isolates in our study. We found two lineages, with enough sequenced isolates for comparison, GPSC3 and GPSC7, which expressed multiple serotypes (Fig. 3j and Supplementary Data 2). We noted higher growth rates for serotype 33F than serotype 8 (P=0.0007; Kruskal-Wallis test).

In contrast, we found no statistically significant differences in the growth rates between serotypes 23A and 23F in lineage GPSC7 (P=0.2568; Kruskal-Wallis test). We speculated that the absence of statistically significant differences between serotypes 23A and 23F reflected their presence in the same serogroup, therefore, had similar capsule properties. Interestingly, we observed similar patterns for the maximum growth density, maximum growth change, and lag phase duration features between serotypes on the GPSC3 and GPSC7 lineages (Fig. 3k-m and Supplementary Data 2). Next, we investigated the growth rates of the same serotypes expressed in distinct genetic backgrounds. Although some serotypes showed similar growth rates across distinct genetic backgrounds, such as serotype 3 (P=0.4599; Kruskal-Wallis test), others showed different growth rates across different lineages, as seen with serotypes 14 (P=0.0006; Kruskal-Wallis test) and 19F (P=0.0105; Kruskal-Wallis test). We noted similar patterns for the other growth features: maximum growth density, maximum growth change, and lag phase duration (Fig. 3n-q, Supplementary Data 2).

These findings demonstrate that the capsular serotype and genetic background significantly impact pneumococcal growth kinetics. The resultant growth features likely depend on the interaction between these factors, whereby some combinations may result in more pronounced growth effects than others.

Pneumococcal growth kinetics show heritability and correlation with the phylogeny of the isolates

Based on the findings that both capsular serotype and genetic background are associated with growth kinetics, we hypothesised that the estimated pneumococcal growth features exhibit a high phylogenetic signal. We assessed the phylogenetic signal for each growth feature using Pagel’s λ statistic³⁸, which evaluates the correlation between given traits and the phylogenetic tree of the isolates, with values close to zero and one indicating weaker and stronger phylogenetic signals, respectively. By reconstructing the estimated growth features on the internal nodes of the phylogeny, we found high but variable associations between the mapped features and different genetic backgrounds, which suggested potentially strong phylogenetic signals (Fig. 4, Supplementary Fig. 2 and 3, and Supplementary Data 3).

Fig. 4.

The amount of the phylogenetic signals based on the normal transformed growth features was highest for the average growth rate (Pagel’s λ=0.72, 95% CI: 0.61 to 0.81; P<2.08×10^-64), maximum growth change (Pagel’s λ=0.60, 95% CI: 0.40 to 0.76; P<1.53×10^-16), and maximum growth density (Pagel’s λ=0.73, 95% CI: 0.57 to 0.84; P<4.18×10^-26), while the lag phase duration showed a lower phylogenetic signal (Pagel’s λ=0.30, 95% CI: 0.17 to 0.47; P<5.42×10^-16) (Fig. 5a). These findings suggest that there is a heterogeneous correlation of the growth features with the phylogeny.

To quantify the variability in the pneumococcal growth kinetics explained by the genotype, we calculated the narrow-sense heritability (h²) for the four normal-transformed estimated growth features and the maximum growth rate from Arends et al¹³. We found strong, though not equal, evidence for a genetic basis for all the estimated growth features: average growth rate (h²=0.76, 95% CI: 0.68 to 0.85), maximum growth change (h²=0.80, 95% CI: 0.71 to 0.89), maximum growth density (h²=0.88, 95% CI: 0.83 to 0.94), and lag phase duration (h²=0.29, 95% CI: 0.14 to 0.43) (Fig. 5b, Supplementary Data 3). To further understand the impact of the capsular serotype and genetic background, we fitted a linear mixed effects model for each growth parameter with serotype and lineage (defined by GPSCs³⁴) as random effects to estimate the amount of variance of the growth features explained by the serotype (h²_serotype) and lineage or GPSC (h²_GPSC) calculated as the ratio of the variance of the serotype and lineage components relative to the total variance, including the residual or environmental effects. We calculated h²_serotype=0.47 and h²_GPSC=0.11 for the average growth rate, h²_serotype=0.13 and h²_GPSC=0.28 for the maximum growth change, h²_serotype=0.29 and h²_GPSC=0.46 for the maximum growth, h²_serotype=0.08 and h²_GPSC=0.39 for the lag phase duration feature.

Fig. 5.

Altogether, these findings demonstrate a substantial impact of the serotype and lineage on the variation of the pneumococcal growth features, whereby most serotypes tend to have similar characteristics regardless of the genetic background. Similarly, some genetic backgrounds tend to show similar phenotypic patterns irrespective of the serotypes, but the serotype appears to have a more significant influence on the growth characteristics.

No individual genomic loci are associated with pneumococcal growth features

Having found moderate to high narrow-sense heritability of the growth features, we then undertook a GWAS analysis to identify specific genomic loci associated with these features independent of the strains’ genetic background using a linear mixed-effects model. We performed separate GWAS analyses of the maximum growth change, maximum growth density, maximum growth change, and lag phase duration phenotypes derived from the growth curve. We conducted the GWAS using variable-length k-mer sequences or unitigs to capture genetic variation in the pneumococcal genomes at different levels, including SNPs, multi-allelic sites, and insertions and deletions^39,40. We found no individual unitigs associated with the derived growth feature phenotypes (Fig. 6a, Supplementary Data 4). Furthermore, repeating the GWAS analyses adjusting for the serotype revealed no statistically significant hits associated with the derived growth phenotypes, suggesting that the variants in the pneumococcal strains are correlated with the growth kinetics independent of the serotype (Fig. 6b, Supplementary Data 4). Overall, the observed and expected statistical significance values in the QQ-plots for our study revealed no apparent issues with the population structure of the isolates, which is a significant confounder in similar studies⁴¹ (Fig. 6c, d). Overall, these findings suggest that no individual genomic variation influences the pneumococcal growth rates, suggesting that multiple genetic loci and lineage effects modulate pneumococcal growth kinetics.

Fig. 6.

Average in vitro growth rates are not strongly correlated with a serotype’s invasiveness or its prevalence in carriage or disease

Previous studies have elucidated that bacterial growth rate may influence colonisation and disease¹⁶, including systemic dissemination during infection²⁶ and disease manifestations^13,27. Pneumococcal strains that are prolific growers may display better survival and replication in systemic tissues than slow growers without vaccine-induced non-serotype-specific immunity. We used the isolates collected before the introduction of the PCV7 vaccine in the Netherlands⁴² to assess the association between in vitro growth rate of the isolates and their population-level characteristics: prevalence in invasive disease and nasopharyngeal carriage and invasiveness, i.e., odds of being isolated in invasive disease compared to the carriage state. Our first hypothesis was that the average in vitro growth rate for each of the pneumococcal serotypes is positively associated with the prevalence of serotypes in invasive diseases. Our analysis revealed no evidence of an association between the prevalence of the serotypes in invasive disease and the in vitro growth rate (Spearman correlation: -0.30, P=0.0856) (Fig. 7a, Supplementary Data 4).

Fig. 7.

Second, we hypothesized that the average in vitro growth rate of the serotypes is positively associated with the prevalence of the serotypes in the nasopharyngeal carriage as fitter strains would grow faster than less fit strains. However, contrary to our hypothesis, we found a negative association between the prevalence of the serotypes in carriage and the in vitro growth rate (Spearman correlation=-0.53, P=0.0098) (Fig. 7b, Supplementary Data 4). We repeated the analysis after excluding serotype 3 and found no evidence of an association between growth rate and carriage prevalence (Spearman correlation=0.38, P=0.0925) which suggests that serotype 3 alone was driving the previously observed association. Finally, we tested the hypothesis that the growth rate is positively associated with the invasiveness of the serotypes. We found no association between the invasiveness and the in vitro growth rates of the serotypes (Spearman correlation: 0.35, P=0.1024) (Fig. 7c and Supplementary Data 4). These findings suggest that growth rates alone may not be sufficient to explain the prevalence and invasiveness of pneumococcal serotypes, which suggests the influence of other pneumococcal, microbiome (e.g., viral infection), and host factors in determining pneumococcal disease susceptibility.

Discussion

Bacterial growth kinetics impact several traits, including gene and protein expression^17,18, antibiotic accumulation and efficacy^23–25, cell surface adhesion¹⁹, intra- and inter-specific competition²⁰, response to environmental conditions^14–16, plasmid replication²¹, bacteriophage dynamics²², colonisation and disease pathogenesis¹², including clinical manifestation^12–14. Whilst previous studies have attempted to understand bacterial, host, and environmental factors influencing bacterial growth, a major drawback has been a narrow focus on a small subset of strains⁴³. Therefore, the insights gained from these studies have mostly been for specific strains, as they ignored potential variability across different genetic backgrounds. These population-level studies are critical for genetically diverse bacterial pathogens with highly structured populations³⁴, substantial antigenic diversity⁴⁴, and recombination rates, such as the pneumococcus^8,36,45. Here, we examine the association between pneumococcal genetics and growth kinetics using a combination of high-throughput quantification of in vitro pneumococcal growth, whole-genome sequencing, and population genomic analyses employing a bacterial GWAS approach. Notably, our findings reveal a strong influence of the capsular serotype compared to the genetic background, and the absence of specific individual genomic loci modulating pneumococcal growth kinetics, particularly the in vitro growth rate. These findings may partially explain why some serotypes and lineages colonise and cause invasive diseases better than others and may have implications for developing measures to prevent and control pneumococcal infections.

Our findings suggest that the capsular serotype and genetic background influence pneumococcal growth kinetics. The findings are generally consistent with results from experimental studies by others^13–15, which used capsule-switched strains to exclude the impact of the genetic background and showed an alternation of growth dynamics due to the capsule, although there may be potential differences in the growth media and conditions. Specifically, encapsulated strains showed higher growth rates than unencapsulated strains in nutrient-rich conditions¹⁴. These findings suggest that the capsule is crucial for pneumococcal growth, but its effect depends on nutrient availability – promoting and interfering with growth when nutrients are abundant and limited, respectively^14,27. A recent study by Tóthpál et al.¹⁵ suggested that the serotype partially explains pneumococcal growth and suggested the genetic backgrounds, which vary in their ability to colonise the nasopharynx and cause invasive diseases^4,46–48, may also play a role. However, the authors did not fully quantify the impact of the genetic background since the isolates were not sequenced¹⁵. Our study confirms the strong influence of both the pneumococcal capsular serotype and genetic background or lineage on the growth kinetics of the isolates. The influence of the serotype appeared to be higher than the genetic background for the average growth rate, while the other growth features appeared to show an opposite trend. These findings suggest that the respective influence of the capsular serotype and lineage varies by specific derived pneumococcal growth features.

Interestingly, although our findings revealed that the capsular serotype and genetic background influence pneumococcal growth kinetics, statistical analysis showed that the variance attributed to the capsular serotype was higher than the genetic background. Specifically, the capsular serotype exhibits a more pronounced growth rate, maximum growth change, and maximum growth density than the genetic background. In contrast, the genetic background had a higher effect on the lag phase duration than the serotype. We acknowledge that our comparison of growth features in the same genetic background could also mean a large amount of genetic variation for some of the lineages, some of which could be responsible for changes in growth. These findings illustrate the importance of the capsular serotype and genetic background on pneumococcal growth kinetics, which may partially explain their variable prevalence in carriage and disease.

Previous studies have revealed individual bacterial genomic loci that impact other phenotypes, including pathogenicity^49–57, virulence^58–61, host, niche, and environmental adaptation^52,54,62, conjugation⁶³, and antimicrobial resistance^64–71, independent of the genetic background. Our GWAS findings revealed no elevated signals in specific genomic loci, suggesting that individual variants may not significantly influence pneumococcal growth kinetics independent of the genetic background and serotype. Previous studies have shown that mutations in other capsule biosynthesis genes (wchA or cpsE) encoding a protein that links activated glucose phosphates to lipid carriers⁴⁴ also influence the growth dynamics of the pneumococcus⁴³. However, our results are contrary to our speculation that genetic variation in the capsule polysaccharide synthesis genes would be strongly associated with pneumococcal growth features. These findings suggest that specific genomic loci in the pneumococcal genome do not strongly influence pneumococcal growth kinetics, hinting at the potential impact of the overall genetic background and several loci possibly with small effects individually. While it’s possible that such follow-up studies are needed to provide additional support for the importance of these variants on pneumococcal growth kinetics.

In contrast to findings in a previous study¹⁴, we found no robust evidence for the association between the in vitro growth rates and the population-level prevalence of the serotypes in nasopharyngeal carriage. However, this may have occurred as a result of a limitation of our study, since we studied a collection of randomly sampled invasive isolates, and the growth experiments were conducted in conditions similar to that of blood, while the Hathaway et al. study¹⁴ used a non-random selection of both carriage and disease isolates. It is likely that isolates have adaptations to the environmental conditions in which they were isolated, and so in vitro growth characteristics for any given isolate may be different under experimental conditions simulating different niches. In addition, the in vitro growth environment cannot fully represent in vivo conditions and fails to consider the impact of selective pressures coming from the serotype-specific responses (antibody) or non-specific responses (phagocytic) and potential interactions with other pathogens. Presumably, the lack of association between the invasiveness and prevalence of serotypes in invasive disease and carriage, and the growth rates could reflect that the conditions in the in vitro growth environment do not fully represent in vivo conditions. Collectively, these findings may suggest that these and other factors may substantially influence the complex pneumococcal disease pathogenesis in vivo.

The strength of this study is the analysis of an extensive collection of well-phenotyped natural isolates and the combination of phenotyping, and estimation of growth features using statistical models and statistical genetics to determine the association between the growth features and genetic variation across the entire genome. As previously noted by Tothpal et al.¹⁵, the culturing history of the isolates may impact the growth kinetics. In our study, we consistently cultured all the isolates after isolation from patients using standardised culturing procedures (single assay, identical initial concentration, environmental conditions or temperature, and equipment to measure the optical density and repeat measurements to spot technical variability¹³). Compared to the present study, future studies would benefit from larger sample sizes and may find more results as we only explained some of the heritability of the growth features. We also acknowledge that the methods used to estimate heritability may also sometimes perform poorlyon some bacteria species, therefore, more robust methods are needed to further validate our findings⁷². Because we conducted our experiments for a relatively short duration (∼15 hours), we could not robustly characterise the kinetics at the stationary phase due to incomplete data or the observed high variability of its kinetics. Therefore, we recommend increasing the culturing time from 15 to >24 hours in follow-up studies, as done elsewhere^14,15, to characterise the stationary phase kinetics of the cultured isolates robustly. Furthermore, since the pneumococcus invades different niches and systemic tissues with varying conditions, including oxygen concentration and temperature^12,14–16, future studies should investigate the effect of these environmental conditions on pneumococcal growth and identify pneumococcal genomic variants associated with growth kinetics for each condition. Considering that our dataset contained only a third of the >100 known pneumococcal serotypes globally³; due to the geographical variation in the distribution of pneumococcal serotypes⁷³, future studies should investigate all serotypes to provide further insights on the impact of different serotypes in pneumococcal growth kinetics. Additionally, considering that we generated the phylogeny of the entire pneumococcal species rather than isolates belonging to a single strain or lineage, as defined by several approaches^34,74,75, we did not remove recombinations using widely used methods designed for clonal bacterial populations, such as Gubbins⁷⁶. Therefore, there is a possibility that some branch lengths in the phylogeny used for the ancestral trait reconstruction analysis may be less accurate although the topology is usually robust to recombination⁷⁷. Furthermore, there may be other variants that potentially influence the pneumococcal growth kinetics which are regarded as not statistically significant in this study based on our conservative genome-wide significance threshold. As more and more similar studies are conducted, particularly those utilising much larger datasets, potential variants associated with pneumococcal growth kinetics, may be discovered.

Our analysis of population-level genomic and high-throughput in vitro pneumococcal growth data suggests that the capsular serotype and lineage influence pneumococcal growth. However, although some growth features are highly heritable, there appears to be no specific genetic variants which strongly influence these phenotypes suggesting that a combination of genomic loci with small effects individually or the entire genetic background may play a greater role in pneumococcal growth kinetics. Overall, our study provides a proof-of-concept for combining phenotyping and population genomics to understand the genetic basis for pneumococcal growth, an approach which can be applied to study other bacteria of clinical significance and those used in bioengineering, genetics, and the food industry^78–80. Ultimately, we anticipate that applying our approach will have implications for devising strategies to reduce life-threatening bacterial infections and improve food production and safety.

Methods

Ethics statement

The isolates used in this work were collected through an observational cohort study Pneumococcal Bacteraemia Collection Nijmegen (PBCN)³³. The non-applicability of the Medical Research Involving Human Subjects Act (WMO) on the study design was confirmed by the Regional Ethics Committee, and the local medical ethics committees of the participating hospitals in the Netherlands approved the study procedures.

Sample characteristics and whole-genome sequencing

Three hundred and forty-eight (n=348) pneumococcal isolates collected from patients with invasive pneumococcal diseases through the PBCN observational cohort study³³ were available for the study. We included a single isolate per patient as invasive pneumococcal infections are typically clonal, so a single isolate is sufficient. The consecutive blood isolates were collected at two hospitals in Nijmegen, the Netherlands, before (2000-2006) and after (2007-2011) the introduction of the PCV7 vaccine into the Dutch paediatric immunisation programme³³. Most of the bacteremia cases (97%) concerned adults. Additional information on the PBCN study and collected clinical variables have been previously described^33,81. To generate whole genome sequences, we sequenced a single isolate per patient using an Illumina HiSeq sequencer (Illumina, San Diego, CA, USA). We then assembled the generated reads using SPAdes genome assembler (version 3.14.0)⁸², and evaluated the draft assemblies using assembly-stats (version: 1.0.1) (https://github.com/sanger-pathogens/assembly-stats). In addition, we used a combination of metrics, including genome size and the number of contigs, to assess the quality of the draft assemblies for inclusion in the analysis.

Measurement of in vitro pneumococcal growth kinetics and derivation of growth features

We measured the pneumococcal growth kinetics for slightly over 15 hours as described by Arends et al.¹³. In short, we incubated a frozen aliquot of standardised pneumococcal inoculum in 1.5 ml of blood-like medium (50% M17, 50% casamino acids tryptone [CAT] medium, 0.25% glucose). We measured growth kinetics at 37°C and 5% CO₂ every 10 min for over 15 hours at OD₆₂₀ with use of a humidity cassette by microplate reader (Spark 10M, Tecan, Switzerland). We added 15 μL of inoculum was added to 1.5 mL of prewarmed rich growth medium (45% M17, 45% CAT, 0.225% glucose, 10% fetal calf serum [Greiner Bio-one], 26 U/mL catalase [Sigma-Aldrich C1345]) in each well of a sterile flat-bottomed 48-well plate (Nunclon Surface, Nunc, Denmark). The growth medium was supplemented with catalase. We chose these nutrient-rich conditions to bring forward differences in the growth features of the pneumococcus. We performed six repeat measurements on three separate days for each isolate. To understand the growth kinetics of the invasive Dutch pneumococcal isolates, we fitted a modified logistic growth curve function, i.e., L₀ + L_max/(1 + e^−r(t−t₀)), separately to the growth data for each isolate for each replicate, using the “optim” function in stats (version 4.3.2) R package (https://www.R-project.org/). We then averaged the inferred values for each replicate to obtain average estimates for each parameter. The parameter L₀ captures the initial cell density measured as absorbance or optical density at a wavelength of 620 nm (OD₆₂₀); L_max is the maximum cell density; t₀ is the initial time, and r is the average growth rate. We fitted two curves in total for each isolate based on the averaged data from the six replicates using the limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS-B) optimisation algorithm by minimising the residual sum of squares⁸³. We fitted the first curve on all the data points and the second curve on the data measured until the start of the end of the stationary phase or the beginning of the death phase, whichever came earlier. We defined the end of the stationary phase or the start of the death phase as the time when the optical cell density decreased by >5% from the maximum cell density. The curve fitted based on the data measurements collected up to the end of the stationary phase or death phase, or the end of the experiment, if none of these phases had been reached, captured the growth dynamics better than the curve fitted using all the data points to accurately estimate the growth rate and the maximum cell density corresponding to the peak of the pneumococcal growth curve, especially for growth curves with observably short stationary phases and rapid autolysis. Therefore, parameters for the downstream analyses were inferred from the curves fitted to the reduced dataset.

We determined four growth phenotypes for each pneumococcal isolate estimated from the fitted growth curve, namely the lag-phase duration, maximum growth density (H_max), i.e., captured by parameter L_max, maximum growth change (ΔH) corresponding to the difference of the maximum and minimum density (L_max-L₀), and average growth rate (r). We defined the lag-phase duration as the time when the optical cell density was at least 1.5-fold higher than the starting initial cell density, as described by Arends and colleagues¹³. Due to the challenges of characterising the kinetics at the stationary or death phase, partly due to the high strain-to-strain variability or short study period, we did not derive the growth features capturing these phases, including the duration of the stationary and death phases. The logistic growth curves showed an excellent fit for the data, which provided an opportunity to investigate the impact of capsular serotype and genetic background or lineages and specific genomic loci, independent of the genetic background, on the estimated growth parameters or features.

Molecular typing of sequenced isolates

Molecular typing of the isolates to determine capsular serotypes and sequence types (STs) was done in silico using Pathogenwatch (https://pathogen.watch/), a global platform for genomic surveillance⁸⁴. We determined the capsular serotypes and STs using SeroBA (version v1.0.1)⁸⁵ and the multilocus sequence typing (MLST) scheme for the pneumococcus³⁵, respectively. Pneumococcal lineages were defined based on the global pneumococcal sequence cluster (GPSC) nomenclature³⁴ using PopPUNK (version 1.1.0)⁸⁶. In addition, we checked the species of each sequenced isolate using Speciator (version 3.0.1) in Pathogenwatch.

Construction of multiple sequence alignments and phylogenetic tree

The generated draft assemblies for each isolate were mapped to the pneumococcal ATCC 700669 reference genome (GenBank accession: NC_011900) using Snippy (version 4.6.0) (https://github.com/tseemann/snippy). We identified 155,305 SNPs in the generated alignment using SNP-sites (version 2.3.2)⁸⁷ and then used them to construct a maximum-likelihood phylogeny using IQ-TREE (version 2.0.3)⁸⁸. We used the alignment of variable sites by mapping to a reference genome to generate a maximum-likelihood phylogeny using IQ-TREE as it has been to produce more accurate phylogenies⁸⁹. The general time reversible (GTR) nucleotide substitution model and Gamma heterogeneity of the substitution rates among sites. We assessed branch support with ultrafast bootstrap approximation using “-B 1000” option⁹⁰. We optimised bootstrapped trees by nearest neighbour interchange on bootstrap alignment by specifying the “--bnni” option. We rooted the phylogenetic tree of the pneumococcal isolates using Streptococcus oralis as an outgroup (not shown in the phylogeny). An annotated and interactive version of the phylogeny is available on the Microreact web tool⁹¹ (https://microreact.org/project/wxmdblfgprepwugvbbtaq7).

Quantifying the phylogenetic signal

To assess the correlation between the estimated growth features and pneumococcal phylogeny, we used models of continuous character evolution to estimate the phylogenetic signal using Pagel’s λ³⁸ implemented in the “phylosig” function in phytools package (version 2.1.1)⁹². Pagel’s λ transforms the internal branches relative to the terminal branches and generates values ranging from zero to one, with smaller values indicating no phylogenetic signal or independent evolution of the phenotype from the phylogeny. In comparison, higher values close to one show a high phylogenetic signal or a decreased phenotype variability among genetically similar taxa at the tips of the phylogeny. First, we transformed the estimated growth features to a normal distribution using a rank-based inverse normal transform implemented in the RNOmni package (version 1.0.1.2) (https://cran.r-project.org/web/packages/RNOmni/). Next, we mapped the estimated continuous growth feature values onto the phylogeny using fast estimation of maximum likelihood for ancestral states⁹².

Bacterial GWAS

We identified the presence and absence patterns of unitig sequences To capture genetic variation in the pneumococcal genomes included in this study for the GWAS analysis. We identified unitig sequences from 31bp k-mers using Bifrost (version 1.0.6)⁹³. First, we inferred unitig sequences from a compacted De Bruijn graph constructed using k-mers of length 31 using draft assemblies of all the pneumococcal isolates included in this study. We then generated De Bruijn graphs for each isolate and queried them with the unitig sequences found in the entire dataset to identify the presence and absence of each unitig found in each sequenced genome. We considered a unitig present when 100% of the k-mer was presented in the De Bruijn graph. We combined the presence and absence patterns of the unitig sequences from all the isolates to construct a single matrix.

Since the linear mixed effects models require the tested continuous phenotype to be normally distributed, we used the normalised growth features after applying the rank-based inverse normal transformation using the RNOmni package (version 1.0.1.2) (https://cran.r-project.org/web/packages/RNOmni/). We performed GWAS analysis of the pneumococcal growth features using robust linear mixed effects models implemented in Pyseer (version 1.3.7-dev)⁹⁴. A kinship matrix specifying the variances and covariances of the random effects for each isolate was generated using the “similarity_pyseer” script in Pyseer. Using this kinship matrix, Pyseer fits an unobserved random effect for each isolate to account for the clonal population structure in the GWAS. We excluded unitigs with <5% minimum allele frequency and >95% maximum allele frequency. The inferred statistical significance for each variant was adjusted to correct for multiple testing using the Bonferroni correction. We estimated the total number of unique unitig presence and absence patterns in the dataset using a script included in Pyseer (https://github.com/mgalardini/pyseer/blob/master/scripts/count_patterns.py), and used the obtained value when adjusting for multiple testing using the Bonferroni correction approach. Therefore, we defined the genome-wide statistical significance threshold as P <2.58×10^-07, i.e., α/G, where α=0.05 and G=194,053 is the number of unique unitig presence and absence patterns while the suggestive threshold as P<5.15×10^-06, whereby α=1. To check potential issues arising from the population structure, we generated quantile-quantile (QQ) plots from the output of GWAS analysis for each growth feature to compare the expected and observed P-values using qqman (version 0.1.7)⁹⁵. The overall proportion of phenotypic variability explained by variation in the genome (narrow-sense heritability [h²]) was estimated using GEMMA (version 0.98.1)⁹⁶ and Pyseer. The input files for GEMMA were prepared based on the unitig presence and absence data using PLINK (version 1.90b4) ⁹⁷. To identify variants associated with the growth features independently of the serotype, we repeated the GWAS by including the serotype as a covariate.

Annotation of variants identified by the GWAS analysis

To annotate the unitig sequences with statistical significance values below the suggestive and genome-wide significance threshold, we used a custom BioPython ⁹⁸ wrapper script for BLASTN (version 2.5.0+)⁹⁹ to identify genomic features in the ATCC 700669 reference genome (GenBank accession: NC_011900) and other complete pneumococcal genomes downloaded from GenBank⁹⁹. The GWAS results were summarised and visualised using Manhattan plots using R (version 4.0.3) [https://www.R-project.org/].

Other statistical analysis

We compared the growth features between serotypes and lineages using the Kruskal-Wallis non-parametric test. The distribution of the growth features was plotted using violin plots using the vioplot package [version 0.4.0] (https://CRAN.R-project.org/package=vioplot). To assess the contribution of serotype and lineage to the variability in the pneumococcal growth features, we used linear mixed-effects models implemented in the lme4 package (version 1.1.35.1)¹⁰⁰ were used, with each growth feature specified as the outcome and the serotype and lineage as random effects. The statistical significance for the linear mixed model was inferred using the lmerTest package [version 3.1.3] (https://CRAN.R-project.org/package=lmerTest). To assess the contribution of the serotype and lineage to the variability of each growth feature, we calculated the proportion of variance for the serotype and lineage components to the total variance. We calculated the average invasiveness of the serotypes based on a meta-analysis of data collected in different countries in infants before the introduction of the PCVs⁴⁷ using the random-effects model implemented in the metafor package [version 2.4.0] (https://CRAN.R-project.org/package=metafor We defined invasiveness as the odds ratio of detecting a pneumococcal serotype among isolates causing invasive disease compared to nasopharyngeal carriage. We added a value of one to cells with zero counts. Due to the unavailability of the nasopharyngeal carriage data, we imputed the prevalence of the serotypes in the carriage by dividing the prevalence in invasive disease by the average invasiveness of each serotype. Next, we used a Spearman correlation test to assess the association between the mean prevalence of serotypes in the nasopharyngeal carriage and invasive disease and average invasiveness with the average in vitro pneumococcal growth rate.

Acknowledgements

The authors would like to thank the study participants, guardians, and the clinical and laboratory staff who collected and processed the samples at various laboratories in the Netherlands. We would also like to acknowledge the support of the sequencing, core, parasites and microbes teams, the Bentley lab at the Wellcome Sanger Institute, and the Weinberger and Pitzer labs at the Yale School of Public Health for their feedback on the analysis. We acknowledge the participating hospitals and affiliated researchers for their contributions and support.

Additional information

Funding information

This study was supported by the TARGET project (grant number: JPIAMR2019-087) funded by JPI-AMR–ZonMW, The Netherlands, and by the SNSF Swiss Postdoctoral Fellowship (SNSF grant number: TMPFP3_209768). Activities at the Wellcome Sanger Institute were funded by the Bill and Melinda Gates Foundation (grant number: OPP1034556) and Wellcome Trust (2016-2021 core award grant number: 206194). The funders had no role in study design, data collection, analysis, the decision to publish, and manuscript preparation. The findings do not necessarily reflect the official views and policies of the author’s institutions and funders. For the purpose of open access, the authors have applied a CC BY public copyright licence to any Author Accepted Manuscript version arising from this submission.

Author contributions

C.C., M.I.dJ., and A.J.H.C. conceived and designed the study. M.Id.J., D.W.A., I.H.P., and A.J.H.C. collected and processed the pneumococcal isolates. S.W.L. and S.D.B. performed whole-genome sequencing, quality control and genomic analysis. C.C. performed statistical and bioinformatic analyses. D.W.A. contributed to data analysis and interpretation. M.I.dJ. and A.J.H.C. contributed to data interpretation. C.C., D.W.A., M.I.dJ., and A.J.H.C. wrote the first draft of the manuscript. C.C., D.W.A., S.W.L., I.H.P., A.Y., J.A.L., A.L.W., D.M.W., S.D.B., M.I.dJ., and A.J.H.C. reviewed and approved the manuscript.

Additional files

Supplementary Information.

Supplementary Data 1. Characteristics of the pneumococcal isolates used in this study and summary of the derived growth features or parameters.

Supplementary Data 2. Summary of the growth features across serotypes and lineages.

Supplementary Data 3. Summary of the phylogenetic signal and narrow-sense heritability of the growth features.

Supplementary Data 4. Summary of the serotype prevalence and growth rates used.