IMPORTANCEThe capsule determines >90 different pneumococcal serotypes, which vary in capsule size, virulence, duration, and prevalence of carriage. Current serotype-specific vaccines elicit anticapsule antibodies. Pneumococcus can take up exogenous DNA by transformation and insert it into its chromosome by homologous recombination. This mechanism has disseminated drug resistance and generated vaccine escape variants. It is hence crucial to pneumococcal evolutionary response to interventions, but there has been no systematic study quantifying whether serotypes vary in recombination and whether this is associated with serotype-specific properties such as capsule size or carriage duration. Larger capsules could physically inhibit DNA uptake, or given the longer carriage duration for larger capsules, this may promote recombination. We find that recombination varies among capsules and is associated with capsule size, carriage duration, and carriage prevalence and negatively associated with invasiveness. The consequence of this work is that serotypes with different capsules may respond differently to selective pressures like vaccines.
Streptococcus pneumoniae causes a high burden of invasive pneumococcal disease (IPD) globally, especially in children from resource-poor settings. Like many bacteria, the pneumococcus can import DNA from other strains or even species by transformation and homologous recombination, which has allowed the pneumococcus to evade clinical interventions such as antibiotics and pneumococcal conjugate vaccines (PCVs). Pneumococci are enclosed in a complex polysaccharide capsule that determines the serotype; the capsule varies in size and is associated with properties including carriage prevalence and virulence. We determined and quantified the association between capsule and recombination events using genomic data from a diverse collection of serotypes sampled in Malawi. We determined both the amount of variation introduced by recombination relative to mutation (the relative rate) and how many individual recombination events occur per isolate (the frequency). Using univariate analyses, we found an association between both recombination measures and multiple factors associated with the capsule, including duration and prevalence of carriage. Because many capsular factors are correlated, we used multivariate analysis to correct for collinearity. Capsule size and carriage duration remained positively associated with recombination, although with a reduced P value, and this effect may be mediated through some unassayed additional property associated with larger capsules. This work describes an important impact of serotype on recombination that has been previously overlooked. While the details of how this effect is achieved remain to be determined, it may have important consequences for the serotype-specific response to vaccines and other interventions.
INTRODUCTIONThe global mortality due to invasive pneumococcal disease (IPD) has been estimated at more than 1 million per annum, with the majority of deaths occurring in children less than 5 years old (1). Use of a 7-valent pneumococcal conjugate vaccine (PCV7) and subsequent 13-valent (PCV13) formulations has been followed by a precipitate reduction in IPD due to vaccine serotypes (VTs) in the United States, including a substantial herd effect in unvaccinated individuals (2, 3). PCVs are plainly an effective selective pressure on the pneumococcal population. The response to that pressure has been notable for the increase in prevalence of nonvaccine serotypes (NVTs). More than 90 distinct serotypes have been described (4−10), so even PCV13 protects against a small minority of antigenic diversity in this organism. Following vaccination, the prevalence of NVTs in carriage, and to a lesser extent in IPD, has increased as the NVTs took advantage of the removal of their competitors in a process termed serotype replacement (11−13). An important component of serotype replacement has been the success of “capsule switch” vaccine escape variants (14). These variants are lineages that were previously associated with a vaccine serotype but have persisted in the postvaccine environment by the introduction of genes for an NVT capsule through homologous recombination (15). The known ability of the pneumococcus to undergo homologous recombination has also been associated with the acquisition of drug resistance (16, 17) and numerous examples of capsule switch variants (15, 18−20). As a result of this recombination-mediated shuffling of serotype and genotype, when PCVs were introduced, the population already contained potential vaccine escape variants, some of which have subsequently become common and important causes of disease.The polysaccharide capsule is known to influence multiple aspects of pneumococcal biology. Serotypes vary in their nasopharyngeal carriage rates (21), and their propensity to cause IPD was usually expressed as the “invasive potential” in order to account for the various exposures to different serotypes (22−24). Moreover, these properties have been linked to basic biochemical features of the capsule: those capsules with fewer carbons in the repeat unit of the polysaccharide tend to be thicker and exhibit a more negative surface charge (25) and are associated with a longer duration of carriage and a higher carriage prevalence—an outcome that may be linked to enhanced resistance to opsonophagocytic killing (21). Those serotypes more associated with virulence meanwhile, tend to have smaller capsules and shorter durations of carriage (21, 22, 24).Despite the importance of both recombination and capsule to pneumococcal biology, whether serotypes vary in their recombination rates in nature and what drives it have not been extensively studied. Genomic studies have shown that different pneumococcal lineages plainly vary in the amount that recombination has contributed to their diversification (20, 26) but have not linked this to the different serotypes associated with each lineage. These studies have also suggested that lineages that completely lack capsule may have a higher relative rate of recombination, suggesting that the capsule can impede the uptake of DNA (26). Support for this comes from the observation that mechanisms such as biofilm formation and downregulation of capsule biosynthesis have been reported to facilitate DNA uptake and chromosomal integration in vitro (27). This finding suggested that the polysaccharide capsule may inhibit genetic exchange in encapsulated isolates either physically or structurally (28).There are hence potentially antagonistic roles for serotype in recombination. A larger capsule is associated with a higher prevalence and longer duration of carriage, therefore offering more opportunities to undergo recombination, but conversely, experimental and observational data suggest that the absence of capsule increases transformation and recombination (26, 27). As yet, there has been no attempt to quantify the association between serotype properties, such as capsule size and carriage duration, and the amount of recombination that those serotypes experience. Understanding this is important because as noted, recombination is crucial to the response of pneumococci to medical interventions (19), and if serotypes vary in recombination in a predictable fashion, we might be able to predict those that are a higher risk of developing drug resistance or vaccine escape.In the present study, we used whole genomes of 439 pneumococcal isolates collected in Malawi, a low-income country in sub-Saharan Africa with high pneumococcal carriage rates (29) and which contains a wide range of serotypes, including serotypes with very low prevalence or that are absent in industrialized countries such as serotype 1 and 5 (18). We quantified the extent to which recombination has contributed to the history of serotype-specific lineages and combined this with data on carriage duration, invasive potential, carriage prevalence, and capsule size to quantitatively investigate the association. Our results demonstrate significant associations between the recombination rate and frequency with the carriage duration, invasive potential, and size of the outer cell wall polysaccharide capsule in pneumococcal isolates.RESULTSWe sequenced the genomes of 439 Streptococcus pneumoniae clinical isolates collected from the Queen Elizabeth Central Hospital (QECH) in Blantyre, Malawi (southeastern Africa), between 2002 and 2010. The set of samples comprised 364 invasive isolates and 75 carriage isolates and representatives of 48 distinct serotypes. The most common serotypes were serotype 1 with 83 isolates (52 isolates from blood, 28 isolates from cerebral spinal fluid [CSF], and 3 isolates from carriers), serogroup 6 with 46 isolates (25 isolates from blood, 11 isolates from CSF, and 10 isolates from carriers), serotype 5 with 25 isolates (17 isolates from blood, 6 isolates from CSF, and 2 isolates from carriers), and other serotypes (285 isolates; 126 isolates from blood, 99 isolates from CSF, and 60 isolates from carriers) (see Table S1 in the supplemental material). Whole-genome sequencing was done using the HiSeq platform (Illumina, CA, USA) as described in Materials and Methods. The raw sequence reads had an average length of 72.23 bp (95% confidence interval [95% CI], 70.58 to 73.87 bp), average read quality of 33.48 (95% CI, 33.22 to 33.75), and average insert size of 263.30 kb (95% CI, 263.3 to 270.10 kb) per genome (Fig. S1). Consensus de novo sequence assemblies were generated (30, 31) with an average of 66.43 contigs (95% CI, 63.65 to 69.22) and an average chromosome size of 2.09 Mb (95% CI, 2.09 to 2.10 Mb) and average contig size of 37.31 kb (95% CI, 35.91 to 38.72 kb) per genome. The phylogeny of all the isolates was constructed using a 0.79-Mb multiple-sequence alignment of 852 concatenated core genes containing 51,389 single nucleotide polymorphisms (SNPs) (32). For the analysis of serotype-specific evolution of lineages identified as described below, SNPs were called against reference genomes (Table S2).Population structure.Analysis of the population using BAPS (33, 34) identified 14 primary sequence clusters (SCs) (Fig. 1). Of these 14 SCs, SC1 to SC13 were monophyletic clades. The remaining SC, SC14, contained sequences that could not be assigned to a monophyletic SC, and the streptococci in this polyphyletic SC are a diverse group of rare lineages, the grouping of which does not necessarily reflect common evolutionary history, and hence, we undertake no further analysis of this group. Most SCs contained isolates with a single serotype and multiple closely related sequence types (STs) as determined by multilocus sequence typing (MLST) (Fig. 1; see Table S1 in the supplemental material). In some cases, SCs contained multiple serotypes separated by long branches. A notable example is SC9, which consisted of clearly distinct subclades of serotype 6A and 35B (Fig. 1), which we split into two subclades 6A-SC9 and 35B-SC9 for the serotype-specific analyses. The subclade with serotype 35B in SC9 contained additional serotypes with the same ST that were extremely closely related at the whole-genome level (Fig. 1), which suggests that this subclade was a single lineage of serotype 35B that recently acquired other capsule types.FIG 1 Population structure of S. pneumoniae strains from carriers and patients with invasive disease. The phylogeny was constructed using a 0.79-Mb multiple-sequence alignment with 51,389 SNPs from 852 universally conserved (core) genes present in single copies. The 14 sequence clusters (SCs) identified by BAPS are labeled on the edge of the phylogeny. Clades corresponding to the 13 different monophyletic SCs (SC1 to SC13) are shown in different colors, while clades for all the strains in the polyphyletic clade (SC14) comprising of the “unclustered” sequences are not colored. The outermost ring around the phylogeny shows whether a serotype is a vaccine type (VT) included in the PCV13 pneumococcal vaccine formulation (black) or nonvaccine type (NVT) not included (white). Phylogeny was rooted using an outgroup method on a branch containing nontypeable (NT) pneumococci, which are genetically distinct from all pneumococcal strains.Variation in recombination among serotypes.Pneumococci are known to undergo homologous recombination, and all SCs containing capsule switch variants reflect a history of homologous recombination events. To determine the specific recombination events responsible for these changes, as well as those that could have occurred elsewhere in the genome, we used Gubbins (35), which identifies regions with an atypically high density of SNPs using a sliding window approach (Fig. 2a and b; see Fig. S3a to S3o in the supplemental material). We distinguish between two separate recombination parameters. The relative recombination rate (μr/m) is the ratio of the numbers of SNPs introduced by recombination to the number introduced by mutation. It hence measures the total amount of diversity accumulated by a lineage. In contrast, the recombination frequency (μre) is defined as the average number of recombination events inferred per isolate, regardless of how much variation they have introduced. Hence, a large amount of recombination with closely related strains could lead to a high μre but a low μr/m. Furthermore, higher presence of shared ancestral recombination events acquired through clonal descent than strain-specific events introduced independently may lead to a high μre. The averages of both these parameters were calculated for each SC. While SC1 and SC2 both contained isolates with a single serotype (serotype 1 and 5, respectively), in the case of SCs containing more than one serotype, the statistics for each were calculated independently (Table 1). Separate estimates were obtained for the two distinct 6A clades in SC3 and SC9, which we refer to as 6A-SC3 and 6A-SC9, respectively, on the basis of the predominant STs in each.FIG 2 Recombination inferred by Gubbins. Recombination events mapped onto internal nodes (red blocks) and terminal branches (blue blocks) are illustrated for a typical low-recombination serotype (serotype 1 [SC1]) (a) and a high-recombination serotype (serotype 6A [SC3]) (b). (i) Maximum likelihood phylogenetic trees of the different serotypes. Branches and tips of the phylogeny are colored according to the MLST sequence type (ST). (ii) Schematic representation of the reference S. pneumoniae genome showing all the genetic annotations and locations of some well-known genes. (iii) Matrix showing the tracks representing each genome. The locations and distribution of regions, which have acquired exogenous DNA through recombination, are colored depending on the number of strains that contain them. Recombination events in internal branches (red) were present in multiple isolates and were shared through clonal descent rather than independent acquisitions, while those in the terminal branches (blue) were isolate specific and represent independent recent acquisitions. The recombination rate (μr/m), i.e., mean number of the inferred distinct recombination events per isolate (each shared ancestral recombination event that occurred once and spread in the clone via clonal descent was counted once). The recombination frequency (μre), i.e., the mean number of SNPs introduced through recombination to those introduced through mutation, is shown. A high presence of recent rather than shared recombination events implies high μre.TABLE 1 Summary statistics for the genetic recombination events identified in each serotypeSerotypeSequence cluster (SC)No. of isolates (n)Mean no. of recombination events/isolate (μre)Mean recombination to mutation (μr/m)Recombination size (bp)Pherotypea
Mean95% CIbMean95% CICSP1CSP2Other5SC12522.820.29−5.357,6423,877−11,407100001SC2830.732.610.74−4.488,7275,753−11,70298.7801.226ASC3118.9123.451.54−45.4411,5798,872−14,2860100019ASC4125.5811.264.64−17.899,1436,308−11,9790010016FSC590.110.440−1.3936,6160 001003SC6108.19.450−23.538,0206,201−9,8390010010BSC7142.074.150−8.709,8554,805−14,90576.4717.655.8823FSC81912.2612.537.33−17.728,0946,737−9,45194.7405.2635BSC91320.0817.025.31−28.749,0227,987−10,056100006ASC91211.9212.854.55−21.157,1155,697−8,5321000012FSC10195.168.172.38−13.977,7096,016−9,40194.7405.267A/FSC11206.48.743.09−14.386,9245,547−8,302950514SC1215413.934.70−23.166,7454,926−8,5631000018B/CSC1386.673.490−11.159,4696,215−12,72483.7816.2204SC1360.751.670−3.679,2080−20,17538.7532.5628.29aProportion of the pneumococcal pherotypes or competence-stimulating peptide (CSP) encoded by the comC gene and its variants in each SC. Detailed information about the CSPs is shown in Fig. S2c and S2d in the supplemental material.b95% CI, 95% confidence interval.Estimates of μr/m ranged from 0.44 to 23.45 with a mean of 8.89 (95% CI, 5.243 to 12.43) (Table 1), and the only estimate less than 1 was found in the 16F lineage SC5. The highest values of μr/m were 12.53 for serotypes 23F-SC8, 23.45 for 6A-SC3, 12.85 for 6A-SC9, and 17.02 for 35B-SC9, while serotypes 1-SC2 (μr/m = 2.61), 16F-SC5 (μr/m = 0.44), 4-SC13 (μr/m = 1.67), and 5-SC1 (μr/m = 2.82) showed the lowest values (Table 1). The distribution of μr/m values deviated from a normal distribution (P < 0.0001) on the basis of three normality tests (see Materials and Methods). In serotypes of the same SC, such as 35B, 6A, and 18B/C in SC4, SC9, and SC13, respectively, the distribution of r/m ratios per branch were not different with P values of 0.7431 and 0.9350, respectively. The average estimated μre was 6.316 (95% CI, 3.316 to 9.316) and ranged from 0.11 to 20.08. The highest μre estimates were 12.26 for 23F-SC8, 20.08 for 35B-SC9, 8.91 for 6A-SC3, and 11.92 for 6A-SC9, while serotypes 1-SC2 (μre = 0.73), 16F-SC5 (μre = 0.11), 4-SC13 (μre = 0.75), and 5-SC1 (μre = 2) showed the lowest μre values (Table 1). Significant differences (P = 0.0005) in the average numbers of recombination events in each serotype were also observed by using the same test (see Fig. S2b in the supplemental material). The Kruskal-Wallis test was used because of significant deviations of the distributions of both μr/m per branch in each serotype and the average number of recombination events for each serotype from the normal distribution. Interestingly, both distinct lineages with 6A capsule type consistently showed high μr/m (23.45 in SC3 and 12.85 in SC9) and μre (8.91 in SC3 and 11.92 in SC9) despite having a dissimilar genetic backbone. In addition, μr/m positively correlated with μre (R2 = 0.5168; P = 0.0025) (Fig. 3). The pneumococcal pherotypes or competence-stimulating peptides (CSP) encoded by the comC gene were diverse and variably distributed in the SCs but showed no obvious association with the recombination parameters (Table 1).FIG 3 Relationship between rate and frequency of genetic recombination estimates in different S. pneumoniae serotypes. The solid red circles show the observed values, while the solid white circles show the estimated values by the univariate linear regression model. The dashed lines connecting the red circles to the white circles represent the residual, the difference between the observed and predicted values by the model. The gray band surrounding the regression line shows the 95% confidence interval (95% CI) for prediction of each data point in the horizontal axis. The two serotype 6A symbols in the figure originated from distinct clades, SC3 and SC9, as shown in Fig. 1 and Table 1.Serotype-specific carriage duration, prevalence, and invasive potential correlate with recombination.Different pneumococcal lineages vary in the extent of recombination that has influenced their genomes (20, 26). Serotypes are known to vary in their duration of carriage and propensity to cause invasive disease per exposure. A longer duration of carriage could offer more opportunities for cocolonization with different strains and possible recombination. We hence tested the association for each serotype between μr/m and μre and previously published serotype-specific estimates of carriage duration (36), carriage prevalence data from this study, capsule size (21), and invasive potential (defined as the odds ratio for causing invasive disease to carriage [ICOR] [22, 37]) using univariate linear regression. We observed a positive association between both μr/m (R2 = 0.4664; P = 0.0071) and μre (R2 = 0.3393; P = 0.0367) with carriage duration (Fig. 4a and Table 2; see Fig. S4a in the supplemental material), but only μr/m showed a significant association (R2 = 0.4173; P = 0.0093) with carriage prevalence of the serotypes (Table 2; see Fig. S6a and S6b in the supplemental material). When modeling the relationship between carriage duration and μre, serotype 35B-SC9 was excluded as an outlier, but its inclusion decreased the P value (R2 = 0.2406; P = 0.0749). Both μr/m (R2 = 0.4549; P = 0.0228) and μre (R2 = 0.7155; P = 0.0010) showed negative association with the invasive potential of the studied serotypes (Fig. 4b and Table 2; Fig. S4b and S4c). In multivariable regression, the association of carriage duration, prevalence, and invasive potential disappeared (Table 2). Repeated univariate and multivariate analyses using only those serotypes for which all response variables were available found similar results (Table S3). A multivariate regression to account for this removed the significant association with serotype invasive potential, while the association with carriage duration and capsule size remained (Table 2).FIG 4 Statistical relationships between recombination rate (μr/m) and frequency (μre) with carriage duration and serotype invasive potential (ICOR). Univariate regression between μr/m and carriage duration (a) and μre and ICOR (invasive to carriage odds ratio) (b). The solid red circles show the observed values, while the solid white circles show the estimated values by the univariate linear regression model. The dashed lines connecting the red circles to the white circles represent the residual, the difference between the observed and predicted values by the model. The gray band surrounding the regression line shows the 95% confidence interval (CI) for prediction of each data point in the horizontal axis. The two serotype 6A symbols in the figure originated from distinct clades, SC3 and SC9, as shown in Fig. 1 and Table 1.TABLE 2 Summarized estimates from the univariate, multivariable, and multivariate multiple regressionResponse variablePredictor variableUnivariateMultivariableMultivariateEstimateSEP valueaEstimateSEP valueEstimatebApprox. F
P valueaRecombination rate (μr/m)Capsule size0.02050.01860.3000−0.01850.02420.4870.904414.1830.0296*Carriage duration0.33360.10300.0071*0.21490.33010.5500.908514.8840.0277*Invasive potential−0.79300.28940.0228*−0.88470.79270.3270.40321.01340.4611Carriage prevalence1.00300.32870.0093*0.10810.94680.9150.54571.80210.3062Recombination frequency (μre)Capsule size0.02630.00780.0095*0.00570.00810.5160.904414.1830.0296*Carriage duration0.17070.07180.0367*0.25950.11010.0780.908514.8840.0277*Invasive potential−0.61420.12910.0010*−0.41270.26450.1940.40321.01340.4612Carriage prevalence0.54430.32600.1189−0.49670.31590.1910.54571.80210.3062aThe asterisks in the P value columns indicate significant values.bAll test statistics, i.e., Pillai-Bartlett trace criterion, Wilks’ lambda, Hotelling-Lawley trace, and Roy’s greatest root, were used to summarize the results from the multivariate analysis of variance (MANOVA) of the multivariate multiple regression model. All the tests showed identical P values. For illustration, only estimates using the Pillai-Bartlett trace are shown.Recombination correlates with larger polysaccharide capsules.Capsule size has been associated with the factors described above—the serotype-specific duration of carriage and prevalence. The polysaccharide capsule size data for different pneumococcal serotypes, determined by the fluorescein isothiocyanate-labeled dextran (FITC-dextran) assay, were obtained from Weinberger and colleagues (21). By using univariate regression, we observed positive association between the capsule size and μre (R2 = 0.5892; P = 0.0095), but not μr/m (R2 = 0.1186; P = 0.300) (Fig. 5; Table 2). Excluding serotype 35-SC9, which was deemed to be an outlier observation, the linear association between capsule size and μr/m increased (R2 = 0.3519; P = 0.0544). This was robust to multivariate analysis (P = 0.0296). We also did a multivariable regression to account for collinearity among capsule size, duration of carriage, carriage prevalence, and invasive potential and identify the independent contribution of each, which showed no association between μr/m and μre (Table 2; see Table S3 in the supplemental material).FIG 5 Relationship between frequency of recombination (μre) and polysaccharide capsule size in pneumococcal isolates. Serotype 35-SC9 was excluded from the regression, as it was deemed an outlier observation. The solid red circles show the observed values, while the solid white circles show the estimated values by the univariate linear regression model. The dashed lines connecting the red circles to the white circles representthe residual, the difference between the observed and predicted values by the model. The gray band surrounding the regression line shows the 95% confidence interval (CI) for prediction of each data point in the horizontal axis. The two serotype 6A symbols in the figure originated from distinct clades, SC3 and SC9, as shown in Fig. 1 and Table 1.DISCUSSIONRecombination is a fundamental process in the evolution and adaptation of many pathogens. It is capable of shuffling existing variation into more fit combinations, and in the case of pneumococcus, it has been an important force in its response to medical interventions like vaccines and antibiotics (19). Despite the obvious benefits to the pathogen that recombination offers in this regard, the rates with which it occurs vary greatly among pneumococcal lineages and serotypes (26) for reasons that are obscure. Understanding the factors that contribute to this variation is important as we attempt to limit the potential of pneumococcal evolution to erode the benefits of clinical innovation.Previous studies have commented on the low presence of recombination in certain serotypes, such as the “highly invasive” serotypes 1 and 7F (20, 38, 39). A systematic comparison of capsular types is challenging because some serotypes are comparatively rare or absent from samples taken in industrialized countries, and while it might be possible to combine studies from different sites, it would be important to account for any bias arising from variation in recombination in different regions—for instance, previous work with MLST data found more evidence of recombination in samples from African carriers and speculated that this might be the consequence of a higher carriage rate (40).We have used a set of pneumococcal samples that includes an unusually wide range of serotypes with diverse properties, and we found a clear association between serotype and μr/m and μre. Because of lack of data in resource-poor settings, we used previously published data on the properties of serotypes, such as capsule size, carriage duration, prevalence, and invasive capacity (all of which vary relatively little among sample sites); we found that all were associated with one recombination parameter or the other recombination parameter or both. In univariate analysis, a negative association with invasive capacity mirrored a positive association of μr/m with carriage duration and prevalence—likely reflecting the fact that the majority of highly invasive serotypes are also infrequently carried. However, it should be noted that we did not have specific invasiveness data from pneumococci infecting the Malawi population, and there are few estimates for this quantity from resource-poor settings, so this is an important focus for future work.We have also found more frequent recombination in serotypes with larger outer cell wall polysaccharide capsules, which shows that despite the potential physical barrier offered by the capsule, the presence of a larger capsule itself does not simply reduce recombination. While this could be considered surprising, it is important to note that serotypes with enlarged capsules are typically carried for longer duration and resist neutrophil phagocytosis (21) and complement-mediated immune responses (41, 42) than serotypes with smaller capsules. In addition, this may be because some of the serotypes with larger capsules have lower metabolic costs (21) associated with capsule expression and may have a higher ability to form biofilms which may prevent elimination (27). This suggests that longer duration of carriage and lower clearance rates by the host’s immune system in serotypes with enlarged capsules leads to sustained exposure to cocolonizers (potential donors of exogenous DNA), which has been previously described to facilitate recombination (40). This in turn results in higher recombination than in serotypes with thinner capsules, which are presented with limited opportunities for recombination due to rapid clearance. The intrinsic recombination rates in vitro in the absence of capsule-induced immune clearance may be different than observed in the natural population where various factors such as strain competition and immune responses can directly affect the recombination rate. Hence, the results of in vitro studies may be misleading when applied to the nasopharynx. In addition, there was no obvious association between the pneumococcal pherotypes such as CSP1, CSP2, and other variants encoded by the comC gene and its allelic variants, which suggests that the observed differences in recombination may not be due to the differences in competence. The evidence presented here certainly strongly suggests that capsule does not prevent recombination in any straightforward fashion. However, while capsule size was associated with a higher frequency of recombination (μre), the association with μr/m was not significant in the univariate regression. Nevertheless, the effect of capsule size and carriage duration in a multivariate analysis accounts for the high correlation between the two recombination parameters, i.e., μr/m and μre, indicating that these two factors are key contributors to the observed variation. The capsule data used in this study were obtained in vitro, and the expressed capsule sizes during carriage and transmission remain unclear and should be investigated in follow-up studies.The properties of the capsule that we have studied here are also highly correlated (Pearson correlation coefficients are shown in Fig. S7 in the supplemental material), and accounting for this removed the association (Table 2). As a result, we have not been able to precisely define the underlying causes. A plausible mechanism is that serotypes with larger capsules are more exposed to other cocolonizing strains because they are carried for a longer period and at a higher prevalence. One way to test this would be to use serotype switches as natural experiments and compare μr/m and μre between switched lineages where the two serotypes involved vary in their size. However, this is limited by the fact that serotype switches identified in Fig. 1 and other studies typically involve serotypes that are effective colonizers such as serogroups 6, 15, 18, 19, and 23 (15). While this means we cannot use the present data set to do the above analysis, we should note that this observation is consistent with the suggestion that more commonly carried types are more likely to undergo recombination with each other.In conclusion, we have shown that in a set of samples from Malawi containing many different serotypes, there is a linear relationship between μre and μr/m and properties of the capsule. In general, there is evidence that recombination increases with carriage duration and polysaccharide capsule size. We have not been able to pick out the causal roots of this relationship due to the complexity of the relationships between other important properties also associated with capsule, but the observation shows conclusively that all serotypes are not alike and that serotypes differ in how recombination has affected their evolution history, which has been associated with many different responses from drug resistance to vaccine escape. The serotypes observed in settings with a high burden of carriage and disease, such as sub-Saharan Africa, are often different from those found in industrialized nations, and there has been concern that lineages with serotypes targeted by vaccine such as serotypes 1, 5, and 7F that are relatively invasive, with a short duration of carriage and small capsules, might survive a serotype switch (21, 22, 36). The results of this study suggest that such concerns may be misplaced.MATERIALS AND METHODSIsolate collection and genome sequencing.We collected 439 S. pneumoniae isolates from the Malawi-Liverpool-Wellcome Trust Clinical Research Programme (MLW) pneumococcal isolate archive in Blantyre, Malawi. The isolates were routinely collected from patients at the Queen Elizabeth Central Hospital, a 1,250-bed referral hospital and the largest hospital located in the Blantyre District in southern Malawi from 2002 to 2011. All isolates were collected before the introduction of the 13-valent pneumococcal conjugate vaccine (PCV13) in routine childhood immunization programs in Malawi in 2011. Isolates were collected from carriers (n = 75) and patients with invasive disease (n = 364) from 2002 to 2011 (see Table S1 in the supplemental material). Isolates were cultured in Todd-Hewitt broth, and genomic DNA was extracted using the Promega Wizard DNA genomic DNA purification kit (Promega, USA). DNA sequencing was done using the Illumina genome analyzer II (Illumina, CA, USA) platform at the Wellcome Trust Sanger Institute. A robust bespoke pipeline (43) that was composed of the Velvet v1.2.09 De Bruijn graph-based DNA sequence assembler (30), Velvet Optimizer (31), SSPACE Basic v2.0 contig-scaffolding software (44), GapFiller (45), and SMALT v0.7.4 short-read aligner (http://sourceforge.net/projects/smalt/) was used to generate sequence assemblies from the paired-end Illumina sequence reads for the study isolates. Serotypes were identified from the whole-genome sequence data using an in silico approach implemented with Perl and BioPerl (46) as described by Croucher et al. (47).Population structure and phylogenetic analysis.Phylogenetic trees were constructed using a 0.79-Mb sequence alignment of 51,389 single nucleotide polymorphisms (SNPs) in 852 concatenated core genes present in single copies (no paralogs). The conserved (core) and nonconserved (accessory) genes were identified using the Roary bacterial genome analysis pipeline (32). The population structure of the isolates was inferred using the hierBAPS module implemented in BAPS v6.0 (33, 34). Maximum likelihood phylogenies of the concatenated core gene alignments for all the study isolates and serotype-specific isolates in the inferred sequence clusters (SCs) from BAPS were constructed using RAxML v7.0.4 (48). We used a generalized time reversible (GTR) (49) model with gamma (Υ) heterogeneity across nucleotide sites and 100 bootstrap replicates. The tree was rooted with an outgroup clade of nontypeable (NT) isolates, which form a distinct clade from other pneumococci (50). Phylogenetic trees and associated metadata, namely, sequence types (STs), serotypes, and inclusion of the sample’s serotype in the PCV13 vaccine formulation, were overlaid on the phylogenies and visualized in iToL (51) and BioPython scripts (52).Recombination detection and phylogeny construction.Phylogenies for specific serotypes were generated for downstream phylogenetic analysis through mapping paired-end short sequence reads against different published reference whole-genome sequences (see Table S2 in the supplemental material) for different pneumococcal serotypes using SMALT v0.7.4. The output binary alignment map (BAM) files were realigned using the Genome Analysis Tool Kit (GATK) v3.3.0 (53). The locations of recombination events in each chromosome in the alignments were detected using Gubbins v1.1.1 (35). The serotype-specific phylogenies with recombination removed were constructed using RAxML v7.0.4 (48) as above and rooted in the middle of the branch, separating the two most divergent isolates. Visualization of the phylogenies was done as described in “Population structure and phylogenetic analysis” above.Carriage duration, serotype invasive potential, and capsule size.We used previously published data from Kilifi in Kenya, data on carriage duration of pneumococcal serotypes in the human nasopharynx from Abdullahi and colleagues (36) and data on serotype invasive potential (odds ratio for IPD to carriage [ICOR]) from Brueggemann and colleagues (22). Additional data for serotype 5 invasive potential were obtained from Smith and colleagues (37). The polysaccharide capsule sizes for the serotypes were obtained from Weinberger and colleagues (21) where the degree of encapsulation was determined by measuring the zone of exclusion using the fluorescein isothiocyanate-labeled dextran (FITC-dextran) assay. The specific methods used to obtain the data sets are described in the respective publications. We excluded serotypes with estimated ICOR values of 0 (22) to prevent overfitting as a consequence of the smaller number of isolates used from invasive disease sample collection. Genomic data on the relative recombination rate and frequency of pneumococcal lineages and prevalence of serotypes in carriage were obtained from Malawi (this study).Statistical analysis.We used univariate, multivariable, and multivariate multiple regression to examine the association between known properties of serotypes and the relative recombination rate (μr/m) and frequency (μre). The univariate regression model were separately formulated as Yi = β0 + β1xi + εi for each dependent variable Yi(μr/m and μre) against each independent variable xi where i designates the polysaccharide capsule size, carriage duration, serotype invasive potential, and carriage prevalence. The intercept and regression coefficient for the independent variables are represented by β0 and β1, respectively, while εi represents the residuals. We formulated the multivariable model similarly but included multiple independent variables to account for collinearity between them as Yi = β0 + β1xi + β2xi + β3xi + β4xi + εi, where β0 represents the intercept and the terms β1, β2, β3, and β4 are the regression coefficients for each independent variable (xi) where i is polysaccharide capsule size, carriage duration, serotype invasive potential and carriage prevalence. We also modeled the relationship between both dependent variables (μr/m and μre) and the independent variables in the multivariable regression this time with a multivariate multipl e regression model. The multivariate model was formulated as Ynxm = Xnx(p+1)β(p+1)xm + εnxm, where n designates each serotype and lineage combination as follows 1-SC2, 5-SC1, 6A-SC3, …, and 4-SC13 (Table 1). The m in the matrix Ynxm takes values 1 and 2 that designate dependent variables μr/m and μre, respectively. The model matrix is the term Xnx(p+1) and consists of p+1 regressors where p represents polysaccharide capsule size, carriage duration, serotype invasive potential, and carriage prevalence of serotypes. The first column in the model matrix Xnx(p+1) consists of 1s, which represents the regression constants. The matrix of regression coefficients for all the independent variables in Xnx(p+1) is represented by parameter β. The residual term in the model is represented by the matrix εnxm. For comparison with the multivariable and multivariate analyses, we repeated the univariate analyses with the same sample size.We calculated the pairwise Pearson correlation coefficients for independent and response variables and plotted them using the ggplot2 package. The regression analysis, normality tests, namely, D’Agostino and Pearson, Shapiro-Wilk, and Kolmogorov-Smirnov tests, Kruskal-Wallis and analysis of variance (ANOVA) tests and summary statistics, such as means and ranges, were calculated using R v3.1.2 (R Core Team, 2014) and GraphPad Prism v6.0 (GraphPad Software, Inc., CA, USA).SUPPLEMENTAL MATERIALFigure S1 Assembly statistics of the pneumococcal study isolates. (a) Length of the assembled genomes; (b) mean number of contigs per genome; (c) mean contig length; (d) N50 assembly statistic; (e) mean number of contigs in N50; (f) mean read length; (g) percentage of total mapped sequence reads; (h) percentage of the total bases mapped. Download Figure S1, EPS file, 2.1 MB
Figure S2 Characteristics of the genetic recombination events identified in the pneumococcal isolates. (a) Distribution of recombination rates on branches of the serotype-specific phylogenies; (b) sizes of recombination blocks identified in the genomes of each pneumococcal serotype; (c) distribution of the competence-stimulating peptide (CSP) variants across the sequence clusters (SCs); (d) protein multiple-sequence alignment of the CSPs encoded by the comC gene. Other variants of the comC gene in additional to the CSP1 and CSP2 were arbitrarily named CPS3 to CSP10. Download Figure S2, TIF file, 1 MB
Figure S3 Recombination inferred by Gubbins. Pneumococcal serotypes 5-SC1 (a), 2-SC1 (b), 6A-SC3 (c), 19A-SC4 (d), 16F-SC5 (e), 3-SC6 (f), 10B-SC7 (g), 23F-SC8 (h), 6A-SC9 (i), 35B-SC9 (j), 12F-SC10 (k), 7AF-SC11 (l), 14-SC12 (m), 4-SC13 (n), and 18 A/F-SC13 (o). (i) Maximum likelihood phylogenetic trees of the different serotypes. Branches and tips of the phylogeny are colored according to the MLST sequence type (ST). (ii) Schematic representation of the reference S. pneumoniae genome showing all the genetic annotations and locations of some well-known genes. (iii) Matrix showing the tracks representing each genome. The locations and distribution of regions, which have acquired exogenous DNA through recombination, are colored depending on the number of strains that contain them. Recombination events in internal branches (red) were present in multiple isolates and were shared through clonal descent, rather than independent acquisitions, while those in the terminal branches (blue) were isolate specific and represent independent recent acquisitions. The recombination rate (μr/m), i.e., mean number of the inferred distinct recombination events per isolate (each shared ancestral recombination event that occurred once and spread in the clone via clonal descent was counted once) and the recombination frequency (μre), i.e., the mean number of SNPs introduced through recombination to those introduced through mutation are shown. A high presence of recent recombination events rather than shared recombination events implies a high μre. Download Figure S3, TIF file, 1.6 MB
Figure S4 Relationships between recombination rate (μr/m) and frequency (μre) with carriage duration and invasive potential in pneumococcal isolates. Univariate regression between μre and carriage duration (a) and μr/m and ICOR (invasive to carriage odds ratio) (b). The solid red circles show the observed values, while the solid white circles show the estimated values by the univariate linear regression model. The dashed lines connecting the red circles to the white circles represent the residual, the difference between the observed and predicted values by the model. The gray band surrounding the regression line shows the 95% confidence interval (CI) for prediction of each data point in the horizontal axis. The two serotype 6A symbols in the figure originated from distinct clades, SC3 and SC9, as shown in Fig. 1 and Table 1. Download Figure S4, EPS file, 1.9 MB
Figure S5 Recombination dynamics and capsule polysaccharide capsule sizes in pneumococcal serotypes. (a) Recombination frequency (μre) and capsule size. Serotype 35-SC9 (light pink) was excluded, as it was deemed an outlier. (b) Recombination rate (μr/m) and capsule sizes. The solid red circles show the observed values, while the solid white circles show the estimated values by the univariate linear regression model. The dashed lines connecting the red circles to the white circles represent the residual, the difference between the observed and predicted values by the model. The gray band surrounding the regression line shows the 95% confidence interval (CI) for prediction of each data point in the horizontal axis. The two serotype 6A symbols in the figure originated from distinct clades, SC3 and SC9, as shown in Fig. 1 and Table 1. Download Figure S5, EPS file, 1 MB
Figure S6 Relationship between frequency (μre) of recombination and carriage prevalence in pneumococcal isolates. (a) μr/m and carriage prevalence; (b) μre and carriage prevalence. The solid red circles show the observed values, while the solid white circles show the estimated values by the univariate linear regression model. The dashed lines connecting the red circles to the white circles represent the residual, the difference between the observed and predicted values by the model. The gray band surrounding the regression line shows the 95% confidence interval (CI) for prediction of each data point in the horizontal axis. The two serotype 6A symbols in the figure originated from distinct clades, SC3 and SC9, as shown in Fig. 1 and Table 1. Download Figure S6, EPS file, 1.9 MB
Figure S7 Composite Pearson correlation matrix derived for the variables. The independent (predictor) variables, namely, carriage duration and capsule size, are indicated by two asterisks, while the different dependent (response) variables, namely, recombination rate (μr/m), recombination frequency (μre), serotype invasive potential (ICOR), and prevalence in carriage and invasive disease, do not have an asterisk. Download Figure S7, EPS file, 1.2 MB
Table S1 Streptococcus pneumoniae isolates used in this study.Table S1, DOCX file, 0.2 MB
Table S2 Reference pneumococcal whole-genome sequences used in this study.Table S2, DOCX file, 0.1 MB
Table S3 Summarized estimates from the univariate and multivariable regression. Only serotypes with all the different independent variables present, i.e., capsule size, carriage duration, invasive potential, and carriage prevalence were used.Table S3, DOCX file, 0.1 MB