Distinct genomic trajectory among invasive Salmonella Typhimurium ST313 infections

Domain: chemistry
Published: 2026-03-09
Editorial outcome: major_revision

Abstract
Invasive non-typhoidal Salmonella (iNTS) disease, primarily manifesting as bloodstream infections and dominated by Salmonella Typhimurium sequence type 313 (ST313), represents a major global health challenge, particularly in sub-Saharan Africa. However, the transmission landscape and evolutionary driver of ST313 remain incompletely understood. Here, we compile the global dataset to date, comprising 11,361 human-derived iNTS genomes and 3,115 ST313 genomes from diverse sources, enabling identification of a highly antimicrobial-resistant sublineage 2.4 and an emerging lineage 4. Phylogenomic analyses suggest patterns consistent with intercontinental dissemination of ST313, which may be associated with modernized human biomass movement. Sublineage 2.4, an antimicrobial-resistant offshoot that displaced lineage 1, harbors extended-spectrum β-lactamase genes, including blaOXA-1, blaCTX-M-15, often linked to mobilome, suggesting a potential trend of increasing resistance within lineage 2. Furthermore, mobilome-driven alterations in antimicrobial-resistant gene carriage correlate with the post-2006 decline of lineage 1 in Africa. Other genetic factors likely enhance invasive potential, as evidenced by enrichment of pseudogenes and positively selected functional genes, potentially driving the replacement by the newly identified ST313 clone. Our findings highlight lineage-specific selection pressures influenced by anthropogenic factors, such as antimicrobial stewardship, that modulate ST313 evolution and underscore the need for enhanced surveillance.

Introduction
Invasive bacterial bloodstream infections pose a significant threat to global health, particularly in resource-limited settings. The World Health Organization (WHO) prioritizes Salmonella species among the top pathogens posing significant risks for gastrointestinal and systemic disorders in humans1–3. Invasive non-typhoidal Salmonella (iNTS) diseases, defined as infections that invade normally sterile sites such as blood or cerebrospinal fluid (excluding fecal sources), are especially alarming due to their mimicry of typhoid fever and high case fatality rates4,5. Recent estimates from the Global Burden of Disease 2019 study indicate that iNTS causes approximately 594,000 cases (95% uncertainty interval [UI] 486,000–718,000) annually worldwide, leading to 79,000 deaths (95% UI 43,000–124,000) and 6.11 million disability-adjusted life years (95% UI 3.32–9.71)6. Notably, incidence rates often exceed 100 cases per 100,000 person-years among high-risk populations, including people living with HIV, those with malnutrition, and individuals at the extremes of age (<5 years and ≥70 years)4,7. Despite this substantial burden, the overall global landscape of iNTS remains incompletely understood, hindered by limited surveillance across regions.

Salmonella enterica serovar Typhimurium (S. Typhimurium) emerges as a primary etiological agent of iNTS disease8. In sub-Saharan Africa, the S. Typhimurium sequence type 313 (ST313) clade predominates, accounting for a substantial proportion of cases and being associated with severe manifestations such as septicemia9,10. Genome-based studies have identified at least four major lineages within ST313—lineage 1 (L1), lineage 2 (L2) with sublineages 2.0–2.3, lineage 3 (L3), and a UK-associated lineage—highlighting ongoing diversification11–14. Nevertheless, the complete global population structure, intercontinental dissemination patterns, and the potential for zoonotic or cross-host transmission remain only partially known.

Emergence of dominant bacterial variants, including within ST313, occurs through mutational processes, horizontal gene transfer, and gene loss, often under selective pressures favoring traits like increased transmissibility, virulence, or immune evasion15. Prior research underscores the pivotal role of antimicrobial resistance (AMR), mediated by acquired resistance genes (ARGs), in lineage dynamics and replacement within ST31316. For instance, chloramphenicol resistance, largely restricted to L2, likely facilitated L1 displacement by L2 around 200117. Since 2005, amid rising use of third-generation cephalosporins and azithromycin, ST313 isolates in Africa—particularly L2 strains—have exhibited adaptive resistance signatures to first-line antibiotics18,19. Beyond AMR, progressive genome degradation—evidenced by the accumulation of pseudogenes and alterations in functional gene expression—has been implicated in enhancing the invasive potential of ST31320–22. These findings highlight existing knowledge gaps in constructing a comprehensive genomic landscape of ST313 to understand its evolutionary trajectory, the driving forces behind its adaptation and the design of precise disease interventions.

To address these gaps, we curated the largest global iNTS genomic dataset to date, comprising 11,361 human-derived genomes from non-fecal clinical sources, and in parallel compiled a partially overlapping collection of 3,115 ST313 isolates collected between 1966 and 2023. Leveraging this unique resource, we aimed to: (1) delineate the global landscape of iNTS, with a particular focus on ST313; (2) characterize its population structure, patterns of intercontinental dissemination, and evidence for cross-host transmission at single-nucleotide polymorphism (SNP) resolution; and (3) investigate the roles of AMR, mobilome dynamics, and other genomic factors in driving lineage replacement within ST313.

Results
The global landscape of the iNTS collection.
Analysis of 11,361 publicly available iNTS genomes (Supplementary Data 1) identifies Africa as the predominant region of iNTS sampling (n = 4,741; 41.7%), followed by North America (n = 4,351; 38.3%). Although Asia accounts for fewer samples (n = 1,350; 11.9%), these were primarily from China (1,173 of 1,350; 86.9%) (Fig. S1a). Calculations of the sampling rates revealed an overall upward trend across countries, with the largest increase observed in South Africa (Fig. S1b). Serovar classification underscored the predominance of S. Typhimurium (n = 3,260, 28.7%) and S. Enteritidis (n = 3,785, 33.3%). Interestingly, in our dataset, S. Typhimurium samples predominated before 2015, but after 2015, S. Enteritidis became predominant (Fig. S1c). Next, to examine the details of this shift within S. Typhimurium, we used multilocus sequence typing (MLST) to identify three major circulating sequence types (STs): ST313 (n = 1,763), ST19 (n = 1,140), and ST34 (n = 237). Among these, the temporal changes in the available ST313 genomes from public databases closely mirrored the post-2015 decline in S. Typhimurium sampling (Fig. S1d–e).

Next, by assembling all publicly available ST313 genomic data, we constructed the most extensive global dataset of ST313 genomes to date, comprising 3,115 isolates (Supplementary Data 2). Our results show that, geographically, ST313 predominated among isolates sampled in Africa (n = 2,266; 72.7%), with a growing but still limited sampling presence in Europe (n = 260; 8.3%) (Fig. 1a–b). From a host perspective, most isolates were human-derived (n = 2,027, 65.1%), with rare detections in animals (n = 15, 0.5%) and food (n = 8, 0.3%) (Fig. 1c). Notably, among human isolates, blood samples were the primary source (1,702/2,027, 84%), far exceeding faecal isolates (96/2,027, 4.7%) (Fig. 1d). Only a small number of ST313 isolates contained sex information (n = 440), of which more than 99% were reported to have been collected at Queen Elizabeth Central Hospital (Blantyre, Malawi)14. Overall, the sex distribution was balanced, with slightly more male cases (n = 226) than female (n = 214) (Fig. 1e).

Global population diversity of ST313 uncovers new lineages
Lineage assignment, based on previous studies13,14, identified five primary lineages: L1 (296/3,115), L2 (2,273/3,115), L3 (164/3,115), L4 (230/3,115), and L5 (152/3,115) (Fig. 2a). Among these, L1 and L2 predominated in Africa. Within L2, phylogenetic analysis revealed a novel sublineage beyond 2.0–2.3, designated 2.4 (61/3,115). The majority of isolates in sublineage 2.4 originated from Kenya (58/61, 95.1%), with two isolates also from the UK (Fig. 2b). Nearly half of L3 genomes (79/164, 48%) were from the Americas. L4, a previously unreported and newly defined lineage here, was concentrated in Africa (159/230, 69.1%). Interestingly, L5, previously referred to as the UK-isolated lineage, included isolates from Australia and Japan. To further assess whether the newly identified sublineage 2.4 and L4 represent single clonal expansion events, we calculated their within-lineage genetic diversity. Pairwise SNP distances within sublineage 2.4 ranged from 0 to 50 nucleotides (mean: 19.1), whereas L4 showed substantially greater divergence, ranging from 0 to 510 nucleotides (mean: 173.4). These findings indicate that neither sublineage 2.4 nor L4 arose from a single clonal expansion event (Fig. S2).

Inference of historical changes in the effective population sizes (Ne) showed distinct declines for L1, L2, L3, and L5. Specifically, L1 exhibited a sharp reduction around the year 2000 (Fig. S3b), while declines in L2 and L3 occurred later, approximately around 2010 (Fig. S4b and Fig. S5b), and L5 around 2017. By contrast, the newly characterized sublineage 2.4 has expanded continuously since its emergence, while L4 has maintained relatively high Ne since 2000 (Fig. 2e).

Phylogeographic analysis suggested that the most recent common ancestor (MRCA) of L1 and L2 originated in Africa, approximately around 1950 (95% highest posterior density [HPD]: 1922-1966) and 1960 (95% HPD: 1955-1965), respectively (Fig. S3a and Fig. S4a). The earliest inferred potential transcontinental transmission event for L2 occurred around 1996 (95% HPD: 1995–1998), while the earliest inferred potential transcontinental transmission event for L1 occurred slightly later, around 2006 (95% HPD: 2003–2009). The earliest transcontinental transmission events of both lineages were potentially associated with Africa and Europe (Fig. S3c and Fig. S4c).

The MRCA of L3, L4, and L5 originated outside Africa. Specifically, the MRCA of L4 is estimated to have appeared around 1731 (95% HPD: 1575–1859) in Asia (Fig. 2c). In contrast, the MRCA of L3 emerged slightly later, around 1874 (95% HPD: 1787–1932) in Asia (Fig. S5a). The MRCA of L5 likely appeared around 1965 (95% HPD: 1943–1982) in Europe (Fig. 2d). The earliest potential transcontinental transmission events of these lineages appear to have occurred around 1899 (95% HPD: 1839–1941) in L3 strains, between Europe and Asia (Fig. S5c). The earliest potential transcontinental transmission events for the newly defined L4 occurred around 1908 (95% HPD: 1864–1946), possibly between Asia and Africa geographically. Finally, the earliest potential transcontinental transmission events for L5 occurred the latest, around 1993 (95% HPD: 1984–2002) (Fig. 2f).

Human-borne ST313 transmission predominates in lineage 2
To identify potential recent transmission events of ST313, we conducted an SNP distance analysis using a 3-SNP threshold as the optimal cutoff (see the Methods section; Fig. 3a and Fig. S6). This identified 1,599 clusters with pairwise distances below the cutoff; 454 clusters (1,970/3,115 isolates) with ≥2 isolates (Fig. 3b). Four clusters (Clusters 357, 342, 131, and 49) showed potential intercontinental transmission events (Fig. 3c). The primary potential recent intercontinental transmission events were concentrated within cluster 357, predominantly associated with lineage 2. These potential transmission events (n = 8) occurred between the Americas (United States) and Africa (Congo), with the estimated timeframe spanning 2019 to 2020. These potential events involved six isolates in total: two from the United States and four from Congo. The geographic information for the Congo isolates was cross-validated against published literature, whereas the isolates from the United States have not yet been reported (Supplementary Data 3). Fewer potential transmission events occurred between Europe (UK) and Africa (Nigeria, n = 3) (Fig. 3d), involving four isolates with well-documented sampling backgrounds (Supplementary Data 3). Considering that these isolates belong to the same cluster and were collected within a relatively short time frame (less than three years), we speculate that the putative recent transmission events represented by these ST313 strains may be associated with international travel; however, alternative exposure routes, such as foodborne transmission through consumption of contaminated food products, cannot be excluded.

Although ST313 has been considered a lineage restricted to human-to-human transmission23, analysis of 454 clusters identified a potential cross-host transmission event within a single cluster (Cluster 453), possibly involving human–canine transmission in the UK (Fig. S7). Despite this observation, potential cross-host transmission events remain extremely rare, accounting for only a single case in our dataset.

Mobilome differences distinguishing ST313 lineages
To characterize genomic differences—particularly within the accessory genome—among the newly identified ST313 lineages, we examined mobile genetic elements (MGEs) previously reported to be carried by ST313, including plasmids (pBT1–pBT3, pSLT-BT, pSLT, pAnK), prophages (BTP1, BTP3, BTP5, BTP6, Gifsy-1, Gifsy-2), transposons (Tn21), and Salmonella pathogenicity islands (SPIs; see Methods). The pSLT was highly conserved across all lineages, consistent with previous reports14. In the newly discovered L4, we observed high carriage rates of pSLT (221/230, 96.1%) and pBT2 (186/230, 80.9%), similar to those in L5. Additionally, sublineage 2.4 exhibited high carriage rates of pSLT (57/61, 93.4%), pBT2 (60/61, 98.4%), and pBT3 (61/61, 100%) (Fig. 4a).

Prophage BTP1, BTP3, BTP5, BTP6, Gifsy-1, and Gifsy-2 are prevalent in the newly discovered L4 and sublineage 2.4. Furthermore, these prophages appear to be lineage-specific within ST313. For example, the carriage rate of BTP3 in the L1 population is only 11.8% (35/296), whereas L3 and L5 rarely carry BTP1 and BTP5. Regarding SPIs, only SPI-5 and SPI-6 are lineage-specific within the ST313 population: SPI-5 is highly prevalent in L1, whereas the carriage rate of SPI-6 in L1 is relatively low (48.3%, 143/296).

Previous studies have shown that Tn21 confers the ability to acquire additional antimicrobial resistance genes (ARGs) upon ST313, and our results indicate that Tn21 is primarily carried by L1 (93%, 275/296) and L2 (72%, 1635/2273), which aligns with prior reports14 (Fig. 4a).

To further investigate the presence of other potentially prevalent MGEs among ST313, we performed sequence alignments between the ST313 genome and a larger-scale MGE database (Fig. S8-S10; Supplementary Data 4). This analysis revealed additional lineage-specific MGEs that were previously unreported (Fig. S11). Specifically, Tn6205 (92%, 2091/2273), Tn6292 (90.1%, 2047/2273), and Tn9 (84.4%, 1919/2273) were highly prevalent in L2, whereas sublineage 2.4 additionally carried Tn602 (100%, 61/61), Tn2012 (96.7%, 59/61), Tn5057 (96.7%, 59/61), and Tn5718 (63.9%, 39/61). Integron In498 exhibited a carriage rate exceeding 60% exclusively in L1 and sublineage 2.4. Notably, prophage 118970_sal3 was highly prevalent across all ST313 lineages, and plasmid classification by replicon type indicated that IncHI2A/rep_cluster_1088-type plasmids were highly prevalent in sublineage 2.4.

Changes in AMR levels are associated with the ST313 lineage shift
AMR could influence the dynamics of bacterial population structure. Our identification of AMR determinants in ST313 revealed that both ARGs and AMR-associated point mutations predominantly occurred in L1 and L2 (Fig. S12; Supplementary Data 5). Sublineage 2.1 is reported as a multidrug-resistance (MDR) branch, conferring ARGs with extended-spectrum β-lactamase (ESBL) production and azithromycin resistance. Unlike sublineage 2.1, sublineage 2.4 additionally carries blaOXA-1, blaCTX-M-15, and tet(A), indicating ongoing AMR evolution in L2. Moreover, a lineage-specific preference for certain ARGs was observed. For instance, the dfrA23 gene conferring trimethoprim resistance was more frequently carried by L1 than L224,25. In contrast, sul2 (sulfonamide resistance gene) and catA1 (chloramphenicol resistance gene) had higher carriage rates in L226,27. The occurrence of AMR-associated point mutations in ST313 is relatively low, with gyrA being slightly more common than gyrB, but overall occurring in only 4.9% (n = 153) of isolates.

The overall MDR rate among ST313 is 63.76% (1,986/3,115), primarily in L1 (34/296, 11.49%) and L2 (1,949/2,273, 85.75%), with the highest prevalence in Africa (1,542/1,986, 77.64%). Notably, the emergence of MDR varies significantly. Unlike L2, which exhibited MDR almost concurrently with its initial emergence, L1 appeared to acquire MDR only post-2000 in Africa (Fig. S3a and Fig. S4a).

A prior study reported that the cat gene was one of the key drivers behind the replacement of L1 by L2 around 2001, which led to the subsequent successful expansion of L2 in Africa17. We also observed a strong concordance between the prevalence of ARGs within L1 and its population dynamics (Spearman ρ = 0.74, p < 0.05, Fig. 5a). To investigate whether the decline of L1 after 2006 was also associated with changes in ARG carriage, we analyzed the temporal trends of 16 types of ARGs carried by L1. Each ARG had a carriage count exceeding 15 (Fig. S13). Further analysis identified six key ARGs, including mcr-9, tet(D), aph(3’)-Ia, sul2, aph(6)-Id, and aph(3”)-Ib, correlating with post-2006 decline (correlation coefficients > 0.6, p < 0.05, Fig. 5, c and d).

To further test our hypothesis that the decline in AMR may have contributed to the disappearance of L1 post-2006, we focused on Malawi, the central source region for L1, and modeled the frequency of ARG transfer. Horizontal ARG transfer frequency modelling showed high rates pre-2000 (~5.4), rising to 2004 amid L1 expansion, then declining post-2005 to half pre-2000 levels by 2010 (Fig. 5b).

Mobilome links with AMR changes in L1 and L2
Our recent study has demonstrated that the mobilome drivers play a critical role in ARG dynamics in Salmonella28. To explore this in ST313, we found that the capacity of the mobilome to contain ARGs varies significantly. Specifically, plasmids exhibit the highest average number of ARGs (n = 5.77), followed by transposons (n = 5.32), integrons (n = 4.23), and prophages (n = 0.02) (Fig. 4b). L2 mobilomes carried more ARGs than L1; sublineage 2.4 plasmids averaged 12.3 ARGs versus 7.3 in other L2 (Fig. 4, c and d; Supplementary Data 6).

The three highly prevalent ARGs—blaTEM-1B, sul1, and ant(3'')-Ia—detected in both L1 and L2, blaTEM-1B was found to be primarily associated with IncFIB/IncFII-type and IncFIB/IncFII/IncQ1-type plasmids. The sul1 is distributed across diverse MGEs, including In1012, In240, Tn6205, IncFIB/IncFII/IncQ1-type plasmids, IncFIB/IncFII-type plasmids, and various prophages (e.g., PWH3a_P1, BcepMu, Escher_520873). Interestingly, ant(3″)-Ia is co-located on In1012, In240, and Tn6205, highlighting the pivotal roles of these MGEs in shaping AMR dynamics in ST313 (Fig. 4 e-h).

Additionally, we investigated the six key ARGs (mcr-9, tet(D), aph(3’)-Ia, sul2, aph(6)-Id, and aph(3”)-Ib) associated with the post-2006 decline of L1 and three key ARGs (blaOXA-1, blaCTX-M-15, and tet(A)) potentially linked to sublineage 2.4 expansions, highlighting IncHI2A/rep_cluster_1088-type resistance plasmids as the key driver (Fig. S14a). Meanwhile, the carriage rates of IncHI2A/rep_cluster_1088-type resistance plasmids in L1 and sublineage 2.4 exhibited consistency with population dynamics (Fig. S14b). By comparing these plasmids from Kenyan isolates with previously reported IncHI2A plasmids in ST313, we found that the IncHI2A/rep_cluster_1088-type resistance plasmids from Kenyan isolates additionally carried an insertion sequence—likely mediated by Tn21—containing blaOXA-1, blaCTX-M-15, and tet(A) (Fig. S15).

Gene gain-or-loss might reshape distinct pathways in ST313 lineages
Although AMR correlated with L1 and L2 shifts, this was less evident in lineages L3-L5, which harbor few ARGs. Bacterial virulence, calculated as invasiveness based on genome degradation, may be involved in this process. To evaluate this, we examined the invasiveness index and found a progressive increase in scores from L1 to L3 (L1: 0.128, L2: 0.132, L3: 0.184), then decreased to L4 (0.131) and L5 (0.129) (Fig. 6a; Supplementary Data 7), corroborated by pseudogene-associated SNPs across different lineages (Fig. S16; Supplementary Data 8). Invasiveness correlated with annual prevalence in L4 (Spearman ρ = 0.67, p < 0.05) and L5 (Spearman ρ = 0.81, p < 0.05) (Fig. 6b).

To explore other genetic determinants for the association study, we conducted a pan-genome analysis and found that genes located in the accessory regions exhibited different patterns (Fig. 6c). Specifically, we identified the largest number of differentially associated genes in L2 (n = 1,526), followed by L1 (n = 982), L4 (n = 710), L3 (n = 689), and L5 (n = 623). To confirm the association between these differential genes and the observed population dynamics, we performed a correlation analysis on individual genes. Our results revealed that most significantly correlated genes exhibit positive correlations (Spearman ρ > 0.6, p < 0.05, Fig. 6d). To elucidate the biological functions of these positively correlated genes, we performed functional enrichment analysis Fig. 6e). The findings indicated that the positive correlation genes enriched in L2 predominantly belonged to pathways associated with replication, recombination, and repair (Supplementary Data 9). In particular, intragenic segments associated with the transposase gene tnpR, which is related to transposons such as Tn2129, have been identified. Additionally, the rusA gene, located within a genomic island, has been shown to enhance Salmonella survival through parenteral transmission30. These determinants might collectively drive ST313 lineage shifts.

Methods
Assembly of the global iNTS genomic dataset
In this genomic epidemiology study, iNTS was defined as non-typhoidal and non-paratyphoidal Salmonella isolates obtained from human specimens of non-fecal origin, including blood, urine, abscesses, and other sterile sites, as previously defined in the literature4. We subsequently retrieved all available Salmonella genomes from EnteroBase and GenBank (data cutoff: March 1, 2024; total number: 488,232). For metadata, when BioSample information was available in EnteroBase, we obtained the corresponding metadata directly from the associated BioSample record, unless the isolate was deposited exclusively in EnteroBase without a linked BioSample entry. We excluded all strains with ambiguous isolation sources—for example, those labeled only as “clinical” or “human”—and retained only those isolates with clearly defined sources that met the above criteria for iNTS. Information on the country and date of isolation was also collected. This process resulted in a final dataset of 11,361 iNTS genomes (Supplementary Data 1).

For the global iNTS genomic dataset constructed above, we further confirmed the species identity as Salmonella using KmerFinder v3.0.234, and serovar determination was performed with SISTR v1.1.235 and SeqSero2 v1.3.136. In terms of isolation sources, most isolates were recovered from blood (n = 7,922; 69.7%) and urine (n = 2,971; 26.2%). Additional isolates primarily originated from abscesses (n = 116), wounds (n = 81), cerebrospinal fluid (n = 74), the respiratory tract (n = 36), and synovial fluid (n = 30) (Fig. S18).

Assembly of the ST313 genomic dataset
Next, we constructed a global dataset comprising 3,115 ST313 genomes (Supplementary Data 2). We not only included human-derived ST313 data but also incorporated isolates from animals, plants, and those with an "Unknown" source label. Through extensive literature searches, we confirmed that 1,584 of the 3,115 ST313 isolates have been explicitly reported in 20 previously published studies (Fig. S19). Genome quality was assessed using the European Reference Laboratory’s Salmonella genomic quality control standards37. A Python3 script was employed to exclude genomes with over 500 contigs and an N50 value less than 30,000, as well as those with a GC content outside the 51%-53% range or a genome size between 4.2 and 5.5 MB.

Furthermore, we verified the completeness and coverage of each ST313 genome. Specifically, genome coverage was assessed using QUAST v5.3.038. The tool was run with the fully sequenced D23580 genome as the reference to compute alignment-based metrics, including genome fraction, which quantifies the proportion of the reference genome aligned by the assembly contigs. Genome completeness was evaluated using BUSCO v6.039, which assesses the presence of near-universal single-copy orthologs in lineage-specific datasets. Analyses were performed in genome mode for prokaryotic sequences, using Prodigal for gene prediction. The enterobacteriaceae_odb dataset, comprising 874 marker genes, was selected to ensure taxon-specific accuracy for the bacterial group under study. To ensure the reliability of downstream analyses, including those involving MGE detection, we retained only genomes with both coverage and completeness exceeding 95% (Fig. S20).

Next, we used SISTR v1.1.2 and SeqSero2 v1.3.1 to determine the serovar of all 3,115 genomes, ensuring they were classified as S. Typhimurium. Additionally, MLST was performed using MLST v2.22 with the senterica_achtman_2 scheme40. Only strains identified as ST313 were included in subsequent analyses. Geographic visualization was performed using the online platform MapChart (https://www.mapchart.net).

Phylogenetic and population structure analyses
Core-genome single-nucleotide polymorphisms (cgSNPs) were identified among 3,115 ST313 isolates using Snippy v4.4.5. To maintain consistency with established practices in ST313 population structure definitions, the D23580 strain (ERS2705687) was utilized as the reference genome14. Then, Gubbins v2.3.441 was used to remove cgSNPs located within recombination regions from the aligned cgSNPs file, which contained a total of 9,238 loci. After filtering, 8,902 high-quality cgSNPs were retained for phylogenetic analysis. The optimal model (TVM+F+ASC+R2) was determined and used to construct the maximum-likelihood phylogenetic tree using IQ-TREE v1.6.242 with 1,000 bootstraps. The final phylogenetic tree was annotated and visualized using iTOL v6.043.

A total of 3,115 global assemblies of ST313 underwent population structure analysis using RhierBAPS v1.1.444, based on the Hierarchical Bayesian clustering algorithm and employing aligned cgSNPs as input. The calculation of population structure utilized R package v4.3.1, supplemented by the R packages "phytools" v2.0.3, "ape" v5.7.1, and "rhierbaps" v1.1.4. Specifically, the aligned ST313 cgSNP sequences were processed with the “hierBAPS” function, with the parameters set to a maximum depth of 2 (max.depth = 2) and an initial number of populations of 20 (n.pops = 20). This hierarchical approach partitions the dataset into nested clusters, with level 1 representing major lineages and level 2 defining sublineages.

Temporal signal and randomization test
We assessed the presence of significant temporal signals in the different ST313 lineages using date-randomization tests (DRT). First, recombination regions were removed from each lineage's genome using Gubbins v2.3. To identify and exclude isolates with inconsistent sampling dates, we applied TreeTime v0.1145 to perform root-to-tip regression analysis (strains whose sampling dates deviated substantially from the regression line, indicating discrepancies between genetic divergence and sampling time). These outliers were excluded before conducting the formal Bayesian Evolutionary Analysis of Sampled Trees (BEAST)46 to ensure robustness (Fig. S21). Subsequently, we used BEAUti to generate XML-formatted input files for DRT and the R package TipDatingBeast47 to randomize tip dates, producing 20 randomized replicates alongside the original dataset. Each of the 21 datasets was analyzed in BEAST2 using Markov chain Monte Carlo (MCMC) with 1,000,000 iterations and sampling every 1,000 steps. Substitution rates were estimated from all runs using BactDating v1.148. A significant temporal signal was inferred if the substitution rate from the original dataset did not overlap with those from the randomized datasets. Significant temporal signals were detected in all lineages.

Population dynamics and phylogeographical analysis of ST313
The strains in the dataset were subjected to phylogeographical reconstruction and inference of effective population sizes using BEAST v2.5. To determine the optimal model, we tested a total of six combinations in the initial phase of our study. These combinations included the strict clock, relaxed lognormal clock, and three population models (Bayesian SkyGrid, Bayesian Skyline, and Constant Size). Before conducting the full BEAST analysis, we evaluated each combination using a MCMC analysis with a total chain length of 100 million and sampling every 10,000 iterations. We then summarized the results using NSLogAnalyser and determined the optimal model based on the marginal likelihood value for each combination. The results indicated that the model incorporating the Bayesian Skyline and the relaxed lognormal clock yielded the highest marginal likelihood value in our sample. Consequently, we performed a time-calibrated Bayesian phylogenetic inference analysis for each lineage. The following settings were configured: the "GTR" substitution model, “4 gamma categories”, the "Relaxed Clock Log Normal" model, and the "Coalescent Bayesian Skyline" tree prior. Convergence was evaluated using Tracer v1.7.2, confirming that all parameters achieved effective sample sizes (ESS) greater than 200. Maximum clade credibility trees were generated using TreeAnnotator v2.6.7. Finally, key divergence time points (with 95% HPD) were estimated and the phylogenetic tree was visualized using Figtree v1.4.3. To mitigate the artefactual downward bias often observed in Bayesian skyline plots toward the end of the inferred demographic trajectory, particularly when estimating effective population size49, we limited our skyline analysis to the period preceding 2020.

SNP-distance-based transmission tracing
To minimize potential bias introduced by the choice of reference during SNP-distance calculation, we selected one complete ST313 genome assembly from each of the five lineages as the corresponding lineage-specific reference. Specifically, strain A130 was used as the reference for lineage 1, D23580 for lineage 2, BKQZM9 for lineage 3, B3589 for lineage 4, and U2 for lineage 514,50. Using the alignment method described above, we obtained cgSNP alignments, and then we calculated pairwise SNP distances between ST313 isolates with SNP-dists v0.8.2.

Furthermore, we performed complete-linkage hierarchical clustering under the optimal SNP threshold (SNP threshold = 3, as described in the Methods section below). This step was implemented in SciPy v1.14.1, ensuring that the SNP distances among isolates within each cluster did not exceed three SNPs. We also incorporated sampling time as additional epidemiological information by restricting the isolation dates of potentially related strains to within three years, a timeframe consistent with the estimated period required for the accumulation of three SNPs in ST313. Two ST313 isolates within the same cluster and collected within a three-year interval were considered to represent a potential transmission event. However, it is important to note that, given the epidemiological resolution achievable by this method, it is used as a pragmatic heuristic to summarize genetic proximity. To facilitate this analysis, we developed a custom script to count transmission events.

Determination of the optimal SNP distance threshold
We first estimated the evolutionary rate of ST313 using the methodology described by Arthur W. Pightling et al.51. Specifically, The number of new SNPs arising per year was determined by multiplying the evolutionary rate estimated with BEAST by the number of core SNP sites identified in the alignments . The results indicated that the overall evolutionary rate of ST313 is approximately 0.97 SNPs per genome per year (95% HPD: 0.75–1.53). By incorporating evolutionary rates, we estimated the number of new SNPs that may arise per year for each lineage: L1: 1.85 (95% HPD: 1.62–2.09), L2: 1.37 (95% HPD: 1.26–1.56), L3: 0.9 (95% HPD: 0.6–1.19), L4: 0.95 (95% HPD: 0.69–1.22), L5: 1.79 (95% HPD: 1.46–2.09). As shown, the annual number of new SNPs in ST313 is highly consistent across lineages, with differences of less than one SNP per year.

Furthermore, we evaluated the source-tracing performance under different SNP thresholds. Considering that an appropriate SNP distance threshold is generally assumed to correspond to the number of SNPs accumulated over approximately 2–3 years, we conducted a more detailed comparison of source-tracing performance using SNP thresholds of 1–5. The underlying assumptions were as follows: (i) the optimal threshold should approximate the number of SNPs expected to accumulate over 2–3 years in ST313; and (ii) when the threshold is set to 1 SNP, inference of potential transmission between isolates is maximally stringent, and the number of inferred potential transmission events should not increase substantially when moving to the optimal SNP threshold.

We found that as the SNP threshold increased, the number of clusters decreased, while the total number of inferred potential transmission events increased. We next compared, for thresholds of 1–5 SNPs, the mean number of potential transmission events per cluster (Fig. S6a). Notably, potential transmission events can be identified only in clusters containing 2 or more isolates; therefore, clusters comprising a single isolate were excluded. Our results show that when the threshold was set to 4 SNPs, the mean number of potential transmission events per cluster was significantly higher than with a threshold of 1 SNP (Mann–Whitney U test, FDR-corrected; p = 0.002).

We also compared, across different SNP thresholds, the mean number of within-country and between-country potential transmission events per cluster. In both cases, the mean number of potential transmission events per cluster increased with increasing SNP threshold. When we examined within-country potential transmission events across clusters, the mean number of events per cluster was significantly higher at a threshold of 4 SNPs than at 1 SNP (Mann–Whitney U test, FDR-corrected; p = 0.03, Fig. S6b). In contrast, for between-country potential transmission events, the mean number of events per cluster did not differ significantly between thresholds of 4 and 1 SNP (Mann–Whitney U test, FDR-corrected; p = 0.83). Taken together, these observations support the use of a 3-SNP threshold for the final source attribution analyses of ST313.

Construction of the SPI and plasmid database carried by ST313
A literature review was conducted to identify known plasmids and Salmonella pathogenicity Island (SPI) carried by ST31314,52–54. The plasmid sequences pBT1–pBT3 and pSLT-BT were derived from the ST313 reference strain D23580, with GenBank accession numbers LS997975.1, LS997976.1, LS997977.1, and LS997974.1. The lineage-specific pSLT plasmid reference sequence was obtained from GenBank (accession number CP060170), and the reference sequence for pAnK was NC_010896.1. All SPI reference sequences (SPI-1, SPI-2, SPI-3, SPI-4, SPI-5, SPI-6, SPI-9, SPI-11, SPI-12, SPI-13, SPI-14, SPI-16) were extracted from the D23580 reference genome (GCA_900538085.1).

Mobile genetic element detecting
Four classes of mobile genetic elements were further identified. Specifically, plasmid sequences carried by each ST313 genome were reconstructed and typed using MOB-suite v3.1.955. Then, BacAnt v3.4.056 was used to detect integrons and transposons in the genomes. Only integrons or transposons with a similarity of more than 60% and a coverage greater than 60% will be identified. The prophages were detected using the Phastest pipeline57 with default parameters.

Detection and co-localization of ARG determinants with MGE
We employed Abricate v1.0.1 with the ResFinder database58 to detect ARGs, setting minimum sequence similarity and coverage thresholds of 80% and 90%, respectively. AMR-associated point mutations were identified using Staramr v0.10.0 with default parameters59. Based on a previous study60, MDR strains are defined as those harboring AMR determinants for chloramphenicol, ampicillin, and co-trimoxazole. The resistance profiles associated with each AMR determinant were retrieved from the Card Database61. Furthermore, ARGs were considered co-localized with plasmids or prophages if they were located within their respective sequence boundaries. For integrons and transposons, given their structural characteristics and detection methods, ARGs located within 5 kb of integrase or transposase genes were also considered potentially co-localized with these MGEs28,62.

Calculation of horizontal ARG transfer frequency
To detect ARGs that may have undergone HGT, we first defined candidate HGT events as ARGs exhibiting perfect nucleotide identity (100 % coverage and 100 % identity) that were located on MGEs present in distinct Salmonella strains. We then applied the HGTphyloDetect pipeline63 to control for ARGs’ vertical inheritance by querying each ARG’s translated amino‑acid sequence against the NCBI non‑redundant protein database and calculating the Alien Index (AI) score and outgroup percentage (out_perc); the AI score was computed as follows:
AI score = In (bbhG+1⋅10−20) - In (bbhO+1⋅10−20) （1）
The bbhG and bbhO scores represent the E-values of the best BLAST hits in the ingroup and outgroup lineages, respectively; the outgroup comprises all taxa outside the focal kingdom, whereas the ingroup includes taxa within the kingdom but outside the subphylum. The out_perc was calculated as follows:
out_pect = 𝑛 / 𝑛 （2）
outside kingdom total hits
The ARGs with an AI score ≥ 45 and an out_perc ≥ 90 % are deemed putative candidates for HGT from evolutionarily distant species, and the frequency of HGT was subsequently calculated to quantify the extent of ARG dissemination as follows:
Horizontal ARG transfer Count
HGT frequency = （3）
Total dissemination count
Specific codes and examples were described in our previous study64.

Invasiveness index calculation and degradation gene analysis
The invasiveness index was calculated following previously established methods65. Specifically, the ability of ST313 to cause intestinal or extraintestinal infections in hosts correlates with the degree of degradation of the top 196 predictor genes. To identify these genes, we employed a machine-learning approach that leveraged a pre-trained random forest classifier and delta-bitscore functional variant-calling. To further validate the presence of these 196 predictor genes, which were confirmed to be associated with invasiveness, we performed BLAST analyses for each gene. This analysis produced a presence/absence matrix derived from the alignment results. Concurrently, predictor genes identified as having undergone mutations were annotated as “degradation”.

Identification of pseudogene-associated SNPs
We employed Snippy v4.4.5 to identify single-nucleotide polymorphisms in all ST313 strains. SNPs predicted to cause “Stop Gained” or “Frameshift Variant” mutations were initially classified as pseudogene-associated SNPs. To minimize potential overestimation, we further examined the positional context of each mutation within its corresponding gene. Following the conservative criterion described by Ochman et al., mutations located within the first 80% of a gene were considered likely to disrupt function, whereas those occurring in the terminal 20% were not classified as pseudogene-associated66.

Pan-genome analysis with functional annotation
For the pan-genome analysis, genome annotation was performed using Prokka v1.14.667 with a minimum contig size of 200bp and the kingdom parameter set to Bacteria. Pan-genome construction was then carried out using Panaroo v1.5.168 with default parameters. The resulting binary presence/absence matrix, indicating gene presence across strains, was analyzed using Scoary v1.6.1669 to identify lineage-specific differentially present genes. Functional annotation of key differentially present genes was conducted using eggNOG-mapper v2.1.11 with the eggNOG 5 database70. We identified a total of 30,061 coding sequences (CDS) annotated within the ST313 population. Based on gene prevalence, these CDS were categorized as core (3,608 genes), soft-core (562 genes), shell (745 genes), and cloud genes (25,146 genes) (Fig. S22).

Statistical analysis
Statistical analyses were performed using SciPy v1.14.1 to evaluate various genomic features. To assess the correlation between the data, we first determined whether the data conformed to a normal distribution using the Shapiro-Wilk test. The results indicated that none of the pairwise comparisons in this study followed a normal distribution. For pairwise comparisons, the Mann–Whitney U test was applied to determine the significance of differences between groups or conditions, and the resulting p-values were adjusted using the Benjamini–Hochberg false discovery rate (FDR) procedure. Consequently, the Spearman correlation test was applied to examine the association between the prevalence of specific genes and population shifts across ST313 lineages. Additionally, in the pan-genome analysis, Benjamini-Hochberg adjusted p-values were employed by Scoary to identify genes that were overrepresented in specific lineages. In all methods, statistical significance was determined based on a p-value threshold of 0.05.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This was a thorough analysis of a global collection of NTS genomes from different ecological and geographical regions and settings, culminating in an excellent attempt at understanding the diversity and transmission complexity and dynamics of invasive NTS disease. Authors leveraged on this unique resource to pursue three objectives: (1) to project the epidemiological landscape of iNTS, emphasising ST313; (2) to define population structure, intercontinental spread, and potential zoonotic transmission at single nucleotide resolution; and (3) to evaluate AMR and other influences on lineage replacement within ST313. The key novel findings and which have a huge implication for evolving NTS epidemiology and genomics were:\n1. NTS L1 and L2 predominated in Africa and within L2, phylogenetic analysis revealed a novel sublineage beyond 2.0–2.3, designated 2.4. Notably, most isolates in sublineage 2.4 originated from Kenya (51/61, 83.6%).\n2. Although ST313 has been considered a lineage restricted to human-to-human transmission, analysis of the 146 major clusters revealed potential cross-host transmission events in clusters 329 and 344, including human-canine transmission in cluster 329 in the UK\n3. This analysis revealed additional previously unreported lineage-specific MGEs that drove evolution of systemic NTS disease\n4.Although AMR correlated with L1 and L2 shifts, this was less evident in lineages L3-L5, which harbor few ARGs\nThe findings in these analyses clearly illustrate that NTS as a species is still evolving and the disease outcomes will largely depnd on what the bacterium picks or sheds along the evolutionary pathway. Clearly, our understanding of the forces that drive the evolution will be crucial in developing tools for mitigation against severe disease.\nA couple of issues that require author response for clarification:\nINTRODUCTION\nPage 3, Line 79: Salmonella enterica serovar Typhimurium (S. Typhimurium) is a major driver of iNTS disease6. By contrast, in sub-Saharan Africa, Malawi in particular, the S.Typhimurium ST313 (ST313) clade predominates, accounting for about 80% of cases and linking to severe outcomes such as septicemia7,8.\nComment: This is not in contrast as Malawi report S. Typhimurium as a major driver of iNTS disease.\nRESULTS\nPage 11, Line 282: Our recent study has demonstrated that the mobilome drivers of ARG dynamics in Salmonella21.\nComment: Please check on this statement as it appears incomplete\nDiscussion\nPage 14, Line 381: Our examination of AMR determinants revealed a robust correlation between L1's post-2006 decline and reduced ARG carriage. In Malawi—the epicentre of L1—the frequency of ARG horizontal transfer decreased markedly during this period.\nComment: What was the role of withdrawal of chloramphenicol as treatment option for iNTS to the loss/reduced ARGs and eventually to ST change during this period? Did the team observe any loss of other ARGs unrelated to chloramphenicol resistance? Was there a case for co-resistance and subsequent co-loss of ARGs?\nPage 14, Line 386: This coincided with Malawi's shift after 2005 to diversified antibiotics, including ceftriaxone as first-line therapy\nComment: Why would introduction of cetriaxone not cause similar selection pressure as chloramphenicol use? Is this a matter of time?\nPage 14, Line 390: In contrast, L4 and L5 harboured minimal ARGs. For these lineages, heightened invasiveness—quantified via genome degradation indices—is potentially correlated with population dynamics.\nComment: Are there specific vir genes associated with genome degradation in L4 and L5 that directly impact their survival and spread?", "reviewer_2": "The manuscript presents a large-scale genomic analysis of invasive Salmonella Typhimurium ST313, focusing on the emergence of sublineage 2.4, an emerging lineage 4, and their antimicrobial-resistance and evolutionary dynamics. This work is mostly descriptive. Clearer methodological justifications, more detailed phylogenetic analyses, and fuller discussion of limitations would improve the clarity, rigor, and interpretability of the study.\n1. L100-103: Most of these identified gaps in ST313 evolution are barely addressed or untouched by the study.\n2. L106-108: What do you mean by One Health framework? It appears the authors simply took a large collection of publicly available genomes, which are heavily biased toward human isolates. With only 0.5% animal isolates and 0.3% food isolates, how does this data set and study feature a One Health approach?\n3. The manuscript reports a novel, highly antimicrobial-resistant sublineage 2.4 and an emerging lineage 4, but provides limited phylogenetic detail. To strengthen these important findings, the authors should include more detailed phylogenetic trees (e.g., high-resolution subtrees or supplementary figures) and present metrics quantifying genetic diversity within each sublineage/lineage—for example, distributions of pairwise SNP distances or other within-group diversity measures. Such analyses would not only better support the distinction and characterization of these newly described lineages but also help confirm that sublineage 2.4 is not merely a single clonal expansion. Demonstrating sufficient genetic diversity is essential to substantiate the reported association between sublineage 2.4 and its key antimicrobial-resistance genes.\n4. The description of the dataset assembly raises some ambiguity regarding the precise definition of invasive non-typhoidal Salmonella (iNTS). The authors state that genomes were “selected from human clinical sources (excluding fecal samples),” but do not specify the exact isolation sources included. Although the authors state that fecal samples were excluded, the term “human clinical sources” is broad and may encompass specimens such as urine, or even incompletely annotated fecal specimens. Given the variability in GenBank and EnteroBase metadata, please clarify how potential ambiguous sources (e.g., samples labeled only as “clinical” or “human”) were identified and excluded. Besides, please provide detailed inclusion and exclusion criteria for sample types, and clarify whether isolates were limited to blood and other sterile sites.\n5. The methodological description of SNP analysis is a little ambiguous. It is not clear whether the core-genome SNP (cgSNP) analysis was restricted to Salmonella enterica ST313 isolates or applied to a broader iNTS dataset. In addition, the statement “a total of 9,238 cgSNPs were identified” does not specify whether this number refers to SNPs remaining after Gubbins recombination masking or to the unfiltered set produced by Snippy.\n6. The description of the temporal signal test is ambiguous regarding the criteria used to declare a significant temporal signal. While the authors state that substitution rates from the original dataset were compared with those from randomized datasets, it is unclear what statistical threshold was applied (e.g., minimum R² from root-to-tip regression, p-value from permutation tests). Without explicit criteria, readers cannot assess the strength of the evidence. The authors should specify the exact decision rule and statistical measures used to conclude that a significant temporal signal was present.\n7. The method used for transmission tracing infers putative transmission events solely from a fixed ≤3 cgSNP threshold. Without any epidemiological support/evidence, this is a questionable exercise that deviates from the established practices of genomic epidemiology. While SNP distance is an important metric, transmission inference based only on a hard cut-off is not best practice, given strain-specific evolutionary rate variation, sequencing artefacts, and potential reference-genome bias. Furthermore, The authors should provide a more detailed description of their simulations used to select this threshold, clarify how “optimal performance” was defined, and, ideally, incorporate additional epidemiological or phylogenetic evidence to corroborate genetic clustering and more reliably identify transmission event.\nAlthough some details appear in the figure legend, the explanation remains incomplete. The authors did not clearly state their justification for using a 3-SNP threshold in the main text. While they explored thresholds from 1 to 50 and observed a plateau in inferred transmission events at 3 SNPs, they did not provide formal statistical criteria for identifying transmission events or independent epidemiological validation. The authors should explicitly present their rationale in the Methods or Results and supply additional statistical, phylogenetic, or epidemiological evidence.\n8. The manuscript infers intercontinental transmission of ST313 based on genomic clustering but does not discuss the possibility that such patterns may reflect human travel or migration rather than direct bacterial transmission. This limitation should be acknowledged, and the authors should clarify that genome-based links alone cannot definitively distinguish person-to-person transmission from human-mediated dispersal across continents.\n9. References are missing or incomplete in several places. Please ensure that all methods (e.g., line 510), software tools and their versions (e.g., lines 522 and 590–597), and key claims (e.g., lines 536–539) are consistently supported by appropriate citations throughout the manuscript.\n10. Figure visualization could be improved to better highlight key findings. In particular, the annotation of important lineages and clusters is insufficient. For example, lineage 2.4 and the clusters emphasized in the text (e.g., Clusters 850, 695, 888, 579, 648) are not clearly or consistently labeled in Fig. 3.", "reviewer_3": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_4": "In this manuscript the authors leverage a large data set including >11,000 isolates from non-fecal clinical samples, as well as an additional dataset of >3,115 ST313. The authors compared isolates from different geographic regions, identifying a novel lineage (Lineage 4) and sublineage (2.4) of ST313, which has important implications for understanding the evolution of ST313, an important cause of invasive disease. The authors also use a robust series of analyses to infer relative dates of emergence and effective population sizes for lineages of ST313. Closer examination of the gene content of these various lineages also suggested a potential difference in invasiveness index for L3 isolates – although pseudogene formation as described in this manuscript, is a bit more nuanced than presence of a ‘premature stop codon,’ and I have some suggestions that could strengthen the robustness of these results. The other key finding of the study is the transmission of ARGs on mobile elements, demonstrating that plasmids and transposons are the primary drivers of ARGs in these isolates. While the study provides novel insights into the population structure of ST313 with an appreciable large dataset, I have concerns about the accuracy of the metadata, and the potential implications of errors in the metadata that would greatly change the conclusions drawn. I think overall, the authors should be careful about the interpretation of the results considering the potential inaccuracies in metadata, which unfortunately cannot be confirmed given their large dataset. I have some recommendation about interpretations of the data and some additional analyses that may help improve the robustness of the authors' findings. There are some cases where additional details would be helpful for understanding what was done and ensuring reproducibility of results.\nMajor Recommendations:\nMy biggest concern lies with the accuracy of the metadata, and the sampling biases that have the potential to greatly influence many of the results reported. The implications of this are not trivial – the authors are making claims about directionality/introduction/spread of a lineage to/from a geographic region. In some cases, an isolate may be collected in one country and then sequenced and the genomic assembly uploaded from another country; it is possible that the ‘source’ country listed in the metadata would be the latter (country where the researcher uploaded data, rather than where the isolate was actually found). Unfortunately, with this dataset, there is no way for the authors to know how accurate the metadata are without consulting the depositors (and it is unreasonable to reach out to public health agencies given that they have >3,000 isolates). The assertions made in the discussion section about ‘an imported, highly AMR branch now circulating in Africa, primarily via US-Kenya transmissions between 2009 and 2016’ would change if the metadata are not accurate, and I worry about the potential harm that could result from such assertions. I think the authors need to carefully consider ways that they can validate the accuracy of the metadata, particularly for any claims made about directionality. The other important thing to consider is that there is undoubtedly bias in relying on WGS data rather than reported case data. It is unfair to assume that cases in each country have the same likelihood of having isolates characterized with WGS. I appreciate that the authors are trying to use the genomic data to uncover helpful information about the spread, but I find the data presentation and conclusions drawn potentially misleading. For example, the temporal analyses in panel 1b and the supplemental data, require equal sampling across all time points. We know this is not true – there were several large, well-funded studies that specifically sought to understand iNTS burden in select countries in Africa. I believe that what is shown in the data in Figure 1B reflects a reduction in focused efforts, rather than a true signal of reductions in ST313 infections (which is what is being assumed by the authors). To address this, the authors should carefully consider reframing their conclusions about decline in ST313 genome sequences over time as this likely reflects projects ending rather than a true decline in infections.\nMinor Recommendations:\nLines 116-131: I think the way that iNTS was defined from the metadata could be problematic and is a key limitation that needs to be discussed; were isolates that had an ‘unknown’ isolation source included? Why did the authors not limit their analyses to just isolates from blood? Most clinical samples do not list the isolation source. How might variations in clinical diagnosis and reporting bias the results in the genomic dataset? In the materials and methods ‘iNTS’ is defined as ‘not from feces’ – in this case most of the dataset should have included Typhi and Paratyphi strains as well but none are reported? Was additional filtering done aside from just ‘excluding feces as a source?’ I think the details about the way the iNTS dataset are missing to allow complete understanding of how these data were collected.\nI see the potential for the reporting of these cases to be biased by public health reporting differences by country. For example, countries with greater public health infrastructure and resources will appear to have higher incidences than those with lower resources. This is an important source of bias that is partially addressed in the discussion.\nLines 129-131: In figure 1B the authors show a marked decline in sequences from Africa after 2014 – given that the authors’ assertion that ST313 is overrepresented among clinical cases in Africa, doesn’t this instead suggest a decrease in genomic sequencing efforts rather than a true decline in ST313 infections? Logically, it does not make sense that ST313 would magically decline; rather, this again suggests sampling bias. This is a limitation in the dataset that needs to be discussed. This is important because these data are from WGS (not just reported cases).\nLines 123-124: How were serovars determined? In silico methods? Serological? In the Materials and Methods it states that SISTR and SeqSero2 were used for the ST313 dataset. Why were these methods not used for the larger 11,631 publicly available dataset as well? Enterobase already runs these analyses on all of their data – it would be nice to ensure that the first genomic dataset also included SISTR and SeqSero2 serovar classification.\nLines 133-141: It is appreciably challenging to gather accurate data from voluntarily reported metadata and I would recommend that the authors account for this in their manuscript. For example, panel 1e: sex data is only provided for 1/5th of the human isolates. Without understanding how these data were collected and by whom, this is potentially misleading.\nLine 152: How was 2.4 identified and what cut-offs were used to distinguish it from the other sub-lineages? This is important for ensuring reproducibility of the results.\nLines 166 and throughout: If BEAST was used, shouldn’t these statistics be HPDs rather than 95% confidence intervals as listed throughout the manuscript?\nLines 191-192: Couldn’t this just be due to international travel?\nLines 198-199: Yes, but there is also bias in the sequencing. We don’t often look for or sequence Salmonella from sick animals unless there is a specific research project looking at this. I think it is safer for the authors to qualify this in their discussion.\nLine 206: As written, it suggests that only pSLT was used in query, but then in line 210 the authors mention that pBT2 and pBT3 were also identified. Please rephrase line 206 so that it is clearer what was queried.\nLines 226-230: Did the authors control for low-sequenced regions by using thresholds for coverage across the genome or in a particular region of the genome? Mobile elements are inherently challenging to sequence and assemble from short-read sequencing data. Did the authors check to see that these were indeed absent? Or is it possible that these regions of the genome had low coverage at these sites and this could explain why they were not detected?\nLine 355: I think caution should be exercised here when claiming that these are ‘spreading around the world.’ Without a more complete dataset, it is not possible to assert this. All that can be concluded is that isolation of ST313 happens outside of Africa as well.\nLine 382: Using terms such as ‘the epicenter of L1’ could be misleading. The authors need to remember that there were targeted efforts to sample and sequence isolates from Malawi (over other African countries).\nLines 398-400: Similar to my comments above, I think it is quite the opposite. There has been an intentional focus on sampling iNTS in African countries. I think this bias is apparent in the data and should be discussed in more detail here. This also contradicts what the authors state in lines 403-404.\nLine 435: Please justify the use of this strain as reference for all genomes. It appears this is most suitable for Lineage 2 ST313 strains but its use as a reference for other lineages may not be valid and should be justified. The easiest way to address this would be to use the most complete genome as reference for each of the other lineages as reference for mapping SNPs among isolates in that lineage.\nLines 574-575: Details about the ‘pre-trained random forest classifier’ need to be provided. Did the authors use an already developed classifier from another paper? Please cite that paper to be fair.\nLines 582-587: This approach may identify PMSCs that would not render a gene a pseudogene. Additional filtering that considers the position of the frameshift is important and necessary for identification of pseudogenes. For example, PMSCs occurring at the beginning of a gene are more likely to result in loss of functon(pseudogene formation) but PMSCs at the very end of a gene may not influence function. Additional filtering based on position is a way for reducing the likelihood (but not removing the possibility) that a gene would be incorrectly identified as a pseudogene.\nLine 840: Sex (male/female) is a more appropriate term here than gender (man/woman)\nFigure 1d: Unknown is spelled wrong; I also question the inclusion of these strains in this dataset. Given that many Salmonella infections are diagnosed from fecal samples, it is likely that at least some of these strains indeed came from feces. Also, the authors stated that ‘feces’ was excluded, but their second most common known source is ‘faeces’ which is just a different spelling of ‘feces.’ Therefore, this dataset likely contains non-invasive infections as well. It is also unclear to me why fecal isolates are included in this dataset, whereas in the larger iNTS dataset they are not?\nFigure 4a: In the pseudo heat map, it looks like there is perfect 0 or 100% correlation between the isolates having a given element (BTP3 for example) or not. Yet in the text percentages are listed (11.8% of L1 carry BTP3). As presented in the figure, it looks like none of the isolates have BTP3. Additional gradients should be used as presented the data in the text does not match the figure. The same applies to Figure S9.\nIn panel 4c, it looks like the comparison is with other subspecies (subsp.)? Please confirm that this should indeed be labeled ‘other subsp.’ and not other lineages.\nFigure 5c: It looks like there was no variance in mcr-9, tet(D), or aph(3’-Ia) in average number of ARGs over time? Is this true? Or the CI overlaps with the line?\nFigure S10: AMR determinants not Determins\nFigure S14: Was correction for multiple comparisons testing done? I suspect not as the SE bars for L4 overlap with L3, yet the comparison suggests that the difference in pseudogene SNPs is significant."}, "author_response": "Response: We sincerely thank you for your thorough and highly positive evaluation of our manuscript. We are encouraged by your recognition that you point out this analysis as “thorough” and the implications of our findings as “huge” for iNTS epidemiology and genomics. In response to your insightful comments, we have carefully addressed your comments point-by-point, and corresponding revisions are highlighted in green in the manuscript.\nINTRODUCTION\nPage 3, Line 79: Salmonella enterica serovar Typhimurium (S. Typhimurium) is a major driver of iNTS disease6. By contrast, in sub-Saharan Africa, Malawi in particular, the S. Typhimurium ST313 (ST313) clade predominates, accounting for about 80% of cases and linking to severe outcomes such as septicemia7,8.\nComment: This is not in contrast as Malawi report S. Typhimurium as a major driver of iNTS disease.\nResponse: Thank you for this helpful comment. We have revised this sentence to ensure for accuracy and clarity. The updated version now reads as follows:\nNew lines 79-83: Salmonella enterica serovar Typhimurium (S. Typhimurium) emerges as a primary etiological agent of iNTS disease6. In sub-Saharan Africa, the S. Typhimurium sequence type 313 (ST313) clade predominates, accounting for a substantial proportion of cases and being associated with severe manifestations such as septicemia7,8.\nRESULTS\nPage 11, Line 282: Our recent study has demonstrated that the mobilome drivers of ARG dynamics in Salmonella21.\nComment: Please check on this statement as it appears incomplete\nResponse: We have revised the sentence accordingly; it now reads as follows:\nNew lines 312-313: Our recent study has demonstrated that the mobilome drivers play a critical role in ARG dynamics in Salmonella24.\nDiscussion\nPage 14, Line 381: Our examination of AMR determinants revealed a robust correlation between L1's post-2006 decline and reduced ARG carriage. In Malawi—the epicentre of L1—the frequency of ARG horizontal transfer decreased markedly during this period.\nComment: What was the role of withdrawal of chloramphenicol as treatment option for iNTS to the loss/reduced ARGs and eventually to ST change during this period?\nResponse: Thank you for raising this point. Chloramphenicol was withdrawn as the empirical treatment for ST313 around 2001; and it was subsequently replaced by the oral fluoroquinolone ciprofloxacin (Nat Microbiol. 2021; DOI: 10.1038/s41564-020-00836-1). In this study, we observed two declines in the L1 population: the first around 2001 and the second around 2006.\nChloramphenicol was withdrawn as the empirical treatment for ST313 around 2001, primarily due to the emergence of a well-documented chloramphenicol-resistant L2 epidemic (Nat. Genet. 2012; DOI: 10.1038/ng.2423). As shown in Response Figure 1, the catA1 gene, which confers chloramphenicol resistance, was present in 88% of L2 isolates.\nResponse Figure 1: The carriage rates of chloramphenicol-related and ceftriaxone-related ARGs among different ST313 lineages are shown. The numbers in the squares indicate the specific carriage rates, with deeper red colors representing higher prevalence.\nTo address your query regarding the impact of chloramphenicol withdrawal—particularly after 2001 (although in Line 381 we discuss the decline of ST313 L1 after 2006)—on ARGs in ST313 L1, we stratified the 287 L1 isolates into two groups (pre-2001 and post-2001). Because some ARGs had expected counts below 5, we used Fisher’s exact test (instead of a Chi-square test) to obtain accurate p-values for small samples. We identified ARGs in L1 that showed significant declines (p < 0.05) after chloramphenicol was discontinued, with p-values adjusted using the Benjamini–Hochberg false discovery rate (FDR) procedure. However, no chloramphenicol resistance gene in L1 showed a significant reduction in carriage. Instead, the carriage rates of dfrA23 (conferring trimethoprim resistance) and ant(2\")-Ia (conferring aminoglycoside resistance) declined significantly (Response Table 1).\nResponse Table 1: Significantly lost antibiotic resistance genes in ST313 L1 before and after chloramphenicol was used as an empirical treatment.\nCarriage rate Carriage rate\nARGs p_adj_Value\n(before 2001) (after 2001)\ndfrA23 0.865 0.627 0.000041\nant(2'')-Ia 0.825 0.689 0.039569\nNotably, the carriage rate of catA2 (a gene conferring chloramphenicol resistance) increased significantly in L1 after 2001, suggesting that the initial decline of L1 around this time was unlikely to be driven by the loss of ARGs. Instead, it may have resulted from the expansion of the L2 lineage, which exhibits high-level chloramphenicol resistance, thereby contributing to the first disappearance of L1 around 2001.\nDid the team observe any loss of other ARGs unrelated to chloramphenicol resistance?\nResponse: Yes. We examined how changes in antimicrobial resistance might relate to the second decline of ST313 after 2006. We found that six key ARGs were lost from ST313 L1, including mcr-9, tet(D), aph(3')-Ia, sul2, aph(6)-Id, and aph(3\")-Ib, which correlated with the post-2006 decline of L1. We hypothesize that these losses are linked to shifts in antimicrobial usage (for example, mcr-9 confers resistance to colistin, a \"last-resort\" polymyxin antibiotic). This finding suggests that ST313 continues to evolve, likely in response to antimicrobial stewardship strategies.\nWas there a case for co-resistance and subsequent co-loss of ARGs?\nResponse: Previous studies have shown that resistance plasmids can drive co-resistance and co-loss of ARGs in Salmonella (mSystems. 2023; DOI: 10.1128/msystems.00883-23). To determine if similar phenomena occur in L1, we examined the temporal dynamics of 54 plasmid types in L1 over the years. We used Fisher’s exact tests (p < 0.05) to identify significant changes in the carriage rate of each plasmid between years. Our analysis revealed 19 plasmid types (out of 54) with significant changes in carriage. Notably, the IncHI2A/rep_cluster_1088 plasmid increased in prevalence between 2000 and 2005, then showed a significant decline around 2006 (Response Figure 2).\nResponse Figure 2: Temporal trends in the gain or loss of 19 plasmids showing significant changes in carriage rates within ST313 L1. In the heatmap, blue indicates plasmids that were significantly lost in ST313 L1 between the compared years, red denotes plasmids that were significantly acquired, and grey represents plasmids with no significant change in carriage rate across the corresponding years.\nAs noted in our manuscript, IncHI2A/rep_cluster_1088 plasmid carries several key ARGs, tet(A), bla , dfrA14, bla , and bla (Supplementary Figure S13).\nTEM CTX-M-15 OXA\nSupplementary Figure S13: The genetic composition of the IncHI2A/rep_cluster_1088-type resistance plasmid. Genome annotation was performed using Bakta, and the IncHI2A/rep_cluster_1088-type resistance plasmid was further annotated with specific color codes: ARGs were highlighted in red, transposons in blue, and prophage-related elements in green.\nTogether, these findings suggest that the ARGs on the IncHI2A/rep_cluster_1088 plasmid underwent co-selection and subsequent co-loss, consistent with your comment.\nPage 14, Line 386: This coincided with Malawi's shift after 2005 to diversified antibiotics, including ceftriaxone as first-line therapy\nComment: Why would introduction of cetriaxone not cause similar selection pressure as chloramphenicol use? Is this a matter of time?\nResponse: Thank you for these insightful questions. Chloramphenicol was the primary empirical treatment before 2001, and ST313 Lineage 2 carried a much higher prevalence of chloramphenicol resistance genes (especially catA1) than Lineage 1 (Response Figure 1). Under chloramphenicol selection pressure, this likely allowed L2 to replace L1 as the predominant lineage.\nResponse Figure 1: The carriage rates of chloramphenicol-related and ceftriaxone-related ARGs among different ST313 lineages are shown. The numbers in the squares indicate the specific carriage rates, with deeper red colors representing higher prevalence.\nIn contrast, after ceftriaxone was introduced (post-2005), no similar replacement pattern was observed. Lineages L3, L4, and L5 carry relatively few ARGs (Supplementary Figure S10). While L2 has an extremely high carriage (~95%) of bla (which confers ceftriaxone resistance; Response Figure 2), the selective pressure of ceftriaxone on L2 is much weaker than chloramphenicol's effect on L1. We infer that the lack of lineage replacement after ceftriaxone's introduction is not simply a matter of insufficient time but also due to L2's intrinsic ceftriaxone resistance, which diminishes any additional selective advantage.\nSupplementary Figure S10: ARG determinants in 3,115 Salmonella ST313 isolates. A phylogenetic tree was constructed based on core-genome single-nucleotide polymorphisms from these isolates, with left branches categorized into distinct lineages using BAPs analysis. Each lineage is indicated by a unique color. The heatmap on the right illustrates the geographic distribution of ST313 isolates by continent, the presence of resistance-associated point mutations in genes (gyrA, gyrB), and the corresponding ARGs.\nPage 14, Line 390: In contrast, L4 and L5 harboured minimal ARGs. For these lineages, heightened invasiveness—quantified via genome degradation indices—is potentially correlated with population dynamics.\nComment: Are there specific vir genes associated with genome degradation in L4 and L5 that directly impact their survival and spread?\nResponse: Thank you for this insightful question. Given the close relationship between gene degradation—especially the loss of virulence genes—and the enhanced invasiveness of ST313, we fully agree that further investigation of virulence gene degradation is necessary and valuable.\nIn the original manuscript, genes containing “Stop Gained” or “Frameshift Variant” SNPs were classified as potential pseudogenes. However, in response to Reviewer #4’s comment — “Additional filtering that considers the position of the frameshift is important and necessary for identification of pseudogenes.” — we refined our detection approach to address this limitation.\n(i) Previous research by Howard Ochman et al. showed that mutations occurring within the first 80% of a gene are more likely to result in pseudogene formation, a conservative indicator of gene inactivation (Nucleic Acids Res. 2005; DOI: 10.1093/nar/gki631). We incorporated this positional filtering strategy into our pseudogene identification pipeline, and we have added information on the relative position of each “Stop Gained” or “Frameshift” mutation in Supplementary Table 8. Specifically, mutations within the last 20% of a gene are not classified as pseudogene-associated SNPs. Using this stricter criterion, we removed 7,344 SNPs from the original Supplementary Table 8 that did not meet the pseudogene definition. In the revised version, we also report the exact position of each SNP within its gene (e.g., “533/722” indicates a 722 bp gene with the SNP at position 533).\n(ii) The Virulence Factor Database (VFDB) is a comprehensive repository of bacterial virulence genes. By comparing our identified pseudogenes with the VFDB, we identified 26 virulence-associated pseudogenes. The carriage rate of virulence-related pseudogenes was 9.87% (15/152) in Lineage 5 and 23.9% (55/230) in Lineage 4. Notably, steC harbored pseudogene-associated SNPs more frequently in Lineage 4, whereas sseG harbored them more frequently in Lineage 5 (Response Figure 3).\nResponse Figure 3: The potential virulence-associated pseudogenes carried by Lineage 4 and Lineage 5. The names of the strains harboring these virulence-related pseudogenes, along with their corresponding lineages, are shown on the left side of the heatmap. Dark blue blocks indicate the presence of virulence-associated pseudogenes, while light blue blocks indicate their absence.\n(iii) Interestingly, both steC and sseG are associated with Salmonella Pathogenicity Island 2 (SPI-2). Previous studies report that SteC is an SPI-2–encoded effector (the only kinase effector) that helps Salmonella breach the intestinal vascular barrier and disseminate to organs (EMBO J. 2024; DOI: 10.1038/s44318-024-00076-7). SseG, another SPI-2 effector, aids the establishment of a replicative niche by inhibiting host autophagy (J. Biol. Chem. 2018; DOI: 10.1074/jbc.M117.811737). Therefore, we speculate that mutations in these two genes may influence the survival and spread of Lineage 4 and Lineage 5.\nReviewer #2 (Remarks to the Author):\nThe manuscript presents a large-scale genomic analysis of invasive Salmonella Typhimurium ST313, focusing on the emergence of sublineage 2.4, an emerging lineage 4, and their antimicrobial-resistance and evolutionary dynamics. This work is mostly descriptive. Clearer methodological justifications, more detailed phylogenetic analyses, and fuller discussion of limitations would improve the clarity, rigor, and interpretability of the study.\nResponse: Many thanks for your helpful comments. We appreciate your suggestions to improve the clarity, rigor, and interpretability of the study. In response, we have implemented several important revisions. Specifically, we have added clearer methodological justifications to better explain our approach. Additionally, we have expanded the phylogenetic analyses to provide more in-depth insights into the evolutionary dynamics of S. Typhimurium ST313. We have also provided a more thorough discussion of our study's limitations to improve transparency. We hope these changes address your concerns and enhance the overall quality of the manuscript. Please see our point-by-point responses below. All corresponding revisions in the manuscript are highlighted in blue.\n1. L100-103: Most of these identified gaps in ST313 evolution are barely addressed or untouched by the study.\nResponse: We apologize for the inappropriate description in the Introduction. The corresponding statement has been completely revised, and the revised paragraph now reads:\nNew lines 90-106: Emergence of dominant bacterial variants, including within ST313, occurs through mutational processes, horizontal gene transfer, and gene loss, often under selective pressures favoring traits like increased transmissibility, virulence, or immune evasion11. Prior research underscores the pivotal role of antimicrobial resistance (AMR), mediated by acquired resistance genes (ARGs), in lineage dynamics and replacement within ST31312. For instance, chloramphenicol resistance, largely restricted to L2, likely facilitated L1 displacement by L2 around 200113. Since 2005, amid rising use of third-generation cephalosporins and azithromycin, ST313 isolates in Africa—particularly L2 strains—have exhibited adaptive signatures of resistance to first-line antibiotics14,15. Beyond AMR, progressive genome degradation—evidenced by the accumulation of pseudogenes and alterations in functional gene expression—has been implicated in enhancing the invasive potential of ST31316–18. These findings highlight existing knowledge gaps in constructing a comprehensive genomic landscape of ST313 to understand its evolutionary trajectory, the driving forces behind its adaptation and the design of precise disease interventions.\n2. L106-108: What do you mean by One Health framework? It appears the authors simply took a large collection of publicly available genomes, which are heavily biased toward human isolates. With only 0.5% animal isolates and 0.3% food isolates, how does this data set and study feature a One Health approach?\nResponse: Thank you for your valuable comment, and we apologize once again for the overstatement in the manuscript. In the revised manuscript, we have removed the overgeneralization about a \"One Health framework.\" We agree that the limited number of animal samples (n = 15, 0.5%) and food samples (n = 8, 0.3%) in our ST313 dataset is insufficient to support a true One Health perspective. This is consistent with the current limited surveillance of ST313, which remains primarily associated with human infections. Given the scarce genomic data available, it is challenging to determine whether cross-host transmission occurs in ST313. In our revision, we removed the description of the ST313 genomic dataset as a “One Health framework” Our aim is to assemble as diverse an ST313 strain collection as possible to investigate potential cross-host transmission, although such events appear to be exceedingly rare in ST313.\n3. The manuscript reports a novel, highly antimicrobial-resistant sublineage 2.4 and an emerging lineage 4, but provides limited phylogenetic detail. To strengthen these important findings, the authors should include more detailed phylogenetic trees (e.g., high-resolution subtrees or supplementary figures) and present metrics quantifying genetic diversity within each sublineage/lineage—for example, distributions of pairwise SNP distances or other within-group diversity measures. Such analyses would not only better support the distinction and characterization of these newly described lineages but also help confirm that sublineage 2.4 is not merely a single clonal expansion. Demonstrating sufficient genetic diversity is essential to substantiate the reported association between sublineage 2.4 and its key antimicrobial-resistance genes.\nResponse: Thank you for your comment. In the revised manuscript, we have included additional high-resolution subtrees for sublineage 2.4 (the high-resolution subtree for L4 is provided in Fig. 2c). We have also quantitatively assessed the within-group genetic diversity of both sublineage 2.4 and L4 using SNP distance as a metric.\nTo address potential reference bias (your Comment 7 below), we used D23580 as the reference strain for sublineage 2.4 and B3589 for L4 (both are complete ST313 genome assemblies from the respective lineages), respectively. We then employed snp-dists to calculate pairwise SNP distances from recombination-filtered core-genome SNP alignments for both lineages. Our results show that pairwise SNP distances within sublineage 2.4 range from 0 to 50 nucleotides (mean = 19.1), whereas L4 shows greater diversity, with distances from 0 to 510 nucleotides (mean = 173.4) (New Fig. S2).\nGiven that L4 is a distinct lineage while sublineage 2.4 is within Lineage 2, the higher genetic diversity in L4 is expected. These results support that neither sublineage 2.4 nor L4 arose from a single clonal expansion.\nNew Fig. S2. High-resolution subtree of the newly identified sublineage 2.4 and comparison of within-lineage genetic diversity between sublineage 2.4 and L4. (a) A phylogenetic tree of sublineage 2.4 was constructed using D23580 as the reference. The heatmaps below the tree indicate the continent of origin and the source of different strains, respectively. (b) Within-lineage genetic diversity was compared between sublineage 2.4 and L4 using SNP distance metrics. The x-axis represents SNP distance, while the y-axis represents the frequency of pairwise genetic distances within the corresponding SNP distance range.\nWe have revised the manuscript accordingly.\nNew lines 163-169: To further assess whether the newly identified sublineage 2.4 and L4 represent single clonal expansion events, we calculated their within-lineage genetic diversity. Pairwise SNP distances within sublineage 2.4 ranged from 0 to 50 nucleotides (mean: 19.1), whereas L4 showed substantially greater divergence, ranging from 0 to 510 nucleotides (mean: 173.4). These findings indicate that neither sublineage 2.4 nor L4 arose from a single clonal expansion event (Fig. S2).\n4. The description of the dataset assembly raises some ambiguity regarding the precise definition of invasive non-typhoidal Salmonella (iNTS). The authors state that genomes were “selected from human clinical sources (excluding fecal samples),” but do not specify the exact isolation sources included. Although the authors state that fecal samples were excluded, the term “human clinical sources” is broad and may encompass specimens such as urine, or even incompletely annotated fecal specimens. Given the variability in GenBank and EnteroBase metadata, please clarify how potential ambiguous sources (e.g., samples labeled only as “clinical” or “human”) were identified and excluded. Besides, please provide detailed inclusion and exclusion criteria for sample types, and clarify whether isolates were limited to blood and other sterile sites.\nResponse: We apologize for the insufficient description in the original manuscript. The definition of invasive non-typhoidal Salmonella (iNTS) used in this study follows the literature (Lancet Infect Dis. 2019; DOI: 10.1016/S1473-3099(19)30418-9). According to current literature, unlike NTS, which cause self-limiting gastrointestinal illness, iNTS refers to Salmonella isolates from normally sterile sites (blood, cerebrospinal fluid, etc.), leading to bacteremia, meningitis, and other invasive infections. Based on this, in our study we defined iNTS as non-typhoidal, non-paratyphoidal Salmonella isolates from human specimens of non-fecal origin (e.g., blood, urine, abscesses, CSF, and other sterile sites).\nWe acknowledge variability in metadata between EnteroBase and GenBank entries. Before assembling our dataset, we consulted Dr. Zhemin Zhou (Previous EnteroBase co-creator, now a professor at the Soochow University). Our approach was: if BioSample information was available in EnteroBase, we retrieved metadata from the BioSample record, unless an isolate was only in EnteroBase without a BioSample. This strategy helped avoid inconsistencies between GenBank and EnteroBase metadata. Importantly, we excluded all isolates with ambiguous sources (e.g., labeled only “clinical” or “human”), as noted by the reviewer.\nWe have expanded Supplementary Table 1 to list each strain’s isolation source and have updated the collation. Among the 11,361 iNTS isolates in our global dataset, 23 distinct isolation sites were identified. The majority were obtained from blood (n = 7,922; 69.7%) and urine (n = 2,971; 26.2%). Additional isolates came from abscesses (n = 116), wounds (n = 81), cerebrospinal fluid (n = 74), the respiratory tract (n = 36), synovial fluid (n = 30), etc. (New Fig. S18).\nNew Fig. S18: Sample isolation sites for the global iNTS genome dataset comprising 11,361 isolates. (a) All strains in the dataset were isolated from humans. (b) Specific isolation sites for the strains, with the majority recovered from blood (n = 7,922; 69.7%) and urine (n = 2,971; 26.2%). (c) A breakdown of isolates from other sources, as indicated in panel b.\nWe have thoroughly revised the “Database assembly” section of Methods. It now reads:\nNew lines 478-500:\nAssembly of the global iNTS genomic dataset\nIn this genomic epidemiology study, iNTS was defined as non-typhoidal and non-paratyphoidal Salmonella isolates obtained from human specimens of non-fecal origin, including blood, urine, abscesses, and other sterile sites, as described in previous literature2. We subsequently retrieved all available Salmonella genomes from EnteroBase and GenBank (data cutoff: 1 March 2024; total number: 488,232). For metadata, when BioSample information was available in EnteroBase, we obtained the corresponding metadata directly from the associated BioSample record, unless the isolate was deposited exclusively in EnteroBase without a linked BioSample entry. We excluded all strains with ambiguous isolation sources—for example, those labeled only as “clinical” or “human”—and retained only those isolates with clearly defined sources that met the above criteria for iNTS. Information on the country and date of isolation was also collected. This process resulted in a final dataset of 11,361 iNTS genomes (Supplementary Table 1).\nFor the global iNTS genomic dataset constructed above, we further confirmed the species identity as Salmonella using KmerFinder v.3.0.230, and serovar determination was performed with SISTR v.1.1.231 and SeqSero2 v.1.3.132. In terms of isolation sources, most isolates were recovered from blood (n = 7,922; 69.7%) and urine (n = 2,971; 26.2%). Additional isolates primarily originated from abscesses (n = 116), wounds (n = 81), cerebrospinal fluid (n = 74), the respiratory tract (n = 36), and synovial fluid (n = 30) (Fig. S18).\n5. The methodological description of SNP analysis is a little ambiguous. It is not clear whether the core-genome SNP (cgSNP) analysis was restricted to Salmonella enterica ST313 isolates or applied to a broader iNTS dataset.\nResponse: Thank you for your comment. In our study, the cgSNP analysis was restricted to 3,115 ST313 isolates. We have clarified this in the revised manuscript, which now reads:\nNew lines 535-538: Core-genome single-nucleotide polymorphisms (cgSNPs) were identified among 3,115 ST313 isolates using Snippy v4.4.5. To maintain consistency with established practices in ST313 population structure definitions, the D23580 strain (ERS2705687) was utilized as the reference genome10.\nIn addition, the statement “a total of 9,238 cgSNPs were identified” does not specify whether this number refers to SNPs remaining after Gubbins recombination masking or to the unfiltered set produced by Snippy.\nResponse: In this study, we first used Snippy to identify cgSNPs among 3,115 ST313 isolates, using D23580 as the reference. The aligned cgSNP set contained 9,238 raw cgSNP loci. Gubbins v2.3.4 was then used to remove cgSNPs in recombination regions, resulting in 8,902 high-quality cgSNPs retained for phylogenetic analysis. We clarified this in the revised manuscript, which now reads:\nNew lines 538-541: Then, Gubbins v2.3.437 was used to remove cgSNPs located within recombination regions from the aligned cgSNPs file, which contained a total of 9,238 loci. After filtering, 8,902 high-quality cgSNPs were retained for phylogenetic analysis.\n6. The description of the temporal signal test is ambiguous regarding the criteria used to declare a significant temporal signal. While the authors state that substitution rates from the original dataset were compared with those from randomized datasets, it is unclear what statistical threshold was applied (e.g., minimum R² from root-to-tip regression, p-value from permutation tests). Without explicit criteria, readers cannot assess the strength of the evidence. The authors should specify the exact decision rule and statistical measures used to conclude that a significant temporal signal was present\nResponse: We apologize for the insufficient detail in the original manuscript. In the revision, we have provided a more complete explanation. We assessed the temporal signal in each dataset using date-randomization tests (DRT), a sensitive approach that accounts for rate variation. Specifically, for each lineage, we employed the optimized BEAST2 model parameters, including the GTR substitution model, four gamma categories, a Relaxed Clock Log Normal model, and the Coalescent Bayesian Skyline tree prior. Using BEAUti and the R package TipDatingBeast, we generated 20 replicate datasets with randomized tip dates, along with the original dataset. Each was run in BEAST2 for 1,000,000 MCMC iterations (sampling every 1,000). We then used BactDating to estimate substitution rates from all 21 runs. We considered a temporal signal significant if the substitution rate from the real dataset fell outside the range of rates from all 20 randomized datasets. By this criterion, all six lineages showed a significant temporal signal (Response Figure 4).\nResponse Figure 4. Comparison of the substitution rates. We estimated the substitution rates for the real dataset (red line) and 20 randomized datasets using BEAST2 and BactDating. The log₁₀-transformed substitution rates from each analysis were visualized for comparison. A significant temporal signal was considered present when the substitution rate estimated from the real dataset did not overlap with those obtained from any of the 20 randomized datasets.\nWe have revised the “Temporal signal and randomization test” section in the Methods accordingly.\nNew lines 558-575:\nTemporal signal and randomization test\nWe assessed the presence of significant temporal signals in the different ST313 lineages using date-randomization tests (DRT). First, recombination regions were removed from each lineage's genome using Gubbins v2.3. To identify and exclude isolates with inconsistent sampling dates, we applied TreeTime v0.1141 to perform root-to-tip regression analysis (strains whose sampling dates deviated substantially from the regression line, indicating discrepancies between genetic divergence and sampling time). These outliers were excluded before conducting the formal Bayesian Evolutionary Analysis of Sampled Trees (BEAST)42 to ensure robustness (Fig. S21). Subsequently, we used BEAUti to generate XML-formatted input files for DRT and the R package TipDatingBeast43 to randomize tip dates, producing 20 randomized replicates alongside the original dataset. Each of the 21 datasets was analyzed in BEAST2 using Markov chain Monte Carlo (MCMC) with 1,000,000 iterations and sampling every 1,000 steps. Substitution rates were estimated from all runs using BactDating v1.144. A significant temporal signal was inferred if the substitution rate from the original dataset did not overlap with those from the randomized datasets. Significant temporal signals were detected in all lineages.\n7. The method used for transmission tracing infers putative transmission events solely from a fixed ≤3 cgSNP threshold. Without any epidemiological support/evidence, this is a questionable exercise that deviates from the established practices of genomic epidemiology. While SNP distance is an important metric, transmission inference based only on a hard cut-off is not best practice, given strain-specific evolutionary rate variation, sequencing artefacts, and potential reference-genome bias. Furthermore, the authors should provide a more detailed description of their simulations used to select this threshold, clarify how “optimal performance” was defined, and, ideally, incorporate additional epidemiological or phylogenetic evidence to corroborate genetic clustering and more reliably identify transmission event. Although some details appear in the figure legend, the explanation remains incomplete. The authors did not clearly state their justification for using a 3-SNP threshold in the main text. While they explored thresholds from 1 to 50 and observed a plateau in inferred transmission events at 3 SNPs, they did not provide formal statistical criteria for identifying transmission events or independent epidemiological validation. The authors should explicitly present their rationale in the Methods or Results and supply additional statistical, phylogenetic, or epidemiological evidence.\nResponse: We agree that selecting an appropriate SNP threshold requires careful consideration of sequencing artifacts, strain-specific rate variation, and reference bias. We also agree that incorporating epidemiological data is essential for reliable transmission inference. In general, a lower SNP-distance threshold is more stringent for defining transmission, but there is no universal consensus on the optimal cut-off (Mol. Biol. Evol. 2019; 36:587–603). For context, several recent genomic epidemiology studies have used fixed SNP thresholds to infer Salmonella transmission, for example: 3 SNPs (J. Infect. 2025; DOI: 10.1016/j.jinf.2025.106469), 4 SNPs (Nat. Med. 2025; DOI: 10.1038/s41591-025-03644-4), and 5 SNPs (NAR Genom Bioinform 2022; DOI: 10.1093/nargab/lqac047).\nTo fully address this concern, we made several revisions:\n(1) Genome sequencing quality assessment: We performed rigorous quality control on all 3,115 ST313 genomes (assessing GC%, N50, size, coverage, completeness). We established thresholds based on European Reference Laboratory standards and prior genomic epidemiology practice (Response Table 2). All ST313 genomes met these criteria (New Fig. S20), ensuring that sequencing quality did not skew SNP distances.\nResponse Table 2. Genome assembly quality-control thresholds used in this study.\nQuality Control Parameters Threshold\nGC (%) 51 - 53\nN50 > 30,000\nGenome size (MB) 4.2 - 5.5\nGenome coverage (%) > 95\nGenome completeness (%) > 95\nNew Fig. S20: Quality control results of 3,115 ST313 genomes. (a) GC content; (b) N50; (c) genome size; (d) genome coverage; and (e) genome completeness.\n(2) Reference genome selection: To minimize potential bias introduced by reference selection during SNP-distance calculation, we used one complete ST313 genome assembly from each of the five lineages as the corresponding lineage-specific reference. In choosing the lineage-specific references, we adhered closely to the established epidemiological practices for ST313(Nat Microbiol. 2021; DOI: 10.1038/s41564-020-00836-1, Infect Genet Evol. 2019: DOI: 10.1016/j.meegid.2019.05.023). Specifically, strain A130 was used as the reference for lineage 1, D23580 for lineage 2, BKQZM9 for lineage 3, B3589 for lineage 4, and U2 for lineage 5.\n(3) SNP threshold based on evolutionary rates: We fully agree that threshold selection should consider the pathogen’s evolutionary rate. We estimated ST313’s rate following Pightling et al. (Front. Microbiol. 2022; DOI: 10.3389/fmicb.2022.797997). The overall ST313 rate is ~0.97 SNPs per genome per year (95% HPD: 0.75–1.53). Using this, we estimated SNPs/year for each lineage (L1: 1.85; L2: 1.37; L3: 0.90; L4: 0.95; L5: 1.79; all within ~1 SNP difference per year). Given this consistency, applying a uniform threshold across ST313 lineages is epidemiologically reasonable. A threshold of ~3 SNPs corresponds to ~2–3 years of evolution, so we considered 3 SNPs a reasonable choice.\n(4) Incorporating isolation time: We agree that SNP-based source tracing should include epidemiological context. In the revised manuscript, when identifying transmission events within clusters, we no longer assume isolates within ≤3 SNP clusters are directly linked. We now incorporate sampling dates: we only count a potential transmission if isolates differ by ≤3 SNPs and were collected within a 3-year interval (approximately the time for 3 SNPs to accumulate in ST313). This excludes unlikely “transmission” pairs that meet the SNP cutoff but are >10 years apart. We developed a custom script for this refined analysis.\n(5) Simulation of source attribution vs. SNP threshold: Using the above refined approach (lineage-specific references and 3-year date criterion), we evaluated transmission inference across SNP thresholds from 1 to 20. We assumed that the optimal threshold is roughly the SNPs accruing over 2–3 years, and that if using 1 SNP is maximally strict, then increasing the threshold toward the optimum should not dramatically increase inferred events.\nOur results showed that as the threshold increased, the number of UPGMA clusters decreased, while total inferred transmission events increased (New Fig. 3a). For thresholds 1–5 SNPs, we compared the mean transmission events per cluster (Mann–Whitney U with FDR correction). The mean per-cluster transmission rate increased with the threshold. At 4 SNPs, the mean per-cluster transmissions were significantly higher than at 1 SNP (p = 0.002, FDR-corrected) (New Fig. S6a).\nNew Fig. 3a: The total number of transmission events was inferred under SNP thresholds ranging from 1 to 20. The x-axis represents SNP thresholds (1–20), and the y-axis shows the corresponding total number of transmission events identified at each threshold.\nNew Fig. S6a: Comparison of source-attribution performance when the SNP threshold is set to 1–5. We performed multiple hypothesis testing using the Mann–Whitney U test to compare the mean number of transmission events per cluster across different SNP thresholds. The resulting p-values were adjusted using the Benjamini–Hochberg false discovery rate (FDR) procedure. Differences were considered statistically significant at p < 0.05. (a) The mean number of transmission events per cluster was calculated under SNP thresholds of 1–5.\nWe also compared, across different SNP thresholds, the mean number of within-country and between-country transmission events per cluster. In both cases, the mean number of transmission events per cluster increased with increasing SNP threshold (New Fig. 3a).\nNew Fig. 3a: The mean numbers of within-country and between-country transmission events per cluster were inferred under SNP thresholds ranging from 1 to 20. The x-axis represents SNP thresholds (1–20), and the y-axis displays the corresponding mean number of within-country and between-country transmission events per cluster identified at each threshold. Red and purple bars indicate between-country and within-country transmission events, respectively.\nWhen we examined within-country transmission events across clusters, the mean number of events per cluster was significantly higher at a threshold of 4 SNPs than at 1 SNP (Mann–Whitney U test, FDR-corrected; p = 0.03) (New Fig. S6b). In contrast, for between-country transmission events, the mean number of events per cluster did not differ significantly between thresholds of 4 and 1 SNP (Mann–Whitney U test, FDR-corrected; p = 0.83).\nNew Fig. S6b: Comparison of source-attribution performance when the SNP threshold is set to 1–5. We performed multiple hypothesis testing using the Mann–Whitney U test to compare the mean number of transmission events per cluster across different SNP thresholds. The resulting p-values were adjusted using the Benjamini–Hochberg false discovery rate (FDR) procedure. Differences were considered statistically significant at p < 0.05. (b) The mean number of transmission events per cluster occurring within the same country was calculated under SNP thresholds of 1–5.\nTaken together, these collective results support using a 3-SNP threshold for our final transmission analyses of ST313. In light of these improvements, we have fully revised the Methods to clarify our assumptions about the optimal SNP threshold and to detail the statistical criteria applied. The methods now state:\nNew lines 603-623:\nSNP distance-based transmission tracing\nTo minimize potential bias introduced by the choice of reference during SNP-distance calculation, we selected one complete ST313 genome assembly from each of the five lineages as the corresponding lineage-specific reference. Specifically, strain A130 was used as the reference for lineage 1, D23580 for lineage 2, BKQZM9 for lineage 3, B3589 for lineage 4, and U2 for lineage 510,46. Using the alignment method described above, we obtained the cgSNP alignments, after which we calculated pairwise SNP distances between ST313 isolates using the SNP-dists tool.\nFurthermore, we performed clustering using the unweighted pair-group method with arithmetic means (UPGMA) under the optimal SNP threshold (SNP threshold = 3, as described in the Methods section below). This step was implemented in SciPy, ensuring that the SNP distances among isolates within each cluster did not exceed three SNPs. We also incorporated sampling time as additional epidemiological information by restricting the isolation dates of potentially related strains to within three years, a timeframe consistent with the estimated period required for the accumulation of three SNPs in ST313. Two ST313 isolates within the same cluster and collected within a three-year interval were considered to represent a single transmission event. To facilitate this analysis, we developed a custom script to compute the number of transmission events.\nNew lines 625-670:\nDetermination of the optimal SNP distance threshold\nWe first estimated the evolutionary rate of ST313 using the methodology described by Arthur W. Pightling et al.47. The number of new SNPs arising per year was determined by multiplying the evolutionary rate estimated with BEAST by the number of core SNP sites identified in the alignments. The results indicated that the overall evolutionary rate of ST313 is approximately 0.97 SNPs (95% HPD: 0.75–1.53) per genome per year. By incorporating evolutionary rates, we estimated the number of new SNPs that may arise per year for each lineage: L1: 1.85 (95% HPD: 1.62–2.09), L2: 1.37 (95% HPD: 1.26–1.56), L3: 0.9 (95% HPD: 0.6–1.19), L4: 0.95 (95% HPD: 0.69–1.22), L5: 1.79 (95% HPD: 1.46–2.09). As shown, the annual number of new SNPs in ST313 is highly consistent across lineages, with differences of less than one SNP per year.\nFurthermore, we evaluated the source-tracing performance under different SNP thresholds. Considering that an appropriate SNP distance threshold is generally assumed to correspond to the number of SNPs accumulated over approximately 2–3 years, we conducted a more detailed comparison of source-tracing performance using SNP thresholds of 1–5. The underlying assumptions were as follows: (i) the optimal threshold should approximate the number of SNPs expected to accumulate over 2–3 years in ST313; and (ii) when the threshold is set to 1 SNP, inference of transmission between isolates is maximally stringent, and the number of inferred transmission events should not increase substantially when moving to the optimal SNP threshold.\nWe found that as the SNP threshold increased, the number of clusters generated by UPGMA decreased, whereas the total number of inferred transmission events increased. We next compared, for thresholds of 1–5 SNPs, the mean number of transmission events per cluster (Fig. S6a). Notably, transmission events can only be identified in clusters containing two or more isolates; therefore, clusters comprising a single isolate were excluded. Our results show that, when the threshold was set to 4 SNPs, the mean number of transmission events per cluster was significantly higher than that obtained with a threshold of 1 SNP (Mann–Whitney U test, FDR-corrected; p = 0.002).\nWe also compared, across different SNP thresholds, the mean number of within-country and between-country transmission events per cluster. In both cases, the mean number of transmission events per cluster increased as the SNP threshold increased. When we examined within-country transmission events across clusters, the mean number of events per cluster was significantly higher at a threshold of 4 SNPs than at 1 SNP (Mann–Whitney U test, FDR-corrected; p = 0.03, Fig. S6b). In contrast, for between-country transmission events, the mean number of events per cluster did not differ significantly between thresholds of 4 and 1 SNP (Mann–Whitney U test, FDR-corrected; p = 0.83). Taken together, these observations support the use of a 3-SNP threshold for the final source attribution analyses of ST313.\nBased on the revised methodology described above, we rewrote the Results section of the manuscript. Because new reference genomes were used, the clusters were renumbered accordingly. Notably, the updated, more stringent source-attribution approach did not alter our principal conclusions regarding the intercontinental transmission of ST313 or its rare cross-host transmission events. The revised Results section now reads as follows:\nNew lines 202-229:\nHuman-borne ST313 transmission predominates in lineage 2\nTo identify recent transmission events of ST313, we conducted an SNP distance analysis using a 3-SNP threshold as the optimal cutoff (see the Methods section; Fig. 3a and Fig. S6). This identified 1,599 clusters with pairwise distances below the cutoff; 454 clusters (1,970/3,115 isolates) with ≥2 isolates (Fig. 3b). Four clusters (Clusters 357, 342, 131, and 49) showed intercontinental transmission events (Fig. 3c). The primary recent intercontinental transmission events were concentrated within cluster 357, predominantly associated with lineage 2. These transmission events (n = 8) occurred between the Americas (United States) and Africa (Congo), with the estimated timeframe spanning 2019 to 2020. These events involved six isolates in total: two from the United States and four from Congo. The geographic information for the Congo isolates was cross-validated against published literature, whereas the isolates from the United States have not yet been reported (Supplementary Table 3). Fewer transmission events occurred between Europe (UK) and Africa (Nigeria, n = 3) (Fig. 3d), involving four isolates with well-documented sampling backgrounds (Supplementary Table 3). Considering that these isolates belong to the same cluster and were collected within a relatively short time frame (less than three years), we speculate that the recent transmission events represented by these ST313 strains may be potentially associated with international travel.\nAlthough ST313 has been considered a lineage restricted to human-to-human transmission19, analysis of 454 clusters identified a potential cross-host transmission event within a single cluster (Cluster 453), involving human–canine transmission in the UK (Fig. S7). Despite this observation, cross-host transmission events remain exceedingly rare, accounting for only a single transmission event.\nFigure 3. The recent transmission events of Salmonella ST313 around the world. (a) The upper bar chart displays the average number of transmission events within each cluster using SNP thresholds ranging from 1 to 20. The lower bar chart illustrates global travel-associated within-country and between-country transmission events of ST313, inferred using average SNP thresholds between 1 and 20. The proportions of within-country and between-country transmission events for each threshold are represented by blue and red bars, respectively. (b) Using an average SNP threshold of 3, UPGMA clustering was performed on the SNP distances of 3,115 ST313 strains. A total of 454 clusters, each containing at least 2 ST313 strains, were selected from 1,599 clusters for further analysis. The bar chart on the right illustrates the number of strains contained within each of the 454 key clusters, where each cluster comprises more than two strains. (c) The phylogenetic tree on the left was constructed using cgSNPs from the 3,115 ST313 strains, with different lineages classified based on BAPs. The middle heatmap highlights the positions of the top 10 clusters from (b) in terms of ST313 numbers. On the far right, the geographical origins of the 454 clusters (each containing at least 2 ST313 strains) were analyzed, with the x-axis representing the proportion of strains from different geographical locations, and different colors indicating different continents. (d) Transmission events of the four clusters with cross-continental transmissions, incorporating sampling time as additional epidemiological evidence, are shown. Different colored lines represent each cluster. The earliest inferred travel-associated events are annotated along these lines. The bar chart on the right depicts the lineage composition of the four clusters with cross-continental transmission events. Transmission events inferred partly from public database metadata are indicated with dashed lines.\nResponse: Thank you for your valuable comment. We completely agree with your point. In the revised version of the manuscript, we have added the limitation you mentioned, along with appropriate references, in the “Discussion” section. The added text now reads as follows:\nNew lines 465-469: Finally, genomic inference of intercontinental ST313 transmission may reflect human travel or migration rather than direct bacterial spread29. Genome-based links alone cannot definitively distinguish person-to-person transmission from human-mediated dispersal across continents.\n9. References are missing or incomplete in several places. Please ensure that all methods (e.g., line 510), software tools and their versions (e.g., lines 522 and 590–597), and key claims (e.g., lines 536–539) are consistently supported by appropriate citations throughout the manuscript.\nResponse: Thank you for your valuable feedback. We have thoroughly checked the manuscript and addressed the areas where references were lacking or where software version numbers were not specified. Relevant citations have been added, and version details for the software used have been included.\n10. Figure visualization could be improved to better highlight key findings. In particular, the annotation of important lineages and clusters is insufficient. For example, lineage 2.4 and the clusters emphasized in the text (e.g., Clusters 850, 695, 888, 579, 648) are not clearly or consistently labeled in Fig. 3.\nResponse: We apologize for the incorrect labeling of clusters in Figure 3d. To address your comment #7, we have thoroughly revised Figure 3. Additionally, we have carefully reviewed all figure visualizations in the manuscript to ensure that similar errors do not occur.\nFigure 3. The recent transmission events of Salmonella ST313 around the world. (a) The upper bar chart displays the average number of transmission events within each cluster using SNP thresholds ranging from 1 to 20. The lower bar chart illustrates global travel-associated within-country and between-country transmission events of ST313, inferred using average SNP thresholds between 1 and 20. The proportions of within-country and between-country transmission events for each threshold are represented by blue and red bars, respectively. (b) Using an average SNP threshold of 3, UPGMA clustering was performed on the SNP distances of 3,115 ST313 strains. A total of 454 clusters, each containing at least 2 ST313 strains, were selected from 1,599 clusters for further analysis. The bar chart on the right illustrates the number of strains contained within each of the 454 key clusters, where each cluster comprises more than two strains. (c) The phylogenetic tree on the left was constructed using cgSNPs from the 3,115 ST313 strains, with different lineages classified based on BAPs. The middle heatmap highlights the positions of the top 10 clusters from (b) in terms of ST313 numbers. On the far right, the geographical origins of the 454 clusters (each containing at least 2 ST313 strains) were analyzed, with the x-axis representing the proportion of strains from different geographical locations, and different colors indicating different continents. (d) Transmission events of the four clusters with cross-continental transmissions, incorporating sampling time as additional epidemiological evidence, are shown. Different colored lines represent each cluster. The earliest inferred travel-associated events are annotated along these lines. The bar chart on the right depicts the lineage composition of the four clusters with cross-continental transmission events. Transmission events inferred partly from public database metadata are indicated with dashed lines.\nReviewer #3 (Remarks to the Author):\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.\nResponse: We would like to thank the reviewer and the Early Career Researcher who co-reviewed this manuscript. We sincerely appreciate their time, effort, and insights, which have helped improve our work. All the requested changes are highlighted in blue.\nReviewer #4 (Remarks to the Author):\nIn this manuscript the authors leverage a large data set including >11,000 isolates from non-fecal clinical samples, as well as an additional dataset of >3,115 ST313. The authors compared isolates from different geographic regions, identifying a novel lineage (Lineage 4) and sublineage (2.4) of ST313, which has important implications for understanding the evolution of ST313, an important cause of invasive disease. The authors also use a robust series of analyses to infer relative dates of emergence and effective population sizes for lineages of ST313. Closer examination of the gene content of these various lineages also suggested a potential difference in invasiveness index for L3 isolates – although pseudogene formation as described in this manuscript, is a bit more nuanced than presence of a ‘premature stop codon,’ and I have some suggestions that could strengthen the robustness of these results. The other key finding of the study is the transmission of ARGs on mobile elements, demonstrating that plasmids and transposons are the primary drivers of ARGs in these isolates.\nResponse: Thank you for your thoughtful and constructive review of our manuscript. We greatly appreciate your positive assessment of the significance of our study, which is highly encouraging. We acknowledge that pseudogene formation involves complexities beyond “stop gained” or “frameshift variant” mutations, and that, as you noted, the specific locations of these mutations should also be considered. In the revised manuscript, we have made extensive revisions in accordance with your suggestions, addressing all of the concerns you raised—including those related to pseudogene detection—to strengthen the robustness of our findings.\nWhile the study provides novel insights into the population structure of ST313 with an appreciable large dataset, I have concerns about the accuracy of the metadata, and the potential implications of errors in the metadata that would greatly change the conclusions drawn. I think overall, the authors should be careful about the interpretation of the results considering the potential inaccuracies in metadata, which unfortunately cannot be confirmed given their large dataset. I have some recommendation about interpretations of the data and some additional analyses that may help improve the robustness of the authors' findings. There are some cases where additional details would be helpful for understanding what was done and ensuring reproducibility of results.\nResponse: Thank you for your constructive feedback. In response, we conducted a thorough validation of the metadata, which confirmed the robustness of our core findings while allowing us to refine certain interpretations with greater caution. Additionally, we implemented the analyses you recommended to further strengthen the reliability of our results. To ensure reproducibility, we have added comprehensive methodological details in the revised manuscript. We believe these revisions address your concerns and enhance the overall quality of the work. Please find our point-by-point responses below. All corresponding changes in the manuscript are highlighted in yellow.\nMajor Recommendations:\nMy biggest concern lies with the accuracy of the metadata, and the sampling biases that have the potential to greatly influence many of the results reported. The implications of this are not trivial – the authors are making claims about directionality/introduction/spread of a lineage to/from a geographic region. In some cases, an isolate may be collected in one country and then sequenced and the genomic assembly uploaded from another country; it is possible that the ‘source’ country listed in the metadata would be the latter (country where the researcher uploaded data, rather than where the isolate was actually found).\nResponse: We appreciate you raising this critical issue. In the revised manuscript, we have taken the following steps to ensure the accuracy and robustness of the metadata:\n1. Comprehensive cross-check of metadata for all 3,115 ST313 isolates\nTo conduct a comprehensive cross-check of metadata for all 3,115 ST313 isolates, we surveyed all published studies on ST313 and cross-verified the metadata for every ST313 isolate in our dataset. Specifically, we updated Supplementary Table 2 to include the corresponding source references and the collectors of each isolate for which published information is available, confirmed by the literature (n = 1,584, 50.8%).\n2. Verification of isolates lacking published reports\nTo verify isolates lacking published reports, we systematically re-evaluated all ST313 isolates in our collection that lacked any published documentation (hereafter referred to as “unpublished” isolates, n = 1,531). We first sought additional guidance from Dr. Zhemin Zhou, one of the developers of EnteroBase, to refine our approach. Our adopted strategy was as follows: when BioSample information was available in EnteroBase, we retrieved the corresponding metadata directly from the BioSample record, unless an isolate was deposited exclusively in EnteroBase without a linked BioSample entry. Given that the metadata for most “unpublished” isolates were primarily obtained from BioSample records, we further evaluated the accuracy of metadata derived from BioSample records. Specifically, we compared the metadata of 1,584 isolates with clear, literature-supported background information to the metadata retrieved solely from their corresponding BioSample links. We found a 100% match for all key metadata fields relevant to this study, including sampling dates, locations, and host sources.\nWe paid particular attention to whether the country listed in the “geographic location” field represented the actual country of isolation rather than the country of data submission. We believe that the metadata retrieval strategy described above adequately addresses your concern. As an illustrative example, BioSample SAMN03098808 was submitted by the U.S. CDC PulseNet program; however, the BioSample record clearly indicates that the isolate was collected in Kenya. As shown in the figure, this demonstrates that the BioSample entry reliably records the true geographic origin of the isolate. In the revised manuscript, we manually examined the geographic information for all 1,531 “unpublished” isolates and corrected the erroneous geographic assignment for sublineage 2.4 reported in the original version.\n3. Refined geographic conclusions and verified provenance of key “unpublished” isolates with collaborators\nWe acknowledge that some key isolates associated with geographic transmission remain “unpublished.” For these key unpublished isolates, we undertook additional verification steps. Using the contact information provided in EnteroBase, we directly contacted the corresponding data holders to further confirm and clarify their metadata; for example, we consulted François-Xavier Weill regarding several of these isolates.\nThrough these procedures, we made every effort to ensure the accuracy of the metadata for these strains. To further guarantee the robustness of our conclusions, we extensively revised the geographic transmission analysis in light of the reviewer’s comments. Specifically, we softened the tone of the statements, more clearly indicated in the figures—using dashed lines—the uncertain transmission branches associated with these key unpublished isolates, and rewrote the corresponding conclusions accordingly. The detailed revisions are described in the responses below.\nUnfortunately, with this dataset, there is no way for the authors to know how accurate the metadata are without consulting the depositors (and it is unreasonable to reach out to public health agencies given that they have >3,000 isolates).\nResponse: The reviewer is correct that it is not feasible to individually consult public health agencies to verify the metadata of all 3,115 isolates. However, we have used a systematic approach, as outlined above, to assess the accuracy of the available metadata. Through extensive literature searches, we confirmed that 1,584 of the 3,115 ST313 isolates have been explicitly reported in published studies. Metadata for these isolates were obtained from the supplementary tables of the corresponding publications, and we identified the original collectors for each isolate. Cross-verification showed that these 1,584 isolates were documented across 20 independent studies (New Fig. S19). Among these, 586 isolates were from Caisey V. Pulford et al. (Nat. Microbiol 2021; DOI: 10.1038/s41564-020-00836-1) and 429 isolates were from Van Puyvelde et al. (Nat. Commun 2023; DOI: 10.1038/s41467-023-41152-6); we have added this source information to Supplementary Table 2. Additionally, our consistency check found no discrepancies between BioSample and published literature with the strain metadata for key fields (sampling dates, locations, and host sources).\nNew Fig. S19. Literature sources of the 1,584 ST313 isolates with explicit published documentation and their corresponding collectors.\nThe assertions made in the discussion section about ‘an imported, highly AMR branch now circulating in Africa, primarily via US-Kenya transmissions between 2009 and 2016’ would change if the metadata are not accurate, and I worry about the potential harm that could result from such assertions. I think the authors need to carefully consider ways that they can validate the accuracy of the metadata, particularly for any claims made about directionality.\nResponse: We sincerely thank the reviewer for raising this important comment. The reviewer is correct that several isolates previously assigned to sublineage 2.4 as being from the United States were, in fact, sampled in Kenya but sequenced by the U.S. Centers for Disease Control and Prevention (CDC). For example, BioSample SAMN03098808 indicates the isolate was collected in Kenya, whereas the submitter is PulseNet USA. This discrepancy explains our earlier metadata misclassification. Following a thorough re-evaluation, we have corrected the isolates previously misassigned as U.S.-origin in sublineage 2.4 of the clade, and the manuscript has been updated accordingly.\nTo further strengthen our transmission-related conclusions, we also incorporated the following revisions:\n1. Lineage-specific reference genomes for SNP analysis. To address the potential impact of reference choice, we now use a complete ST313 genome from each lineage as that lineage’s reference. Specifically, we selected strain A130 for Lineage 1, D23580 for Lineage 2, BKQZM9 for Lineage 3, B3589 for Lineage 4, and U2 for Lineage 5, and re-conducted the phylogenetic and SNP distance analyses with these lineage-specific references.\n2. Additional isolation time criterion for transmission events. In the revised analysis, we identify putative transmission clusters using a 3-SNP cutoff and require that isolates were collected within a three-year window. This time frame is consistent with the estimated period for ~3 SNPs to accumulate in ST313. Using this approach allows us to exclude unlikely transmission events in which isolates exhibit fewer than three SNP differences, but were collected even more than a decade apart, where genetically close isolates were collected many years apart. We developed a custom script to perform this refined analysis. (The Methods section has been updated to describe this approach.)\nNew lines 603-623:\nSNP distance-based transmission tracing\nTo minimize potential bias introduced by the choice of reference during SNP-distance calculation, we selected one complete ST313 genome assembly from each of the five lineages as the corresponding lineage-specific reference. Specifically, strain A130 was used as the reference for lineage 1, D23580 for lineage 2, BKQZM9 for lineage 3, B3589 for lineage 4, and U2 for lineage 510,46. Using the alignment method described above, we obtained the cgSNP alignments, after which we calculated pairwise SNP distances between ST313 isolates using the SNP-dists tool.\nFurthermore, we performed clustering using the unweighted pair-group method with arithmetic means (UPGMA) under the optimal SNP threshold (SNP threshold = 3, as described in the Methods section below). This step was implemented in SciPy, ensuring that the SNP distances among isolates within each cluster did not exceed three SNPs. We also incorporated sampling time as additional epidemiological information by restricting the isolation dates of potentially related strains to within three years, a timeframe consistent with the estimated period required for the accumulation of three SNPs in ST313. Two ST313 isolates within the same cluster and collected within a three-year interval were considered to represent a single transmission event. To facilitate this analysis, we developed a custom script to compute the number of transmission events.\n3. More cautious interpretation of geographic spread. We revised the Results section to temper our conclusions regarding transmission “directionality.” At the same time, we updated Supplementary Table 3 to provide the metadata sources for all isolates inferred to be geographically linked based on the SNP distance approach (for example, whether an isolate has a clearly documented literature source or whether its metadata were confirmed directly with the data owner). Finally, for transmission inferences supported only by metadata from public databases, we modified Figure 3 so that these routes are depicted with dashed lines (primarily clusters 357 and 131). The revised Results section and corresponding figure are shown below.\nNew lines 202-229:\nHuman-borne ST313 transmission predominates in lineage 2\nTo identify recent transmission events of ST313, we conducted an SNP distance analysis using a 3-SNP threshold as the optimal cutoff (see the Methods section; Fig. 3a and Fig. S6). This identified 1,599 clusters with pairwise distances below the cutoff; 454 clusters (1,970/3,115 isolates) with ≥2 isolates (Fig. 3b). Four clusters (Clusters 357, 342, 131, and 49) showed intercontinental transmission events (Fig. 3c). The primary recent intercontinental transmission events were concentrated within cluster 357, predominantly associated with lineage 2. These transmission events (n = 8) occurred between the Americas (United States) and Africa (Congo), with the estimated timeframe spanning 2019 to 2020. These events involved six isolates in total: two from the United States and four from Congo. The geographic information for the Congo isolates was cross-validated against published literature, whereas the isolates from the United States have not yet been reported (Supplementary Table 3). Fewer transmission events occurred between Europe (UK) and Africa (Nigeria, n = 3) (Fig. 3d), involving four isolates with well-documented sampling backgrounds (Supplementary Table 3). Considering that these isolates belong to the same cluster and were collected within a relatively short time frame (less than three years), we speculate that the recent transmission events represented by these ST313 strains may be potentially associated with international travel.\nAlthough ST313 has been considered a lineage restricted to human-to-human transmission16, analysis of 454 clusters identified a potential cross-host transmission event within a single cluster (Cluster 353), involving human–canine transmission in the UK (Fig. S5). Despite this observation, cross-host transmission events remain exceedingly rare, accounting for only a single transmission event.\nFigure 3. Recent Salmonella ST313 transmission events worldwide. (a) The upper bar chart displays the average number of transmission events within each cluster using SNP thresholds ranging from 1 to 20. The lower bar chart illustrates global travel-associated within-country and between-country transmission events of ST313, inferred using average SNP thresholds between 1 and 20. The proportions of within-country and between-country transmission events for each threshold are represented by blue and red bars, respectively. (b) Using an average SNP threshold of 3, UPGMA clustering was performed on the SNP distances of 3,115 ST313 strains. A total of 454 clusters, each containing at least 2 ST313 strains, were selected from 1,599 clusters for further analysis. The bar chart on the right illustrates the number of strains contained within each of the 454 key clusters, where each cluster comprises more than two strains. (c) The phylogenetic tree on the left was constructed using cgSNPs from the 3,115 ST313 strains, with different lineages classified based on BAPs. The middle heatmap highlights the positions of the top 10 clusters from (b) in terms of ST313 numbers. On the far right, the geographical origins of the 454 clusters (each containing at least 2 ST313 strains) were analyzed, with the x-axis representing the proportion of strains from different geographical locations, and different colors indicating different continents. (d) Transmission events of the four clusters with cross-continental transmissions, incorporating sampling time as additional epidemiological evidence, are shown. Different colored lines represent each cluster. The earliest inferred travel-associated events are annotated along these lines. The bar chart on the right depicts the lineage composition of the four clusters with cross-continental transmission events. Transmission events inferred partly from public database metadata are indicated with dashed lines.\nFigure S5. Global host-associated transmission patterns of ST313. (a) The phylogenetic tree on the left was constructed using core-genome Single Nucleotide Polymorphisms (cgSNPs) from the 3,115 ST313 strains, with different lineages classified based on BAPs. The heatmap on the right indicates the positions of clusters 437 and 453, associated with cross-host transmission events. (b) A statistical analysis was performed on the host origins of Salmonella within 454 clusters, each containing at least two ST313 strains. The x-axis represents the proportion of ST313 strains originating from different host sources. Blue bars denote animal-derived strains, yellow bars indicate food-derived strains, and red bars represent human-derived strains. (c) A transmission event between a human and a canine was identified within cluster 453. Transmission events inferred partly from public database metadata are indicated with dashed lines.\nThe other important thing to consider is that there is undoubtedly bias in relying on WGS data rather than reported case data. It is unfair to assume that cases in each country have the same likelihood of having isolates characterized with WGS.\nResponse: You are correct that relying on sequenced genomes can introduce sampling bias. First, we acknowledge this limitation. Second, we attempted to quantify this sampling bias within our global dataset of 11,361 iNTS genomes to enable a more informed discussion. Specifically, we examined the relationship between each country’s economic level (measured as GDP per capita) and its corresponding sample count (New Fig. S17).\nNew Fig. S17: Source-bias analysis of available iNTS WGS data reported by each country. Each circle represents a country (circle size indicates the number of iNTS isolates). The x-axis is national GDP per capita, and the y-axis is the number of iNTS isolates from that country.\nThe analysis confirms this bias: countries with greater public health infrastructure and resources (e.g., the USA, Canada, China) appear to have higher recorded numbers of sequenced isolates. At the same time, regions with fewer resources—such as Malawi, Congo, and South Africa—also show high isolate counts. This pattern likely reflects large, well-funded studies conducted specifically to investigate iNTS in select African countries.\nI appreciate that the authors are trying to use the genomic data to uncover helpful information about the spread, but I find the data presentation and conclusions drawn potentially misleading. For example, the temporal analyses in panel 1b and the supplemental data, require equal sampling across all time points. We know this is not true – there were several large, well-funded studies that specifically sought to understand iNTS burden in select countries in Africa. I believe that what is shown in the data in Figure 1B reflects a reduction in focused efforts, rather than a true signal of reductions in ST313 infections (which is what is being assumed by the authors). To address this, the authors should carefully consider reframing their conclusions about decline in ST313 genome sequences over time as this likely reflects projects ending rather than a true decline in infections.\nResponse: We fully agree with your point, and we apologize for the lack of precision in our earlier description. In the revised manuscript, we have thoroughly rewritten the section previously titled “Global landscape of iNTS reveals an epidemiological transition,” with a clearer emphasis on the shift in global sampling efforts. We now attribute the post-2014 decline in ST313 genomes to a reduction in sequencing efforts (as major projects wound down) rather than a true epidemiological decline. For example, we highlight that Africa was the predominant source of iNTS isolates and that the apparent drop in ST313 after 2014 coincides with decreased sequencing in Africa, not necessarily fewer infections. The revised paragraph is as follows：\nNew lines 120-151:\nThe global landscape of the iNTS collection\nAnalysis of 11,361 publicly available iNTS genomes (Supplementary Table 1) identifies Africa as the predominant region of iNTS sampling (n = 4,741; 41.7%), followed by North America (n = 4,351; 38.3%). Although Asia accounts for fewer samples (n = 1,350; 11.9%), these were primarily from China (1,173 of 1,350; 86.9%) (Fig. S1a). Calculations of the sampling rates revealed an overall upward trend across countries, with the largest increase observed in South Africa (Fig. S1b). Serovar classification underscored the predominance of S. Typhimurium (n = 3,260, 28.7%) and S. Enteritidis (n = 3,785, 33.3%). Interestingly, in our dataset, S. Typhimurium samples predominated before 2015, but this shifted to S. Enteritidis after 2015 (Fig. S1c). Next, to examine the details of this shift within S. Typhimurium, we used multilocus sequence typing (MLST) to identify three major circulating sequence types (STs): ST313 (n = 1,763), ST19 (n = 1,140), and ST34 (n = 237). Among these, the temporal changes in the available ST313 genomes from public databases closely mirrored the post-2015 decline in S. Typhimurium sampling (Fig. S1d–e).\nFig S1. The global landscape of the iNTS collection (n=11,361). (a) Geographic distribution of iNTS within the dataset, with red circles indicating the number of iNTS genomes in each country; larger circles represent higher genome counts. (b) The annual iNTS sampling rates across different countries were calculated by performing linear fits on the yearly number of iNTS isolates reported by each country. Red areas indicate an increasing trend, while blue areas represent a decreasing trend. (c) Temporal distribution of the two predominant iNTS serovars, S. Typhimurium and S. Enteritidis, represented by red and blue, respectively. (d) Temporal distribution of the three major S. Typhimurium STs: ST313 (red), ST19 (blue), and ST34 (green). (e) Linear regression was used to analyze the post-2015 growth trends of the major S. Typhimurium STs, confirming that ST313 in the dataset was the primary driver of the overall decline in S. Typhimurium.\nNext, by assembling all publicly available ST313 genomic data, we constructed the most extensive global dataset of ST313 genomes to date, comprising 3,115 isolates (Supplementary Table 2). Our results show that, geographically, ST313 predominated among isolates sampled in Africa (n = 2,266; 72.7%), with a growing but still limited sampling presence in Europe (n = 260; 8.3%) (Fig. 1a–b). From a host perspective, most isolates were human-derived (n = 2,027, 65.1%), with rare detections in animals (n = 15, 0.5%) and food (n = 8, 0.3%) (Fig. 1c). Notably, among human isolates, blood samples were the primary source (1,702/2,027, 84%), far exceeding faecal isolates (96/2,027, 4.7%) (Fig. 1d). Only a small number of ST313 isolates contained sex information (n = 440), of which more than 99% were reported to have been collected at Queen Elizabeth Central Hospital (Blantyre, Malawi)10. Overall, the sex distribution was balanced, with slightly more male cases (n = 226) than female (n = 214) (Fig. 1e).\nFigure 1. Geographic, source, and temporal distribution of the ST313 genomic dataset (n = 3,115). (a) Geographic distribution of ST313 in the constructed dataset. Blue circles denote the number of ST313 genomes in different countries, with larger circles representing a larger number of genomes. (b) The number of ST313 isolations in the constructed dataset across six continents over multiple years. The y-axis represents the number of genomes, and different colors indicate different continents. (c) Statistics on the location and sampling time of ST313 genomes isolated from \"Human,\" \"Food,\" and \"Animal\" sources. (d) Statistics of the source of human-isolated ST313 in the dataset. (e) Statistics of the sex of patients from whom ST313 was isolated. The y-axis represents the number of ST313 genomes from patients of each sex.\nIn addition, based on New Fig. S17, we further discuss how several large, well-funded studies conducted in multiple African countries have contributed to sampling biases in the number of available isolates:\nNew lines 451-464:\nThis study has limitations. First, disparities in surveillance and reporting are likely to bias prevalence estimates toward specific regions. Because the available data are derived from WGS rather than reported case counts, this issue becomes even more apparent when assembling a global iNTS genomic dataset comprising 11,361 isolates. To quantify this sampling bias, we examined the relationship between a country’s economic level (measured as GDP per capita) and its corresponding sample count (Fig. S17). The results indicate that countries with greater public health infrastructure and resources—particularly the USA, Canada, and China—appear to have higher recorded incidences than those with more limited resources. At the same time, regions with fewer resources—such as Malawi, Congo, and South Africa—also seem to show higher incidences. This pattern reflects existing large, well-funded studies conducted specifically to investigate the burden of iNTS in select African countries.\nMinor Recommendations:\nLines 116-131: I think the way that iNTS was defined from the metadata could be problematic and is a key limitation that needs to be discussed; were isolates that had an ‘unknown’ isolation source included? Why did the authors not limit their analyses to just isolates from blood?\nResponse: We apologize for the insufficient description in the original manuscript. First, regarding the 11,361 global invasive non-typhoidal Salmonella (iNTS) isolates in Supplementary Table 1: we did not include any isolate with an “unknown” source. We have expanded Supplementary Table 1 to list the isolation source for each strain. Among the 11,361 iNTS isolates, 23 distinct isolation sites were identified. The majority were from blood (n = 7,922; 69.7%) and urine (n = 2,971; 26.2%), with additional isolates from abscesses (n = 116), wounds (n = 81), cerebrospinal fluid (n = 74), the respiratory tract (n = 36), synovial fluid (n = 30), etc. (New Fig. S18).\nNew Fig. S18: Sample isolation sites for the global iNTS genome dataset comprising 11,361 isolates. (a) All strains in the dataset were isolated from human. (b) Specific isolation sites for the strains, with the majority recovered from blood (n = 7,922; 69.7%) and urine (n = 2,971; 26.2%). (c) A breakdown of isolates from other sources, as indicated in panel b.\nSecondly, the definition of iNTS used in our study follows the established literature (Lancet Infect Dis. 2019; DOI: 10.1016/S1473-3099(19)30418-9). In brief, invasive NTS refers to Salmonella (excluding S. Typhi and S. Paratyphi) isolated from normally sterile sites in the human body (e.g., blood, cerebrospinal fluid, etc.). Based on this definition, our dataset includes not only blood isolates but also those from other sterile sites (urine, abscesses, CSF, etc.). We have now explicitly clarified this in the Methods section.\nNew lines 478-500:\nAssembly of the global iNTS genomic dataset\nIn this genomic epidemiology study, iNTS was defined as non-typhoidal and non-paratyphoidal Salmonella isolates obtained from human specimens of non-fecal origin, including blood, urine, abscesses, and other sterile sites, as described in previous literature2. We subsequently retrieved all available Salmonella genomes from EnteroBase and GenBank (data cutoff: 1 March 2024; total number: 488,232). For metadata, when BioSample information was available in EnteroBase, we obtained the corresponding metadata directly from the associated BioSample record, unless the isolate was deposited exclusively in EnteroBase without a linked BioSample entry. We excluded all strains with ambiguous isolation sources—for example, those labeled only as “clinical” or “human”—and retained only those isolates with clearly defined sources that met the above criteria for iNTS. Information on the country and date of isolation was also collected. This process resulted in a final dataset of 11,361 iNTS genomes (Supplementary Table 1).\nFor the global iNTS genomic dataset constructed above, we further confirmed the species identity as Salmonella using KmerFinder v.3.0.230, and serovar determination was performed with SISTR v.1.1.231 and SeqSero2 v.1.3.132. In terms of isolation sources, most isolates were recovered from blood (n = 7,922; 69.7%) and urine (n = 2,971; 26.2%). Additional isolates primarily originated from abscesses (n = 116), wounds (n = 81), cerebrospinal fluid (n = 74), the respiratory tract (n = 36), and synovial fluid (n = 30) (Fig. S18).\nWe also acknowledge that identifying iNTS solely based on metadata may introduce bias. The revised sentence is as follows:\nNew lines 375-381:\nIn this study, we assembled the largest dataset of iNTS genomes to date (n = 11,361), using isolation-source information and established epidemiological definitions of iNTS to delineate global patterns. However, it is important to note that our global iNTS genomic collection included only isolates with clearly defined sources; thus, some potential iNTS strains that were not explicitly reported may have been overlooked. This likely reflects variations in clinical reporting and sampling practices across different regions.\nMost clinical samples do not list the isolation source. How might variations in clinical diagnosis and reporting bias the results in the genomic dataset?\nResponse: You are correct that this is a limitation of our study. By including only isolates with clearly defined sources, we may have missed some invasive Salmonella cases that were not explicitly documented as such, owing to regional differences in clinical reporting and sampling. We have added a statement to the Discussion acknowledging this potential bias.\nNew lines 375-381:\nIn this study, we assembled the largest dataset of iNTS genomes to date (n = 11,361), using isolation-source information and established epidemiological definitions of iNTS to delineate global patterns. However, it is important to note that our global iNTS genomic collection included only isolates with clearly defined sources; thus, some potential iNTS strains that were not explicitly reported may have been overlooked. This likely reflects variations in clinical reporting and sampling practices across different regions.\nIn the materials and methods ‘iNTS’ is defined as ‘not from feces’ – in this case most of the dataset should have included Typhi and Paratyphi strains as well but none are reported?\nResponse: This appears to be a matter of conceptual clarification. As detailed above, the primary distinction between invasive non-typhoidal Salmonella (iNTS) and non-typhoidal Salmonella (NTS) is iNTS's capacity to cause bloodstream and other invasive infections from non-fecal sterile sites. Genomic epidemiology studies of Salmonella routinely provide explicit definitions of NTS, for example:\n1) Lancet. 2012; DOI: 10.1016/S0140-6736(11)61752-2:\n“With regard to human disease, salmonellae are divided into typhoidal serotypes (Salmonella enterica var Typhi [S. Typhi] and Salmonella enterica var Paratyphi [S. Paratyphi]) and thousands of non-typhoidal Salmonella serotypes (frequently referred to as NTS serotypes).”\n2) Pediatrics. 2023; DOI: 10.1542/peds.2023-062357:\n“NTS was defined as Salmonella spp. that were not of serotype typhi or paratyphi.”\nTherefore, no S. Typhi or S. Paratyphi strains are included in the global iNTS genomic dataset, as these were intentionally excluded. We believe that the additional clarification of the iNTS definition provided in the revised manuscript will effectively address this concern.\nNew lines 479-482:\nIn this genomic epidemiology study, iNTS was defined as non-typhoidal and non-paratyphoidal Salmonella isolates obtained from human specimens of non-fecal origin, including blood, urine, abscesses, and other sterile sites, as described in previous literature2.\nWas additional filtering done aside from just ‘excluding feces as a source?’ I think the details about the way the iNTS dataset are missing to allow complete understanding of how these data were collected.\nResponse: We apologize for the insufficient description in the original manuscript. We have now provided a detailed description of how the iNTS dataset was assembled. In summary:\n(1) Data collection: We retrieved all available Salmonella genome data from EnteroBase and GenBank (data cutoff 1 March 2024; ~488,000 total entries). Because metadata can vary between databases, we consulted Dr. Zhemin Zhou (a previous EnteroBase co-creator) to devise a consistent strategy. When BioSample information was available, we obtained metadata from the BioSample record. We included only isolates explicitly annotated as human-derived and from sterile sites (blood, urine, abscess, etc.), excluding any isolates from non-sterile sources (e.g., fecal samples).\n(2) In silico validation: We confirmed each genome’s species as Salmonella using KmerFinder and determined each isolate’s serovar using SISTR and SeqSero2. During this process, any S. Typhi or S. Paratyphi entries were removed, yielding a final global iNTS dataset of 11,361 genomes (Supplementary Table 1). We also assigned sequence types (STs) using MLST (Achtman scheme) for context, but we did not exclude any isolates based on ST.\nWe have thoroughly revised the Methods section to include all of these details, ensuring transparency in how the data were collected and filtered.\nNew lines 478-500:\nAssembly of the global iNTS genomic dataset\nIn this genomic epidemiology study, iNTS was defined as non-typhoidal and non-paratyphoidal Salmonella isolates obtained from human specimens of non-fecal origin, including blood, urine, abscesses, and other sterile sites, as described in previous literature2. We subsequently retrieved all available Salmonella genomes from EnteroBase and GenBank (data cutoff: 1 March 2024; total number: 488,232). For metadata, when BioSample information was available in EnteroBase, we obtained the corresponding metadata directly from the associated BioSample record, unless the isolate was deposited exclusively in EnteroBase without a linked BioSample entry. We excluded all strains with ambiguous isolation sources—for example, those labeled only as “clinical” or “human”—and retained only those isolates with clearly defined sources that met the above criteria for iNTS. Information on the country and date of isolation was also collected. This process resulted in a final dataset of 11,361 iNTS genomes (Supplementary Table 1).\nFor the global iNTS genomic dataset constructed above, we further confirmed the species identity as Salmonella using KmerFinder v.3.0.230, and serovar determination was performed with SISTR v.1.1.231 and SeqSero2 v.1.3.132. In terms of isolation sources, most isolates were recovered from blood (n = 7,922; 69.7%) and urine (n = 2,971; 26.2%). Additional isolates primarily originated from abscesses (n = 116), wounds (n = 81), cerebrospinal fluid (n = 74), the respiratory tract (n = 36), and synovial fluid (n = 30) (Fig. S18).\nI see the potential for the reporting of these cases to be biased by public health reporting differences by country. For example, countries with greater public health infrastructure and resources will appear to have higher incidences than those with lower resources. This is an important source of bias that is partially addressed in the discussion\nResponse: We fully agree. We analyzed the relationship between each country’s economic level (GDP per capita) and its number of sequenced isolates (New Fig. S17, described earlier). The results confirmed that wealthier countries often contribute more WGS data, whereas high contributions from some lower-resource countries were driven by targeted studies.\nNew Fig. S17: Source-bias analysis of available iNTS WGS data reported by each country. Each circle represents a country (circle size indicates the number of iNTS isolates). The x-axis is national GDP per capita, and the y-axis is the number of iNTS isolates from that country.\nWe have expanded our discussion to explicitly acknowledge this sampling bias and its impact on our findings, as follows:\nNew lines 451-464:\nThis study has limitations. First, disparities in surveillance and reporting are likely to bias prevalence estimates toward specific regions. Because the available data are derived from WGS rather than reported case counts, this issue becomes even more apparent when assembling a global iNTS genomic dataset comprising 11,361 isolates. To quantify this sampling bias, we examined the relationship between a country’s economic level (measured as GDP per capita) and its corresponding sample count (Fig. S17). The results indicate that countries with greater public health infrastructure and resources—particularly the USA, Canada, and China—appear to have higher recorded incidences than those with more limited resources. At the same time, regions with fewer resources—such as Malawi, Congo, and South Africa—also seem to show higher incidences. This pattern reflects existing large, well-funded studies conducted specifically to investigate the burden of iNTS in select African countries.\nLines 129-131: In figure 1B the authors show a marked decline in sequences from Africa after 2014 – given that the authors’ assertion that ST313 is overrepresented among clinical cases in Africa, doesn’t this instead suggest a decrease in genomic sequencing efforts rather than a true decline in ST313 infections? Logically, it does not make sense that ST313 would magically decline; rather, this again suggests sampling bias. This is a limitation in the dataset that needs to be discussed. This is important because these data are from WGS (not just reported cases).\nResponse: The newly added analysis (New Fig. S17) further confirms that this decline is due to reduced sequencing efforts. In the revised manuscript, we have thoroughly rewritten this section to emphasize the reduction in sequencing efforts for ST313 after 2014, rather than a true drop in infections. We also explicitly discuss this limitation in the Discussion, noting that the apparent decline in ST313 genomes corresponds to the conclusion of specific sequencing projects (and we cite this as a bias in our study rather than a genuine epidemiological trend).\nNew lines 120-151:\nThe global landscape of the iNTS collection\nAnalysis of 11,361 publicly available iNTS genomes (Supplementary Table 1) identifies Africa as the predominant region of iNTS sampling (n = 4,741; 41.7%), followed by North America (n = 4,351; 38.3%). Although Asia accounts for fewer samples (n = 1,350; 11.9%), these were primarily from China (1,173 of 1,350; 86.9%) (Fig. S1a). Calculations of the sampling rates revealed an overall upward trend across countries, with the largest increase observed in South Africa (Fig. S1b). Serovar classification underscored the predominance of S. Typhimurium (n = 3,260, 28.7%) and S. Enteritidis (n = 3,785, 33.3%). Interestingly, in our dataset, S. Typhimurium samples predominated before 2015, but this shifted to S. Enteritidis after 2015 (Fig. S1c). Next, to examine the details of this shift within S. Typhimurium, we used multilocus sequence typing (MLST) to identify three major circulating sequence types (STs): ST313 (n = 1,763), ST19 (n = 1,140), and ST34 (n = 237). Among these, the temporal changes in the available ST313 genomes from public databases closely mirrored the post-2015 decline in S. Typhimurium sampling (Fig. S1d–e).\nFig S1. The global landscape of the iNTS collection (n=11,361). (a) Geographic distribution of iNTS within the dataset, with red circles indicating the number of iNTS genomes in each country; larger circles represent higher genome counts. (b) The annual iNTS sampling rates across different countries were calculated by performing linear fits on the yearly number of iNTS isolates reported by each country. Red areas indicate an increasing trend, while blue areas represent a decreasing trend. (c) Temporal distribution of the two predominant iNTS serovars, S. Typhimurium and S. Enteritidis, represented by red and blue, respectively. (d) Temporal distribution of the three major S. Typhimurium STs: ST313 (red), ST19 (blue), and ST34 (green). (e) Linear regression was used to analyze the post-2015 growth trends of the major S. Typhimurium STs, confirming that ST313 in the dataset was the primary driver of the overall decline in S. Typhimurium.\nNext, by assembling all publicly available ST313 genomic data, we constructed the most extensive global dataset of ST313 genomes to date, comprising 3,115 isolates (Supplementary Table 2). Our results show that, geographically, ST313 predominated among isolates sampled in Africa (n = 2,266; 72.7%), with a growing but still limited sampling presence in Europe (n = 260; 8.3%) (Fig. 1a–b). From a host perspective, most isolates were human-derived (n = 2,027, 65.1%), with rare detections in animals (n = 15, 0.5%) and food (n = 8, 0.3%) (Fig. 1c). Notably, among human isolates, blood samples were the primary source (1,702/2,027, 84%), far exceeding faecal isolates (96/2,027, 4.7%) (Fig. 1d). Only a small number of ST313 isolates contained sex information (n = 440), of which more than 99% were reported to have been collected at Queen Elizabeth Central Hospital (Blantyre, Malawi)10. Overall, the sex distribution was balanced, with slightly more male cases (n = 226) than female (n = 214) (Fig. 1e).\nFigure 1. Geographic, source, and temporal distribution of the ST313 genomic dataset (n = 3,115). (a) Geographic distribution of ST313 in the constructed dataset. Blue circles denote the number of ST313 genomes in different countries, with larger circles representing a larger number of genomes. (b) The number of ST313 isolations in the constructed dataset across six continents over multiple years. The y-axis represents the number of genomes, and different colors indicate different continents. (c) Statistics on the location and sampling time of ST313 genomes isolated from \"Human,\" \"Food,\" and \"Animal\" sources. (d) Statistics of the source of human-isolated ST313 in the dataset. (e) Statistics of the sex of patients from whom ST313 was isolated. The y-axis represents the number of ST313 genomes from patients of each sex.\nWe have also further discussed this limitation in the revised manuscript.\nNew lines 451-464:\nThis study has limitations. First, disparities in surveillance and reporting are likely to bias prevalence estimates toward specific regions. Because the available data are derived from WGS rather than reported case counts, this issue becomes even more apparent when assembling a global iNTS genomic dataset comprising 11,361 isolates. To quantify this sampling bias, we examined the relationship between a country’s economic level (measured as GDP per capita) and its corresponding sample count (Fig. S17). The results indicate that countries with greater public health infrastructure and resources—particularly the USA, Canada, and China—appear to have higher recorded incidences than those with more limited resources. At the same time, regions with fewer resources—such as Malawi, Congo, and South Africa—also seem to show higher incidences. This pattern reflects existing large, well-funded studies conducted specifically to investigate the burden of iNTS in select African countries.\nLines 123-124: How were serovars determined? In silico methods? Serological? In the Materials and Methods it states that SISTR and SeqSero2 were used for the ST313 dataset. Why were these methods not used for the larger 11,631 publicly available dataset as well? Enterobase already runs these analyses on all of their data – it would be nice to ensure that the first genomic dataset also included SISTR and SeqSero2 serovar classification.\nResponse: To clarify, serovars were determined in silico using SISTR and SeqSero2 for both the ST313 dataset (n = 3,115) and the larger 11,361-isolate iNTS dataset. We applied these tools uniformly across all genomes to ensure methodological consistency. We have explicitly stated this in the Methods section and have added the relevant details (including version numbers and citations for SISTR and SeqSero2) for the global dataset as well.\nNew lines 494-496:\nFor the global iNTS genomic dataset constructed above, we further confirmed the species identity as Salmonella using KmerFinder v.3.0.230, and serovar determination was performed with SISTR v.1.1.231 and SeqSero2 v.1.3.132.\nLines 133-141: It is appreciably challenging to gather accurate data from voluntarily reported metadata and I would recommend that the authors account for this in their manuscript. For example, panel 1e: sex data is only provided for 1/5th of the human isolates. Without understanding how these data were collected and by whom, this is potentially misleading.\nResponse: Only a very small fraction of human ST313 isolates have patient sex information available (440 out of 3,115, ~14%). Moreover, >99% of those sex data (439/440) come from a single source: the large Blantyre, Malawi study by Pulford et al. (Nat. Microbiol. 2021; DOI: 10.1038/s41564-020-00836-1), which recorded patient sex for thousands of invasive Salmonella cases. The remaining one isolate with sex information (female) was retrieved via a GenBank BioSample entry (accession SAMN11866370). We have revised the manuscript to make this clear.\nNew lines 147-151:\nOnly a small number of ST313 isolates contained sex information (n = 440), of which more than 99% were reported to have been collected at Queen Elizabeth Central Hospital (Blantyre, Malawi)10. Overall, the sex distribution was balanced, with slightly more male cases (n = 226) than female (n = 214) (Fig. 1e).\nLine 152: How was 2.4 identified and what cut-offs were used to distinguish it from the other sub-lineages? This is important for ensuring reproducibility of the results.\nResponse: Thank you for your comments. Sublineage 2.4 was identified through a population structure analysis of Salmonella ST313 assemblies using the rhierbaps package in R, which implements the hierarchical Bayesian clustering algorithm described in the fastbaps framework (Nucleic Acids Res. 2019; DOI: 10.1093/nar/gkz361). Specifically, the aligned ST313 cgSNP sequences were processed with the “hierBAPS” function, with the parameters set to a maximum depth of 2 (max.depth = 2) and an initial number of populations of 20 (n.pops = 20). This hierarchical approach partitions the dataset into nested clusters, where level 1 represents major lineages and level 2 defines sublineages. Sublineage 2.4 was identified as one such cluster based on optimized log marginal likelihood values.\nA key advantage of this phylogeny-free approach is that it does not rely on explicit cut-offs—such as predefined SNP distance thresholds—to differentiate sublineage 2.4 from others. Instead, cluster boundaries are determined through Bayesian model optimization, which maximizes the marginal likelihood under a uniform prior for the number of clusters. This methodology is consistent with previous Salmonella epidemiological studies, including those focusing on ST313 (Nat Microbiol. 2021; DOI: 10.1038/s41564-020-00836-1, Nat Commun. 2016; DOI: 10.1038/ncomms12827).\nWe have further revised the Methods section, it now reads:\nNew lines 547-556:\nA total of 3,115 global assemblies of ST313 underwent a population structure analysis using the RhierBAPS v.1.1.440, based on the Hierarchical Bayesian clustering algorithm and employing aligned cgSNPs as input. The calculation of population structure utilized R package v.4.3.1, supplemented by the R packages \"phytools\" v.2.0.3, \"ape\" v.5.7.1, and \"rhierbaps\" v.1.1.4. Specifically, the aligned ST313 cgSNP sequences were processed with the “hierBAPS” function, with the parameters set to a maximum depth of 2 (max.depth = 2) and an initial number of populations of 20 (n.pops = 20). This hierarchical approach partitions the dataset into nested clusters, where level 1 represents major lineages and level 2 defines sublineages.\nLines 166 and throughout: If BEAST was used, shouldn’t these statistics be HPDs rather than 95% confidence intervals as listed throughout the manuscript?\nResponse: You are correct—these ranges should be described as 95% highest posterior density (HPD) intervals rather than 95% confidence intervals. We have corrected this terminology throughout the manuscript to accurately reflect that BEAST produces HPD intervals. All instances of “95% CI” have been replaced with “95% HPD interval” in the revised text.\nLines 191-192: Couldn’t this just be due to international travel?\nResponse: We agree that international travel is a likely explanation for the observed genomic links. However, due to the absence of patient-specific travel histories, we cannot conclusively attribute recent ST313 Lineage 2 transmissions to travel. Nevertheless, we have now acknowledged in the manuscript that these intercontinental transmission events could very well be travel-related. For example, we note that the UK and African isolates in question were collected within a relatively short interval (~3 years) and likely reflect introduction via human travel. We also emphasize that genome-based links alone cannot definitively distinguish direct person-to-person transmission from importation via travelers. Accordingly, we have added the following sentences:\nNew lines 219-222:\nConsidering that these isolates belong to the same cluster and were collected within a relatively short time frame (less than three years), we speculate that the recent transmission events represented by these ST313 strains may be potentially associated with international travel.\nNew lines 465-469:\nFinally, genomic inference of intercontinental ST313 transmission may reflect human travel or migration rather than direct bacterial spread29. Genome-based links alone cannot definitively distinguish person-to-person transmission from human-mediated dispersal across continents.\nLines 198-199: Yes, but there is also bias in the sequencing. We don’t often look for or sequence Salmonella from sick animals unless there is a specific research project looking at this. I think it is safer for the authors to qualify this in their discussion.\nResponse: We recognize this bias. As described in our revised analysis above, we refined our SNP-based transmission tracing by using lineage-specific references and adding a time window criterion. This refinement did not alter our main conclusions (for instance, we still detected the rare human–canine transmission event in Lineage 5). We have now expanded the Discussion to address the potential sampling bias regarding animal infections. In particular, we note that current ST313 isolates are predominantly from human clinical samples and that the prevalence of ST313 in animals is likely under-investigated. This raises the possibility that ST313 has a broader host range than presently documented, underscoring the need for further One Health surveillance efforts:\nNew lines 399-408:\nRefined SNP-distance–based tracing, incorporating isolation time as epidemiological support (≤3 SNPs and ≤3 years apart), revealed evidence of cross-host transmission within UK-associated cluster 353. These findings suggest infrequent reverse zoonosis, including possible human–canine transmission. However, potential sampling bias should be considered, as current ST313 collections are predominantly derived from human clinical samples, and its prevalence in sick animals has received insufficient attention. This raises the possibility that ST313 may possess a broader host range than previously recognized, underscoring the need for further investigation from a One Health perspective\nLine 206: As written, it suggests that only pSLT was used in query, but then in line 210 the authors mention that pBT2 and pBT3 were also identified. Please rephrase line 206 so that it is clearer what was queried.\nResponse: This was indeed a labeling error on our part. We have corrected the phrasing to clarify that we queried all known plasmids associated with ST313, not just pSLT. The revised sentence now explicitly lists the plasmids we examined (pBT1–pBT3, pSLT-BT, pSLT, pAnK), along with the prophages, transposons, and Salmonella pathogenicity islands included in our analysis. We have also updated the Methods section (“Construction of the SPI and plasmid database carried by ST313”) to enumerate the reference sequences for each plasmid and SPI, ensuring that the scope of our query is clear.\nNew lines 232-237:\nTo characterize genomic differences—particularly within the accessory genome—among the newly identified ST313 lineages, we examined mobile genetic elements (MGEs) previously reported to be carried by ST313, including plasmids (pBT1–pBT3, pSLT-BT, pSLT, pAnK), prophages (BTP1, BTP3, BTP5, BTP6, Gifsy-1, Gifsy-2), transposons (Tn21), and Salmonella pathogenicity islands (SPIs; see Methods).\nNew lines 672-682:\nConstruction of the SPI and plasmid database carried by ST313\nA literature review was conducted to identify known plasmid and SPI carried by ST31310,48–50. The plasmid sequences pBT1–pBT3 and pSLT-BT were derived from the ST313 reference strain D23580, with GenBank accession numbers LS997975.1, LS997976.1, LS997977.1, and LS997974.1. The lineage-specific pSLT plasmid reference sequence was obtained from GenBank (accession number CP060170), and the reference sequence for pAnK was NC_010896.1. All SPI reference sequences (SPI-1, SPI-2, SPI-3, SPI-4, SPI-5, SPI-6, SPI-9, SPI-11, SPI-12, SPI-13, SPI-14, SPI-16) were extracted from the D23580 reference genome (GCA_900538085.1).\nLines 226-230: Did the authors control for low-sequenced regions by using thresholds for coverage across the genome or in a particular region of the genome?\nResponse: We conducted rigorous quality control on the 3,115 ST313 strains included in the genomic analysis, evaluating parameters such as GC content, N50, genome size, coverage, and completeness. Specifically, genome GC content, N50, genome size, and genome coverage were assessed using the Python 3 script provided by QUAST. The tool was executed with the reference sequence (the completely sequenced D23580 genome) as input to compute alignment-based metrics, including genome fraction, a measure of coverage across the genome that quantifies the percentage of the reference genome aligned by the assembly contigs. Genome completeness was evaluated using BUSCO, which assesses the presence of near-universal single-copy orthologs from lineage-specific datasets. The analysis was conducted in genome mode for prokaryotic sequences, employing Prodigal as the gene predictor. The enterobacteriaceae_odb lineage dataset was selected, comprising 874 marker genes, to ensure specificity for the bacterial taxon under study.\nWe established distinct genome quality control thresholds based on the Salmonella genome quality standards developed by the European Reference Laboratory, as well as insights from prior genomic epidemiology practices. Details are provided in Response Table 2.\nResponse Table 2. Genome assembly quality-control thresholds used in this study.\nQuality Control Parameters Threshold\nGC (%) 51 - 53\nN50 > 30,000\nGenome size (MB) 4.2 - 5.5\nGenome coverage (%) > 95\nGenome completeness (%) > 95\nThe results demonstrate that the ST313 genome dataset used in our analysis meets the quality criteria defined above (New Fig. S20). Specifically, the GC content of all strains ranged between 51% and 53%, the N50 values of all genomes exceeded 30,000, and the genome sizes were within the range of 4.2–5.5 Mb. In addition, the genome coverage relative to the reference sequence was greater than 95% for all strains, and the genome completeness of each strain also exceeded 95%.\nNew Fig. S20 Quality control results of 3,115 ST313 genomes. (a) GC content; (b) N50; (c) genome size; (d) genome coverage; and (e) genome completeness.\nIn addition, we expanded the Methods section in the manuscript, which now reads as follows:\nNew lines 514-526:\nFurthermore, we verified the completeness and coverage of each ST313 genome. Specifically, genome coverage was assessed using the QUAST v5.3.034. The tool was run with the fully sequenced D23580 genome as the reference to compute alignment-based metrics, including genome fraction, which quantifies the proportion of the reference genome aligned by the assembly contigs. Genome completeness was evaluated with BUSCO v6.035, which measures the presence of near-universal single-copy orthologs from lineage-specific datasets. Analyses were performed in genome mode for prokaryotic sequences, using Prodigal for gene prediction. The enterobacteriaceae_odb dataset, comprising 874 marker genes, was selected to ensure taxon-specific accuracy for the bacterial group under study. To ensure the reliability of downstream analyses, including those involving MGE detection, we retained only genomes with both coverage and completeness exceeding 95% (Fig. S20).\nMobile elements are inherently challenging to sequence and assemble from short-read sequencing data. Did the authors check to see that these were indeed absent? Or is it possible that these regions of the genome had low coverage at these sites and this could explain why they were not detected?\nResponse: We considered this possibility. However, our assemblies have >95% reference genome coverage and >95% completeness (see the aforementioned New Fig. S20), which gives us confidence that major mobile genetic elements (MGEs) would be captured in the data. Additionally, we examined the MGEs using long-read-sequenced reference genomes for each ST313 lineage: L1 (strain name: A130), L2 (strain name: D23580), L3 (strain name: BKQZM9), L4 (strain name: B3589), and L5 (strain name: U2) (sublineage 2.4 currently lacks a complete long-read assembly). In these fully sequenced genomes, we identified the key transposons and integrons present in each lineage and compared them to our short-read findings. The presence/absence patterns of MGEs in the long-read assemblies matched our short-read results, indicating that our short-read assemblies did not miss those elements due to coverage issues (Response Figure 5). Based on these analyses, we believe that if a particular integron or transposon was not detected in our data, it is unlikely to be solely due to insufficient short-read coverage.\nResponse Figure 5: Detection of transposons and integrons carried by the genome assemblies of five ST313 lineages. The maximum-likelihood phylogeny on the left illustrates the approximate positions of the five ST313 lineages. The dashed line indicates that L2.4 lacks an available genome assembly. The map on the right shows the dominant transposon and integron types within each lineage. The gray shaded areas represent the presence of the corresponding transposons or integrons, while the white areas indicate their absence.\nLine 355: I think caution should be exercised here when claiming that these are ‘spreading around the world.’ Without a more complete dataset, it is not possible to assert this. All that can be concluded is that isolation of ST313 happens outside of Africa as well.\nResponse: Thank you—we have revised this wording. In the Results, instead of suggesting that ST313 lineages are “spreading around the world,” we now state that our focused dataset of 3,115 ST313 genomes is still largely dominated by sub-Saharan African isolates, with only sporadic detections on other continents. We attribute those non-African cases to human-mediated dispersal (for example, travel-related introductions) without implying widespread global establishment. This rephrasing conveys that ST313 infections occur outside of Africa, but avoids over-asserting the extent of global spread:\nNew lines 390-394:\nNext, we assembled a focused dataset of 3,115 rigorously quality-controlled ST313 genomes spanning diverse geographic and temporal origins. In contrast to the wider iNTS landscape, ST313 infections appear predominantly restricted to sub-Saharan Africa, with sporadic detections elsewhere potentially attributable to human-mediated dispersal.\nLine 382: Using terms such as ‘the epicenter of L1’ could be misleading. The authors need to remember that there were targeted efforts to sample and sequence isolates from Malawi (over other African countries).\nResponse: We agree and have rewritten that part of the Discussion. We no longer describe Malawi as “the epicenter of L1.” Instead, we provide context about sampling and drug pressure:\nNew lines 425-439:\nAMR exerts a profound selective pressure on ST313 evolution, as evidenced by historical lineage replacements10. Chloramphenicol, which served as an empirical therapy for ST313 in Africa until its replacement around 200113, promoted resistance acquisition—predominantly in L2 via ARGs such as catA1 (88% carriage rate)—facilitating lineage displacement from L1 to L2. Our AMR profiling, stratified by year of isolation, demonstrates a robust correlation between the post-2006 contraction of L1 and reduced carriage of ARGs such as mcr-9, tet(D), aph(3')-Ia, sul2, aph(6)-Id, and aph(3\")-Ib. This pattern is consistent with the decline in ARG horizontal transfer frequencies observed among L1 isolates from Malawi. We further show that loss of key AMR drivers, such as Tn6205 and IncHI2A/rep_cluster_1088-type plasmids, is a major contributor to this phenomenon. Interestingly, this temporal alignment with Malawi’s post-2005 shift in clinical antimicrobial policies supports AMR modulation as a driver of lineage turnover, extending the evidence for resistance-fueled expansions.\nLines 398-400: Similar to my comments above, I think it is quite the opposite. There has been an intentional focus on sampling iNTS in African countries. I think this bias is apparent in the data and should be discussed in more detail here. This also contradicts what the authors state in lines 403-404.\nResponse: We have expanded our discussion of this sampling bias. As noted earlier, we now explicitly acknowledge that surveillance and sequencing efforts have been concentrated in certain African countries, which biases the available data. We also emphasize that our analysis is based on genomes rather than case incidence, so it inherently mirrors where sequencing has been done. Additionally, we reiterate that intercontinental genomic links (e.g., between Africa and Europe) may be explained by human travel or migration rather than direct local transmission. We have thoroughly rewritten this section, and it now reads as follows:\nNew lines 451-475:\nThis study has limitations. First, disparities in surveillance and reporting are likely to bias prevalence estimates toward specific regions. Because the available data are derived from WGS rather than reported case counts, this issue becomes even more apparent when assembling a global iNTS genomic dataset comprising 11,361 isolates. To quantify this sampling bias, we examined the relationship between a country’s economic level (measured as GDP per capita) and its corresponding sample count (Fig. S17). The results indicate that countries with greater public health infrastructure and resources—particularly the USA, Canada, and China—appear to have higher recorded incidences than those with more limited resources. At the same time, regions with fewer resources—such as Malawi, Congo, and South Africa—also seem to show higher incidences. This pattern reflects existing large, well-funded studies conducted specifically to investigate the burden of iNTS in select African countries. Second, our focus on genome-based transmission features excluded host factors such as HIV status or other comorbidities. Finally, genomic inference of intercontinental ST313 transmission may reflect human travel or migration rather than direct bacterial spread29. Genome-based links alone cannot definitively distinguish person-to-person transmission from human-mediated dispersal across continents. Nevertheless, by harnessing the most comprehensive iNTS and ST313 genomic collections to date, we deliver an integrated spatiotemporal framework, uncover novel lineages, and illuminate evolutionary pressures shaping this pathogen. These findings advocate for expanded, multisectoral surveillance encompassing non-African regions and animal reservoirs to capture emerging threats. Ultimately, such insights pave the way for refined global strategies to mitigate iNTS morbidity and mortality.\nLine 435: Please justify the use of this strain as reference for all genomes. It appears this is most suitable for Lineage 2 ST313 strains but its use as a reference for other lineages may not be valid and should be justified. The easiest way to address this would be to use the most complete genome as reference for each of the other lineages as reference for mapping SNPs among isolates in that lineage.\nResponse: You are correct, and we have implemented this change. In the revised analyses, we assigned lineage-specific reference genomes for variant calling and BEAST dating. Specifically, we selected five published complete ST313 genomes corresponding to Lineages 1–5 (A130 for L1, D23580 for L2, BKQZM9 for L3, B3589 for L4, and U2 for L5). We then re-ran the BEAST evolutionary analyses for Lineages 1, 3, 4 and 5 using these lineage-appropriate references.\nUsing strains from each lineage as the mapping reference did not alter the overall patterns of effective population size (N ) change—Lineages 1, 3, 5 still show sharp Ne reductions at their respective historical inflection points, and Lineage 4 still maintains a high Ne in Africa over time. However, switching to more refined references did slightly shift the inferred divergence times for some lineages. We have updated the Results to reflect these changes. For example, we now report that Lineage 4 (the newly defined lineage) has maintained a relatively high N in Africa since around 2000, whereas Lineages 1, 3, and 5 each experienced significant N declines at different times (L1 ~2000, L3 ~2010, L5 ~2017). We also provide revised estimates for the most recent common ancestors and earliest spread events of each lineage using the appropriate references for each. These adjustments address the reviewer’s concern by ensuring that each lineage’s analysis is internally consistent and not biased by using a Lineage 2 reference across the board.\nNew lines 153-200:\nGlobal population diversity of ST313 uncovers new lineages\nLineage assignment, based on previous studies9,10, identified five primary lineages: L1 (296/3,115), L2 (2,273/3,115), L3 (164/3,115), L4 (230/3,115), and L5 (152/3,115) (Fig. 2a). Among these, L1 and L2 predominated in Africa. Within L2, phylogenetic analysis revealed a novel sublineage beyond 2.0–2.3, designated 2.4 (61/3,115). The majority of isolates in sublineage 2.4 originated from Kenya (58/61, 95.1%), with two isolates also from the UK (Fig. 2b). Nearly half of L3 genomes (79/164, 48%) were from the Americas. L4, a previously unreported and newly defined lineage here, was concentrated in Africa (159/230, 69.1%). Interestingly, L5, previously referred to as the UK-isolated lineage, included isolates from Australia and Japan. To further assess whether the newly identified sublineage 2.4 and L4 represent single clonal expansion events, we calculated their within-lineage genetic diversity. Pairwise SNP distances within sublineage 2.4 ranged from 0 to 50 nucleotides (mean: 19.1), whereas L4 showed substantially greater divergence, ranging from 0 to 510 nucleotides (mean: 173.4). These findings indicate that neither sublineage 2.4 nor L4 arose from a single clonal expansion event (Fig. S2).\nInference of historical changes in the effective population sizes (N ) showed distinct declines for L1, L2, L3, and L5. Specifically, L1 exhibited a sharp reduction around the year 2000 (Fig. S3b), while declines in L2 and L3 occurred later, approximately around 2010 (Fig. S4b and Fig. S5b), and L5 around 2017. By contrast, the newly characterized sublineage 2.4 expanded continuously since emergence, while L4 maintained relatively high Ne since 2000 (Fig. 2e).\nPhylogeographic analysis suggested that the most recent common ancestor (MRCA) of L1 and L2 originated in Africa, approximately around 1950 (95% highest posterior density [HPD]: 1922-1966) and 1960 (95% HPD: 1955-1965), respectively (Fig. S3a and Fig. S4a). The earliest inferred transcontinental transmission event for L2 occurred around 1996 (95% HPD: 1995–1998), while the earliest inferred transcontinental transmission event for L1 occurred slightly later, around 2006 (95% HPD: 2003–2009). The earliest transcontinental transmission events of both lineages were potentially associated with Africa and Europe (Fig. S3c and Fig. S4c).\nThe MRCA of L3, L4, and L5 originated outside Africa. Specifically, the MRCA of L4 is estimated to have appeared around 1731 (95% HPD: 1575–1859) in Asia (Fig. 2c). In contrast, the MRCA of L3 emerged slightly later, around 1874 (95% HPD: 1787–1932) in Asia (Fig. S5a). The MRCA of L5 likely appeared around 1965 (95% HPD: 1943–1982) in Europe (Fig. 2d). The earliest transcontinental transmission events of these lineages appear to have occurred around 1899 (95% HPD: 1839–1941) in L3 strains, between Europe and Asia (Fig. S5c). The earliest transcontinental transmission events for the newly defined L4 occurred around 1908 (95% HPD: 1864–1946), possibly between Asia and Africa geographically. Finally, the earliest transcontinental transmission events for L5 occurred the latest, around 1993 (95% HPD: 1984–2002) (Fig. 2f).\nLines 574-575: Details about the ‘pre-trained random forest classifier’ need to be provided. Did the authors use an already developed classifier from another paper? Please cite that paper to be fair.\nResponse: We have added the requested details. The invasiveness index in our study follows the method (and uses the same pre-trained model) described by Langridge et al. (PLoS Genet. 2018; DOI: 10.1371/journal.pgen.1007333). Specifically, this approach employs 196 top predictor genes and integrates a machine-learning framework based on a random forest classifier combined with delta-bitscore functional variant calling. This method quantifies the predicted invasiveness of ST313 toward the host. The code used to calculate the invasiveness index is available at: https://github.com/Gardner-BinfLab/invasive_salmonella. In the revised manuscript, we have added the relevant citations and included the link to the pre-trained model used in our analysis in the Code Availability section.\nLines 582-587: This approach may identify PMSCs that would not render a gene a pseudogene. Additional filtering that considers the position of the frameshift is important and necessary for identification of pseudogenes. For example, PMSCs occurring at the beginning of a gene are more likely to result in loss of function (pseudogene formation) but PMSCs at the very end of a gene may not influence function. Additional filtering based on position is a way for reducing the likelihood (but not removing the possibility) that a gene would be incorrectly identified as a pseudogene.\nResponse: Thank you for this insightful suggestion. We agree that position-based filtering is necessary to avoid overestimating pseudogene formation. We have now examined the relative position of each “stop gained” or “frameshift” SNP within its gene. Following the conservative criterion of Ochman et al. (Nucleic Acids Res. 2005; DOI: 10.1093/nar/gki631), we do not classify mutations occurring in the last 20% of a gene as causing pseudogenization. We updated Supplementary Table 8 to report the gene length and the SNP position for each putative pseudogene mutation (for example, “533/722” indicates the SNP is at position 533 in a 722 bp gene). Applying this stricter filter removed 7,344 SNPs from our original pseudogene list that were located in the terminal 20% of genes. We have revised the Methods (“Identification of pseudogene-associated SNPs”) to describe this change. In the revised analysis, we only consider nonsense or frameshift mutations in the first 80% of a coding sequence as likely to result in pseudogene formation.\nNew lines 747-755:\nIdentification of pseudogene-associated SNPs\nWe employed Snippy v.4.4.5 to identify single-nucleotide polymorphisms in all ST313 strains. SNPs predicted to cause “Stop Gained” or “Frameshift Variant” mutations were initially classified as pseudogene-associated SNPs. To minimize potential overestimation, we further examined the positional context of each mutation within its corresponding gene. Following the conservative criterion described by Ochman et al., mutations located within the first 80% of a gene were considered likely to disrupt function, whereas those occurring in the terminal 20% were not classified as pseudogene-associated62.\nLine 840: Sex (male/female) is a more appropriate term here than gender (man/woman)\nResponse: Thank you for your comment. The corresponding sentence has been revised and now reads as follows:\nNew lines 1071-1073:\n(e) Statistics of the sex of patients from whom ST313 was isolated. The y-axis represents the number of ST313 genomes from patients of each sex.\nFigure 1d: Unknown is spelled wrong;\nResponse: The label “Unknwon” in Figure 1d has been fixed to “Unknown” in the revised figure. All figure labels have been carefully proofread to prevent such errors.\nFigure 1. Geographic, source, and temporal distribution of the ST313 genomic dataset (n = 3,115). (a) Geographic distribution of ST313 in the constructed dataset. Blue circles denote the number of ST313 genomes in different countries, with larger circles representing a larger number of genomes. (b) The number of ST313 isolations in the constructed dataset across six continents over multiple years. The y-axis represents the number of genomes, and different colors indicate different continents. (c) Statistics on the location and sampling time of ST313 genomes isolated from \"Human,\" \"Food,\" and \"Animal\" sources. (d) Statistics of the source of human-isolated ST313 in the dataset. (e) Statistics of the sex of patients from whom ST313 was isolated. The y-axis represents the number of ST313 genomes from patients of each sex.\nI also question the inclusion of these strains in this dataset. Given that many Salmonella infections are diagnosed from fecal samples, it is likely that at least some of these strains indeed came from feces. Also, the authors stated that ‘feces’ was excluded, but their second most common known source is ‘faeces’ which is just a different spelling of ‘feces.’ Therefore, this dataset likely contains non-invasive infections as well. It is also unclear to me why fecal isolates are included in this dataset, whereas in the larger iNTS dataset they are not?\nResponse: In this study, we constructed two genomic datasets: a global iNTS dataset (11,361 isolates) and an ST313-specific dataset (3,115 isolates). These two datasets were assembled using distinct criteria. The global iNTS dataset, by definition, excluded isolates from fecal samples (focusing only on invasive isolates from sterile sites). By contrast, our ST313 dataset was intended to be a comprehensive collection of ST313 genomes, regardless of source. Thus, the ST313 dataset does include isolates from fecal sources and from non-human hosts (e.g., the “faeces” entries), provided they passed our quality filters. We chose not to restrict the ST313 dataset to only invasive cases because our goal was to analyze ST313’s overall population structure, evolution, and transmission potential. This inclusive approach to ST313 data has precedent (e.g., Pulford et al., 2021 analyzed ST313 from various sources).\nTo avoid confusion, we have substantially rewritten the Methods section to separately describe the assembly of the global iNTS dataset versus the ST313 dataset. The corresponding revised text in the Methods section now reads as follows:\nNew lines 478-532:\nAssembly of the global iNTS genomic dataset\nIn this genomic epidemiology study, iNTS was defined as non-typhoidal and non-paratyphoidal Salmonella isolates obtained from human specimens of non-fecal origin, including blood, urine, abscesses, and other sterile sites, as described in previous literature2. We subsequently retrieved all available Salmonella genomes from EnteroBase and GenBank (data cutoff: 1 March 2024; total number: 488,232). For metadata, when BioSample information was available in EnteroBase, we obtained the corresponding metadata directly from the associated BioSample record, unless the isolate was deposited exclusively in EnteroBase without a linked BioSample entry. We excluded all strains with ambiguous isolation sources—for example, those labeled only as “clinical” or “human”—and retained only those isolates with clearly defined sources that met the above criteria for iNTS. Information on the country and date of isolation was also collected. This process resulted in a final dataset of 11,361 iNTS genomes (Supplementary Table 1).\nFor the global iNTS genomic dataset constructed above, we further confirmed the species identity as Salmonella using KmerFinder v.3.0.230, and serovar determination was performed with SISTR v.1.1.231 and SeqSero2 v.1.3.132. In terms of isolation sources, most isolates were recovered from blood (n = 7,922; 69.7%) and urine (n = 2,971; 26.2%). Additional isolates primarily originated from abscesses (n = 116), wounds (n = 81), cerebrospinal fluid (n = 74), the respiratory tract (n = 36), and synovial fluid (n = 30) (Fig. S18).\nAssembly of the ST313 genomic dataset\nNext, we constructed a global dataset comprising 3,115 ST313 genomes (Supplementary Table 2). We not only included human-derived ST313 data but also incorporated isolates from animals, plants, and those with an \"Unknown\" source label. Through extensive literature searches, we confirmed that 1,584 of the 3,115 ST313 isolates have been explicitly reported in 20 previously published studies (Fig. S19). Genome quality was assessed using the European Reference Laboratory’s Salmonella genomic quality control standards33. A Python3 script was employed to exclude genomes with over 500 contigs and an N50 value less than 30,000, as well as those with a GC content outside the 51%-53% range or a genome size between 4.2 and 5.5 MB. Furthermore, we verified the completeness and coverage of each ST313 genome. Specifically, genome coverage was assessed using the QUAST v5.3.034. The tool was run with the fully sequenced D23580 genome as the reference to compute alignment-based metrics, including genome fraction, which quantifies the proportion of the reference genome aligned by the assembly contigs. Genome completeness was evaluated with BUSCO v6.035, which measures the presence of near-universal single-copy orthologs from lineage-specific datasets. Analyses were performed in genome mode for prokaryotic sequences, using Prodigal for gene prediction. The enterobacteriaceae_odb dataset, comprising 874 marker genes, was selected to ensure taxon-specific accuracy for the bacterial group under study. To ensure the reliability of downstream analyses, including those involving MGE detection, we retained only genomes with both coverage and completeness exceeding 95% (Fig. S20).\nNext, we used SISTR v.1.1.2 and SeqSero2 v.1.3.1 to determine the serovar of all 3,115 genomes, ensuring they were classified as S. Typhimurium. Additionally, MLST was performed using MLST v.2.22 with the senterica_achtman_2 scheme36. Only strains identified as ST313 were included in subsequent analyses.\nFigure 4a: In the pseudo heat map, it looks like there is perfect 0 or 100% correlation between the isolates having a given element (BTP3 for example) or not. Yet in the text percentages are listed (11.8% of L1 carry BTP3). As presented in the figure, it looks like none of the isolates have BTP3. Additional gradients should be used as presented the data in the text does not match the figure. The same applies to Figure S9.\nResponse: Thank you very much for your comment. In the pseudo–heat map shown in Figure 4a, we used two colors—white and gray—to indicate the carriage rates of different MGEs within each lineage. Specifically, when the carriage rate of a given MGE in a lineage is below 60%, we represent it with white; when the carriage rate exceeds 60%, we represent it with gray. Therefore, although the carriage rate of BTP3 in L1 is 11.8%, the corresponding block appears white. To improve the readability of Figure 4a and Figure S9, we have revised both figures. The updated versions of Figure 4 and Fig. S9 (now corresponding to Fig. S11) are provided below:\nFigure 4. Different lineages have distinct mobilome repertoires, which are associated with the resistome. Statistical analysis was performed using the Mann-Whitney test, with significance defined as P < 0.05. (a) The maximum-likelihood phylogeny on the left illustrates the approximate positions of the ST313 lineages. The map on the right illustrates the virulence plasmid (pSLT), prophages, Salmonella pathogenicity islands (SPIs), and Tn-21 elements compiled from the literature and associated with different ST313 lineages. The color blocks on the right use gray and black to indicate whether a specific mobilome is widely carried by the corresponding ST313 lineage, with \"widely carried\" defined as a carriage rate exceeding 60% within the lineage. (b) Using a more comprehensive mobilome database, we performed mapping of integrons, transposons, prophages, and plasmids present in ST313 genomes. In parallel, the average number of antimicrobial resistance genes (ARGs) associated with each of the four mobilome types was calculated. (c) The average number of ARGs carried by the integrons, transposons, prophages, and plasmids in ST313 lineages 1 and 2. (d) Comparison of the number of ARGs carried by plasmids in ST313 sublineage 2.4 versus sublineages 2.0–2.3. (e–h) The distribution of resistome across mobilomes in the 3,115 ST313 strains, with e-h representing integrons, transposons, plasmids, and prophages, respectively. Red font is used to highlight the categories of resistome, while black font indicates the names of the various mobile genomic elements. The connecting lines between the resistome and the mobilome illustrate their potential association.\nFig. S11 Mobilome distribution across ST313 lineages. The maximum-likelihood phylogeny on the left indicates the approximate positions of the five ST313 lineages. The heatmap on the right shows dominant integron, transposon, prophage, and plasmid types within each lineage. Gray color blocks indicate whether a given mobilome is carried by >60% of strains in that lineage.\nIn panel 4c, it looks like the comparison is with other subspecies (subsp.)? Please confirm that this should indeed be labeled ‘other subsp.’ and not other lineages.\nResponse: Thank you for your comment. We assume that the panel you are referring to is panel 4d, which compares sublineage 2.4 with the other sublineages (2.0–2.3). We have redrawn the labels in the corresponding panel of Figure 4. The revised Figure 4 is shown below:\nFigure 4. Different lineages have distinct mobilome repertoires, which are associated with the resistome. Statistical analysis was performed using the Mann-Whitney test, with significance defined as P < 0.05. (a) The maximum-likelihood phylogeny on the left illustrates the approximate positions of the ST313 lineages. The map on the right illustrates the virulence plasmid (pSLT), prophages, Salmonella pathogenicity islands (SPIs), and Tn-21 elements compiled from the literature and associated with different ST313 lineages. The color blocks on the right use gray and black to indicate whether a specific mobilome is widely carried by the corresponding ST313 lineage, with \"widely carried\" defined as a carriage rate exceeding 60% within the lineage. (b) Using a more comprehensive mobilome database, we performed mapping of integrons, transposons, prophages, and plasmids present in ST313 genomes. In parallel, the average number of antimicrobial resistance genes (ARGs) associated with each of the four mobilome types was calculated. (c) The average number of ARGs carried by the integrons, transposons, prophages, and plasmids in ST313 lineages 1 and 2. (d) Comparison of the number of ARGs carried by plasmids in ST313 sublineage 2.4 versus sublineages 2.0–2.3. (e–h) The distribution of resistome across mobilomes in the 3,115 ST313 strains, with e-h representing integrons, transposons, plasmids, and prophages, respectively. Red font is used to highlight the categories of resistome, while black font indicates the names of the various mobile genomic elements. The connecting lines between the resistome and the mobilome illustrate their potential association.\nFigure 5c: It looks like there was no variance in mcr-9, tet(D), or aph(3’-Ia) in average number of ARGs over time? Is this true? Or the CI overlaps with the line?\nResponse: We thank the reviewer for highlighting this aspect of Figure 5c. The apparent lack of variance in the mcr-9, tet(D), and aph(3')-Ia data is due to the loss of these ARGs in ST313 Lineage 1 since 2006, particularly since 2008, when the prevalence of these ARGs in Lineage 1 consistently remained at its lowest levels. We have also provided the annual carriage rate of mcr-9, tet(D), and aph(3')-Ia in Lineage 1 since 2006 in Response Table 3 below.\nResponse Table 3: The annual carriage rates of the mcr-9, tet(D), and aph(3')-Ia resistance genes in ST313 Lineage 1 since 2006.\nFigure S10: AMR determinants not Determins\nResponse: Thank you for your comment. The revised Fig. S10 is shown below:\nFigure S10. ARG determinants in 3,115 Salmonella ST313 isolates. A phylogenetic tree was constructed based on core-genome single-nucleotide polymorphisms from these isolates, with left branches categorized into distinct lineages using BAPs analysis. Each lineage is indicated by uniquely colored branches. The heatmap on the right illustrates the geographic distribution of ST313 isolates by continent, the presence of resistance-associated point mutations in genes (gyrA, gyrB), and the corresponding antimicrobial resistance genes.\nFigure S14: Was correction for multiple comparisons testing done? I suspect not as the SE bars for L4 overlap with L3, yet the comparison suggests that the difference in pseudogene SNPs is significant.\nResponse: Thank you for your valuable comment. Following your suggestion, we incorporated the positional information of “Stop Gained” and “Frameshift Variant” mutations into the revised manuscript. Across 3,115 ST313 strains, we identified 24,926 potential pseudogene-associated SNPs. We quantified the number of pseudogene-associated SNPs carried by isolates from Lineages 1–5 and subsequently re-evaluated our statistical approach.\nTo assess whether the data followed a normal distribution, we first applied the Shapiro–Wilk test. The results indicated that the data did not conform to a normal distribution (p < 0.05), suggesting that non-parametric methods would be more appropriate for subsequent statistical analyses. Accordingly, we performed multiple hypothesis testing using the Mann–Whitney U test to compare the number of pseudogene-associated SNPs among different lineages, and adjusted the resulting p-values using the Benjamini–Hochberg false discovery rate (FDR) procedure. The results showed that the difference between Lineage 3 and Lineage 4 remained significant after FDR correction (New Fig. S16). Moreover, after adjustment, the standard error (SE) bars for Lineage 3 and Lineage 4 no longer overlapped, further supporting the robustness of this difference. At the same time, we have also revised the corresponding figure legend, which now reads as follows:\nNew Fig. S16. Pseudogene-associated SNPs in ST313. The number of pseudogene-associated SNPs (including nonsense and frameshift mutations) identified across different lineages. Multiple comparisons were performed using the Mann–Whitney U test, and the resulting p-values were adjusted using the Benjamini–Hochberg false discovery rate (FDR) procedure. p-values less than 0.05 were considered statistically significant (*)."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The responses given by authors as rebuttal to comments and queries from my review are well detailed and adequate to address my concerns. The additional explanations in form of accompanying figures was very helpful", "reviewer_2": "I appreciate the authors’ detailed response and the substantial revisions addressing several concerns. These additions improve transparency and strengthen the methodological description.\nHowever, I remain concerned that the clustering strategy, as described, may not fully support the transmission-related inferences. Specifically, the manuscript states that clustering was performed using UPGMA (average linkage) under a 3-SNP threshold and was used to ensure that “the SNP distances among isolates within each cluster did not exceed three SNPs” (Lines 613–617). For my understanding, UPGMA/average-linkage hierarchical clustering does not necessarily enforce a strict maximum pairwise SNP-distance constraint within each resulting cluster. It is therefore unclear whether all isolate pairs within the reported “≤3 SNP clusters” truly satisfy a hard cutoff. Because downstream identification and counting of “transmission events” depends directly on these clusters and the chosen threshold, this ambiguity could influence cluster membership and bias transmission estimates. More generally, even with an added temporal window, interpreting genetic proximity under a fixed SNP threshold as discrete “transmission events” risks overstating the epidemiological resolution achievable from SNP-distance clustering alone, especially under incomplete sampling.\nThese issues also bear on the manuscript’s broader interpretation that “phylogenomic analyses revealed patterns consistent with intercontinental dissemination of ST313, predominantly associated with human movement.” While the observed intercontinental mixing of lineages may be compatible with human-mediated dissemination, the phrasing “predominantly associated with human movement” implies a level of directionality and causal attribution that is difficult to support using genomic data alone—particularly if between-country patterns are inferred largely from SNP-threshold-based event counting.\nAlternative explanations (e.g., repeated independent introductions via common travel corridors, outbreak-driven sampling, and unsampled intermediate transmission chains) could plausibly generate similar phylogenetic and geographic patterns.\nCritically, the transmission analyses depend jointly on SNP thresholds and metadata-based constraints (sampling time and geographic origin). Although the authors have made commendable efforts in the revision to cross-check and refine the metadata, residual uncertainty in metadata accuracy and sampling representativeness could still substantially affect inferred between-region links. Because these interpretations depend strongly on threshold-based assumptions and the accuracy/representativeness of the underlying metadata, they may be misleading without independent epidemiological corroboration. While I appreciate the authors acknowledging these limitations, the key conclusions and terminology in the main text still appear stronger than warranted.", "reviewer_3": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_4": "The authors have put in a substantial amount of work revising and strengthening their manuscript. I thank the authors for addressing my initial suggestions as well as those of the other reviewers; this was not a light task.\nWith the additional attention to the metadata, and the authors’ findings that the data were indeed different from their original analysis, I think it is safer to phrase it as ‘potential’ or ‘putative’ transmission events (such as in line 1101) figure 3 and throughout the manuscript.\nLine 222: could also be foodborne (eating food contaminated with the ST313 strain)\nLine 326: It would be helpful to clarify if the ant(3”)-la was confirmed to be present on the MGEs rather than just correlated with. Judging by the materials and methods section, it was indeed co-located with the MGE and not just correlated with it, which is stronger and the manuscript text should be revised to explicitly state this.\nLine 467: Please add ‘global distribution of food supply’ or something similar as well here as this is an important source of infections as well\nFigure 2: Since none of the strains are MDR, is there a reason for including this column in the two figures 2c and 2d?\nFigure 3: I still find it safer to suggest that these are ‘potential’ transmission events. Without strong epidemiological evidence, I find it challenging to report these as fact."}, "author_response": "Response: We would like to thank the reviewer and the Early Career Researcher who co-reviewed this manuscript. We sincerely appreciate their time, effort, and insights, which have helped improve our work. All the requested changes are highlighted in blue."}, {"round": 3, "reviewer_comments": {"reviewer_2": "All of my comments and concerns have been fully addressed by the authors.", "reviewer_3": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": ""}]

Decision letter