Pneumococcal genetic variability influences age-dependent bacterial carriage

The pneumococcal conjugate vaccine (PCV) primarily reduces disease burden in adults through a reduction in carriage prevalence of invasive serotypes in children. Current vaccine formulations are the same for both adults and children, but tailoring these formulations to age category could optimize vaccine efficacy. Identification of specific pneumococcal genetic factors associated with carriage in younger or older age groups may suggest alternative formulations and contribute to a better mechanistic understanding of immunity. Here, we used whole genome sequencing to dissect pneumococcal variation associated with age. We performed genome sequencing in a large carriage cohort, and conducted a meta-analysis with an existing carriage study. We compiled a dictionary of pathogen genetic variation including serotype, sequence cluster, sequence elements, SNPs, burden combined rare variants, and clusters of orthologous genes (COGs) for each cohort – all of which used in a genome-wide association with host age. Age-dependent colonization had some heritability, though this varied between cohorts (h2 = 0.10, 0.00 – 0.69 95% CI in the first; h2 = 0.46, 0.33 – 0.60 95% CI in the second cohort). We found that serotypes and genetic background (strain) explained most of the heritability in each cohort (h2serotype = 0.06 and h2GPSC = 0.04 in the first; h2 serotype= 0.20 and h2 GPSC = 0.23 in the second cohort). When looking for genetic effects independent of genetic background and serotype, no individual genetic variants were robustly associated with carriage age. Overall, association with age was highly cohort and strain dependent, supporting proposals for a future vaccination strategy which is primarily targeted using serotypes rather than proteins, and is tailored towards specific pathogen populations.


Introduction
Streptococcus pneumoniae is a common commensal of the human upper respiratory tract and nasopharynx, but can also cause invasive diseases such as pneumonia, sepsis or meningitis.(1) Invasive pneumococcal disease (IPD) has a high mortality, and the overall mortality rate from IPD is higher in extreme age ranges, such as infants and the elderly. (2,3) In the Netherlands, pneumococcal carriage rates are higher in children than in adults, with a prevalence of up to 80% at 2 years of age. (4) Host age is known to affect carriage prevalence and carriage duration of different serotypes (5,6), which is suggested to be driven by differences in immunity. (7) Studies in mice and humans showed evidence for agedependent host-pathogen interactions involving interleukin (IL)-1 response in reaction to the pore-forming pneumoplysin (ply) toxin. (8) IgA secretion is important in clearing S. pneumoniae from host upper respiratory tract mucosa and this secretion more effective in previously exposed individuals, the adults.(9) Bacterial genetics has shown to explain over 60% of the variability in carriage duration, and specifically that presence of a bacteriophage inserted in a mediator of genomic competence was associated with a decreased carriage duration. (10) Pneumococci are highly genetically variable, displaying over 100 diverse capsular serotypes (11), which are a major antigen and the strongest single predictor of carriage prevalence. (12) Pneumococcal conjugate vaccines, targeting up to thirteen capsule serotypes with high burden of invasive disease, cause decreased the rate of nasopharyngeal carriage and invasive disease. (13,14) Besides a direct effect of vaccination with an pneumococcal conjugate vaccine (PCV) on the disease burden in the target population, i.e. young children, it also reduces the disease burden caused by pneumococci with vaccine serotypes in the population not eligible for vaccination through indirect protection from colonization -reducing carriage rates in children reduces overall transmission of the most invasive serotypes. (12,15,16) However, the introduction of PCV has resulted in the replacement of serotypes not covered by the vaccine (17,18), which in some countries reaches levels of invasive disease return towards pre-vaccine levels. (19,20) As not all serotypes can be included in a conjugate vaccine, three perspectives will lead to improved pneumococcal vaccination have been proposed: whole-cell vaccines (21,22), protein vaccines (23), or changing components in the conjugate vaccine in response to the circulating population. (24) Whole-cell vaccination trials are ongoing, but efficacy remains unproven in human populations. (25) Protein vaccines contain antigens which illicit a strong mucosal immune response, with their targets chosen to be common or conserved in the target population, and ideally reducing onward transmission. (26) In their current form, protein vaccines are not thought to be effective on their own, but if administered with serotype conjugates they may help to reduce serotype replacement. Detailed modelling of the dynamics of pneumococcal population genetics has shown that targeting these vaccines towards serotypes prevalent in specific populations would likely be a superior strategy. This work further shows that providing age-specific vaccine design, using complementary adult-administered vaccines (CAVs) is predicted to have the greatest effect on total IPD burden. (24) For a future pneumococcal vaccination strategy based on age, we should understand the differences between infant and adult carriage. Differences between host niches have been found, some with a potential effect on onward transmission. (27)(28)(29) Treating age as a niche in a systematic analysis of genetic variation would therefore be a powerful way to discover useful vaccine targets. Here, we aim to determine how pathogen variation affects colonization. The identification of age-associated genetic variation, could provide further targets for protein vaccination, whereas ruling these out could provide confidence in age-targeted vaccine formulations based on serotype differences alone.
We performed a pathogen genome-wide association study on pneumococci isolated from nasopharyngeal swabs of 4320 infants and adults from the Netherlands and Myanmar. To dissect pneumococcal variation associated with age we compare prevalence of pneumococcal strains and serotypes between infants and adults in different settings, calculate the contribution of pathogen genetic variation towards predilection for host age, and search for genetic regions associated with host age.

Cohort collection
The Dutch cohort consists of carriage samples from individuals obtained from three prospective carriage surveillance studies. (30)(31)(32) In these studies, carriage was assessed by conventional culture of nasopharyngeal or oropharyngeal swabs of vaccinated children (11 and 24 months of age) and their parents in 2009, in 2010/2011, in 2012 and 2013.(30) All children were vaccinated with PCV-7 or PHiD-CV10 according to the Dutch national immunization program at 2, 3, 4 and 11 months of age. Vaccination status of the parents was unknown. Exclusion criteria are described elsewhere. (30,31) Nasopharyngeal swabs were collected from all individuals and oropharyngeal swabs were collected from all adult subjects by trained study personnel using flexible, sterile swabs according to the standard procedures described by the World Health Organization. (33) After sampling, swabs were immediately placed in liquid Amies transport medium and transported to the microbiology laboratory at room temperature and cultured within 12 hours. Pneumococcal isolates were identified using conventional methods, as described previously.

Host age distribution in sequenced carriage cohorts
In the Dutch cohort, children had a median age of 23 months (interquartile range (IQR) 10 -24 months).
Adults had a median age of 35 (IQR 32 -38) years. In the Maela cohort, the median age of children was 13 months (IQR 6 -19 months), and for mothers (women of childbearing age) the exact age was unknown (Supplementary Figure S1). (6,36) In the Dutch cohort, all children were vaccinated with PCV-7 or PHiD-CV10.
None of the members of the Maela cohort had received PCV.

DNA extraction and whole genome sequencing
For the Dutch cohort, DNA extraction was performed with the Gentra Puregene Isolation Kit (Qiagen), and quality control procedures were performed to determine yield and purity. Sequencing was performed using multiplexed libraries on the Illumina HiSeq platform to produce paired end reads of 100 nucleotides in length (Illumina, San Diego, CA, USA). Quality control involved analysis of contamination with Kraken (version 1.1.1)(37), number and length of contigs, GC content and N50 parameter. Sequences for which one or more of these quality control parameters deviated by more than 3 standard deviations from the mean were excluded.
Sequences were assembled using a standard assembly pipeline. (38) Assembly statistics can be found in the Supplementary Table S1. Genome sequences were annotated with PROKKA, version 1.11.(39) For the Maela cohort, DNA extraction, quality control and whole genome sequencing have been described elsewhere. (40) Serotypes were determined from the whole-genome sequence by in-house scripts.(41) Sequence clusters (strains) were defined as Global Pneumococcal Sequence Clusters (GPSC) using PopPUNK (version 2.2.0), using a previously published reference database. (42,43) For 114 and 401 sequences in the Dutch and Maela cohorts respectively, the GPSC couldn't be inferred due to low sequence quality.

Data availability
Fastq sequences of bacterial isolates from the Dutch cohort were deposited in the European Nucleotide Archive (ENA, study and accession numbers in Supplementary Table S2). Sequences of bacterial isolates in the Maela cohort are available at ENA under study numbers ERP000435, ERP000483, ERP000485, ERP000487, ERP000598 and ERP000599 (Supplementary Table S3).

Phylogenetic tree
A core genome for sequences from both cohorts together was generated with Roary (version 3.5.0, default parameters), using a 95% sequence identity threshold.(44) A maximum likelihood phylogeny of singlenucleotide polymorphisms (SNPs) in the core genome of all sequenced isolates from both cohorts together was produced with iqtree (version 1.6.5, including fast stochastic tree search algorithm, GTR+I+G) assuming a general time reversible model of nucleotide substitution with a discrete γ-distributed rate heterogeneity and the allowance of invariable sites. (45) Heritability analysis

Determining bacterial genetic variation -unitigs, SNPs and COGs
Using the whole-genome sequence reads from both cohorts, we called SNPs, small insertions and deletions and SNPs clustered as rare variants (deleterious variants at an allele frequency < 0.01) based on the S.

Genome-wide association study
The association analysis for SNPs, unitigs, rare variants and COGs was run as a linear mixed model in Pyseer (version 1.1.1), with a minimum minor allele frequency of 0.05. (48) To correct for population structure, the model included a kinship matrix as covariates, which was calculated from the midpoint rooted phylogenetic tree. An association analysis not corrected for population structure was run with unitigs as sequence elements using a simple fixed effects model in Pyseer. Rare variants were clustered in their corresponding gene and analyzed in a burden test. Meta-analysis was performed on summary statistics from the Pyseer results files with METAL (version released on August 28 2018, default parameters) for each variant.(51) A threshold for association of the phenotype with meta-analyzed variants was determined using a Bonferroni correction with alpha < 0.05 and the number of independent tests in the Dutch cohort, giving: p < 1.0x10 -7 for unitigs, p < 1.0x10 -6 for SNPs, p < 2.0x10 -5 for COGs and p < 1.0x10 -4 for rare variants. Unitigs were mapped to the S. pneumoniae D39V reference genome with bowtie-2 (version 2.2.3, with equal quality values and length of seed substrings 7 nucleotides). In accordance with the study populations in both cohorts, the phenotype was dichotomized as host age 0 to 24 months versus adult age (Supplementary Figure S1). Manhattan plots were generated in R, version 3.5.1, with the package ggplot2 (version 3.1.0). Presence or absence of pilus genes was detected by nucleotide BLAST (version 2.6.0, default parameters) analysis. Pilus gene presence association to carriage age was calculated with a likelihood-ratio test in Pyseer (version 1.1.1), corrected for population structure by including a kinship matrix as covariates.
The prediction analysis used the elastic-net mode of pyseer. This fitted an elastic net model with a default mixing parameter (0.0069 L1/L2) to the unitigs counted in each cohort, using the strains from PopPUNK as folds to try and reduce overfitting.(50) ROC curves for each cohort were drawn using the linear link output, with the R package pROC (version 1.16.2) using smoothing. To test inter-cohort prediction, the called unitigs from the other cohorts were used as predictors with the model from the opposing cohort.

Results
We first analyzed the observed distribution of serotypes and strains in each of the two cohorts, to assess overall trends of differences between adults and children, and look at the genetic heterogeneity between the two cohorts. Although our cohorts were broadly matched in the primary phenotype, age, large differences between the pathogen population are expected due to different geographies, social backgrounds, and only children in one cohort being vaccinated.
A phylogenetic tree of pooled sequences from both cohorts, with serotype, sequence cluster, age group and cohort for each sequence, revealed clonal discrimination between cohorts ( Figure 2). Combined with the effects shown in Figure 1, this highlighted a key feature of our analysis of these datasets, which was the genetic heterogeneity between the two cohorts. Individually, each dataset clearly has strains and serotypes with strong signals of host age differences, but the overall makeup of each dataset is very different (twelve common serotypes are shared, but only a single common GPSC), and where there are shared serotypes many have different effect directions between the two cohorts.

Host age is heritable and mostly explained by strain and serotype
To quantify the amount of variability in carriage age explained by variability in the genome, we calculated a heritability estimate (h 2 ) for each cohort. For isolates in the Dutch cohort, we did not find strong evidence that genetic variability in bacteria was related to variance in host age (h 2 = 0.10, 0.00 -0.69 95% CI). In the Maela cohort, we found significant evidence that affinity with host age was heritable (h 2 = 0.46 0.33 -0.60 95% CI) and thus genetic variation in this cohort explained variation in carriage age to a greater degree. In both cohorts pan-genomic variation could be used to predict host age to some degree of accuracy (area under the ROC  To further investigate the association of serotype and sequence cluster to carriage age, we determined the proportion of variation in carriage age explained by serotype and sequence cluster alone. Here, we estimated h 2 serotype = 0.06 and h 2 GPSC = 0.04 for the Dutch cohort and h 2 serotype = 0.20 and h 2 GPSC = 0.23 for the Maela cohort, confirming the larger contribution of serotype and sequence cluster to carriage age in sequences from the Maela cohort. We also performed a genome-wide association analysis, but without controlling for population structure. This reveals genetic variants specific to serotype as determinants for carriage age (pvalues < 5.0x10 -8 ) in both cohorts (Supplementary Table S8 and Supplementary Table S9). Among the genetic variants with the lowest p-values were variants in capsule locus genes (Cps) in both cohorts. This further supports a role of strain and serotype in association with host age, but does not distinguish between the two.

Genome-wide association analysis does not find genetic variants independent of strain
Following these observations that serotype and strain do not explain the full heritability, we performed a pathogen genome-wide association analysis to investigate whether we can detect genetic variants irrespective of the genetic background that are associated with carriage in children or adults. Though the cohorts have little genetic overlap in terms of genetic background, we would be well-powered to detect genetic variation independent of background ('locus' associations). (52,53) In the Dutch cohort, none of the unitigs, SNPs, COGs or rare variants surpassed the threshold for multiple testing correction (Supplementary Figure 1). The burden (sum) of rare variants in a gene for tryptophan synthase, trpB, approach the multiple testing threshold, but was not significant. In the Maela cohort unitigs in the ugpA gene surpassed the threshold for statistical  Figure 2). After meta-analysis, there were two hits which surpass the threshold for statistical significance (Figure 4). The first, a nucleotide sequence marked by multiple unitigs of which the lowest has a p-value of 1.2x10 -9 . This sequence does not map to the S. pneumoniae D39V reference sequence, or to any other S. pneumoniae genome sequences in the NCBI BLAST database (Supplementary Table S10). (54) For this reason, it is not visualized on the Manhattan plot, for which unitigs were mapped to the S. pneumoniae D39V reference sequence. Upon inspection of the individual sequences these unitigs are called from, we find them to map to a sequence outside an open reading frame near to a bacterial protein export gene called accessory Sec (aSec).
The nucleotide sequence near the aSec gene does not map to a region containing an open reading frame, and is a single short sequence, therefore a biologically plausible effect is unlikely. The second hit is a burden of rare variants in a gene for tryptophan synthase, trpB, that surpass the threshold for statistical significance at a pvalue of 5.0x10 -5 . The variants are 2 frameshift variants of very low frequency. These result in a predicted dysfunctional trpB gene in 9 out of 1282 (1%) sequences in the Dutch cohort and in 12 out of 3073 (0.4%) sequences in the Maela cohort. This association of the trpB gene is likely to be an artefact of low allele frequency, as we estimate we are only powered to detect variation in at least 5% of isolates.

Pilus gene presence does not determine carriage age independent of genetic background
Finally, we investigated whether pneumococcal isolates containing a pilus gene preferentially colonize children in the Dutch cohort, as has been previously described in the Maela cohort. (9,55) This study analyzed the Maela cohort and found that 934 out of 2557 (37%) isolates in children versus 95 out of 592 (16%) isolates in adults trpB regulatory divIVA cps transporter had pilus genes present. However, this association of pilus gene presence to carriage age was dependent on lineages within the population. (9) In the Dutch cohort, we found no evidence that host age was dependent on pilus gene presence (22 out of 208 (10%) in adults versus 129 out of 1099 (12%) in children). This was the case whether or not the genetic background was adjusted for (p = 0.35, uncorrected for population structure and p = 0.69, corrected for population structure). Based on these findings, we suggest that the previously reported pilus-IgA1 association is not a universal explanation for difference in colonization between hosts of different ages.

Discussion
The age of the host is known to have an important effect on pneumococcal colonization. (12) Observational studies have demonstrated variation in serotype prevalence and carriage duration between infants and adults.
Mechanistic studies in mice and humans have shown examples of differing immune responses depending both on host factors and pathogen factors. Findings from these studies include the observation that capsular polysaccharides (determinants of serotype) inhibit phagocytic clearance in animal models of upper respiratory tract colonization.(56) A pneumolysin-induced IL-1 response determined colonization persistence in an agedependent manner (8); and pilus expressing strains were found to preferentially colonize children, because of immune exclusion via secretory IgA in non-naïve hosts. (9) Building upon these observations, we sought to investigate and quantify the contribution of pathogen genetic variation to carriage in infant versus adult hosts, using a top-down approach. Through whole genome sequencing and application of statistical genetic methods to two large S. pneumoniae carriage cohorts, we show evidence that bacterial genetic variability influences predilection for host age, though this appears to be highly variable between populations. One important difference between our study cohorts was that children from the Dutch cohort were vaccinated, while children from the Maela cohort were not. While our findings demonstrate that vaccinated versus unvaccinated children were colonized with different bacterial serotypes and different sequence clusters, we observed differences in prevalence beyond just the serotypes included in the vaccine. Another difference between the cohorts was that adults from the Dutch cohort were males and females, while adults from the Maela cohort were female only.
Strain, or genetic background, appears to be the main effect, explaining roughly half of the total heritability in each cohort. We did not find any large genetic effects on host age independent of genetic background, either individually or in a meta-analysis of the two cohorts. This is suggestive of a polygenic architecture of many variants with low effect sizes, along with larger effects between strains. Three reasons can contribute to this: the proportion of the heritability which is caused by lineage effects; rare locus effects which could not be detected with the current sample size; and by sampling from a cohort with vaccinated children and unvaccinated adults and comparing with a cohort of unvaccinated children and adults, we had lower power due to the reduced overlap within and between cohorts in pan-genome content. Although differences in vaccination status between cohorts is a plausible explanation for our findings, we were unable to rule out other factors, for example a population-specific host effect, or the broad effects of different socio-economic status between these cohorts.
In previous bacterial GWAS studies of antimicrobial resistance (such as a single gene which causes antibiotic resistance), large monogenic effects have typically been found to have high heritabilities close to one, and the GWAS identify the causal variant precisely. (50,52,57) When applied to virulence and carriage duration phenotypes, heritable effects have also been found, but these only explained some of the variation in the phenotype. These appeared to be caused by weaker polygenic effects, not all of which could be detected using the relatively small cohorts available. (10,58) We found similar results for host-age heritability in these two cohorts. It is also notable that these are all host-based phenotypes without a clear lab-based assay to define them, and we would expect either little or no selection for the phenotype.
We also could not distinguish between genetic background or serotype being the primary effect due to their correlation. We did note a difference in effect size of serotype between the two cohorts, which may make it unlikely to be the single largest effect on host age. This difference in cohorts could be explained by strain/GPSC being the main and consistent effect on host age. As strains are different between cohorts and each serotype appears in multiple strains, combining them in different amounts would create different directions of effect for serotype. We did not replicate the association of piliated genomes in infant hosts in our newly sequenced cohort, further demonstrating important differences between populations.
In summary, we found an effect of pneumococcal genetics on carriage in children versus adult hosts, which varies between cohorts, and is likely primarily driven by strain (lineage) effects rather than large populationwide effects of individual genes. An important corollary of our work is on future pneumococcal vaccine optimization efforts. A promising approach for future vaccination strategies is to target the different age groups. (24) Whether these should consist of the dominant disease-causing serotypes overrepresented in carriage by each age group, or whether there are age-specific pathogen proteins that should be included is an open question. Our study suggests that targeting these age groups using serotype makeup alone would be sufficient, and supports previous observational and modelling studies which advise targeting the serotype makeup in the vaccine at specific populations to maximize their effect.