This study involved subjects who participated in our prospective, longitudinal analysis of NP colonization and pneumococcal infections in young children from June 2006 through June 2019 (Casey et al., 2010; Casey et al., 2013; Kaur et al., 2016; Kaur et al., 2017; Kaur et al., 2020). The focus in this study was pneumococcal serotype 35B for 2006–2019.

Under the study protocol, children were enrolled from community-based pediatric practices, mainly from middle-class, suburban communities in the Rochester, New York area (USA) at 6 months of age at seven time points (6, 9, 12, 15, 18, 24, and 30–36 months of age). All enrolled children had received the full primary series of PCV7 or PCV13 before the enrollment. During these well-child check-ups, NP washes (instilling and withdrawing ~2 ml of saline in each nostril with bulb syringe) were prospectively collected and cultured for detection of pneumococci. The details of the study design have been previously described (Casey and Pichichero, 2004; Pichichero and Casey, 2007; Casey et al., 2010; Kaur et al., 2011; Casey et al., 2013; Kaur et al., 2016; Kaur et al., 2017; Pichichero et al., 2018; Kaur et al., 2020). Pneumococcal AOM was used as the sentinel disease to distinguish commensals from virulent strains. For precision, clinical AOM infections were microbiologically confirmed by culture of tympanocentesis fluid. Children were treated with antibiotics, most often with amoxicillin/clavulanate unless the child was allergic to penicillin, in which case cefdinir was typically used.

Written informed consent was obtained from parents prior to enrollment in the study, and the IRB of the Rochester Regional Health System approved the study.

Standard microbiology processing and identification techniques were used for isolation of bacteria including SPN from NP and MEF samples. Serotypes of SPN were determined using pure cultures of SPN by Quellung reaction using Latex pools and serotype-specific pneumococcal antisera (Serum Staten Institute, Denmark). Antibiotic susceptibility testing was conducted for 69% of the (173/250) serotype 35B clinical isolates obtained. When an NP and MEF isolate at onset of AOM were both serotype 35B, only the MEF isolate was included in antibiotic susceptibility analyses as described (Kaur et al., 2020). The susceptibility of SPN to antibiotics was determined with VITEK-2 AST-GP68 or GP74 susceptibility cards (bioMérieux) in the clinical laboratories of Rochester General Hospital. Benzylpenicillin, amoxicillin, ceftriaxone, cefotaxime, meropenem, ertapenem, ofloxacin, levofloxacin, moxifloxacin, erythromycin, telithromycin, vancomycin, linezolid, tetracycline, chloramphenicol, and trimethoprim-sulfamethoxazole (TMP-SMX) susceptibility was tested by automated modified broth microdilution. SPN strain ATCC49619 was used as a control for each batch of testing. Pneumococci were classified as susceptible, intermediate, or resistant based on 2021 CLSI breakpoints, with oral cutoffs used for penicillin (CLSI, 2021).

Multi-locus sequence typing (MLST) was conducted for 62% (154/250) of the serotype 35B isolates using either standard PCR and Sanger sequencing for seven housekeeping genes (n=82 isolates) or by in silico analyses based on their WGS assemblies (n=72 isolates). Sequence types (STs) of SPN by PCR method were determined as previously reported (Casey et al., 2013; Kaur et al., 2016; Pichichero et al., 2018; Kaur et al., 2020). For extraction from WGS assemblies, “Sequence query” on PubMLST was used (Jolley and Maiden, 2010), and STs associated with 35B were extracted using the goeBURST algorithm from PhyloViZ (https://online.phyloviz.net/index) to estimate the clonality of the bacterial population (Feil et al., 2004; Francisco et al., 2012). Strains that shared greater than five out of seven MLST alleles were classified as belonging to a clonal complex (CC) (Tyrrell, 2011). Due to cost constraints and other technical issues (SPN could not be recovered from saved culture stocks), antibiotic susceptibility (69%) and MLST typing (62%) were confined to a subset of available SPN 35B isolates, selected from our repository to be representative of the time frame and sample types of this project.

WGS was mainly performed by Illumina. Subsequent sequencing was conducted using PacBio (Pacific Biosciences) for some isolates to ensure higher contiguity and out of convenience, due to multiplexing requirements of our sequencers. Barcoded Illumina sequencing libraries (insert sizes ~150–600 bp) were produced using Nextera XT kits according to the manufacturer’s protocols, and pooled libraries were sequenced on a NextSeq500 at 2 × 150 nt, targeting 200-fold genomic coverage per library. Barcoded PacBio SMRTbell libraries (insert sizes ~3–20 kb) were produced according to protocol for sequencing on the PacBio Sequel instrument, targeting 100-fold genomic coverage per library. We performed WGS sequencing on a subset of our isolates ( Technical Appendix 1 ) to identify potential genetic factors within ST558 associated with the emergence of serotype 35B isolates in pediatric AOM cases and to validate a putative 15A-to-35B capsular serotype switch. This included ST558 isolates from pre-2015 (13 isolates) and post-2015 (36 isolates), ST14687 (one isolate), ST156 (seven isolates), ST10493 (eight isolates), and the novel 35B/ST14683 (one isolate) associated with 15A/ST63 (five isolates) collected for our patient cohort over the same time period to act as a basis for comparison.

All sequence data and assemblies have been deposited at NCBI under PRJNA734910. Illumina reads had technical sequences, such as barcodes and adapters, removed by Trimmomatic (v 0.33) (Bolger et al., 2014). Paired ends that overlapped were merged into single reads using COPE (v 1.1.2) (Liu et al., 2012). Merged and unmerged reads had errors corrected using ErrorCorrectReads.pl from allpathslg. Corrected reads were assembled using Spades (v 3.7.0) (Bankevich et al., 2012). Heavily contaminated assemblies were removed by dropping those that were not SPN according to Taxator-tk. This was done before removing bad contigs because, in our experience, small amounts of contamination can be removed by contig filtering, but this is more difficult as the level of contamination increases. Poor quality contigs were removed by removing those in which more than 97% of the contig was a single nucleotide, the total length was less than 256 bases, or the coverage per million reads remaining after Trimmomatic was less than 4. Genomes were removed from further analysis if they had more than 130 contigs remaining.

PacBio reads were assembled using Falcon (smrtlink v 7.0.1) (Chin et al., 2016). Contigs were circularized using Circlator (v 1.5.5) (Hunt et al., 2015), and then the resulting assemblies were polished using Arrow (smrtlink v 7.0.1). Contigs were removed if their normalized coverage was less than 0.06 or their quality score was less than 60. The remaining contigs had errors corrected using the corrected Illumina reads (the reads after running ErrorCorrectReads.pl) and Pilon (v 1.23) (Walker et al., 2014). Assemblies for some strains (PP6, PP9, PP11, PP12, PP13, PP14, PP15, PP16, and S589) were not corrected because they were not sequenced on Illumina. Finally, circular contigs were permuted to start at dnaA and linear contigs (only S602) were split into two pieces such that one piece started with dnaA. For subsequent analyses, PacBio assemblies were used if available; otherwise, Illumina assemblies were used.

Before annotation, the species of each assembly was confirmed using Taxator-tk (v 1.2) (Droge et al., 2015). The genes in each assembly were annotated using Prokka (v 1.13) (Seemann, 2014), and then gene clusters were identified using Roary (downloaded from github on January 22, 2019) (Page et al., 2015). Gene clusters with at most one copy per genome were aligned using Prank (v 0.100802) (Loytynoja, 2014) and then concatenated, adding gaps for genes absent from any alignments. This alignment was used as the input alignment to make a phylogenetic tree with RaxML (v 8.2.4) (Stamatakis, 2006). The tree was midpoint rooted using ETE (v 3.0) (Huerta-Cepas et al., 2016). To search for gene presence/absence differences associated with the time periods of interest, we gave the roary gene clusters and the phylogeny to Scoary (v 1.6.16). Both the gene clusters and the phylogeny were filtered so the only contained strains were 35B and ST558.

To find SNPs between the 35B, ST558 strains, we called variants using Snippy (v 4.4.5) (Seemann, 2015). Isolate PP6 was used as the reference strain. For each query strain, we used Illumina reads; otherwise, we used the PacBio assembly for that strain (strains without Illumina reads are listed above). The list of variants was filtered for significant hits using the filtered phylogeny and Scoary (v 1.6.16). These top hits were mapped to genes using Prokka’s gene annotations and SnpEff (v 4.3) (Cingolani et al., 2012). Finally, the gene clusters from Roary were used to map hits between strains. The genes that were not identified by Prokka were further identified by BLAST search.

The best candidate for the parent strain was identified by examining the phylogeny and using BLAST (v 2.2.28) with the sequences flanking the capsule. The candidate for the capsule donor was identified by using BLAST with the sequence of the capsule and the surrounding region. Each strain was then aligned to S558 (the strain with the capsule switch) using Mauve (v 2015-02-13) (Darling et al., 2004). In order to improve the quality of the alignments, contigs were first reordered using Mauve’s contig mover. SNPs and indels were extracted, and sliding windows of each were computed using a window size of 1,000 bp and an increment of 100 bp. Capsular serotypes identified by serotyping methods were further confirmed with PneumoCat (v1.2) (Kapatai et al., 2016).

For pilus-1 operon (Regev-Yochay et al., 2010), genes associated with antibiotic resistance [β-lactam: pbp1a, pbp2x, pbp2b (Li et al., 2017), fluoroquinolones: gyrA, parC (Brueggemann et al., 2002), macrolides: ermB, and tetracycline: TETM (Chochua et al., 2017)] were annotated and extracted with PATRIC. (Brettin et al., 2015; Wattam et al., 2018). PBP typing was classified based on the penicillin binding protein gene type list available at (https://www.cdc.gov/streplab/pneumococcus/mic.html?CDC_AA_refVal=https%3A%2F%2Fwww.cdc.gov%2Fstreplab%2Fmic-tables.html) as of February 2021. For fluoroquinolones, where extracted gene sequences have previously reported point mutations on gyrA and parC were analyzed (Brueggemann et al., 2002). The presence of the ermB gene and TETM were checked using PATRIC.

The difference in the case (AOM) to colonization ratio (the number of isolates from MEF at the time of AOM/the number of isolates from NP at the time of health) was calculated for each respiratory year (June–May). The ratio between pre-2015 (June 2006–May 2015) and post-2015 (June 2015–May 2019) was analyzed by the Mann-Whitney test. The change in ratio among ST558 was analyzed by Fisher’s exact test. The trend in antibiotic MICs was computed using a linear regression model. Comparison of MIC50 values was conducted using the Mann-Whitney test with a Bonferroni correction.

We re-examined 250 serotype 35B isolates (consisted of 11.3% of total SPN isolates for the study period) that had been collected from 140 children (18.6% of enrolled children) during 2006 to 2019. Of these, 172 isolates (12% of SPN isolates) from the NP of healthy subjects, 60 isolates (11.7%) from the NP at the onset of AOM, and 18 isolates (13%) from MEF were collected. Yearly prevalence of 35B isolates as a proportion of all pneumococci is shown in Figure 1 . Before introduction of PCV13 in 2010, serotype 35B was rarely isolated. After 2011, the relative prevalence of colonizing serotype 35B strains (from healthy NP) significantly increased from 2.7 to 6.1% of all SPN (from June 2006 to May 2011) to 12.3–18.8% (from June 2011 to May 2015). Despite colonization prevalence, AOM infections were infrequent. Beginning from June 2015, the prevalence of 35B increased, and it became the most frequently isolated serotype from MEF during AOM bouts. To see the risk of AOM given the frequency of colonization, we calculated the case-to-colonization ratio changes among 35B isolates over time and found AOM cases (MEF isolates during AOM) as a proportion of controls (NP isolates from healthy colonization) significantly increased from 0.08 for the period before June 2015 to 1.68 for the period after May 2015 (p=0.002). Thus, introduction of PCV13 in 2010 was followed by increased prevalence of NP colonization with serotype 35B strains, and an additional shift to higher prevalence and increased AOM disease burden beginning in 2015.

Among 163 35B strains tested, 13 distinct sequence types (STs), including two newly identified STs were identified ( Figure 2A ). Most 35B strains belonged to ST558 (89% of pre-2015 and 70% of post-2015 isolates). Although ST558 isolates were detected at comparable levels in healthy children among pre-2015 and post-2015 isolates (6.0 and 5.6%), their isolation among children causing AOM significantly increased post-2015. When the ratio of AOM disease cases to healthy colonization was calculated only among ST558, values for pre-2015 and post-2015 were 0.09 and 2.18, respectively (p=0.0002), underlining the emergent AOM disease burden caused by ST558 SPN strains.

To search for evidence of 35B capsular switching, we clustered STs into clonal complexes (CCs) using goeBurst ( Figure 2B ). We found three STs in pre-2015 and two STs observed in post-2015 samples were closely related to ST558 (circled “CC558” consisting of single allele variants from ST558, Figure 2B ). A more distantly related cluster of three STs separated by a single allele variant including one newly described ST14687 (CC156) was identified; this likely represents a recently described 35B capsule switch between 35B/ST558 and 9V/ST156 serotype strains by Chochua et al. (2017). We also found evidence of a new potential capsular switch. We identified a post-2015 serotype 35B isolate belonging to newly described ST14683, which shares only one allele (ddl) with CC558 and none with those in the “CC156” strains. The ST14683 isolate was detected from two different children. The new ST14683 is a single allelic step from ST63, which was associated with serotype 15A (“CC63”) in our cohort.

Genome assemblies were used to produce a phylogenetic tree based on protein-coding genes ( Figure 3 ), which found three distinct clades and was consistent with the ST-based results. The largest clade included ST558 and related STs, and strains were isolated from all categories (further discussed below): AOM (MEF and NP) and colonization and pre-2015 and post-2015. The second largest clade contained CC156 isolates that were all collected post-2015, consistent with reports of a recent emergence of this lineage (Chochua et al., 2017). The final branch of serotype 35B consisted of a single ST14683 (S558) strain that was isolated from MEF at an AOM visit that grouped closely with serotype 15A isolates from ST63.

To confirm this novel capsule switch and reconstruct the recombination breakpoints, we produced and aligned nearly finished genomes for the 35B/ST14683 capsule switch strain S558 and two representatives of the parental strain backgrounds (for the donor, 35B/ST558 isolate S548, and for the recipient, 15A/ST63 isolate S534). Comparison of nucleotide divergence across the capsule locus and its flanking sequences shows near-perfect identity between the capsule locus and its immediate flanks with the 35B donor representative, whereas the more distant flanks are nearly identical to the 15A recipient representative ( Figure 4A and Supplementary Figure 1 ). The precise breakpoints were obscured by high divergence from both parents at intermediate distances from the flanks, including the presence of indels, suggesting the involvement of a third donor or more complex history of divergence within 15A strains. Consistent with past recombination at other loci from 35B/ST558 and other strains into ST14683, most MLST alleles in ST14683 were shared with 15A/ST63 strains, but also included two alleles shared with ST558 and three novel alleles, including a novel potential recombinant allele of PBP 1a ( Figure 4B ). This ST isolate was negative for rrgC (PI-type 1) that is a common character of 15A and positive for ermB and tetM that were the same profiles as 15A ST63 detected in our cohort.

In order to identify potential SPN genetic determinants that could be responsible for the increased disease burden associated with the 35B strains post-2015, we focused on genomic variation within ST558, since most of our sequenced 35B strains belonged to this ST and these had been isolated from across the study period. The ST558 gene phylogeny did not group strains by time period, colonization/AOM status, or NP/MEF isolation during AOM, indicating that emergent 35B disease burden post-2015 is not due to a specific ST558 clone ( Figure 3 ). We next sought to identify whether specific genetic variants—either gene presence/absence differences or SNPs—within the ST558 lineage were associated with isolation from disease post-2015. To examine gene presence/absence differences within ST558, 2,072 clusters of orthologous genes were identified, of which 1,804 (87%) were found in >95% of strains. The remaining accessory (or distributed) genes showed no significant associations with time period or health status ( Supplementary Figure 2 ). Because the resurgence of rrgC (PI-type 1) positive strains were observed among non-vaccine type over time and proposed to be selective virulence advantages (Croney et al., 2013), we specifically checked for rrgC. We confirmed that most of ST558 isolates through the entire study period were positive for rrgC.

SNP analyses identified 138- post-2015-occurring SNPs that were not detected from strains isolated prior to 2015 that were significantly associated with the post-2015 era ( Table 1 ). These SNPs were distributed among 34 different genes. Ten percent of these SNPs were found within the nrdG (metal binding) gene; 9% within the metP and metN genes (members of an ABC transporter); 6% within the corA (Mag2+ transporter) gene; 5% within the priA gene (DNA replication); and 3% within the enzymatic gene, ldcB (LD-carboxypeptidase).

The MICs for 35B isolates significantly increased over time for 7 of 16 antibiotics tested ( Figure 5 ). Six of the seven antibiotics were β-lactams, and one was a fluoroquinolone. For amoxicillin, cefotaxime, and meropenem, linear regression analysis showed that the resistance category changed from sensitive to intermediate over the period from 2006 to 2019. The MIC trend for ertapenem and ofloxacin were close to category change (sensitive to intermediate). To identify STs associated with the increase in antibiotic MICs, MIC50s for each ST were also compared ( Table 2 ). In general, STs detected after 2015 showed higher MIC50 in six antibiotics (four different β-lactams, ofloxacin, and TMP/SMX). Ofloxacin and penicillin MIC₅₀ among ST558 showed significant increases between the two time periods. Although only two isolates were detected as ST14683, it is noteworthy that this ST showed multidrug resistance against 6 β-lactams, ofloxacin, tetracycline, and TMP/SMX.

Since we detected increased MICs in β-lactams and ofloxacin among 35B strains over time, we classified β-lactam resistance determinants (PBP types:PBP1a, 2b and 2x) (Li et al., 2017), and fluoroquinolones resistance determinants (gyrA/parC genes) (Brueggemann et al., 2002) in our whole genome sequence data. 35B from pre-2015 only associated with PBP type 4:7:7 except two isolates that had one amino acid mutation on PBP1a and PBP2b. CC558 (ST558, 10493) from post-2015 was also exclusively associated with PBP type 4:7:7. ST156 and ST3280 were associated with PBP type 4-12-7. Our only sequence from ST14683 had a novel PBP1a type that was not found in the list of PBP types provided by the US CDC ( Technical Appendix 1 ).