Introduction

The acquisition of antimicrobial resistance (AMR) by pathogens is well-established as one of the most severe threats of the 21st century (1–3). It is estimated that over 10 million people could die per year due to AMR by 2050 (4). Antimicrobial resistance genes (ARGs), collectively known as the resistome, play a pivotal role in AMR development, progression, and amplification (5). Pathogens typically acquire resistome through mobilomes, like plasmids, via horizontal gene transfer (HGT) (6, 7). From an evolutionary perspective, the dynamic change in the resistome is crucial for enhancing pathogen fitness and enabling it to colonize new ecological niches (8). Key factors contributing to resistome alterations include disease burden, human activities, climate change, and geographical selection forces (9–13). Understanding the evolutionary trajectory of the AMR and related reservoirs is necessary, as disease management strategies will differ substantially (14). However, knowledge regarding the overall profile of pathogen regional adaptation related to the stepwise dynamics of the resistome remains limited.

Before recent advancements in whole-genome sequencing (WGS) technology, little was understood about resistome diversity that affects pathogen endemicity (15). Traditional typing methods generally exhibit low resolution, making monitoring and quantitatively comparing horizontal ARG transfer events challenging. The scarcity of appropriate models also poses a limitation. Salmonella, a pathogen with distinct geographical characteristics, has more than 90% of its serovars frequently categorized as geo-serotypes (16). Among the thousands of geo-serotypes, Salmonella enterica serovar Gallinarum (S. Gallinarum) is an avian-specific pathogen that causes severe mortality, with particularly detrimental effects on the poultry industry in low- and middle-income countries (17, 18). As once a globally prevalent pathogen in the 20th century, S. Gallinarum was listed by the World Organization for Animal Health (WOAH) and gradually became an endemic pathogen with sporadic outbreaks following the implementation of eradication programs in most high-income countries, making it a perfect model for study (19, 20).

Nowadays, antimicrobial therapy remains a priority choice against S. Gallinarum infections (21). The overuse or misuse of antimicrobials has led to increased epidemiological escalation of S. Gallinarum, with a regionally higher risk of AMR, especially for sulfonamides, penicillin, and tetracyclines (22, 23). To fill the gaps in understanding the evolution of S. Gallinarum under regional-associated AMR pressures and its adaptation to endemicity, we collected the most comprehensive set S. Gallinarum isolates, consisting of 580 genomes, spanning the period from 1920 to 2023. Using such unique WGS dataset, we investigated: 1) the population structure and potential geographical transmission history of S. Gallinarum at a single base level; 2) the dynamic resistome and mobilome changes that are associated with S. Gallinarum transition to an endemic variant; and 3) horizontal resistome transfer frequency and pattern between distinct regions.

Results

Global distribution of S. Gallinarum links with distinct sublineages

To understand the global geographic distribution and genetic relationships of S. Gallinarum, we assembled the most comprehensive S. Gallinarum WGS dataset (n=580), comprising 535 publicly available genomes and 45 newly sequenced genomes (24). The core-genome single nucleotide polymorphism sites (cgSNPs) were obtained by using the fully sequenced genome of S. Gallinarum R51 as a reference. Through hierarchical Bayesian analysis of cgSNPs, it was confirmed that S. Gallinarum divides into three biovars: S. Gallinarum biovar Pullorum (bvSP) (n=528/580, 91.03%), S. Gallinarum biovar Gallinarum (bvSG) (n=50/580, 8.6%), and S. Gallinarum biovar Duisburg (bvSD) (n=2/580, 0.34%). An association was identified between the biovar type of S. Gallinarum and its global geographic distribution. Most bvSP isolates were from Asia (436/528), with fewer occurrences in Europe (54/528), South America (12/528), and North America (5/528). In bvSG, South America (18/50, 36%) is the primary source of isolates, followed by North America (9/50, 18%), Europe (9/50, 18%), Africa (7/50, 14%), and Asia (3/50, 6%) (S1 Fig).

From a lineage perspective, bvSP can be further classified into five lineages: L1 (10/528), L2a (45/528), L2b (163/528), L3b (169/528), and L3c (141/528). Importantly, L3a, previously considered distinct, has been revealed as a sub-lineage of L3c. The predominant lineage types vary across different continents. Regarding predominant bvSP sources, Asia stands out, with L3c, L3b, and L2b being the main lineage types. On the other hand, L3b and L2a are more widespread in Europe and America, respectively (Fig 1a). For the bvSP strains from Asia included in our dataset, we found that all originated from China. To further investigate the distribution of bvSP across different regions in China, we categorized them into three distinct regions: southern, eastern, and northern (Supplementary Table 3). Our findings indicate significant variations in bvSP across various regions from geo-temporal aspects. Prevalence was highest in the eastern region (233/436), with lower incidences in the northern (137/436) and southern regions (65/436). Interestingly, the predominant bvSP lineage type also varied. In eastern and southern China, L2b and L3c exhibited as the predominant lineage types, while in northern China, L3b and L3c held nearly equivalent positions (Fig 1c).

Figure 1.

An analysis of the temporal prevalence of bvSP in China revealed a gradual replacement pattern of the lineages. Before the 2000s, the predominant lineage type was L3b. Then, from the 2000s to the 2010s, L3b and L1 were continually replaced by L3c, with L2b becoming the predominant lineage after the 2010s (Fig 1b). However, the replacement pattern varied among different regions. Specifically, the same pattern was observed in eastern China, but for northern and southern China, only replacements from L3b to L2b and L1 to L3c were observed (S2 Fig). These findings might indicate accelerated localized adaptation of bvSP.

Genomics portrait bvSP historic transmission

Considering the historical pandemic status of bvSP, we then investigated its global geographic transmission routes to understand its evolution into an endemic pathogen. L2b and L3b were identified as the dominant global lineages due to their antimicrobial resistance risks, and potential intercontinental transmission events. Other lineages were also analyzed using a Bayesian model with a relaxed molecular clock to infer their historical evolution (S3-S5 Figs). The temporal structure was verified (S6 Fig). Our findings show that the origin of L3b in China can be traced back to as early as 1683 (95% CI: 1608 to 1839). In contrast, the earliest possible origin of L2b in China dates back to 1880 (95% CI: 1838 to 1902) (Fig 2). Then, we specifically estimated the time points for the first intercontinental transmission events for the two major lineages. Our results indicate that L2b, likely underwent two major intercontinental transmission events. The first occurred around 1893 (95% CI: 1870 to 1918), with transmission from China to South America. The second major transmission event occurred in 1923 (95% CI: 1907 to 1940), involving the spread from South America to Europe. In contrast, the transmission pattern of L3b, an S. Gallinarum lineage originating in China, appears relatively more straightforward. It underwent only one intercontinental transmission event, from China to Europe, likely around 1790 (95% CI: 1661 to 1890) (S7 Fig).

Figure 2.

Endemic isolates further support localized dissemination by bvSP

To investigate the dissemination pattern of bvSP, we obtained forty-five newly isolated bvSP from 734 samples (6.1% overall isolation rate) collected from diseased chickens at two farms in Yueqing and Taishun, Zhejiang Province. Those bvSP can be classified into two dominant sequence types (STs): ST3717 and ST92. Isolates from Yueqing were all identified as ST3717, while most Taishun isolates were characterized as ST92. Only the E404 isolate was identified as ST2151 (Fig 3a). Of note, a lineage-preferential association was also observed, with bvSP isolated from Taishun belonging to L3b and those from Yueqing belonging to L2b. To further understand the genetic relationship among isolates in the two regions, we calculated cgSNPs distances between isolates. Isolates with cgSNP distances less than five were determined as the threshold to be involved in the same transmission event. Interestingly, genetic clustering revealed that isolates from the same farm exhibited a significantly close genetic correlation (Fig 3b). No historical transmission events of bvSP were found between Taishun and Yueqing.

Figure 3.

Other genetic characteristics also exhibited a similar pattern. The ARG and plasmid profile disclosed that IncFII(S) and ColpVC plasmids exhibited the highest carriage percentages (100%) among the 45 bvSP isolates. The IncX1 plasmid was exclusively identified in the E404 isolate. Furthermore, E404 harboured two distinctive resistance genes, bla_TEM-1B and sul2, in contrast to the absence of ARGs in the remaining bvSP isolates. We also estimated the invasiveness index to assess the invasiveness of isolated bvSP and showed that bvSP from the same farm had comparable invasive abilities. However, isolates from Taishun demonstrated a higher invasive ability than those from Yueqing (Fig 3c, Supplementary Table 6).

Furthermore, we simulated potential transmission events between the bvSP strains isolated from Zhejiang Province (n=95) and those from China with available provincial information (n=435) using five SNPs as the threshold. As a result, the transmission events showed a strong geographic-preferential distribution. We identified a total of 91 potential transmission events, all of which occurred exclusively within Zhejiang Province. No inter-provincial transmission events were detected, supporting the statement that bvSP in China is a highly localized pathogen (S8 Fig, Supplementary Table 8).

Mobilome and resistome distinct lineage-preferential association in a regional manner

The mobilome drives partitions of the resistome and plays a pivotal role in shaping the ecological niches and local adaptation of bacteria. Therefore, a quantitative evaluation of the resistome and mobilome is necessary to enhance the understanding of the localized distribution of bvSP. For S. Gallinarum, a total of 13 classes of ARGs, which can be further classified into six categories, were identified. The results revealed that bvSP exhibited a significantly greater resistome than bvSG and bvSD (Fig 4a). Among them, sul2 (196/528, 37.1%) has the highest prevalence, followed by bla_TEM-1B (183/528, 34.7%) and tet(A) (104/528, 19.7%). Moreover, the resistome also demonstrated a lineage-preferential association. We observed that L3b is more inclined to carry bla_TEM-1B and sul2 than other lineages, whereas tet(A) exclusively exists within L3 (Supplementary Table 4). The diversity of resistome carried by L3 may elevate the risk of multi-drug resistance (MDR). Notably, AMR risks vary among different regions of China, with the highest risk observed in southern China, followed by the northern and eastern regions (Fig 4c). Interestingly, tet(A), bla_TEM-1B, and sul2 were predominant resistome types prevalent across all regions of China. Meanwhile, aadA5 was typically more commonly found in southern China (Fig 4b).

Figure 4.

Four categories of the mobilome—prophage, plasmid, transposon, and integron—were characterized, revealing strong lineage-specific patterns (S9a Fig, Supplementary Table 5). In silico analyses revealed prophages emerge as the most diverse mobilome type, with Entero_mEp237, Escher_500465_1, and Klebsi_4LV2017 being more likely associated with L1, L2a, and L2b, respectively, while the carrier percentage of Escher_pro483 is significantly higher among L3c. Interestingly, Salmon_SJ46, Gifsy_2, Escher_500465_2, and Shigel_SfIV were consistently observed across all bvSP lineages and were found to be duplicated on the chromosome. For plasmids, IncFII(S) was predominantly carried across all bvSP lineages (520/528, 98.5%) (S9b Fig). Notably, ColpVC had the lowest carriage percentage among L1 (2/10, 20%), whereas IncX1 was explicitly associated with L3c (140/141, 99.3%).

The diversity of transposons and integrons is comparatively lower in bvSP. Specifically, we identified seven types of transposons and four types of integrons. Transposons were found to be more abundant, comprising a total of 255 in 528 bvSP. Among them, Tn801 (185/528, 35%) had the highest carriage rate, followed by Tn1721 (48/528, 9.1%). Interestingly, L3c is more likely to carry Tn801 and Tn1721 transposons than other lineages. In contrast, a total of 42 integrons were identified among the 528 bvSP, all of which belong to class I: In498 (n=22), In1440 (n=14), In473 (n=4), and In282 (n=2). The integrons also exhibited a lineage-preferential distribution, predominantly carried by L3b.

From a geographic distribution perspective, bvSP from northern and southern China carried the most extensive mobilome, followed by those from eastern China. Interestingly, bvSP isolated from northern China typically have more diverse mobilome types (S10 Fig).

Plasmid and transposon guide resistome geo-temporal dissemination

Variations in regional antimicrobial use may result in uneven pressure for selecting AMR. The mobilome is considered the primary reservoir for spreading resistome, and a consistent trend between the resistome and the mobilome has been observed across different lineages, from L1-L3c. We observed an overall gradual rise in the resistome quantity carried by bvSP across various lineages, correlating with the total mobilome content (S11 Fig). Furthermore, we investigated the interplay between particular mobile elements and resistome types in bvSP.

In bvSP, the predominant resistome is mainly associated with certain types of mobilome (Fig 5a). Specifically, regarding bla_{TEM_1B}, tet(A), and sul1, which are highly prevalent among bvSP, we found most bla_{TEM_1B} were carried by the transposon Tn801 and Tn1721, the plasmid IncX1, and the prophage SJ46. Additionally, tet(A) was primarily associated with the transposon Tn1721 and the plasmid IncX1. Integrons also facilitated the dissemination of the resistome. For instance, In498 and In1440 facilitate the dissemination of sul1, leading to quinolone resistance. Our results indicate that plasmids and transposons are the predominant reservoirs for the resistome in bvSP, with prophages and integrons following closely behind (Fig 5b). Interestingly, we observed that the primary reservoir for the resistome shows regional variations. In the southern and northern regions of China, the primary reservoirs for bla_{TEM_1B} are Tn801, Tn1721, IncX1, and prophage SJ46. However, Tn1721 was less frequently found in bvSP, which carries bla_TEM-1B from the eastern region. Similarly, the reservoir for tet(A) exhibits greater diversity in the east region but less diversity in the southern region (Fig 5c-e).

Figure 5.

Horizontal transfer of resistome occurs widely in localized bvSP

Horizontal transfer of the resistome facilitates the acquisition of AMR among bacteria, which may record the distinct acquisition event in the bacterial genome. To compare these events in a geographic manner, we further investigated the HGT frequency of each ARG carried by bvSP isolated from China and explored the HGT frequency of resistome between three defined regions. Potentially horizontally transferred ARGs were defined as those with perfect identity (100% identity and 100% coverage) and were located on MGEs across different strains (Fig. 6a). We first categorized a total of 621 ARGs carried by 436 bvSP isolates collected in China and found that 415 ARGs were located on MGEs (Supplementary Table 9). After excluding the ARGs not associated with MGEs, our findings reveal that horizontal gene transfer of ARGs is widespread among Chinese bvSP isolates, with an overall transfer rate of 92%. Specifically, 50% of the ARGs exhibited an HGT frequency of 100%, indicating that these ARGs might underwent extensive frequent horizontal transfer events (Fig. 6b). It is noteworthy that certain resistance genes, such as tet(A), aph(3’’)-Ib, and aph(6)-Id, appear to be less susceptible to horizontal transfer.

Figure 6.

However, different regions generally exhibited a considerable difference in resistome HGT frequency. Overall, bvSP from the southern areas in China showed the highest HGT frequency (HGT frequency=95%). The HGT frequencies for bvSP within the eastern and northern regions of China are lower, at 92% and 91%, respectively (Fig. 6c). For specifical ARG type, we found tet(A) is more prone to horizontal transfer in the southern region, and this proportion was considerably lower in the eastern region. Interestingly, certain ARGs such as aph(6)-Id, undergo horizontal transfer only within the eastern and northern regions of China (Fig. 6d). Notably, as a localized transmission pathogen, resistome carried by bvSP exhibited a dynamic potential among inter-regional and local demographic transmission, especially from northern region to southern region (HGT frequency=93%) (Fig. 6e, Supplementary Table 7).

Discussion

Over the past decades, the emergence of AMR has garnered significant attention from both public health agencies and the poultry industry (25). Region-associated pressures shaped the dissemination patterns of bacterial pathogens, transitioning from pandemic to endemic, and led to dynamic changes in the resistome. Here, utilizing S. Gallinarum as a prototype for an endemic pathogen that once spread globally, and assembling the most comprehensive S. Gallinarum WGS database, we have discovered that region-associated pressures can influence the AMR risks of S. Gallinarum within lineage- and region-specific patterns. These pressures might also influence lineage-distinct evolutionary structures and transmission histories of S. Gallinarum.

The biovar types of S. Gallinarum have been well-defined as bvSP, bvSG, and bvSD historically (26). Among these, bvSP can be further subdivided into five lineages (L1, L2a, L2b, L3b, and L3c) using hierarchical Bayesian analysis. Different sublineages exhibited preferential geographic distribution, with L2b and L3b of bvSP being predominant global lineage types with a high risk of AMR. The historical geographical transmission was verified using a spatiotemporal Bayesian framework. The result shows that L3b was initially spread from China to Europe in the 18^th-19^th century, which may be associated with the European hen fever event in the mid-19th century (27). L2b, on the other hand, appears to have spread to Europe via South America, potentially contributing to the prevalence of bvSP in the United States. Considering the economic losses caused by bvSP, the United States, Europe, and other industrialized countries implemented eradication programs in the mid-19th century, thus eliminating the risk of pullorum disease. Implementing similar measures is challenging due to China’s vast geographical area and various economic factors, resulting in bvSP becoming an endemic pathogen, mainly in China.

Additionally, forty-five bvSP was isolated from Taishun and Yueqing, to explore the transmission pattern of bvSP in China. The cgSNP distances between bvSP isolates from Taishun and Yueqing, showed significant diversity, with cgSNP distances of less than 100 isolate pairs observed exclusively within a single region. The transmission history of 45 global bvSP isolates revealed that most transmission events occurred within a single region. Additionally, the number of transmission events decreased rapidly with increasing distance, confirming that bvSP in China belongs to a highly localized pathogen.

Humanity appeared powerless against epidemics until the advent of antimicrobials, which are still the primary agents used to combat S. Gallinarum infection (28, 29). As a manifestation of localization-related forces, the emergence of AMR allows pathogens to acquire additional fitness benefits (30).

Our findings further indicate that the mobilome facilitates the acquisition of AMR among bvSP, mainly through plasmids and transposons, such as IncX1, IncQ1, and IncN type plasmids, as well as Tn801, Tn6205, and Tn721 type transposons. A survey conducted in China between 1962 and 2010 documented the association of class 1 integrons with AMR among bvSP (31). Our study further demonstrates that In498, In1440, In473, and In282, all belonging to class 1 integrons, play a significant role in carrying ARGs associated with sulfonamides, and trimethoprim. However, compared to other Salmonella serovars (32–34), the level of resistome in bvSP remains lower.

Considering lineage prospects, recently emerged lineages present a higher risk of AMR, with a lineage-specific distribution of resistome. In L2b and L3c, the most prevalent groups of ARGs include bla_TEM-1b and sul2. L3c, which exhibits the highest resistome richness with sul1, tet(A), aadA5, and dfrA17, has a higher likelihood of MDR. However, the transmission routes of resistome showed a high degree of consistency across continents. Globally, a small proportion of resistome carriers belonging to L2b and L3b were observed in the Americas and Europe, while isolates from Asia exhibited a notable rise in resistome carriage. Therefore, we speculate that the mobilome in bvSP may have been acquired locally in Asia, possibly due to increased antimicrobial usage and an intensive industrialized poultry farming pattern within the regional poultry industry (35). The L3c, being a local Asian lineage and more prone to acquiring additional resistome locally, provided further support.

In China, the endemic region of Asia, an unexpected rising trend has been observed from the 1970s to the 2020s. However, the risk of AMR varies among distinct regions (36–39). Overall, the average number of ARGs carried by bvSP isolated from eastern China was the lowest, while it was higher among bvSP isolated from northern and southern China. This variation may be influenced by differences in the economy, climate, or farming practices among different regions. Interestingly, the widespread horizontal transfer of the resistome across regions seems to indicate that the resistome has a different mode of transmission compared to its pathogenic bacteria. Specifically, as a localized pathogen, our results show a low frequency of cross-regional transmission of bvSP (as shown in S5 Fig), but up to 90% of ARGs are transmitted cross-regionally. This risk has not been previously observed in S. Gallinarum, suggesting that potential intermediate hosts may play a key role.

In summary, the findings of this study highlight that S. Gallinarum remains a significant concern in developing countries, particularly in China. Compared to other regions, S. Gallinarum in China poses a notably higher risk of AMR, necessitating the development of additional therapies, i.e. vaccine, probiotics, bacteriophage therapy in response to the government’s policy aimed at reducing antimicrobial use (40). Furthermore, given the dynamic nature of S. Gallinarum risks across different regions, it is crucial to prioritize continuous monitoring in key areas, particularly in China’s southern regions where the extensive poultry farming is located. Lastly, from a One-Health perspective, controlling AMR in S. Gallinarum should not solely focus on local farming environments, with improved overall welfare on poultry and farming style. The breeding pyramid of industrialized poultry production should be targeted on the top, with enhanced and accurate detection techniques (41). More importantly, comprehensive efforts should be made to reduce antimicrobial usage overall and mitigate potential AMR transmission from environmental sources or other hosts (42–44).

However, the current study has some limitations. Firstly, despite assembling the most comprehensive WGS database for S. Gallinarum from public and laboratory sources, there are still biases in the examined collection. The majority (438/580) of S. Gallinarum samples were collected from China, possibly since the WGS is a technology that only became widely available in the 21st century. This makes it impractical to sequence it on a large scale in the 20th century, when S. Gallinarum caused a global pandemic. So, we suspect that human intervention in the development of this epidemic is the main driving force behind the fact that most of the strains in the data set originated in China. In our future work, we aim to actively gather more data to minimize potential biases within our dataset, thereby improving the robustness and generalizability of our findings. Secondly, in silico analysis relies heavily on sophisticated, continuously updated databases. Although we utilized the most up-to-date databases, the exact prevalence of mobilome genetic elements (MGEs) and ARGs may be underestimated. Nevertheless, by using the most extensive global WGS data from S. Gallinarum, we elucidate a dynamic resistome evolution framework in a single bacterial pathogen from pandemic to endemic, which is hallmarked by a gradual increase in antimicrobial resistance risks and a stepwise resistome adaptation. The insights gained in this study will be helpful in further prevention and surveillance of the AMR risks of localized pathogens.

Material and Methods

Ethics statements

The protocols used for animal studies were approved by the Committee of the Laboratory Animal Center of Zhejiang University (ZJU20190094; ZJU20220295).

Bacterial isolates

All 734 samples of dead chicken embryos aged 19 to 20 days were collected from Taishun and Yueqing in Zhejiang Province, China. After a thorough autopsy, the liver, intestines, and spleen were extracted and added separately into 2 mL centrifuge tubes containing 1 mL PBS. The organs were then homogenized by grinding. In the initial enrichment phase, we utilized Buffered Peptone Water (BPW, Haibo Biotechnology Co, Qingdao, China), employing a 1:9 dilution ratio (sample in PBS: BPW). Subsequently, the composite was incubated at 37°C for 16-18 hours in a rotary incubator set to 180 rpm. For selective enrichment, Tetrathionate Broth Base (TTB, Land Bridge Biotechnology Co, Beijing, China), fortified with iodine solution and brilliant green solution (both from Land Bridge Biotechnology Co, Beijing, China), was employed at a ratio of 1:10 (sample in BPW: TTB). The mixture was subjected to incubation at 42°C for a duration ranging between 22 and 26 hours within a rotary incubator set at 180 rpm. Isolated Salmonella colonies from positive samples were obtained by sub-culturing selectively enriched samples on Xylose Lysine Deoxycholate (XLD, Land Bridge Technology Co, Beijing, China) agar, followed by an 18–22-hour incubation at 37°C. Typical and pure colonies were selected after sub-culturing on XLD agar and transferred into Luria-Bertani (LB) broth. Finally, the transferred bacterial culture underwent an additional incubation for 18-22 hours at 37°C in a rotary incubator operating at 180 rpm.

DNA extraction and genomic assembly

The DNA of forty-five isolates was extracted using the Vazyme Fastpure^® Bacteria DNA Isolation Mini Kit (Vazyme Biotech Co., Ltd.), which was then quantified using the NanoDrop1000 system (Thermo Fisher Scientific, USA). Subsequently, DNA libraries were constructed and subjected to sequencing using the Illumina Novaseq 6000 platform (Beijing Novogene Co., Ltd). Assembly of genome sequences was performed by SPAdes (45) v3.12.0 with default parameters.

Global dataset assembly

To provide global context, we expanded the dataset by including 540 additional WGS genome data, consisting of bvSP (n=483), bvSG (n=50), bvSD (n=2), and Salmonella enterica serovar Enteritidis (n=5). The inclusion of five strains of Salmonella enterica serovar Enteritidis enhanced evolutionary analyses. Notably, in this collection of 540 genomic data from public databases, 325 sequences were previously preserved in our laboratory.

All WGS data passed strict quality control according to criteria set by the European Reference Laboratory (46). Genomic data exceeding 500 contigs and having an N50 of less than 30,000 were excluded. Lastly, the bacterial species associated with each genomic data set was confirmed by utilizing KmerFinder (47) v3.2. Finally, the assembled database includes a total of 585 high-quality sequences.

Genotyping and phylogenetic analyses

Snippy v.4.4.5 was used to identify cgSNPs, with the complete genome of R51 used as a reference. The optimal model (TVM+F) was determined and employed to construct the maximum-likelihood phylogenetic tree using IQ-TREE v.1.6 (48). The final phylogenetic tree was annotated and visualized using iTOL v.6.0 (49).

The global assembly of S. Gallinarum underwent a population structure analysis using the RhierBAPS v1.1.4 (50), based on the Hierarchical Bayesian clustering (BHC) algorithm and employing 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. Default parameters were employed for all population structure analysis.

Temporal and phylogeographical analysis

To examine the emergence and geographical transfers, we first utilized Gubbins v.2.3 (51) to eliminate the recombination regions for each lineags (S12 Fig). Then, TreeTime (52) was used to assess the temporal structure. This was accomplished through regression analysis of the root-to-tip branch distances within the maximum likelihood tree, considering the sampling date as a variable. Then, the strains within the dataset underwent phylogeographical reconstruction via Bayesian Evolutionary Analysis of Sampled Trees (BEAST) (53) v.2.5. To determine the optimal model for running BEAST, 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 Markov Chain Monte Carlo (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 proceeded to perform 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, the “Coalescent Bayesian Skyline” tree prior. Convergence was assessed using Tracer, with all parameter effective sampling sizes (ESS) exceeding 200. Maximum clade credibility trees were generated using TreeAnnotator v.2.6.7. Finally, key divergence time points (with 95% credible intervals) were estimated and the phylogenetic tree was visualized using Figtree v.1.4.3.

SNP distance-based S. Gallinarum geographical tracing

To understand potential transmission events, we calculated the cgSNPs distance between distinct bacterial strains employing the SNP-dists. We estimated the overall evolutionary rate of the S. Gallinarum using BEAST. We applied the methodology described previously (54). The numbers of SNPs per year were determined by multiplying the evolutionary rates estimated with BEAST by the number of core SNP sites identified in the alignments. We hypothesize that a slower evolutionary rate in bacteria typically requires a lower SNP threshold when tracing transmission events using SNP distance analysis. Previous research found an average evolutionary rate of 1.97 SNPs per year (95% HPD, 0.48 to 4.61) across 22 different Salmonella serotypes. Our updated BEAST estimation for the evolutionary rate of S. Gallinarum suggests it is approximately 0.74 SNPs per year (95% HPD, 0.42 to 1.06). Based on these findings, as well as our previous experience with similar studies (55), we set a threshold of 5 SNPs. Additionally, the geographic location was considered to enhance the precision of speculation regarding potential transmission sites and transmission direction.

MGE detecting

Four types of MGEs were detected: plasmids, transposons, integrons, and prophages. Plasmids were detected using Abricate v.1.0.1 with the PlasmidFinder (56) database. Only plasmids with a similarity of more than 90% and a coverage over 95% were identified. BacAnt (57) 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 Phaster pipeline. The genomic data were split into two temporary databases based on the number of contigs, one dataset containing single contig files and the other containing multiple contig files. The two databases were imported into the Phaster pipeline separately with default parameters.

Determination of the potential MGE carried ARG and horizontal ARG transfer frequency

We have devised a pipeline to discern the coexistence of MGEs and ARGs. For plasmids and prophages, only ARGs situated within MGE regions are deemed MGE-carried. Regarding transposons and integrons, the search region was expanded to 5 kilobases upstream and downstream of ARGs, accounting for potential splicing errors. In this study, HGT events involving ARGs were defined as instances where ARGs exhibited perfect identity (100% coverage and 100% identity) and were located on MGEs across different strains. In a previous study (58), horizontal transfer of ARGs was identified based on gene regions with 100% coverage and over 99% identity. We believe that these stricter criteria in our pipeline more accurately reflect the transfer of ARGs within the S. Gallinarum population. HGT frequency was utilized to evaluate the extent of horizontal transfer of ARGs as follows:

The pipeline developed in this study accepts Resfinder or Resistance Gene Identifier (RGI) results as input. Specific code and examples are uploaded to: https://github.com/tjiaa/Cal_HGT_Frequency/tree/main.

Genomic data analysis

The serovars of all WGS genome data were confirmed through SISTR (59) v.1.1 and Seqsero2 (60). Multilocus Sequence Typing was carried out using MLST v.2.22 with the senterica_achtman_2 scheme. ARGs were detected by ResFinder (61). Both similarity and coverage were set to a minimum value of 90. The invasiveness index was calculated using methods previously reported (62). Specifically, Salmonella’s ability to cause intestinal or extraintestinal infections in hosts is related to the degree of genome degradation. We evaluated the potential for extraintestinal infection by 45 newly isolated S. Gallinarum strains (L2b and L3b) using a model that quantitatively assesses genome degradation. We analyzed each sample using the 196 top predictor genes for measuring the invasiveness of S. Gallinarum, employing a machine-learning approach that utilizes a random forest classifier and delta-bitscore functional variant-calling. This method evaluated the invasiveness of S. Gallinarum towards the host, and the distribution of invasiveness index values for each region was statistically tested using unpaired t-test. The code used for calculating the invasiveness index is available at https://github.com/Gardner-BinfLab/invasive_salmonella.

Additional information

Data availability

For the newly isolated 45 strains of Salmonella Gallinarum, genome data have been deposited in NCBI Sequence Read Archive (SRA) database. The “SRA Accession” for each strain are listed in Supplementary Table 1. Additionally, the genome data for the 540 publicly available genomes have been uploaded to Figureshare and can be accessed at: https://figshare.com/articles/journal_contribution/Genomic_data/26028592.

Code availability

The open-source software used in this study includes:

BacAnt (https://github.com/xthua/bacant),

Abricate (https://github.com/tseemann/abricate),

Snippy (https://github.com/tseemann/snippy),

SISTR (https://github.com/phac-nml/sistr_cmd),

MLST (https://cge.food.dtu.dk/services/MLST/)

SeqSero2 (https://github.com/denglab/SeqSero2),

IQtree (https://github.com/iqtree/iqtree2).

BEAST (https://www.beast2.org)

SNP-dists (https://github.com/tseemann/snp-dists)

Fast BAPS (https://github.com/gtonkinhill/fastbaps)

R (https://www.r-project.org)

Phaster (http://phaster.ca)

KmerFinder (https://bitbucket.org/genomicepidemiology/kmerfinder)

Treetime (https://github.com/neherlab/treetime)

Author contributions

M.Y., C.J. designed the study; finalized the manuscript;

C.J. collected isolates and clinical data, integrated the data, performed metadata and genomic analysis, prepared tables and figures; and wrote the draft manuscript;

C.H. performed metadata and genomic analysis; interpreted the data; finalized the manuscript;

X.Z., X.K., Q.C., Y.H. and Z.W. contributed to the data analysis;

C.H., M.Y. prepared tables and figures;

M.Y., C.H., A.S., Y.L., F.H., and C.J. finalized the manuscript;

M.Y. managed the project;

All the authors contributed to the article and approved the submitted version.

Declaration of Interests

None.

Acknowledgements

This work was supported by the National Program on the Key Research Project of China (2022YFC2604201), the Zhejiang Provincial Natural Science Foundation of China (LZ24C180002; LR19C180001), the Hainan Provincial Joint Project of Sanya Yazhou Bay Science and Technology City (2021JJLH0083), the Zhejiang Provincial Key R&D Program of China (2023C03045, 2022C02024) and the Open Project Program of the Jiangsu Key Laboratory of Zoonosis (R1902). We also thank Dr. Peide Li for providing the valuable resource for this investigation.