We created a pangenome of 84 different pneumococcal isolates, spanning both genes and the non-coding IGRs. To put the nature of the pneumococcal IGR pangenome into perspective, we performed the same analysis for S. aureus. Both species may colonize the human upper respiratory tract, both are opportunistic pathogens and both possess open pangenomes, making them prime candidates for comparison (Laux et al., 2019). The analysis shows that the otherwise non-coding IGRs of both species are conserved in a similar manner to genes across the pangenome, although the number of unique genes outnumber the number of unique IGRs in both species (Figure 1).

While the proportion of core genes roughly scales relative to the size of the genome (pneumococcus 2.1 Mbp/S. aureus 2.8 Mbp) the proportion of core IGRs relative to genome size is lower in pneumococcus (Figure 1). However, pneumococcus seemingly compensates for the lower number of unique IGRs by having multiple copies of several core IGRs in each genome. On average, each core IGR is present 1.23 times in each pneumococcal genome compared to 1.03 times in S. aureus. Each core gene is present 1.08 times in each pneumococcal genome and 1.02 times in S. aureus, this indicates that the high copy number of pneumococcal core IGRs is quite unusual.

The average pneumococcal genome has 79% of its genes as core genes and 66% of its IGRs as core IGRs. The average S. aureus genome is comparatively close to that observed in pneumococcus, here core genes constitute 79% and core IGRs 68% of each genome (Table 1).

Despite pneumococcus having fewer unique core IGRs relative to genome size than S. aureus, as stated earlier, their copy number is higher in each genome, thus the percentage of core IGRs per genome is roughly equivalent in the two species (Table 1). The higher copy number of core IGRs in pneumococcus is also illustrated by the fact that 669 unique core IGRs exist in the pneumococcal core genome (Figure 1) but on average each genome has 817 core IGRs (Table 1).

The number of unique IGRs in the pneumococcal pangenome increases with the number of isolates analyzed in a similar manner to the number of unique genes (Figure 2A). Overall, fewer unique IGRs are present in the pangenome than genes, part of this is due to the exclusion of IGRs of <30 bp in length, which are most often intraoperonic (Thorpe et al., 2018).

Most IGRs are either present in almost all pneumococcus isolates or unique to only a few, that is, they are either very common or very rare (Figure 2B). Pneumococcus genes show a similar distribution across the pangenome, though a larger proportion of IGRs are confined to only a few isolates than genes.

Pneumococcus retains more unique genes than IGRs within its pangenome (Figure 2A), and most unique IGRs are only found in single isolates, making them rare (Figure 2B). This scarcity of unique IGRs could indicate that IGRs experience a higher evolutionary selection threshold than genes, thereby lowering the likelihood of a newly acquired IGR of spreading to more isolates through HRT.

IGRs can be categorized according to the orientation of their flanking genes. IGRs that are downstream of two convergently transcribed genes are considered non-regulatory (NR), IGRs that are upstream one gene and downstream another gene are considered single regulatory (SR) and IGRs that are between two divergently transcribed genes are considered double regulatory (DR) (Figure 3).

Looking at the distribution of the IGR types across the pneumococcal pangenome, NR and DR regions are rare compared to SR IGRs (Figure 3). DR regions also constitute a greater relative proportion of the core IGRs than seen in the accessory genome.

Next, we analyzed the degree of linkage between core IGRs and core genes, that is, how often a core IGR is directly upstream a core gene. IGRs and their flanking genes were identified and any IGRs directly upstream the start codon of a gene was selected. The status of the IGR/gene pairs as accessory or core genome was then assessed and the ratio of each combination calculated. On average 81% of core genes in S. pneumoniae are associated with a core IGR, whereas only 74% of accessory genes are linked to accessory IGRs (Table 2). For comparison, the linkage of core genes to core IGRs is greater in S. aureus at 86%, and accessory IGRs are flanking accessory genes 82% of the time.

None of the IGRs of the capsular polysaccharide synthesis (cps) operon were found to be core IGRs, however the highly conserved flanking genes dexB and aliA were both associated with core IGRs (Appendix 1).

Next, we examined the IGR candidates for regulatory switching. Regulatory switching describes when one IGR is replaced by a different non-homologue IGR. The origin of these switched IGRs is not inferred in this analysis, thus they can both originate from within the isolate itself or even from a separate species. For this analysis, only switches where the IGRs share no significant sequence homology with a BLASTN were included.

We detected three switches within the pneumococcus pangenome and only one of these is flanked by core genes. We designated the core switched IGR as csIGR (Table 3). While the two versions of the csIGR are highly conserved on their own, with both having a nucleotide identity of >99% amongst themselves, aligning the two versions with each other results in an insignificant nucleotide identity of 57%. These results were manually confirmed with a blastn and confirmed that all pneumococcal isolates always have one of these two csIGRs but only in a single copy and always between the same flanking genes.

The flanking genes were both single copy core genes and were identified in the common lab strains S. pneumoniae D39 and TIGR4 (Figure 4). These two strains have distinct csIGR types, with D39 having csIGR1 and TIGR4 having csIGR2 (Figure 4). Interestingly, rather than flanking an operon, the IGRs are predicted to sit in the middle of an operon. Little is known about the flanking genes, other than their status as single copy core genes found in this study. SPD_1559/SP_1749 is considered essential in pneumococcus and is a homologue to ygeH, a gene involved in biogenesis of the 30 S ribosome subunit (Pek et al., 2007; Liu et al., 2017). The sequence of both csIGR types is provided in Supplementary Appendix SA2. Interestingly, neither of the csIGR types were confined to a specific phylogenetic cluster of pneumococci (Figure 5). The fact that the csIGR types are spread across the phylogenetic tree indicates that their distribution is due to HRT.

We performed a pangenome wide association study to see if any genes within the accessory genome were significantly co-occurring with the csIGR alleles across the pangenome. However, no genes were exclusively associated with neither of the csIGR types. Both sequences were also screened for promoters, riboswitches and homology to known regulatory RNAs with no significant hits. However, the translated RNA sequence of both sequences was predicted to form significant secondary structures, the significance of which remains to be elucidated. The predicted secondary structures are provided in Supplementary Appendix SA2.

All 84 complete S. pneumoniae genomes available from the National Center for Biotechnology Information, GenBank resource was downloaded in raw FASTA format. Additionally, 84 randomly selected S. aureus genomes were retrieved for comparison with S. pneumoniae (12/5/2021). Genomes were then annotated with Prokka (v 1.14.5), using the standard parameters of the software (Seemann, 2014). The genomes used are listed in Supplementary Appendix SA3.

Initially a pangenome of the coding sequences (CDS) was created using Roary (v3.13.0) (Page et al., 2015). Then a complementary pangenome of the IGRs was created using Piggy (v1.5), an intergenic pangenome analysis tool that emulates Roary (Thorpe et al., 2018). Some steps were taken to ensure comparability between the outputs of the software. Roary was set to cluster CDSs with -e -n (to perform alignments using MAFFT (Katoh et al., 2002)), -i 90 (90% sequence identity cut-off) and -s (to not split paralogs into separate clusters). The settings for running Piggy were set at the standard parameters of the software, except for -len_id 10 (the minimum percentage of length identity to form a cluster). The length identity was reduced for comparability with Roary, as gene clusters generated by Roary only require a sequence length identity of 120 bp for clustering CDSs, thus the len_id of 10 is recommended by the creators of Piggy for Roary consistency as IGRs are not erroneously placed into separate clusters (Koonin et al., 2001; Dagan et al., 2008; Molina and Van Nimwegen, 2008; Ragan and Beiko, 2009; Peters et al., 2011; Matus-Garcia et al., 2012; McCutcheon and Moran, 2012; Somvanshi et al., 2012; Oren et al., 2014; Page et al., 2015; Xiao et al., 2015; Ochman and Caro-Quintero, 2016; Thorpe et al., 2017; Thorpe et al., 2018; Jørgensen et al., 2020). The randomly assigned locus tags provided by Prokka were translated when necessary, by aligning the GFF files of the Genbank annotated files and Prokka output.

The linkage of core genes and core IGRs was quantified using R (v. 4.1.0). The gene_presence_absence file from Roary and the IGR_presence_absence file from Piggy was loaded as dataframes in R. For each gene and IGR cluster in the files their status as a core or accessory gene was identified and assigned. For each genome all IGRs were paired with their upstream gene. Thus, NR regions were removed from the dataset and both flanking genes for DR regions were analyzed separately, if any of the two genes were core, the DR IGR was assigned as flanking a core gene. The R code is provided in Supplementary Appendix SA4.

For identification of switched IGRs, a separate analysis using Piggy was performed with -len_id 90 (the minimum percentage of length identity to form a cluster). This was done to perform a more strict analysis of the IGRs, as the higher threshold for forming a cluster ensured that homologue IGRs were not identified as switched IGRs (Thorpe et al., 2018). IGR switches were identified using the “gene-pair” method of Piggy, here two or more different IGR sequences that occupy the same space between a specific gene pair are analyzed. The candidate IGR sequences are then aligned with BLASTN with low complexity filtering turned off and if there are no significant matches between the IGR they are identified as “switched”. If there is a significant match Piggy aligns the sequences using MAFFT and provides the nucleotide identity of the alignments.

The identified switch was validated manually with a BLASTN against all the genomes (data not shown).

A phylogenetic tree of the strains included in this study was created, based on single nucleotide polymorphisms (SNPs) in the core genes. Roary was run separately with the same settings as previously mentioned with the exception of -e (Core gene alignment with PRANK) (Page et al., 2015). This produced a highly accurate alignment of the core genes within the pangenome. FastTree (v2.1.11) was then run to infer an approximately-maximum-likelihood phylogenetic tree based on SNPs within the core genes (Price et al., 2009). The resulting newick file was then visualized using iTol (v5.7) and exported to Adobe Illustrator (Letunic and Bork, 2007).

To identify whether any accessory genes were significantly associated with any of the csIGR alleles, a pangenome-wide association study was performed using Scoary (v1.6.16) (Brynildsrud et al., 2016). A trait matrix was created as an input for Scoary, indicating which of the two csIGRs alleles were present in which genomes. Scoary then sorted the accessory genome provided by the gene_presence_absence file from Roary, scoring each accessory gene according to their co-occurrence with each csIGR.

To assess homology to existing regulatory RNAs, a BLASTN was performed for each csIGR against the RefSeq RNA database. Both sequences were screened for potential riboswitches using Riboswitch Finder (Bengert and Dandekar, 2004). Any potential promoter or terminator regions were screened for using BPROM and FindTerm (Softberry) (Solovyev et al., 2011). Secondary structures were predicted for both sequences using the RNA structure package available through Mathews lab, at standard parameters (Reuter and Mathews, 2010).