The bacterial genus Salmonella consists of organisms responsible for substantial human and livestock morbidity and mortality by causing a spectrum of diseases manifesting as enteric fever, septicemia, and gastroenteritis (Dekker and Frank, 2015). Salmonella enterica is the main species of medical and veterinary importance and consists of about 2,500 serovars. Based on their degree of host adaptability in both humans and animals, Salmonella are classified into two groups (Gal-Mor et al., 2014; Jajere, 2019). The serovars adapted to specific hosts typically cause systemic, typhoid-like symptoms, e.g., Typhi, Paratyphi, Gallinarum, Abortusovis, and Pullorum, and are distinct in behavior from the generalists which can infect many hosts. The generalists cause only a limited pathology, most commonly a gastroenteritis, and are referred to as non-typhoidal Salmonella (NTS), e.g., Typhimurium and Enteritidis. The NTS represent a major cause of diarrhea globally (Cheng et al., 2019). Although only a small proportion of the Salmonella serovars have been associated with foodborne outbreaks (Khakhria et al., 1997; Centers for Disease Control and Prevention, 2016; European Food Safety Authority, 2022), a single food source may be contaminated with multiple serovars or multiple strains of the same serovar (Whitney et al., 2021; Schwensohn et al., 2022), which makes identifying the source of an infection very difficult and continues to hinder the control of human salmonellosis.

Whole-genome sequencing (WGS) is now the main approach for characterizing isolates of Salmonella for outbreak investigation and surveillance and has almost completely replaced the traditional phenotypic and molecular typing tools such as serotyping, pulsed-field gel electrophoresis (PFGE), and phage typing (PT). Determining the serotype or serovar is desired because of its relevance to Salmonella epidemiology, and the serological agglutination procedure traditionally used is now mostly replaced with bioinformatics analysis of sequence reads as performed by means of two main algorithms, namely, the Salmonella In Silico Typing Resource (SISTR) software (Yoshida et al., 2016) and the SeqSero2 software (Zhang et al., 2015). Delineating strains of the same serovar was traditionally achieved using the PFGE and PT methods but is now accomplished by one of two approaches. The first approach is to determine the sequence type of the organism and is commonly achieved using the multilocus sequence typing (MLST) method which was originally designed to evaluate nucleotide sequences of five to seven housekeeping genes of the organism (Maiden et al., 1998). This procedure can now be more easily performed, e.g., by uploading the genome reads into the web-based EnteroBase platform to provide multiple, hierarchical levels of characterization and strain differentiation (Zhou et al., 2020). The second approach for strain differentiation relies on single-nucleotide polymorphism at variable sites in the genome. Many algorithms have been developed for the purpose of characterizing SNPs in either the core genome or the pan-genome. The SNVPhyl algorithm developed by the Public Health Agency of Canada identifies high-quality SNPs among a set of selected isolates and is useful for generating phylogenetic trees from these SNPs (Petkau et al., 2017). A similar algorithm, the SnapperDB which was developed by the United Kingdom Health Security Agency (UKHSA) also identifies high-quality SNPs useful for evaluating genetic distances among the genomes and for inferring relatedness among strains (Dallman et al., 2018). Parsnp is yet another algorithm which detects core-genome SNP in bacterial genomes and with the aid of an adjunct interactive tool known as Gingr, displays informative overviews for specific sub-clades and genomic regions (Treangen et al., 2014). The kSNP tool detects SNPs in the pan-genome but is uniquely able to carry out comparisons among genomes without a requirement for genome alignment or a reference genome (Gardner et al., 2015).

WGS tools were applied to investigate an outbreak of Salmonella enterica serovar Adjame in England in 2017. The organism is a rare NTS for which the first clinical case was reported in Cote d’Ivoire in 1967 (Le Minor et al., 1967). The first documented clinical case in the western world was observed in the United Kingdom in 1993 and was followed by 13 sporadic cases leading up to the year 2016 (Chattaway et al., 2019). An epidemiological investigation of the latest outbreak identified herbs and spices as the most likely food vehicles (Chattaway et al., 2019), and WGS analysis using both EnteroBase Salmonella core-genome MLST (cgMLST, allele-based method¹) and UKHSA SnapperDB pipeline (Dallman et al., 2018) showed the presence of heterogeneous strains which were grouped into two clusters (Chattaway et al., 2019). The two outbreak clusters were shown to be distinct from a third cluster comprised of S. Adjame strains from sporadic cases present in the UKHSA archive. Strain characterization and identification of multiple strains relied on core-genome analysis but without a good understanding of the population structure of the rare organism, the question arises whether the strain characterization was sufficiently rigorous to provide information that may facilitate an intervention to forestall a spread of a current outbreak or combat a future outbreak. This prompted us to investigate the accessory genome of isolates of S. Adjame recovered before, during, or around the time of the 2017 England outbreak to further characterize the population and ask whether the analysis could lead to a more defined clustering or alternatively confirmed the heterogeneity of strains involved.

A major contributor to the accessory genome of the Salmonella organism are prophages or bacteriophages which have become integrated into the bacterial chromosome (Ogunremi et al., 2014). Other contributors to the accessory genome include a myriad of genetic sequences that enter the bacteria cell by horizontal gene transfer and may end up in the chromosome or in a plasmid and includes the following: antibiotic resistance genes, insertion sequences, genomic islands, integrative conjugative elements, and other transposons (Subedi et al., 2019). Bacterial plasmids as an entity also meet the definition for inclusion in the accessory genome since they are not always present in all or a majority of the strains of a species. Bacteriophages are considered the most abundant “life form” on the planet (Wommack and Colwell, 2000), representing an incredibly diverse gene pool, and have been shown to contribute significantly to host bacterial evolution (Fortier and Sekulovic, 2013). As part of their parasitic life cycle in bacteria, bacteriophages can either be virulent, and follow a lytic life cycle or can integrate into unique attachments sites in the chromosome of the bacteria by following the lysogenic cycle in which case they are referred to as prophages (Kropinski et al., 2007; Switt et al., 2015). Prophages have proven to be suitable markers for differentiating S. enterica subtypes because of the large but stable variation present in phages (Gao et al., 2020). This variation was exploited to develop a highly discriminatory prophage subtyping tool (PST) which not only distinguished serovars but possessed a very high resolution that is useful for defining epidemiologically unrelated Salmonella enterica serovar Enteritidis strains during foodborne outbreaks (Mottawea et al., 2018). The discriminatory ability of the method is based on variations in the composition of prophages among strains of the same organism as well as size and nucleotide differences that may occur in the same prophage in the different strains (Mottawea et al., 2018). Additional observations from the study showed that sufficient diversity in a single phage genome, namely RE-2010 phage has enough variability for distinguishing among Salmonella Enteritidis strains from different outbreaks (Mottawea et al., 2018). Furthermore, Goodridge and colleagues demonstrated a similar discriminatory capacity in the tyrosine integrase gene of 32 enteric prophages for the strain characterization of members of S. enterica (Colavecchio et al., 2017).

In the current study, we used the PST pipeline to analyze 38 strains of S. Adjame which represents the entire population of the organism with genome sequences available in the GenBank at the initiation of this study, providing an opportunity to characterize the strains by focusing on the accessory genome to evaluate differences among strains and to compare outbreak and sporadic strains. The study has provided an early understanding of the genomic diversity and population structure of the rare S. Adjame based on the available genome resources available, albeit limited. An application of core and accessory genomes analytical tools for future S. Adjame strains, either from outbreaks or sporadic sources should help to accurately define relatedness to hitherto observed strains and possibly facilitate a decision in support of an early intervention to prevent a developing outbreak which could have been dismissed in the absence of any knowledge of the genomic structure of the population.

In this study, we have carried out whole-genome analyses of all 38 strains of the rare Salmonella serovar Adjame available in the GenBank at the inception of this study. The strains included those recovered from 14 cases identified in the only known outbreak of S. Adjame, which occurred in the summer of 2017 in England. The small size of the entire population of outbreak and sporadic strains provided an excellent cohort to dissect the variability between the core and accessory genome of the strains and to overlay the analytical inferences on the epidemiology of S. Adjame, especially regarding the time of transmission to individuals. Analysis of the 2017 outbreak by Chattaway et al. (2019) using the Snapper DB high-quality SNP pipeline and EnteroBase Salmonella core-genome multilocus sequence typing scheme showed that the clinical cases were due to heterogeneous strains of S. Adjame rather than a point source exposure. In this study, our goal was to carry out a comparative analysis of the core and accessory genomes of the organisms to evaluate the extent of the heterogeneity of organisms involved in the outbreak, compare the outbreak strains with the sporadic strains including those isolated around the time of the outbreak, in order to provide an insight into the population structure of the rare Salmonella organism. At the onset, we evaluated the variation in the core genome of the Adjame strains by means of two SNP analyses, namely Parsnp and kSNP, and observed a clustering of the outbreak strains similar to that reported by Chattaway et al. (2019). Our results confirmed a considerable genetic variation among the strains identified as part of the outbreak. Our Parsnp analysis (Figure 1) showed that organisms identified as part of the 2017 outbreak (green and red cluster) had a closer ancestry when compared to all the other strains analyzed including organisms belonging to the blue cluster, which is in contrast with the observations by Chattaway et al. (2019). These results indicate that different SNP analysis algorithm or analysis of the core genome may give different clustering outputs, and additional information may be required from further analysis of the genome to clarify strain relatedness. To that end, the accessory genome may shed the required light to evaluate the characteristic of the population including relatedness of the strains and provide evidence-based inferences about evolutionary changes in the population. Prophage genes form the majority of the bacterial accessory genomes and have great potential to be used as highly discriminatory markers for Salmonella (Mottawea et al., 2018) and Listeria monocytogenes (Zamudio et al., 2020). In this study, we analyzed the prophage sequences in the S. Adjame strains previously reported as part of the atypical outbreak (June—July 2017), and of a group of strains recovered just prior to the outbreak, mainly between March and April 2017), as well as other sporadic strains recovered as early as October 2008 or as late as December 2018.

The pronounced diversity in the prophage genomes among our collection of S. Adjame strains verified and reinforced the heterogeneity of the organisms recovered from the atypical outbreak (Figures 3, 4) Twenty-four different prophages were detected while only four of them were present in the majority of the studied strains. Upon further dissection of the individual prophage sequence profile, we found that one nucleotide variant in each of the conserved prophage region seemed sufficiently informative for assigning certain strains to a cluster (Figures 6–8). For instance, the Adenine/Guanine base in the conserved regions of prophage Entero_p88 (Figure 5) and the Cytosine/Thymine base in the conserved regions of the prophage Burkho_phiE255 (Figure 7) were used to separate the March 2017 sporadic strains from the rest of the strains recovered during the period spanning the outbreak, i.e., June to July 2017, highlighting the potential use of a single nucleotide in the accessory genome as a marker to differentiate between sporadic and outbreak strains. Notably, the detected Adenine/Guanine location in the prophage Entero_p88 was not included in any of the predicted ORF, which showed that mutations in a non-coding region could provide adequate information in strain characterization. Other authors have demonstrated the potential role of non-coding regions of Salmonella in pathogenesis in animals and virulence (Gong et al., 2011; Hammarlof et al., 2018). For prophages that showed greater variability within the conserved regions such as prophage Edward (Figure 8A), three different prophage patterns were observed among the 28 strains that harbored the prophage, and these translated into three distinct clusters in the phylogenetic tree (Figure 8). The clustering pattern observed in the phylogenetic tree developed for the prophage Edward sequence (Figure 8) was similar to the SNP analysis clustering using the SnapperDB algorithm (Chattaway et al., 2019). Thus, the inclusion of additional prophages in our analysis predictably led to a demonstration of an even greater heterogeneity among the Adjame strains. Indeed, when we combined the four conserved prophage sequences to generate a phylogenetic tree and overlaid the groups with the previously described cluster designations (i.e., red, green and blue clusters, Figure 9), we were able to clearly show an even greater degree of heterogeneity, which became even enhanced in the phylogenetic tree produced using all 24 different prophages (Figure 4). We have used the PST to show that highly related Salmonella are indistinguishable from one another (Ogunremi et al., 2014); however in this study, only a few Adjame strains appear to be very closely related (e.g., SRR6466751 and SRR6466752). We propose that the 2-year time differences between these strains (2011 and 2013) may have led to mutational changes detectable by the PST since the assay is affected by nucleotide variation in a target prophage if the assay parameters are set at a high stringency as we have done. Our study showed that prophage sequence analysis provided enough granularity to distinguish all the strains of Adjame in the population from one another. It is possible that the biology of the organism or the niche environment occupied by the strains predisposed them to one or all of the following: receptivity of the organisms to phages, integration into the genomes, and rapid changes in the chromosomal integrated bacteriophages. Prophages are known to represent a diverse arsenal that allows bacteria to adapt to their ever-changing and competitive environments (Bobay and Ochman, 2017), and the diversity within and among the many prophages observed in this study is supportive of this claim. However, we also noticed limited size variation in two prophages namely, the prevalent Burkho_phiE255 (35 of 38 strains) and the less prevalent Entero-cdtl (19 out of 38 strains). We speculate that the high conservation seen in these prophages suggests that they may not be a factor in driving the diversity of Adjame. It will be interesting to see if the size conservation of the prophages has arisen from their low susceptibility to recombination with other phages.

We noticed that the outbreak strains obtained during the month of June 2017 and characterized by SNP analysis as the green cluster showed the most closely related prophage profiles (Figure 4). Another isolate not identified as part of the 2017 outbreak, namely, SRR6237069, but nevertheless recovered during the same time as the outbreak (June 2017) showed relatedness in profile prophage with the other green isolates as defined by SNP (Figures 1, 2) and prophage analysis (Figure 4). Furthermore, the four prevalent prophages were present in SRR6237069 and segregated as did the rest of the green isolates. The prophage analysis also placed the ERR2071995 isolate obtained from Ireland and grouped by SNP analysis with the green cluster (Chattaway et al., 2019; Figures 1, 2, 6). The additional observation obtained here from application of the highly discriminatory prophage typing would therefore suggest that the outbreak was international in scope based on close genetic relatedness of the isolates both at the core and prophage genome levels.

One of the isolates from the green cluster, SRR6233939, showed a similar prophage profile to another isolate from the red cluster. This green isolate was exceptional in that it was the only one missing the Entero_p88 isolate and therefore occupied the most distant position among the green group by prophage analysis. Along the same line, SRR6192965, which has a core genome that grouped as a red cluster showed noticeable variation in the prophage profile as indicated by a long branching of the prevalent Burkho_phiE255 prophage (Figure 7). This, and other prophage changes led to a dissociation of SRR6192965 from the other strains that clustered into the red group when all four prevalent phages were analyzed (Figure 9). The prophage analysis showed a demarcation of the outbreak strains that were grouped together into the red cluster based on the core genome into two groups and two singletons, indicating an even greater heterogeneity than previously described by Chattaway et al. (2019). The variation seen among the prevalent prophages may shed light into the commonality among strains, provide evidence of the substructures within the population, and enable the tracking of closely related strains (e.g., green isolates from England and Ireland).

Two strains obtained during the course of outbreak in June 2017 (SRR6190984) and July 2017 (SRR6237093) as well as a third strain obtained months after the outbreak in March 2018 (SRR7351460) were unique in their prophage profile and did not show relatedness with any other isolate (Figure 4). Only one of these strains was grouped as part of the outbreak by clustering with the red group (SRR6190984). While this strain is scored as part of the red cluster based on core-genome level, it is very different from the others at the prophage level. Indeed, analysis of the prophage Entero_p88 in this strain revealed a strong similarity with another group of five, non-outbreak strains obtained between October 2008 and June 2017 (Figure 5) and these five strains formed a loose cluster at the core-genome level according to Chattaway et al. (Figure 5). Thus, it is remarkable that a strain involved in the England outbreak (SRR6190984) and which clustered with other outbreak strains (red cluster) could have at some point acquired a version of the Entero_p88 prophage or mutations that are different from those by other members in the cluster, but similar to another version of the Entero_p88 seen in another group with a different core-genome attribute. This may suggest that although the “genetic backbone” of the SRR61090984 as defined by the core genome is related to the other outbreak strains, it did acquire accessory genomic attributes by occupying a similar niche with the five additional strains that were not involved in the outbreak. In that sense, the evaluation of the core- and accessory-genome properties of the strains helped to show lineage similarities regardless of the time of recovery from patients. Thus, while the atypical outbreak was defined by the time the patients presented with symptoms using valid epidemiological parameters (i.e., between June to July 2017), there were strains obtained during and outside the time period, by a combination of core-genome and accessory-genome analysis, but not the former alone, that demonstrably and convincingly displayed close relatedness. Prophage level relatedness could therefore point to a shared history of strains (e.g., niche) even when there is little similarity at the core-genome level.

For such a rare serovar, we were interested in identifying the closest genetic serovar relatives. Our Neighbor-joining phylogenetic tree consisting of all 38 S. Adjame strains and 9,066 Refseq Salmonella enterica genomes from the NCBI RefSeq database (Figure 10), showed that the closest serovar to S. Adjame was S. Mississippi. Apart from the closer relatedness shown by WGS-based analyses, these two serovars also shared some somatic O (O) antigen and Flagella H phase 2 (H2) antigen, i.e., S. Adjame with O (13, 23) and H2 (1, 6) and as S. Mississippi with O (1, 13, 23) and H2 (1, 5), respectively. In contrast to S. Adjame, the serovar Mississippi has been detected from many parts of the world and from various sources including humans, animals, and environmental sources (Ashbolt and Kirk, 2006; Ford et al., 2019).

Although the prophage genome is very distinct from the core genome as measured by GC content, which is further accentuated by the unique biology of the phage virus, e.g., high recombination, evolutionary changes in the core genome can be matched to changes in the prophage genome suggesting a common but unidentified signaling mechanism contributing to a synchronized evolutionary change in the disparate parts of the bacterial chromosome. The detection of synchronized changes in both genomic fractions, core and accessory, may indicate a considerable period of co-existence in the strain. The observations in this study of identical single-nucleotide mutations in prophages (Entero_p88 and Burkho_phiE25) belonging to the green and red clusters appear to support our inference of the cluster relatedness however the contrast is observed for the third common prophage (Salmon_SEN34) in which identical SNP was observed in the red and blue clusters, indicating that caution may be warranted not to over-interpret the mutational changes in the prophages as they may not always infer a direct relationship.

In conclusion, this study explored the diversity of prophages, the main component of the accessory genomes in Salmonella. Our results indicate that S. Adjame prophages collectively displayed a considerable degree of variation that allowed each and every strain in our cohort to be characterized at a very high resolution in a manner that allowed each to be distinguished from one another. Surprisingly, the pattern of variation observed in one of the most prevalent prophages, Edward_GF, accurately reproduced the clustering pattern seen by core-genome analysis, indicating that correlative mutational changes may be found in the core and accessory genomes. The routine use of core-genome analysis for characterizing Salmonella stands to benefit from the application of our highly discriminatory PST pipeline to develop a very comprehensive definition of genomic attributes and identify similarities or dissimilarities among isolates. Prophages as well as their patterns of variation can shed light on strain relatedness and niche characteristics even when these changes are not imprinted on the core genome.

MAC is affiliated to the National Institute for Health Research Health Protection Research Unit (NIHR HPRU) in Genomics and Enabling Data at University of Warwick in partnership with the UK Health Security Agency (UKHSA). She is based at UKHSA. The views expressed are those of the author(s) and not necessarily those of the NIHR, the Department of Health and Social Care or the UK Health Security Agency. LG is the Leung Family Professor of Food Safety, Department of Food Science, University of Guelph, Guelph, Canada. DO is an Adjunct Professor, Department of Food Science, University of Guelph., Guelph Canada. RG is currently employed at the National Microbiology Laboratory, Public Health Agency of Canada.

The original contributions presented in the study are included in the article/Supplementary material, further inquiries can be directed to the corresponding author. Raw read genome sequences are available in the Short Read Archivs of the National Centre for Biotechnology Information (NCBI) database using the accension number and SRA identification provided in Table 1.

RG analyzed the core-genome and prophage sequences, prepared the figures and tables, and drafted the manuscript. M-OD analyzed the prophages and wrote the scripts. MC participated in study design, provided background data, and reviewed manuscript. LG contributed to study conception and reviewed the manuscript. DO conceived the study and participated in the data analysis, and in manuscript development and review. All the authors revised and approved the manuscript.

This study was funded by grants from Genome Canada to the Salmonella Syst-OMICS project, Genome Research and Development Initiative of the Government of Canada, Canadian Security and Science Program of the Department of National Defense, and the Canadian Food Inspection Agency.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.