Distinct Salmonella Enteritidis lineages associated with enterocolitis in high-income settings and invasive disease in low-income settings.

An epidemiological paradox surrounds Salmonella enterica serovar Enteritidis. In high-income settings, it has been responsible for an epidemic of poultry-associated, self-limiting enterocolitis, whereas in sub-Saharan Africa it is a major cause of invasive nontyphoidal Salmonella disease, associated with high case fatality. By whole-genome sequence analysis of 675 isolates of S. Enteritidis from 45 countries, we show the existence of a global epidemic clade and two new clades of S. Enteritidis that are geographically restricted to distinct regions of Africa. The African isolates display genomic degradation, a novel prophage repertoire, and an expanded multidrug resistance plasmid. S. Enteritidis is a further example of a Salmonella serotype that displays niche plasticity, with distinct clades that enable it to become a prominent cause of gastroenteritis in association with the industrial production of eggs and of multidrug-resistant, bloodstream-invasive infection in Africa.

Introduction

Salmonella enterica serovar Enteritidis (hereafter referred to as S. Enteritidis) has been a global cause of major epidemics of enterocolitis, which have been strongly associated with intensive poultry farming and egg production1. The serovar is usually considered to be a generalist in terms of host range and has a low human invasiveness index, typically causing self-limiting enterocolitis2. Following a number of interventions in the farming industry involving both improved hygiene and poultry vaccination, epidemic S. Enteritidis has been in decline in many countries including the United Kingdom and USA3,4. S. Enteritidis has also been used extensively since the early 1900s as a rodenticide (named the “Danysz virus”), following development at Institut Pasteur, France. Although by the 1960s, Salmonella-based rodenticides had been banned in the US, Germany and the UK, S. Enteritidis is still produced as a rodenticide in Cuba, under the name Biorat®5.

Serovars of Salmonella that cause enterocolitis in industrialised settings are strongly associated with life-threatening invasive nontyphoidal Salmonella (iNTS) disease in sub-Saharan Africa (SSA). S. Enteritidis and Salmonella enterica serovar Typhimurium (S. Typhimurium) are the two leading causes of iNTS disease in SSA6 and both are associated with multidrug resistance (MDR)7. The clinical syndrome iNTS disease is associated with immunosuppression in the human host, particularly malnutrition, severe malaria and advanced HIV in young children and advanced HIV in adults8. It has been estimated to cause 681,000 deaths per year9.

Salmonella is a key example of a bacterial genus in which there is a recognizable genomic signature that distinguishes between a gastrointestinal and an extra-intestinal/invasive lifestyle10, whereby functions required for escalating growth in an inflamed gut are lost when the lineage becomes invasive11. In order to investigate whether there were distinct bacterial characteristics explaining the very different epidemiological and clinical profile of epidemic isolates of serotype S. Typhimurium from SSA and industrialised settings, whole-genome sequence (WGS) investigations of this serovar were previously undertaken. These revealed a novel pathotype of multilocus sequence type (MLST) ST313 from SSA, which differed from clades that cause enterocolitis in industrialised settings, by showing patterns of genomic degradation potentially associated with more invasive disease and differential host adaptation12–17.

In relation to S. Enteritidis, there is a growing body of literature on the evolutionary history, phylogeny and utility of WGS for surveillance of S. Enteritidis outbreaks18–20. The broadest study of the phylogeny to date revealed five major lineages, but contained only two African isolates21. There have also been limited reports of isolates of S. Enteritidis from African patients living in Europe that are MDR and which display a distinct phage type (PT 42)22,23. We therefore hypothesized that there are distinct lineages of S. Enteritidis circulating in both the industrialised and developing world with different origins, likely distinct routes of spread and that are associated with different patterns of disease, which will display the distinct genomic signatures characteristic of differential adaptation. To investigate this we have collected a highly diverse global collection of S. Enteritidis isolates and compared them using whole-genome sequencing, the highest possible resolution typing methodology.

Results

Isolate collection

In total, 675 isolates of S. Enteritidis isolated between 1948 and 2013 were sequenced. The collection originated from 45 countries and six continents (Table 1). 496/675 isolates were from Africa, with 131 from the Republic of South Africa (RSA), a further 353 from the rest of SSA, and 12 from North Africa (Table 1). There were 343 isolates from normally sterile human sites (invasive), 124 non-invasive human isolates (predominantly stool samples) and 40 from animal, food or environmental sources. The full metadata are described in Supplementary Table 1 and have been uploaded to the publically available database Enterobase.

Table 1

Summary of metadata (n) by region in numbers

Region	Total	Site of isolation			Antimicrobial resistance phenotype
		Human Invasive	Human non-invasive	Food/Animal/Environment	Drug susceptible	Resistant to 1-2 1st line	MDR*	Fluoroquinolone	ESBL†
Asia	11	5	5	1	0	0	0	0	0
Europe	61	0	16	24	2	0	0	0	0
South America	27‡	3	6	7	8	0	0	0	0
North Africa	12	9	1	1	9	0	0	2	0
Sub-Saharan Africa	353	269	22	7	99	64	149	0	3
Republic of South Africa	131	57	74	0	83	44	4	0	0

Multidrug resistant: resistant to ≥3 antimicrobials

Extended spectrum beta lactamase producing

Uruguay strains previously characterised by Betancor 44

Phylogeny

675 S. Enteritidis genomes and one Salmonella enterica serovar Gallinarum were mapped to the S. Enteritidis strain P125109 reference sequence, variable regions excluded and the remaining sites were screened for single nucleotide polymorphisms (SNPs). This left an alignment containing a total of 42,373 variable sites, from which a maximum likelihood (ML)-phylogeny was constructed using S. Gallinarum, which is a closely related serovar, as an out-group (Figure 1). HierBAPS was run over two rounds, which provided clear distinction between clades/clusters24. The phylogeny of S. Enteritidis revealed evidence of three clades associated with epidemics, one which we have termed the ‘global epidemic clade’ and includes the reference PT4 isolate P125109 and two African clades: one predominantly composed of West African isolates (labeled the ‘West African clade’) and a second composed of isolates predominantly originating in Central and Eastern Africa, called the ‘Central/Eastern African clade’). Figure 1 also shows the other clades and clusters predicted by HierBAPS, the largest of which is a paraphyletic cluster from which the global epidemic clade emerged (Outlier Cluster in Figure 1), and a further five smaller clades or clusters predicted by HierBAPS.

Figure 1

Maximum likelihood phylogeny of S. Enteritidis based on 675 isolates rooted to S. Gallinarum. There are 3 epidemic clades; 2 African epidemic clades and a global epidemic clade. Scale bar indicates nucleotide substitutions per site.

The global epidemic clade contains isolates of multiple phage types, including 4 and 1, which have been linked to the global epidemic of poultry associated human enterocolitis25. It comprised 250 isolates from 28 countries, including 43 from Malawi and 82 from RSA. They were isolated from across a 63-year period (1948-2013). Antimicrobial susceptibility testing had been performed on 144 isolates and 104 were susceptible to all antimicrobials tested, five were multidrug resistant (MDR: resistant to 3 or more antimicrobial classes), one was nalidixic acid resistant and none were extended-spectrum beta-lactamase (ESBL)-producing isolates. Database comparison of the genomes from this clade revealed that 221 (88%) of them contained no predicted antimicrobial resistance (AMR) genes apart from the cryptic resistance gene aac(6')-Iy26.

The global epidemic clade has emerged from a diverse cluster previously described by Zheng27, which encompassed 131 isolates (Figure 1: ‘Outlier Cluster’). In addition to being paraphyletic, this group was geographically and temporally diverse, and predominantly drug susceptible (59/71 isolates). Whilst the majority of the diversity of phage typed isolates was contained within the global epidemic clade, this cluster alone contained isolates of phage type 14b, which was recently associated with a multi-country outbreak of S. Enteritidis enterocolitis in Europe associated with chicken eggs from Germany28. There were also 41 isolates from RSA in this clade, where it has been a common cause of bloodstream infection, and 39 bloodstream isolates from Malawi. Database comparison of the genomes from this clade revealed that 122 (82%) of these genomes contained no predicted AMR genes apart from the cryptic resistance gene aac(6')-Iy.

There were two related, but phylogenetically and geographically distinct, epidemic clades that largely originated from SSA. The Central/Eastern African clade included 166 isolates, all but two of which (from RSA) came from this region. Of these, 126/155 (82%) were MDR and 148/153 (97%) displayed phenotypic resistance to between one and four antimicrobial classes. All of these genomes contained at least five predicted resistance genes and 128 (77%) contained nine (Table 2 and Supplementary Table 2). 155/165 (94%) of these isolates were cultured from a normally sterile compartment of a human (i.e. blood or cerebrospinal fluid) and were considered to be causing invasive disease (Table 2). The second African epidemic clade was significantly associated with West Africa with 65/66 isolates coming from this region and one isolate from USA. This clade was also associated with drug resistance (62 [94%] resistant to ≥1 antimicrobial class by phenotype and genotype) and human invasive disease (61 [92%]). It also included two isolates that were subtyped as phage type 4.

Table 2

Metadata summarised by clade

Major Clade/cluster	Site of isolation N (%)				Number (%) of antimicrobial resistance genes*
	Human Invasive	Human non-invasive	Food/Animal/Environment	Unknown	1-3	4-6	7-9
West African	61 (92)	1 (2)	0 (0)	4 (6)	22 (33)	9 (14)	35 (66)
Central/Eastern African	155 (93)	7 (4)	0 (0)	5 (3)	0 (0)	11 (7)	156 (93)
Global epidemic	94 (38)	95 (38)	31 (12)	30 (12)	243 (97)	7 (3)	0 (0)
Outlier cluster	51 (38)	36 (27)	27 (20)	20 (15)	128 (96)	3 (2)	3 (2)

All isolates contained cryptic aminoglycoside acetyltransferase gene aac(6’)-ly26

The remaining 58 isolates included in this study were extremely diverse, phylogenetically, temporally and geographically. Only two displayed any phenotypic AMR, one of which was MDR. Inspection of the genome revealed that five had predicted AMR genes in addition to aac(6')-Iy, four of which were isolated in sub-Saharan Africa. Twenty were associated with invasive human disease, and six were recovered from stool. Three isolates were from stocks of rodenticide and these were phylogenetically remote from both global-epidemic and the two African epidemic clades.

To add further context to these findings we screened the entire publically available Public Health England (PHE) sequenced Salmonella routine surveillance collection, which includes 2,986 S. Enteritidis genomes, 265 of which were associated with travel to Africa (Supplementary Figure 1). Within this huge collection, including 61 (2.0%) bloodstream isolates and 2670 (89.4%) stool isolates, only 6 isolates (4 from blood culture, 1 from stool) fell within to the West African clade and 1 (from stool) belonged to the Central/Eastern African clade. Notably, these isolates were all either associated with travel to Africa and/or taken from patients of African origin.

It is apparent from the location of the archetypal reference isolate and archetypal phage types in the phylogeny (Supplementary Figure 2) that the majority of S. Enteritidis studied previously belonged to the global epidemic clade associated with enterocolitis in industrialised countries. Furthermore, its also clear that two additional, previously unrecognized S. Enteritidis lineages have emerged, largely restricted to Africa, that are strongly associated with MDR and invasive disease.

To understand how recently these Africa-associated lineages emerged we used Bayesian Evolutionary Analysis by Sampling Trees (BEAST) to reconstruct the temporal history of the epidemic clades29. These data (Supplementary Figure 3) estimate the most recent common ancestor (MRCA) of the Central/Eastern African clade dates to 1945 (95% Credible Interval [CrI]: 1924-1951) and for the West African clade it was 1933 (95% CrI: 1901-1956). We estimate the MRCA of the global epidemic clade originated around 1918 (95% CrI: 1879-1942 – Supplementary Figure 4), with a modern expansion occurring in 1976 (95% CrI: 1968-1983), whereas the paraphyletic cluster from which it emerged dates to approximately 1711 (95% CrI: 1420-1868).

Contribution of the accessory genome

Prophages have the potential to carry non-essential "cargo" genes, which suggests they confer a level of specialization to their host bacterial species, whilst plasmids may confer a diverse array of virulence factors and AMR 30,31. Therefore it is critical to evaluate the accessory genome in parallel with the core. 622 sequenced genomes were used to determine a pangenome, which yielded a core genome comprising 4,076 predicted genes present in ≥90% isolates, including all 12 recognised Salmonella Pathogenicity Islands as well as all 13 fimbrial operons found in the P125109 reference32. The core gene definition was set to minimize stochastic loss of genes from the core due to errors in individual assemblies across such a large dataset. The accessory genome consisted of 14,015 predicted genes. Of the accessory genes, 324 were highly conserved across the global and two African epidemic clades, as well as the outlier cluster. Almost all were associated with the acquisition or loss of mobile genetic elements (MGEs) such as prophage or plasmids. Prophage regions have been shown to be stable in Salmonella genomes and are potential molecular markers, the presence of which has previously been used to distinguish specific clades13,33.

The lineage-specific whole gene differences of the major clades are summarized in Figure 2 and plotted against the representatives of the four major clades in Supplementary Figure 5. The lineage specific sequence regions include 57 predicted genes found to be unique to the global epidemic clade (Figure 2), all of which were associated with prophage ΦSE20, a region shown to be essential for invasion of chicken ova and mice in one previous study34. There were a further 39 genes conserved in the global epidemic and the paraphyletic outlying cluster, which were absent from both African clades, 26 of which correspond to region of difference (ROD) 2132. The Central/Eastern Africa clade contained 77 predicted genes that were absent in the other clades. 33 were associated with the virulence plasmid and a further 40 chromosomal genes were associated with a novel, Fels-2 like prophage region (ϕfels-BT). The West African clade had only 15 distinct predicted genes, 11 of which were plasmid-associated. The two African clades shared a further 102 genes: 48, including a leucine-rich repeat region, were associated with a novel prophage region closely related to Enterobacter phage P88, 44 were associated with a Gifsy-1 prophage found in S. Bovismorbificans and eight were associated with a Gifsy-2 prophage which has degenerated in the reference P125109.

Figure 2

Differences in accessory genomes of 4 major clades. Approximate position of prophages in chromosome is depicted, although prophages are not drawn to scale

The S. Enteritidis plasmid is the smallest of the generic Salmonella virulence plasmids at 58 kb and is unusual in that it contains an incomplete set of tra genes that are responsible for conjugative gene transfer. The phylogeny of the S. Enteritidis virulence plasmid backbone was reconstructed using reads that mapped to the S. Enteritidis reference virulence plasmid, pSENV. 120/675 (18%) genomes lacked pSENV. The virulence plasmid phylogeny is similar to that of the chromosome, suggesting that they have been stably maintained by each lineage and diversified alongside them (Supplementary Figure 6).

The virulence plasmids from the African clades were much larger than those held in the other clades at ~90 kb. A representative example was extracted from Malawian isolate D7795, sequenced using long read technology to accurately reconstruct it (PacBio; see methods) and denoted pSEN-BT (Accession number LN879484). pSEN-BT is composed of a backbone of pSENV with additional regions that are highly similar to recently sequenced fragments of an novel S. Enteritidis virulence plasmid (pUO-SeVR) isolated from an African patient presenting with MDR invasive S. Enteritidis in Spain22. Plasmid pSEN-BT harbours nine AMR genes (full list in Supplementary Table 2), plus additional genes associated with virulence and a toxin/antitoxin plasmid addiction system. Of note, plasmids from the West African isolates carry resistance gene chloramphenicol acetyl transferase A1 (catA1), whereas the Central/Eastern African strains carry catA2 and tetracycline resistance gene tet(A). Like pSENV, the African virulence plasmid contained an incomplete set of tra genes and so is not self-transmissible. This was confirmed by conjugation experiments and is consistent with previous reports22,23. These observations suggest that the evolution of the S. Enteritidis plasmid mirrors that of the chromosome; it is thus not a ‘novel’ plasmid, but in different SSA locations has acquired different AMR genes.

Multiple signatures of differential host adaptation

It has been observed in multiple serovars of Salmonella including S. Typhi, S. Gallinarum and S. Typhimurium ST313 that the degradation of genes necessary for the utilization of inflammation-derived nutrients is a marker of that lineage having moved from an intestinal to a more invasive lifestyle13,14,32,35. Accordingly, we have looked for similar evidence within a representative example of a MDR, invasive, Central/Eastern African clade isolate, D7795, that was isolated from the blood of a Malawian child in 2000. The draft genome sequence of D7795 closely resembles that of P125109, however, in addition to the novel prophage repertoire and plasmid genes described above, it harbours a number of predicted pseudogenes or hypothetically disrupted genes (HDGs)11.

In total, there were 42 putative HDGs in D7795, many of which are found in genes involved in gut colonisation and fecal shedding as well as various metabolic processes such as cobalamine biosynthesis which is a cofactor for anaerobic catabolism of inflammation-derived nutrients, such as ethanolamine, following infection36. Curation of the SNPs and insertions or deletions (indels) predicted to be responsible for pseudogenisation across the Central/Eastern African clade and West African clade revealed 37/42 predicted HDGs were fixed in other representatives of the Central/East African clade, with 27 of them being present in over 90% of isolates from that clade. Relatively fewer HDGs in D7795 (19/42) were present in representatives of the West African clade, although 13 were present in ≥90% of isolates (Supplementary Table 3).

In addition to this evidence of reductive evolution in D7795, there were 363 genes containing non-synonymous (NS)-SNPs, which change the amino acid sequence and so may have functional consequences37. The two African clades were screened for the presence of these NS-SNPs and 131 were found to be present and completely conserved across both clades, including NS-SNPs in 43 genes encoding predicted membrane proteins, 36 metabolic genes and 23 conserved hypothetical genes (Supplementary Table 4). Furthermore many of these NS-SNPs fall in genes within the same metabolic pathways as the HDGs (see Supplementary Note for detailed description). Supplementary Table 5 provides a list of some of the common traits identified amongst the functions of genes lost independently by D7795, S. Typhi and S. Gallinarum. The disproportionate clustering of mutations in membrane structures observed in the African clades is yet another sign of differential host adaptation analogous to that reported in both S. Typhi35 and S. Gallinarum32.

Biolog™ growth substrate platform profiling

The Biolog™ platform was utilized to generate a substrate growth utilisation profile for selected S. Enteritidis isolates. Corresponding signal values of replicate pairs of a Central/Eastern African isolate (D7795) and a global epidemic isolate (A1636) were compared using principal component analysis and found to be highly consistent. In total, 80 metabolites showed evidence of differential metabolic activity (Figure 3). Evaluation of data from the Central/Eastern African isolate using Pathway Tools software revealed that 14/27 (52%) of pathways with evidence of decreased metabolic activity at 28°C had a corresponding component of genomic degradation. This was also true for 12/30 (40%) of pathways with evidence of decreased metabolic activity at 37°C.

Figure 3

Heat map revealing changes in metabolic activity of Central/Eastern African clade isolate D7795 when compared to global epidemic isolate A1636 at 28 and 37°C. The figure also displays whether there are corresponding mutations in genes related to the affected metabolic pathway. (NSSNP=Non-synonymous single nucleotide polymorphism, HDG = Hypothetically disrupted gene)

Instances of reduced metabolic activity in a Central/Eastern African strain (D7795) compared to a global epidemic strain (A1636) included dulcitol and glycolic acid in the glycerol degradation pathway, propionic acid in the propanediol pathway and ethylamine and ethanolamine. These are all vitamin B12 (cobalamin) dependent reactions, for which there was a corresponding signature of genomic degradation. Also there was reduced activity in response to three forms of butyric acid, alloxan and allantoic acid metabolism. Allantoin can be found in the serum of birds, but not humans and is utilised as a carbon source during S. Enteritidis infection of chickens38, and HDGs relating to allantoin have been noted in S. Typhimurium ST31313. The full list of differences is detailed in Supplementary Table 6 and 7. This is a further sign of decreased metabolism of the Central/Eastern African isolate in the anaerobic environment of the gut.

Chick infection model suggests divide in host range

Given the phenotypic differences observed in the genotypically distinct global and African clades, we hypothesized that these lineages could have differing infection phenotypes in an in vivo challenge model. We compared the infection profile of a member of the Central/Eastern African clade (D7795) to the reference global epidemic strain P125109 in an avian host. The chicken group infected with P125109 showed mild hepatosplenomegaly consistent with infection by this Salmonella serovar and cecal colonization (Figure 4A-C). In contrast, the Central/Eastern African strain displayed significantly reduced invasion at 7 dpi of both liver (p=0.027) and spleen (p=0.007), however cecal colonization was not significantly reduced (p=0.160). This is in marked contrast to the behavior of S. Typhimurium ST313, which is more invasive in a chick infection model12.

Figure 4

Salmonella isolation from a chick infection model demonstrates failure of Central/Eastern African clade isolate to invade chicken spleen (4A) and liver (4B) or to colonize chicken caeca (4C) at 7 days post infection (dpi) (n=24 at this time point) compared to the global epidemic clade. Numbers are expressed as colony forming units (CFU) per gram of tissue

Discussion

S. Enteritidis is an example of a successful Salmonella lineage with the apparent ability to adapt to different hosts and transmission niches as and when opportunities for specialization have presented. Langridge et al recently evaluated the Enteritidis/Gallinarum/Dublin lineage of Salmonella, revealing components of the nature and order of events associated with host-range and restriction39. In the present study, we have highlighted the plasticity of S. Enteritidis, providing evidence of three distinct epidemics of human disease. In addition we show multiple additional clades and clusters that demonstrate the huge reservoir of diversity amongst S. Enteritidis from which future epidemics might emerge.

An important question posed by this study is why have distinct clades of Salmonella emerged to become prominent causes of iNTS disease in Africa, from a serotype normally considered to be weakly invasive? The presence of a highly immunosuppressed population due to the HIV pandemic is clearly a key host factor that facilitates the clinical syndrome iNTS disease40,41. In addition to human host factors, there are two distinct African epidemic lineages that have emerged in the last 90 years. Both lineages are significantly associated with a novel prophage repertoire, an expanded, MDR-augmented virulence plasmid, and patterns of genomic degradation with similarity to other host-restricted invasive Salmonella serotypes including S. Typhi and S. Gallinarum and to clades of S. Typhimurium associated with invasive disease in Africa13,32,35. This pattern of genomic degradation is concentrated in pathways specifically associated with an enteric lifestyle, however it is noteworthy that in the chick infection model, the African S. Enteritidis invaded the chick liver and spleen less well than the global pandemic clade. This raises the possibility that the two clades occupy different ecological niches outside the human host or that they behave differently within the human host and screening of the huge S. Enteritidis collection from routine Salmonella surveillance by PHE supports the assertion that these lineages are geographically restricted to Africa. This study therefore indicates a need to understand what these ecological niches might be, and then to define the transmission pathways of African clades of S. Enteritidis, in order to facilitate public health interventions to prevent iNTS disease.

The evolution of the S. Enteritidis virulence plasmid is intriguing; pSENV is the smallest of the known Salmonella virulence-associated plasmids, but in SSA, the plasmid has nearly doubled in size partly through the acquisition of AMR genes. The absence of tra genes necessary for conjugal transfer either indicates that MDR status has evolved through acquisition of MGEs multiple times or through clonal expansion and vertical transmission of the plasmid to progeny. The available data suggest that the former scenario has happened twice, once in West Africa, and once in Central/Eastern Africa.

Despite S. Enteritidis being reported as a common cause of bloodstream infection (BSI) in Africa6,7 the Global Enteric Multicenter Study (GEMS) found that Salmonella serotypes were an uncommon cause of moderate to severe diarrhoea in African children less than 5-years of age42. Our data associating the African lineages with invasive disease is also consistent with data presented in a recent Kenyan study comparing a limited number and diversity of S. Enteritidis isolates from blood and stool. Applying the lineages defined in this study to the genome data reported from Kenya showed that 20.4% of isolates from that study belonging to the global clade were associated with invasive disease, whereas 63.2% of the isolates in that study belonging to our Central/Eastern African clade were associated with invasive disease43. The remaining isolates were associated with cases of enterocolitis or asymptomatic carriage, confirming that the Central/Eastern African clade can also cause enterocolitis. The association of S. Enteritidis clades circulating in sub-Saharan Africa with iNTS disease may reflect the fact that their geographical distribution permits them to behave as opportunistic invasive pathogens in a setting where advanced immunosuppressive disease is highly prevalent in human populations.

In summary, two clades of S. Enteritidis have emerged in Africa, which have different phenotypes and genotypes to the strains of S. Enteritidis circulating in the industrial world. These strains display evidence of changing host adaptation, different virulence determinants and multi-drug resistance, a parallel situation to the evolutionary history of S. Typhimurium ST313. They may have different ecologies and/or host ranges to global strains and have caused epidemics of BSI in at least three countries in SSA, yet are rarely responsible for disease in South Africa. An investigation into the environmental reservoirs and transmission of these pathogens is warranted and urgently required.

Online Methods

Bacterial Isolates

S. Enteritidis isolates were selected on the basis of six factors; date of original isolation, antimicrobial susceptibility pattern, geographic site of original isolation, source (human [invasive vs stool], animal or environmental), phage type (where available), and multilocus variable number tandem repeat (MLVA) type (where available). S. Enteritidis P125109 (EMBL accession no. AM933172) isolated from a poultry farm from the UK was used as a reference32. The full metadata are in Supplementary Table 1. Isolates have been attributed to region according to United Nations statistical divisions.

Sequencing, SNP-calling, construction of phylogeny and comparative genomics

PCR libraries were prepared from 500 ng of DNA as previously described45. Isolates were sequenced using Illumina GA II, HiSeq 2000 and MiSeq machines (Illumina, San Diego, CA, USA) and 150 bp paired-end reads were generated. The strains were aligned to Salmonella Enteritidis reference genome P125109 using a pipeline developed in-house at the Wellcome Trust Sanger Institute (WTSI). For each isolate sequenced, the raw sequence read pairs were split to reduce the overall memory usage and allow reads to be aligned using more than one CPU. The reads were then aligned using SMALT, a hashing based sequence aligner. The aligned and unmapped reads were combined into a single BAM file. Picard was used to identify and flag optical duplicates generated during the making of a standard Illumina library, which reduces possible effects of PCR bias. All of the alignments were created in a standardized manner, with the commands and parameters stored in the header of each BAM file, allowing for the results to be easily reproduced.

The combined BAM file for each isolate was used as input data in the SAMtools mpileup program to call SNPs and small indels, producing a BCF file describing all of the variant base positions46. A pseudo-genome was constructed by substituting the base call at each variant or non-variant site, defined in the BCF file, in the reference genome. Only base calls with a depth of coverage >4 or quality >50 were considered in this analysis. Base calls in the BCF file failing this quality control filter were replaced with the “N” character in the pseudo-genome sequence.

All of the software developed is freely available for download from GitHub under an open source license, GNU GPL 3.

Phylogenetic modelling was based on the assumption of a single common ancestor, therefore variable regions where horizontal genetic transfer occurs were excluded47,48. A maximum likelihood (ML) phylogenetic tree was then built from the alignments of the isolates using RAxML (version 7.0.4) using a GTR+I+G model49. The maximum-likelihood phylogeny was supported by 100 bootstrap pseudo-replicate analyses of the alignment data. Clades were predicted using Hierarchical Bayesian Analysis of Population Structure (HierBAPS)24. This process was repeated to construct the plasmid phylogeny, using reads that aligned to pSENV.

To ascertain the presence of the clusters defined by HierBAPs in the Public Health England (PHE) routine Salmonella surveillance collection, seventeen isolates representing the diversity of the collection were compared against 2986 S. Enteritidis PHE genomes. Single linkage SNP clustering was performed as previously described50. A maximum-likelihood phylogeny showing the integration of the seventeen isolates with 50-SNP cluster representatives of the PHE S. Enteritidis collection was constructed as above. FASTQ reads from all PHE sequences in this study can be found at the PHE Pathogens BioProject at the National Center for Biotechnology Information (Accession PRJNA248792).

Temporal reconstruction was performed using Bayesian Evolutionary Analysis Sampling Trees (BEAST version 1.8.2)51. A relaxed lognormal clock model was initially employed. The results of this model indicated that a constant clock model was not appropriate, as the posterior of the standard deviation of the clock rate did not include zero. A range of biologically plausible population models (constant, exponential and skyline) was investigated. Skyline models can be biased by non-uniform sampling and we observed a strong similarity between reconstructed skyline population and the histogram of sampling dates and so this model was excluded. The exponential models consistently failed to converge and were excluded. Thus, for all datasets, lognormal clock and constant population size models were used. The computational expense required for this analysis precluded running estimators for model selection. However, we note that Deng et al used the same models in their analysis of 125 S. Enteritidis isolates. Default priors were used except for ucld.mean, Gamma(0.001,1000), initial: 0.0001; exponential.popSize, LogNormal(10,1.5), initial: 121.

Three chains of 100 million states were run in parallel for each clade of the four major HierBAPS clades, as well as a fourth chain without genomic data to examine the influence of the prior, which in all cases was uninformative. The final results, as used here, all had effective sample sizes (ESS) of over 200 and had convergence between all three runs. For the Global and Global Outlier lineages, the datasets were not computationally feasible to analyse. We thus created 3 further random subsets of the data by drawing n isolates from each sampled year where n was sampled from a Poisson distribution where λ=2. The posteriors of all subsets were extremely similar and runs were combined to produce the final most recent common ancestor (MRCA) estimates.

In order to gain a detailed insight into genomic differences, a single high quality sequence from Malawian S. Enteritidis isolate D7795 was aligned against the P125109 using ABACAS and annotated52. Differences were manually curated against the reference using the Artemis Comparison Tool (ACT)53. Sections of contigs which were incorporated into the alignment, but which did not align with P125109 were manually inspected and compared to the public databases using BLASTn. When these regions appeared to be novel prophages, they were annotated using the phage search tool PHAST and manually curated54. In order to investigate whether the SNPs and/or indels that were predicted to be responsible for pseudogene formation in D7795 were distinct to that isolate or conserved across both African epidemic clades, all isolates were aligned to P12509 and the relevant SNPs/indels investigated using in-silico PCR of the aligned sequences. Manual curation was performed to confirm the nature of all pseudogene-associated SNPs/indels. NS-SNPs identified in D7795 were sorted throughout the African clades by extracting and aligning the appropriate gene sequences from P125109 and D7795. The coordinates of the NS-SNPs were then used to identify the relevant sequence and determine the nature of the base.

Accessory genome

The pangenome for the dataset was predicted using ROARY 55. Genes were considered to be core to S. Enteritidis if present in ≥90% of isolates. A relaxed definition of core genome was used as assemblies were used to generate it and the more assemblies one uses, the more likely it is that a core gene will be missed in one sample due to an assembly error. The remaining genes were considered to be core to the clades/clusters predicted by HierBAPS if present in ≥75% if isolates from within each clade/cluster. These genes were then curated manually using ACT to search for their presence and position in P125109 or the improved draft assembly of representative isolates of each of the other clades if not present in P125109. Any large accessory regions identified were blasted against the assembled genomes of the entire collection to confirm they were grossly intact.

Plasmid identification

Plasmid DNA was extracted from isolate D7795 using the Kado & Liu method and separated by gel-electrophoresis alongside plasmids of known size, to estimate the number and size of plasmids present56. Plasmid conjugation was attempted by mixing 100 μL of overnight culture of donor and recipient strains (rifampicin resistant Escherichia coli C600) on Luria-Bertani agar plates and incubating overnight at 26°C and 37°C. The plasmid was sequenced using the PacBio platform to gain long reads and a single improved draft assembly, which was aligned against P125109 plasmid pSENV (Accession Number HG970000). For novel regions of the plasmid from isolate D7795, genes were predicted using GLIMMER and manual annotations applied based on homology searches against the public databases, using both BLASTn and FASTA. The plasmid phylogeny was reconstructed using the same methodology as the chromosome; a maximum likelihood (ML) phylogenetic tree was built from the alignments of the isolates using RAxML (version 7.0.4) using a GTR+I+G model

Identification of AMR genes

A manually curated version of the Resfinder database was used to investigate the isolates for the presence of AMR genes57. To reduce redundancy, the database was clustered using CD-HIT-EST58, with the alignment length of the shorter sequence required to be 90% the length of the longer sequence. All other options were left as the defaults. The representative gene of each cluster was then mapped with SMALT to the assemblies of each isolate to identify and matches with an identity of 90% or greater were considered significant, in line with the default clustering parameters of CD-HIT-EST. Where partial matches were identified at the ends of contigs, having an identity of 90% or greater to the matched region of the gene, potential AMR gene presence was recorded. To confirm presence of these partial matches, raw sequencing reads of the pertinent isolates were mapped using SMALT to these genes to check for 90% identity across the entire gene.

The Biolog™ platform enables the simultaneous quantitative measurement of a number of cellular phenotypes, and therefore the creation of a phenotypic profile of a variety of assay conditions59. Incubation and recording of phenotypic data were performed using an OmniLog® plate reader.

In these experiments, two replicates of D7795 were compared to two replicates of a PT4 like strain at 28°C and 37°C to represent environmental and human temperatures. Biolog™ plates PM1-4 and 9 (Carbon source [PM1,PM2], nitrogen source [PM3] and phosphor and sulphur source [PM4] metabolism and osmotic pressure [PM9]) were used. Each well was inoculated as described below, thereby testing 475 conditions at once (each plate has one negative control well).

The isolates were cultured overnight on LB-agar at 37°C in air to exclude contamination. Colonies were scraped off plates and dispensed into IF-0a solution (Biolog) to a cell density corresponding to 81% transmittance. For each plate used, 880 µL of this cell suspension was added to 10 mL IF-10b GP/GP solution (Biolog) and 120 µL dye mix G (Biolog). This was then supplemented with a 1 mL solution of 7.5 mM D-ribose (Sigma), 2 mM magnesium chloride, 1 mM calcium chloride, 2 mM sodium pyrophosphate (Sigma), 25 µM L-arginine (Sigma), 25 µM L-methionine (Sigma), 25 µM hypoxanthine (Sigma), 10 µM lipoamide (Sigma), 5 µM nicotine adenine dinucleotide (Sigma), 0.25 µM riboflavin (Sigma), 0.005% by mass yeast extract (Fluka) and 0.005% by mass Tween 80 (Sigma). 100µl of this mixture was dispensed into each well on the assay plate. Plates were then allowed to equilibrate in air for 5 min prior to being sealed in airtight bags and loaded into the Omnilog machine (Biolog). Plates were scanned every 15 min for 48 hours while incubated at 28°C and 37°C in air. Culture under anaerobic conditions was unavailable. Two paired replicates were performed for each of the two isolates.

After completion of the run, the signal data were compiled and analysed using the limma package in ‘R’ described previously60. A log-fold change of 0.5 controlling for a 5% false discovery rate was used as a cut-off for investigating a specific metabolite further using Pathway Tools61 and whether the metabolic change was related to pseudogenes and non synonymous(NS)-SNPs in genes in the respective genomes.

In vivo Infection Model

Two isolates were used in the animal models: S. Enteritidis P125109 and D7795. Unvaccinated commercial female egg-layer Lohmann Brown chicks (Domestic Fowl [Gallus gallus]) were obtained from a commercial hatchery and housed in secure floor pens at a temperature of 25°C. Eight chicks per strain per time point were inoculated by gavage at 10 days (d) of age and received a dose of ~108 Salmonella colony forming units (CFU) in a volume of 0.2 mL. Subsequently, four to five birds from each group were humanely killed at 3, 7 or 21 d post-infection (p.i.). At post mortem, the liver, spleen, and caecal contents were removed aseptically, homogenised, serially diluted and dispensed onto Brilliant Green agar (Oxoid) to quantify colony forming units (CFU) as described previously62. Statistical analysis was performed using SPSS, version 20 (IBM). Kruskal-Wallis was used to compare bacterial loads between infected groups.

All work was conducted in accordance with the UK legislation governing experimental animals, Animals (Scientific Procedures) Act 1986, under project licence 40/3652 and was approved by the University of Liverpool ethical review process prior to the award of the project license. The licensing procedure requires power calculations to determine minimal group sizes for each procedure to ensure results are significant. For these experiments a group size of 8 birds per time point was chosen, based on a variation in 1.0 log10 in bacterial count between groups as being significant along with prior experience of Salmonella infection studies. Groups were randomly selected on receipt from the hatchery and investigators conducting animal experiments were not blinded, as the current UK code of practice requires all cages or pens to be fully labeled with experimental details. No animals were excluded from the analysis. All animals were checked a minimum of twice daily to ensure their health and welfare.

Code availability & URLs

Software is referenced and URLS are provided below. All software is open source.

BEAST: http://beast.bio.ed.ac.uk/

Biolog™: http://www.biolog.com

BLASTn: http://blast.ncbi.nlm.nih.gov

limma package: www.bioconductor.org

PacBio platform: http://www.pacificbiosciences.com/

Picard: https://broadinstitute.github.io/picard

‘R’: www.R-project.org

SMALT: www.sanger.ac.uk/science/tools/smalt-0

United Nations statistical divisions: www.unstats.un.org/unsd/methods/m49/m49regin.htm