Introduction

Invasive nontyphoidal Salmonella (iNTS) infections are a major cause of global morbidity and mortality (1). Of the >2,500 NTS serovars, iNTS infections are associated with a relatively limited number of serovars. These include Salmonella enterica serovar Dublin (S. Dublin), S. Panama, S. Virchow, S. Choleraesuis and S. Typhimurium ST313 (2–4). S. Dublin, a bovine-adapted pathogen, is one of the leading iNTS serovars both in Australia and globally (5–7), and is associated with bacteraemia in a high proportion of infections (2,5–7). Symptomatic cattle may develop septicaemia and transplacental transmission (7,8), while asymptomatic infections in cattle enable persistence within the herd population, facilitating disease spread (7,9). S. Dublin causes zoonotic infections predominately through the consumption of contaminated dairy and meat produce (10–12), and as such is considered a ‘One Health’ pathogen. (10–12).

Antimicrobial resistance (AMR) in NTS has increased over the past several decades, leading to the reduction of therapeutic options (6,13). Although S. Dublin isolates are largely susceptible to therapeutics, AMR has become increasingly prevalent in S. Dublin (7,12,14–16). For example, a multidrug-resistant (MDR; defined as resistance to ≥ 3 antimicrobial classes (17)) lineage has recently been described in North America, characterised by varying profiles of AMR to ampicillin, streptomycin, chloramphenicol, sulphonamides, tetracyclines and third generation cephalosporins. These MDR S. Dublin isolates all type as sequence type 10 (ST10), and the AMR determinants have been demonstrated to be carried on an IncC plasmid that has recombined with a virulence plasmid encoding the spvRABCD operon with IncX1 and IncFII replicons, to produce a 329kb megaplasmid (12,16). S. Dublin is closely related to S. Enteritidis and S. Gallinarum (18), with differences in accessory genome content reflective of adaptation to different ecological niches (18). The majority of the STs reported for S. Dublin are ST10 or single-locus variants of this ST, although there is a second lineage, ST74, that is distinct from the main S. Dublin population (19). This uncommon S. Dublin lineage, ST74, may be variably serotyped as S. Dublin or S. Enteritidis dependent upon the strains used to generate typing sera in different laboratories (19), which may have also contributed to it not being identified. The antigenic formula of S. Dublin (1,9,12:g,p:-) is highly similar to S. Enteritidis (1,9,12:g,m:-), with three non-synonymous point mutations in fliC (Ala220Val, Thr315Ilc and Thr318Ala) differentiating the two serovars. Although the Vi flagella antigen (encoded by the viaB locus) is part of the antigenic definition for S. Dublin, it has been shown that the Vi antigen is absent in many S. Dublin isolates (19,20). Other putative virulence determinants in S. Dublin include Salmonella pathogenicity islands (SPIs) 6 and 19, which both encode type VI secretion system (T6SS) machineries (21,22). The presence of two T6SSs in S. Dublin contrasts with other NTS serovars, which encode a single T6SS (21,22). T6SSs are molecular nanomachines that deliver effector proteins into eukaryotic and/or prokaryotic cells to promote pathogen survival within multispecies and complex immune environments such as the gut. Recent studies have demonstrated a role for SPI-6 and SPI-19 in interbacterial competition by specific S. Dublin strains (23,24), however their effect on host immune processes remains largely unknown.

We previously undertook a genomic investigation of another emerging NTS lineage in Australia, Salmonella enterica serovar 4,[5],12:i:-, and identified several bacterial traits that may have contributed to the successful spread of this lineage (25–29). In particular, we hypothesised that flagellin deletion and MDR status may have contributed to host immune evasion and positive selection, respectively, in Salmonella enterica serovar 4,[5],12:i:-(30). Here, we sought to understand the contribution of virulence and resistance to the recent emergence of MDR S. Dublin. We interrogated a global dataset of 1,303 isolates collected from 13 countries on five continents and identified distinct geographical lineages associated with specific AMR and virulence profiles, including a lineage unique with increased capability for intracellular survival and immune evasion in human macrophages.

Results

Distinct S. Dublin populations circulate globally

A total of 1,303 S. Dublin genomes were included in our analysis, spanning four decades (Figure 1A, Supplementary Table 1-2). This dataset included 53 genomes from humans in Australia, and 1,250 genomes from Africa, Europe, North and South America collected from humans (n=462), animals (n=551), food (n=223) and other sources that were not stated (n=67) (Figure 1B, Supplementary Table 1). In total, 1,235 genomes (99.4%) were ST10 and 58 were a single locus variant (SLV) of ST10, with two double locus variant genomes (e.g. ST4100) restricted to isolates from South America. The exceptions were ST74 and SLV ST1545 that were only detected in Australia (Figure 1C-E, Supplementary Table 3). All isolates, including the ST74 lineage, were serotyped in silico as S. Dublin. However, the ST74 lineage and a single ST10 isolate had codon profiles in fliC of Gly60, Leu138, Ala220, Thr315 and Ala318, in contrast to the remaining ST10 and SLV, which both had a profile of Gly60, Leu138, Val220, Ile315 and Ala318 (the previously described ‘Du2’ and ‘Du1’ electrophoretic profiles (20), respectively) (Supplementary Table 4).

Figure 1.

A maximum likelihood (ML) phylogeny was inferred from the 1,303 genomes, revealing two major lineages of S. Dublin, ST74 and ST10 (Figure 1D, Supplementary Figure 1). Within the ST10 lineage, five clades were identified using Bayesian Analysis of Population Structure (BAPS), inferred from the alignment of 2,957 core single nucleotide polymorphisms (SNPs). Four clades were broadly associated with geographical region (Supplementary Figure 1, Supplementary Table 1). The Australian ST10 isolates (n = 39) clustered in clade 1, while isolates from North America were associated with clade 5 (Supplementary Figure 1). The BAPS clade 2 were the isolates that didn’t have sufficient evidence to be split into further clades, hence are spread across the tree spanning a range of geographical regions.

Emergence of S. Dublin lineages in the 20^th century

To investigate the emergence of S. Dublin clades, Bayesian phylodynamic analysis was undertaken on a representative subset of 660 genomes, selected for phylogenetic, geographical and AMR diversity (see methods). We inferred a maximum clade credibility (MCC) phylogeny from an alignment of 3,971 SNPs using the highest supported model of a relaxed lognormal molecular clock and coalescent exponential tree prior. The estimated substitution rate was 1.05 x 10^-7 (95% highest posterior density (HPD): 1.16 x 10^-7– 9.41 x 10^-8), slightly higher than previous estimated mean substitution rates for host-restricted Salmonella serovars (S. Typhi and S. Paratyphi A) but lower than rates reported for host-generalist NTS serovars S. Kentucky and S. Agona (31).

The topology of the MCC tree was consistent with the ML tree, with geographically associated sub-clades were observed in the MCC tree (Figure 2, interactive tree is available in Microreact under https://microreact.org/project/sdublinpaper). The most recent common ancestor (MRCA) for the S. Dublin global population was inferred to be 1901 (HPD 1878-1924), with the MRCA of the different clades and sub-clades largely falling between 1950 and 2000. For example, the most recent common ancestor (MRCA) of clade 5 MDR isolates from North America was estimated to be ∼1964 (HPD 1952 - 1974) (Figure 2). The MRCA of the Australian sub-clade was estimated to be ∼1967 (HPD 1953-1980); with the Australian and North American isolates inferred to have a MRCA ∼1931 (HPD 1914-1947). Moreover, three clades (1, 3, and 4) had European sub-clades associated with Sweden, the UK, and Denmark, respectively that were all largely susceptible to antimicrobials. Some regions, such as Europe and South America, were associated with several sub-clades of S. Dublin, potentially suggesting multiple importations and expansions, while other regions, such as Australia were comprised of a single sub-clade within clade 1.

Figure 2.

Resistance determinants are associated with geographical regions

Different plasmid types co-occurred with different AMR and HMR determinants in the S. Dublin population, and these were also associated with differential geographical distribution (Figure 2, Supplementary Figure 3-4, Supplementary Table 1). Overall, 51.5% (667/1,295) of S. Dublin genomes had at least one AMR determinant, with 632/667 (94.8%) isolates classified as MDR. MDR S. Dublin was predominantly identified from North American isolates, although some AMR determinants were detected in S. Dublin from other geographic regions (Figure 3, Supplementary Figures 2-3). AMR S. Dublin collected from animals was common (n=305/667; 45.7%). AMR in S. Dublin from food sources and humans was similar in proportion (n=174/667; 26.1% and n=173/667; 25.9% respectively). Although these profiles are likely impacted by the sampling frame of different studies, the source of collection was not considered as a selection criterion for this study, and as such, may broadly represent the AMR profiles of S. Dublin from different sources (Figure 3A-B, Supplementary Figure 1A-C).

Figure 3.

The most prevalent AMR determinants were strA-B (n=620), floR (n=556), tet(A) (n=608) and sul2 (n=596) (Supplementary Figure 2). AMR mechanisms to ciprofloxacin (triple point mutations in quinolone resistance determining regions gyrA[83:S-F], gyrA[87:D-N], parC[80:S-I]), macrolides (mph(A)), and colistin (mcr-9.1) were each detected in single genomes from S. Dublin clade 5 at a low frequency of 0.08% (n=1/1,303) in the dataset. The only third generation cephalosporin (3GC) resistance mechanism detected was bla_CMY-2, present in 467/667 (70%) of AMR S. Dublin genomes. Although strA-B, tet(A) and bla_TEM-1 were observed in multiple clades, floR, sul2 and bla_CMY-2 were restricted to North American isolates that fell in clade 5 (Supplementary Figure 2-3). This MDR profile, in addition to strA-strB and tet(A), was mediated by IncC plasmids in 593 genomes (n=509 from the United States of America [USA], n=82 and n=2 from the UK) (Supplementary Figure 3, 4C). Nearly all clade 3 isolates linked to the UK had no AMR determinants, except for four isolates collected between 2010-2015 (Supplementary Figure 1-3, Figure 3, Supplementary Table 1).

In addition to AMR determinants, we characterised biocide and heavy metal resistance (HMR) profiles of the S. Dublin isolates. Although gold resistance determinants golS and golT, were ubiquitous in S. Dublin (n=1,302/1,303; 99.9%), biocide and HMR determinants generally co-occurred with AMR determinants (Figure 3C, Supplementary Figures 2-3). The biocide quaternary ammonium compound resistance determinant, qacEdelta1, was primarily identified in MDR genomes from clade 5 isolates (n=76/614; 12.4%), associated with North America. Further, HMR determinants ter (tellurium metal), pcoS (copper metal) and arsC (arsenate metal) were detected in a single genome, SRR5905515 from North American animal isolate collected in 2016, with AMR mechanisms to seven different antimicrobial classes, including colistin (mcr-9.1) (Supplementary Figure 3, Supplementary Table 1). Mercury resistance genes (merABDEPRT) were found in 341/614 (55.5%) MDR clade 5 genomes, while merCPRT genes were detected in single MDR genomes from Australia and South America (Supplementary Table 1). The mer operon co-occurred in 90.1% (n = 534) of the isolates mediated by the IncC MDR plasmid, suggesting an association of this operon with resistance determinants (Supplementary Figure 3, Supplementary Table 1).

Characterisation of S. Dublin plasmids and identification of a novel hybrid plasmid

Both IncX1 (n=1269/1303; 97.4%) and IncFII(S) (n=1213/1303; 93.1%) plasmid replicons were found at high frequency within the dataset, in line with previous findings (12). Three plasmid replicons, namely IncN, IncX1 and IncFII(S), and a MOBF/MOBP relaxase, were detected in one of the Australian S. Dublin ST10 genomes from clade 1 (AUSMDU00035676). Comparative analysis revealed the AUSMDU00035676 plasmid lacked the spv operon and had low homology to reference IncN and IncX1/IncFII(S) plasmids previously associated with Salmonella (Supplementary Figure 5, Supplementary Table 5). The loss of the spv operon may be associated with a less invasive phenotype (32). This finding suggests the detection of a potentially novel hybrid plasmid with IncN/IncX1/IncFII(S) markers that encodes AMR and HMR determinants flanked by IS26 elements, which may aid gene mobilisation (33). In contrast, the AUSMDU00056868 plasmid had IncX1/IncFII(S) plasmid replicon markers and was found to not encode any mechanisms for resistance. This plasmid was >99% homologous to known virulence reference S. Dublin plasmids pOU1115 and pSDU2_USMARC_69807 and encoded the spv operon (Supplementary Table 5, Supplementary Figure 6). This suggests that the S. Dublin virulence plasmid is conserved, while the spv operon that has been shown to enhance disease severity (14,16,34–36) is maintained on the plasmid.

Diversity in virulence determinants across S. Dublin lineages

To investigate the potential for invasiveness within S. Dublin, we applied a previously described bioinformatic tool to our dataset of 1,303 S. Dublin genomes (37). In brief, this tool infers the likelihood of invasiveness by assigning an index score between 0 to 1 that reflects the level of genome degradation, whereby a higher score suggests genomic features that may be more associated with a narrower host range and great likelihood of invasiveness (37). We performed initial validation work using 36 complete publicly available reference genomes, covering known invasive and non-invasive phenotypes (Supplementary Table 8). Using this validation set, invasive and largely host-restricted serovars, including S. Choleraesuis, S. Paratyphi and S. Typhi, scored >0.5, suggestive of an invasive genomic profile (Figure 4A, Supplementary Table 8). Conversely, non-invasive NTS serovars with broader ecological niches, including S. Anatum, ST19 S. Typhimurium and S. Derby, scored <0.5 (Figure 4A, Supplementary Table 8). When applied to our S. Dublin dataset, only eleven genomes (n=11/1,303; 0.84%) had an index score < 0.5. These included three ST10 genomes and all eight ST74 genomes (Figure 4B, Supplementary Table 6). No association was observed between index scores and region or source attribution (Figure 4B).

Figure 4.

Putative virulence determinants identified in S. Dublin genomes included SPIs, the spv operon (16), prophages and other individual genes (Supplementary Table 6). Some putative virulence factors, such as pagN [cell adhesion and invasion (38,39)] and pgtE, [virulence (26,40)], were found in 1,303 (100%, n=1,303/1,303) and 1,302 (99.9%, n=1,302/1,303) genomes respectively, and no SNPs in the promoter region of pgtE were detected, a mutation previously reported in invasive S. Typhimurium ST313 (26). In contrast, ST313-gi of unknown function (26,40), the Gifsy-2 like prophage [unknown function (26,40)] and ggt gene thought to be involved in gut colonisation (34), were detected in the major ST10 S. Dublin lineages, and absent in ST74 S. Dublin. The viaB operon encoding the Vi flagella antigen was only detected in SLV ST73, a historical isolate from France that displayed the Du3 electrophoretic type that was consistent with a previous study (20). This suggests that the Vi antigen is not conserved in the S. Dublin population or associated with virulence. Further, there was no association between the number of virulence determinants and the AMR and HMR profile of isolates, suggesting a potential lack of fitness cost associated with resistance mechanisms and virulence (Supplementary Table 6).

The presence of SPI-6 and SPI-19 differed across S. Dublin isolates. Both SPIs encode T6SSs that provide a means for injecting effector proteins into adjacent competing bacteria or host cells for immune evasion (22,41). SPI-19 was absent in the ST74 S. Dublin lineage but present in all ST10 S. Dublin (Supplementary Figure 7). The genetic differences between ST74 and ST10 were reflected in the disparity of their index scores of <0.5 and >0.5, respectively. The SPI-6 region was detected in all isolates, however <90% sequence coverage was detected in 422/1,303 genomes (32.4%). Long read sequencing on three genomes of the ST74 lineage and three genomes from the ST10 clade 1 lineage corroborated these findings (Figure 4C-E), demonstrating truncation of SPI-6 with the T6SS machinery (Figure 4D). Long-read sequencing also revealed an absence of SPI-19 in the three ST74 isolates and presence of SPI-19 in the three ST10 isolates (Figure 4C, E). The absence of SPI-19 in ST74 suggests that this deletion event occurred during its divergence from the primary S. Dublin lineage, ST10.

Infection of human macrophages reveals distinct intracellular replication patterns between ST10 and ST74 lineages

Given the distinct genomic profiles of ST10 and ST74 lineages, we compared the phenotypic profile of a selection of Australian isolates from each lineage in vitro. The immortalised human macrophage cell line, THP-1, was used to assess invasion, intracellular replication and host cytotoxicity induced by a subset of clinical isolates from the two lineages. Specifically, we assayed an ST10 population comprising eleven isolates, including an isolate from the highly related ST4293 lineage, and an ST74 population comprising seven isolates including an isolate from the highly related ST1545 lineage (Supplementary Table 10). Although predicted invasive index scores were higher for the ST10 population (Figure 4B), each isolate failed to replicate more than 2-fold over 24 hours in THP-1 cells (Figure 5A), whereas isolates within the ST74 population displayed a sustained increase in fold change of intracellular replication, up to 8-fold, over 24 hours (Figure 5B). Comparing bacterial growth (CFU/ml) across both sequence types, both lineages were comparable for initial bacterial uptake (Figure 5C) however, the ST74 population replicated to significantly higher levels in THP-1 macrophages at 9-and 24-hours post-infection (hpi) compared to the ST10 population (Figure 5D-F). Given that NTS also replicate in epithelial cells, we utilised an immortalised human intestinal epithelial cell line, HT-29, to assess S. Dublin replication. In contrast to the results obtained in macrophages, we found that 6 out of 11 isolates within the ST10 population replicated more efficiently (up to 10-fold) in HT-29 cells, whereas all seven isolates within the ST74 population showed minimal replication (up to 2-fold) over 24 hours (Supplementary Figure 8A-C).

Figure 5.

Increased replication of ST74 isolates in human macrophages does not correlate with increased host cell death

It is well documented that protype S. Typhimurium (ST19) potently induce macrophage cytotoxicity both in vitro and in vivo (42). Here, we utilised a cellular lactate dehydrogenase (LDH) release assay to assess THP-1 cytotoxicity in response to infection with S. Dublin ST74 and ST10 populations. Surprisingly, the substantial intracellular replication of the ST74 population in THP-1 cells induced significantly less host cell cytotoxicity at 3 and 9 hpi compared with the ST10 population suggesting potential immune evasion by ST74 (Figure 6A, B). By 24 hpi however, cytotoxicity induced by both ST74 and ST10 populations was relatively high, with no significant difference observed across the sequence types (Figure 6C). In contrast, there was minimal cytotoxicity observed in infected HT-29 cells across the board (Supplementary Figure 8D), suggesting a less potent activation of host cell death in infected epithelial cells compared with macrophages.

Figure 6.

To understand how the ST74 population may be supressing macrophage death at earlier stages of infection, we measured protein secretion of IL-1ß via cytokine release and immunoblot assays at 9 hpi, respectively. IL-1ß is a potent pro-inflammatory cytokine that is released upon activation of the NLRC4 inflammasome by components of the Salmonella type III secretion apparatus (43), resulting in a lytic cell death process known as pyroptosis (44). Here, we observed significantly less IL-1ß release by THP-1 cells infected with isolates from the ST74 population at 9 hpi compared with those from the ST10 population (Figure 6D), directly reflecting the LDH cytotoxicity data represented in Figure 6B. To further support our observations, we selected single representative isolates of ST74 and ST10 to perform immunoblot analyses of cleaved IL-1ß (activated form) within THP-1 cells. Here, ST74 induced less IL-1ß cleavage at 9 hpi compared to the ST10 (Figure 6E). We included the ubiquitous laboratory NTS reference strain, S. Typhimurium SL1344 in the immunoblot analysis as this is known to potently induce pyroptosis in human macrophages (45). As expected, SL1344 induced significant IL-1ß cleavage during infection of THP-1 macrophages (Figure 6E). To assess the correlation of cellular cytotoxicity as measured by LDH, IL-1ß release and intracellular replication, we generated a PCA plot of all isolates within the ST74 and ST10 populations. Here we can clearly demonstrate that isolates from the ST74 population have increased replicative ability, induce less cytotoxicity, and lower levels of IL-1ß release (Figure 6F) compared with those within the ST10 population. There were however, some ST10 isolates that induced similarly low levels of host cell death, presumably due to lower levels of intracellular replication (Figure 6F and Supplementary Figure 9A, B). These isolates generally clustered away from the ST10 isolates that induced more cell death (Supplementary Figure 9C).

A pangenome analysis of the isolates from both ST10 and ST74 lineages revealed that out of a total of 5,792 genes, there were 473 and 1,143 genes specific to the ST10 and ST74 lineage, respectively (Supplementary Figure 10, Supplementary Table 11). Further, out of 1,143 genes, there were 362 core genes present in all isolates from the ST74 lineage (Supplementary Table 11). While these findings require further investigation, an increased gene content and genome size were consistently observed in all isolates of the ST74 lineage (Supplementary Table 9 and 11). These results suggest an association with unique genome content of the ST74 lineage and may account for altered host responses and replicative ability (Supplementary Table 9).

Discussion

Our data provide comprehensive insights into the emergence and spread of distinct lineages of S. Dublin, with striking geographic associations with AMR and virulence, consistent with previous work (12). Analysing the largest global dataset of S. Dublin genomes to date, we show that highly resistant S. Dublin is almost entirely confined to North America, with emergence in the late 20^th century and ongoing local microevolution. North American S. Dublin isolates in clade 5 were characterised by the presence of an IncC plasmid, often with co-location of AMR and heavy metal resistance determinants. Heavy metals and biocides are both used in the agri-food industry and may contribute to the persistence of AMR via co-selection (46). The co-occurrence of AMR, HMR and biocide resistance determinants in animal, food and human S. Dublin genomes highlights the potential for selection pressure in one ecological niche (e.g. the use of biocides in farming or heavy metals in animal feeds) leading to AMR in another (e.g. AMR in humans).

The phylogeographic separation of S. Dublin lineages is consistent with other studies (7,12,47), and is likely to be influenced by several factors, including local AMR use and patterns of livestock movement. Given the livestock reservoir of S. Dublin is relatively host-restricted to cattle, it is conceivable that both historic and contemporary global cattle movement may have contributed to region-specific differences in S. Dublin. For example, Australian S. Dublin genomes in clade 1 shared a common ancestor with North American isolates in clade 5. During the early 1930’s, Brahman cattle were imported from the U.S into Australia (48), supporting our phylodynamic finding that the common ancestral link between S. Dublin isolates from North America and Australia emerged in the early-mid 20^th century. Following region-specific introductions, it is likely that local antimicrobial and husbandry practices may then contribute to the local evolution of S. Dublin. In contrast to North American MDR isolates, ST10 S. Dublin isolates from Australia were predominantly fully susceptible to antimicrobials. Selective pressure, misdirected treatment, the misuse of antimicrobials in cattle and MDR plasmids in agricultural settings have been associated with the increase of MDR S. Dublin cases in North America (49). Conversely, the use of antimicrobials has been extensively minimised in the Australian dairy industry, hence cattle infections are rarely associated with AMR (50). Consistent with regional differences in AMR profiles were geographic differences in plasmid profiles. IncC plasmids are large (70-250kb), can encode up to 10 AMR determinants and are widespread in North America (51–53), while IncN plasmids are smaller (30-70kb), encode up to seven AMR determinants and are widespread in Europe (51,53,54).

Of note, amongst S. Dublin isolates from Australia, genomic analysis (cf. serotyping) revealed a smaller population of ST74 isolates. S. Dublin ST74 was previously identified as S. Enteritidis ST74, and sits as an intermediate between S. Dublin and S. Enteritidis (19), and has been understudied to date in studies of S. Dublin. The difference between the antigenic formula of S. Dublin and S. Enteritidis is in the H1 flagella region, which is typed as “g,p” for S. Dublin and “g,m” for S. Enteritidis. ST74 S. Dublin genomes were characterised by: (i) a lack of any AMR or HMR determinants; (ii) lack of complete SPI-6 regions, specifically the T6SS, and (iii) a lack of putative virulence determinants including ST313-gi, Gifsy-2 like prophage and SPI-19 (which encodes a T6SS mechanism for intestinal colonisation and survival that interacts with the host by mediating macrophage cytotoxicity (55)). Our data demonstrated that while machine learning predicted the ST74 population to have a lower potential for invasiveness compared to ST10, the genetic differences in ST74 may contribute to increased replication in human macrophages. Collectively, these findings highlight the co-circulation in Australia of two S. Dublin populations, an ST10 and ST74 population, which may be phenotypically characterised as S. Dublin.

The lack of the Vi antigen in 1,302/1,303 (99.9%) S. Dublin genomes has several implications. First, the antigenic formula of S. Dublin, 1,9,12[Vi] may need to be reviewed if most S. Dublin genomes lack the machinery to express the Vi flagella. These limitations with serotyping have been highlighted in a previous study, suggesting a transition to a novel genotyping scheme that would prevent these phenotypic discrepancies (56). Second, the Vi antigen is a vaccine target for S. Typhi (57), and the Vi antigen has been explored as a target for S. Dublin (58). While previous studies have noted that some S. Dublin genomes may be capable of producing the Vi antigen (58), the genomic data from this study shows that these genomes are in the minority of the S. Dublin population and may not be stable targets for S. Dublin specific vaccines.

Our phenotypic data demonstrated a striking difference in replication dynamics between ST10 and ST74 populations in human macrophages, whereby isolates within the ST74 population replicated significantly over 24 hours compared with the rapid clearance of those within the ST10 population after 9 hours of infection. The ST74 population generally induced significantly less host cell death during the early-mid stages of macrophage infection, supported by a limited processing and release of the pro-inflammatory cytokine IL-1ß at 9 hpi. Research to date suggest that Salmonellae are potent activators of inflammasome activation (59), however most of the data that support this theory have been generated using laboratory adapted strains of S. Typhimurium. Our data, however, suggest that ST74 (and highly related) isolates may utilise specific immune evasion mechanisms to avoid host recognition and the activation of host cell death signalling pathways in phagocytes in the early stages of infection. This observation is not unique to NTS. In 2015, Carden et al. (60) demonstrated that S. Typhimurium ST313, a major cause of bacteraemia in humans in sub-Saharan Africa, induced significantly less inflammasome activation compared with S. Typhimurium ST19 isolates during infection of murine bone marrow-derived macrophages. Such events could potentially facilitate increased pathogen replication and dissemination in the host at later stages of infection. Indeed, we observed similarly high levels of host cell cytotoxicity in response to ST10, ST74 and related isolates at 24 hpi, indicating that ST74 likely induce cell death via alternative mechanisms once intracellular bacterial levels become unsustainable within the macrophage. Further investigation is required to identify and validate specific genomic factors that may contribute to the mechanisms of enhanced host cell survival, in particular inhibition of inflammasome activation and subsequent IL-1ß release during infection with ST74 isolates. Indeed, a pangenome analysis of the ST10 and ST74 lineages included in phenotypic assays have revealed multiple SPI and non-SPI-related genes that are unique to the ST74 lineage (Supplementary Figure 11), >900 of which remain uncharacterised. A systematic analysis of these unique genes, including comparison with other invasive NTS serovars, protein function prediction and targeted mutagenesis may reveal genes of interest that contribute to pathogen survival in the hostile immune environment of a macrophage.

To date, SPI-19 has been identified in NTS serovars Dublin, Weltevreden, Agona, Gallinarum and Enteritidis (22), however the lack of SPI-19 in the S. Dublin ST74 lineage (Figure 4C and E, Supplementary Figure 7, Supplementary Table 6) is an interesting observation and may suggest a specific adaptation that supports increased replicative ability in phagocytes. In previous studies, S. Dublin SPI-19 has been associated with colonisation and manipulation of intestinal microbiome in poultry (21,55), as well as contributing to oral (but not systemic) virulence in mice (55). Collectively, these outcomes suggest a role for SPI-19 at the mucosal-luminal interface of the intestinal tract, possibly contributing to colonisation of intestinal epithelial cells. Given that ST74 isolates replicate so efficiently in macrophages, it may be that they simply do not require the added benefit of SPI-19-mediated mucosal colonisation and rely solely on their ability to access phagocytes via M cell or dendritic cell uptake at the apical surface of the intestine. The fact that ST74 isolates have a much larger pangenome than ST10 isolates, may also contribute to successful colonisation, replication and survival within the host and thus compensate for the lack of SPI-19 function. Regardless, this work highlights how little we know about the specific host-pathogen interactions of NTS and what specifically drives the success of potential emerging lineages of iNTS.

This study has both strengths and limitations. It represents the most comprehensive genomic study of S. Dublin to date, systematically characterising antimicrobial, heavy metal and biocide resistance determinants, in addition to exploring the virulome in the global population of S. Dublin. Importantly, our phenotypic data reveals the caution required when making generalised statements about infection outcomes and virulence mechanisms of Salmonella serovars. Our data clearly shows that within the Dublin serovar, there are distinct phenotypic differences observed across sequence types and host cell types, all of which could impact significantly on future development of effective therapeutics or vaccine strategies for NTS infections. The key limitation of this dataset is the geographical distribution bias with an overrepresentation of North American genomes and the lack of data from other regions, such as Asia. However, this is a common limitation of many genomic studies and often results in a bias towards high income countries that have greater resources for whole genome sequencing. This situation may change in the future as many countries build their public health genomics capacity.

In conclusion, we provide genomic insights into the evolutionary, resistance and virulence dynamics in S. Dublin circulating globally. The increasing incidence of iNTS caused by S. Dublin may warrant this pathogen becoming a notifiable disease, similar to typhoid and paratyphoid fever. The distinct virulence and AMR profiles of S. Dublin globally warrant revisitation of the current classification scheme for NTS, in keeping with the rapid and ongoing evolution of these significant pathogens. Importantly, our data further highlight the need for a One Health approach to genomic surveillance.

Methods

Setting and dataset

For Australian isolates, whole genome sequencing (WGS) was performed on all viable S. Dublin isolates from human clinical samples received at the Microbiological Diagnostic Unit Public Health Laboratory (MDU PHL), the primary bacterial reference laboratory for the State of Victoria in Australia. In total, 53 S. Dublin genomes between 2000-2019 were included. For publicly available data, we included S. Dublin sequences that had associated metadata (geographical origin and date of sample collection), providing a total of 1,250 publicly available sequences. Data were collected in accordance with the Victorian Public Health and Wellbeing Act 2008. Ethical approval was received from the University of Melbourne Human Research Ethics Committee (study number 1954615.3).

Genome quality control, assembly and serovar prediction

Genomes were assembled and quality controlled using the Nullarbor pipeline v2.0.20191013 (https://github.com/tseemann/nullarbor). Briefly, Skesa v2.4.0 (61) was used to assemble the genomes before determining the sequence type (ST) using the multilocus sequence typing (MLST) scheme ‘senterica’ within mlst v2.19.0 (https://github.com/tseemann/mlst). Kraken v1.1.1 (62) was then used to confirm the taxonomic classification. Snippy v4.6.0 (https://github.com/tseemann/snippy) with parameters set to a 0.90 variant rate proportion against a minimum coverage of 10X mapped the reads against a well-studied reference genome, S. Dublin CT_02021853 (NC_011205.1). To be included in subsequent analysis, genomes had to pass the following quality control parameters: mapped to ≥ 80% of the reference genome; had a ≥ 80% taxonomic sequence similarity to S. enterica; N50 >25000, and <300 contigs. To verify that each genome was S. Dublin, all sequences were typed using Salmonella in silico Typing Resource (SISTR) v1.1.1 (63). The electrophoretic type of each isolate was determined using Antimicrobial Resistance Identification by Assembly (ARIBA) v2.14.5 (64) against the reference Du 2 (Genbank accession: M84972.1).

Core genome phylogenomic analysis

Variant calls from snippy were used to construct a core genome single nucleotide polymorphism (SNP) alignment using snippy-core v4.6.0 (https://github.com/tseemann/snippy). Phaster v2020-12-22 (65) was used to create a BED file containing prophage regions. The option ‘—mask’ in snippy-core was used to exclude the identified prophage regions from the BED file. Constant sites from the reference genome were obtained using the script iqtree-calc_const_sites.sh (https://github.com/MDU-PHL/mdu-tools/blob/master/bin/iqtree-calc_const_sites.sh). A maximum-likelihood (ML) phylogeny was inferred using IQ-TREE v2.1.0 (66) with parameters set to a GTR+G4 nucleotide substitution model with constant sites and ultrafast bootstrapping (67) to 1000 replicates. The resulting ML tree was then visualised using ggtree v2.3.3 R package (68). Population clades were defined using rhierBAPS v1.1.3 to one level on an alignment of 2,957 core SNPS of 1,295 isolates in the ST10 S. Dublin lineage (69)

In silico prediction of resistance determinants and plasmid replicons

To screen for quinolone point mutations and acquired AMR determinants, ARIBA and AbriTAMR v0.2.2 (https://github.com/MDU-PHL/abritamr) were used, respectively. The AMR database used for this study was the NCBI AMRFinderPlus database v2020-01-22.1. Each isolate’s susceptibility profile was classified as either susceptible; resistant to 1-2 antimicrobials, or multidrug resistant (MDR). Heavy metal and biocide determinants were identified by screening the AMRFinder database v2020-07-16.2 with AMRFinder v3.8.4 (70). Replicon types were determined using ABRicate v1.0.1 (https://github.com/tseemann/abricate) at a threshold of 80% nucleotide identity and 95% coverage, and screening against the Plasmidfinder database (71). A strict threshold was applied to ensure the detection of only unambiguous replicon types. To establish the most common resistant profiles, co-occurrence matrices of resistant determinants were constructed in R using igraph v1.2.4.1 (72).

Virulome analysis

Virulence genes were identified using ABRicate at 90% nucleotide identity/coverage against the virulence factor database (vfdb) v2018-03-18 (73). Databases from previous studies (16,26,34,74–78) were also interrogated and collated. This included a database containing additional virulent determinants and Salmonella pathogenicity islands (SPIs) -1, -2, -6 and -19 (Supplementary Table 7). Using ABRicate, the presence of SPIs -1, -2, -6 and -19 were identified at 80% ID and 10% COV. ABRicate was also used to screen the manually created SPI determinant database for genes relevant to SPI-1, -2, -6 and -19 with a threshold of 90% identity/coverage (Supplementary Table 7). To quantify invasiveness, index scores were calculated in silico for each isolate using Wheeler et al.’s tool developed in 2018 for inferring the likelihood of invasiveness (http://www.github.com/UCanCompBio/invasive_Salmonella), with a higher score for an isolate deemed more invasive (37). Prior to running the tool, a set of reference genomes consisting of invasive and non-invasive Salmonella serovars were used as a validation set (Supplementary Table 8). The validation results were then visualised using ggtree.

Bayesian phylodynamic analysis

A subset of 660 S. Dublin genomes from the main cluster was used for BEAST v1.10.4 (79) analysis. This subset comprised all Australian genomes and a proportion of international S. Dublin genomes. The sampling strategy for selecting representative international genomes involved capturing all geographical regions and sources harbouring distinct susceptibility profiles. When two genomes demonstrated the same AMR profile and were from the same region and source, the oldest genome was preferred. Recombination sites were filtered using Gubbins v2.4.1 (80). A core genome SNP alignment of these genomes was created using snippy-core v4.6.0 and then a ML phylogeny was inferred as above. Temporal signal was assessed using TempEst v1.5.3 (81) by assessing the best-fitting root using the heuristic residual mean squared method. The core alignment was then used to create an xml file in BEAUti v1.10.4 using the GTR+G4 substitution model. To identify the most appropriate model for our dataset, a combination of using either a “relaxed lognormal” or a “strict” molecular clock with either a “coalescent constant” or a “coalescent exponential” population prior was tested. For all models, the Markov chain Monte Carlo chain length was set to 300,000,000 and sampling was set to every 20,000 trees with the marginal likelihood estimation set to “generalised stepping stone (GSS)” to provide statistical support. The analysis was repeated in triplicates with and without tip dates to ensure there was no temporal bias. Bayesian evaluation of temporal signal (BETS) was used to determine whether a temporal signal was present (82). To optimise computational resource and time, BEAST v1.10.4 was used with BEAGLE v3.0.2 (83). A Maximum Clade Credibility (MCC) tree from the highest supported model using GSS was extracted using TreeAnnotator v1.10.4 and median heights. The resulting MCC tree was then visualised using ggtree.

Long read sequencing

The Oxford Nanopore Technologies MinION platform was used to obtain complete genomes of six Australian S. Dublin genomes (Supplementary Table 9). These genomes, comprised of ST10 lineage clade 1 (AUSMDU00056868 and AUSMDU00035676,) and ST74 (AUSMDU00016676 and AUSMDU00021231), were selected based on their AMR profile, virulome and index scores. To prepare new DNA libraries for long read sequencing, genomes were freshly re-cultured. DNA was extracted using the Virus/Pathogen DSP midi kit on the QiaSymphony SP instrument. Sequencing libraries were prepared using a Ligation sequencing kit (SQK-LSK109) with native barcode expansion and sequenced for 48 hours on a GridION Mk2 instrument using a R9.4 flowcell. WGS was performed using the Illumina NextSeq 500/550 platform (Illumina, San Diego, CA). Demultiplexed long reads and corresponding WGS reads were subjected to QC to exclude unreliable reads using fastp v0.20.1 (84) and Filtlong v0.2.1 (https://github.com/rrwick/Filtlong), respectively. Long reads were assembled using Trycycler v0.4.1 (85). The final sequence was polished using Medaka v1.4.3 (https://github.com/nanoporetech/medaka). To ensure small plasmids were not excluded during the initial QC step, hybrid assemblies using both long and WGS reads were generated Unicycler v0.4.4 (86). Trycycler assemblies were checked against Unicycler assemblies with Bandage v0.8.1 (87). The final genome assemblies were confirmed as S. Dublin using SISTR and annotated using Prokka v1.14.6 (88) (Supplementary Table 9).

Chromosomal assemblies were aligned to complete genomes of known S. Dublin references CT_02021853 (CP001144.1, without Vi flagella) and CFSAN000518 (CP074226.1, with Vi flagella) and visualised using Blast Ring Image Generator (BRIG) v0.95 (89). Plasmid assemblies were screened using PlasmidFinder (90) and mob-typer v3.1.0 (91) to identify replicon and relaxase types, respectively. Plasmid profiles were subjected to BLASTn against previously identified IncX1/IncFII(S) virulence plasmids reported in S. Dublin and IncN plasmids associated with serovars of Salmonella (Supplementary Table 5). Genome comparisons of plasmids demonstrating the top four highest BLASTn max scores were visualised using BRIG. Prokka was used to annotate the regions of interest, while Easyfig v2.2.2 (92) was used to compare regions of the genome. Regions of the genome were visualised using the R package gggenomes v0.9.12.9000 (https://github.com/thackl/gggenomes).

Human cell line infections

Human macrophage (THP-1) (ATCC® TIB-202™) and colonic epithelial (HT-29) (ATCC® HTB-38™) cell lines were maintained in Roswell Park Memorial Institute (RPMI) 1640 media + 200 mM GlutaMAX (Life Technologies) supplemented with 10% v/v Fetal Bovine Serum (FBS; Bovogen) and grown in a humidified 5% CO2 37°C incubator. THP-1 monocytes were differentiated into macrophages with 50 ng/mL phorbol 12-myristiate-12 acetate (PMA, Sigma-Aldrich) for 72 hours prior to infection.

Single colonies of each S. Dublin isolate were inoculated into 10 mL Luria Bertani (LB) broth and grown overnight at 37°C at 200 rpm, then sub-cultured for 3 hours in the same conditions in 5 mL LB. Colony forming units/mL were measured by absorbance at 600 nM (OD₆₀₀). THP-1 and HT-29 cells were infected at a multiplicity of infection (MOI) of 10 in RPMI in triplicate wells of a 96-well plate and centrifuged at 525 x g for 5 minutes to synchronize bacterial uptake. Infected cells were incubated at 37°C for 30 minutes, then infective media removed and replaced with 100 μg/mL gentamicin in 10% v/v FBS/RPMI to inhibit extracellular bacterial growth. 0 hours post infection (hpi) samples were collected at this time to measure initial bacterial uptake. Cells were then incubated in high dose gentamicin for 1 hour, then replaced with 10 μg/mL gentamicin in 10% v/v FBS/RPMI for the remainder of the infection. For THP-1 cells, samples were collected at time 0-, 3-, 9- and 24-hours post-infection. For HT-29 cells, samples were collected at 0- and 24-hours post-infection.

Cytotoxicity of infection was measured by lactate dehydrogenase (LDH) release by infected cells before lysis for intracellular bacterial enumeration. Briefly, 35 uL of infected (and uninfected) cell culture supernatants were collected from each well and LDH quantified as per manufacturer’s instructions (Promega), with percent cytotoxicity calculated as a proportion of lysed control cells and blank media. Following this, infected cells were washed three times with warm phosphate-buffered saline (PBS) to remove any lingering extracellular bacteria, then lysed in 0.1% Triton X-100/PBS and plated out onto LB agar after serial dilution. Bacterial enumeration and cell viability assays were performed over 3 biological replicates, each performed independently, with a new passage of host cells, on a separate day. Within each biological replicate, technical triplicates were performed for each isolate. All statistical analyses were performed using Prism software (GraphPad Software v9.0) and determined by nested students t-test, which allowed for inclusion of individual isolates in sub-groups of each sequence type (ST10 vs ST74).

Measurement of cellular responses

THP-1 cells were infected at an MOI of 10 with a representative isolate from each of the ST10 and ST74 lineages as described above, then lysed in radio-immunoprecipitation (RIPA) buffer (20 mM Tris pH 8.0, 0.5 mM EDTA, 150 mM NaCl, 1% v/v Triton X-100, 0.5% w/v sodium deoxycholate, 0.1% w/v SDS) supplemented with Complete protease and phosphatase inhibitor cocktails (Roche) at 0-, 3-, 8- and 24 hpi . Total protein concentration was quantified by Bicinchoninic acid (BCA) assay and 15 μg separated on 4-12% Bolt Bis-Tris polyacrylamide gels (Thermo Scientific), transferred to nitrocellulose (iBlot, Thermo Scientific) and blotted for cleaved IL-1β (clone D3A3Z, Cell Signalling). Cell culture supernatants were also collected from 96-well plate infections at 9 hpi and secreted IL-1β measured by DuoSet© Human IL-1β/IL-1F2 ELISA as per manufacturer’s instructions (R&D Systems). A core genome SNP ML phylogeny was constructed on the Australian ST10 lineage to demonstrate the genetic relationship of representative ST10 isolates that underwent in vitro phenotypic assays.

Pangenome analysis

The pangenomes of the isolates used in the macrophages and epithelial infections were analysed using Panaroo v1.3.4 (93) set to default threshold settings, under strict mode and with the removal of invalid genes on the draft genome assemblies (Supplementary Table 11).

Data availability

Paired end reads of public data were obtained from either the European Nucleotide Archive (ENA) or National Center for Biotechnology Information (NCBI). The Australian short read data and complete genomes for the four Australian S. Dublin isolates are available at the NCBI Sequence Read Archive under BioProject PRJNA319593. Details of all 1,303 genomes included in the study are available in Supplementary Tables 1-11. These supplementary tables also includes i) the resistance, virulence, Vi and fliC profiles of all 1,303 genomes; ii) the genbank accession numbers of all reference genomes, virulence determinants and IncN plasmids used in this study, and iii) all raw CFU, cytotoxicity and ELISA values from the phenotypic analyses. The interactive tree is available in Microreact under https://microreact.org/project/sdublinpaper).

Code availability

The machine learning tool developed by Wheeler et al. for predicting the likelihood of invasiveness is available in a public Github depository http://www.github.com/UCanCompBio/invasive_Salmonella.

Acknowledgements

We thank all staff members of the Enterics and MDDI section at the Microbiological Diagnostic Unit Public Health Laboratory (MDU PHL) for their technical assistance. We would also like to thank Jake Lacey for his assistance uploading the reads to NCBI. The BEAST analysis for this paper was undertaken using the LIEF HPC-GPGPU Facility hosted at the University of Melbourne. This Facility was established with the assistance of LIEF Grant LE170100200.

Additional information

Contributions

Study design and oversight: CMS, RA, DAW, JSP, DJI. Isolate collection and metadata: MV, PA, SAB, BJP. Performed wet lab experiments: RA, JSP. Data analysis: CMS, RA, JSP, DJI. Manuscript writing: CMS, RA, DAW, JSP, DJI. All authors contributed to manuscript editing.

Ethics declarations

Competing interests

The authors declare no competing financial interests.

Funding

DAW was supported by a National Health and Medical Research Council (NHMRC) Investigator Grant (GNT1174555) and a Medical Research Future Fund (MRFF) Grant (FSPGN000045). BPH is supported by NHMRC Investigator Grant (GNT1196103). DJI is supported by an NHMRC Investigator Grant (GNT1195210). JSP was supported by a Sylvia and Charles Viertel Senior Medical Fellowship. MDU PHL is funded by the Victorian Government, Australia. No conflict of interest declared.