Introduction

The acquisition of antimicrobial resistance (AMR) by pathogens is well-established as one of the most serious 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 targeted human interventions, climate change, and geographical selection forces (9, 10). Understanding the evolutionary trajectory of the AMR and related reservoirs is necessary, as disease management strategies will differ substantially (11). 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 (12). 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 strong geographical-related characteristics, has more than 90% serovars, frequently categorized as geo-serotypes (13). Among a few thousand geo-serotypes, Salmonella enterica serovar Gallinarum (S. Gallinarum) is an avian-specific agent, leading to horrific mortality rates that remain a disproportionate burden to the poultry industry in low- and middle-income countries (14, 15). 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 (16, 17).

Nowadays, antimicrobial therapy is still used against S. Gallinarum infection (18). 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 (19, 20). 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 of 580 S. Gallinarum isolates, covering the period from 1920 to 2023. Using WGS dataset we investigated: 1) the population structure and 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 relationship of S. Gallinarum, we assembled the most comprehensive S. Gallinarum WGS dataset (n=580). The core genome SNPs (cgSNPs) were obtained by using the fully sequenced genome of S. Gallinarum R51 as a reference. Recombination regions in the genome were eliminated, and through hierarchical Bayesian analysis of cgSNPs, it was confirmed that S. Gallinarum divides into three biovars: bvSP (n=528/580, 91.03%), bvSG (n=50/580, 8.6%), and 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 (141/528), and L3c (169/528). Importantly, L3a, previously considered distinct, has been revealed as a sub-lineage of L3b. 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, L3c and L2a are more widespread in Europe and America, respectively (Fig 1a). Notably, all bvSP isolates from Asia were exclusively found in China, which can be manually divided into three distinct regions (southern, eastern, and northern). 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 L3b emerged 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 L3c. Then, from the 2000s to the 2010s, L3c and L1 were continually replaced by L3b, 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 L3c to L2b and L1 to L3b 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 L3c were selected as the dominant global lineage types, and a Bayesian model with a relaxed molecular clock was used to infer the evolutionary dynamics and historical transmission routes using the BEAST software (Fig 2). The temporal structure was verified (S3 Fig), and our findings show that L2b emerged in China earlier than L3c, with potential spread from China to Europe as early as around the year 1800. Then, a long period of L2b spread internationally between South America and Europe was observed between 1800 and 1980. Notably, the L2b exhibited a return transmission from Europe to China around 1980, coinciding with the normalization of China-EU relations, as well as the timing of increased economic cooperation between China and Europe. In contrast, the historical transmission route for the L3c branch is more distinct, with L3c introduced to Europe from China around 1890 and reintroduced to China around 1910 (S4 Fig).

Figure 2.

Endemic isolates further support localized dissemination pattern of bvSP

To investigate the dissemination pattern of bvSP in China, 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). Interestingly, a lineage-preferential distribution was also observed, with bvSP isolated from Taishun belonging to L3c 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 100 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).

Furthermore, we simulated potential transmission events between the bvSP strains isolated from Zhejiang Province and those from other provinces of China. As a result, the transmission events showed a strong geographic-preferential distribution. The number of transmission events decreased quickly with increasing distance, which provides additional support for bvSP in China belonging to a highly localized pathogen (S5 Fig).

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

The mobilome drives partitions of the resistome and plays a pivotal role in shaping the ecological niches 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 distribution. We observed that L3b is more inclined to carry bla_TEM-1B and sul2 than other lineages, whereas tet(A) exclusively exists within L3. 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 (S6a Fig). 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 L3b. 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%) (S6b Fig). Notably, ColpVC had the lowest carriage percentage among L1 (2/10, 20%), whereas IncX1 was explicitly associated with L3b (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, L3b 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 L3c. 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 a greater variety of mobilome types (S7 Fig).

Plasmids and transposons guide resistome geo-temporal dissemination

Variations in regional antimicrobial use result in an uneven distribution of AMR risks. 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 (S8 Fig). In more detail, both the resistome and mobilome exhibited a steady decline until the 1980s, followed by a consistent increase from the 1980s to the 2010s. However, after the 2010s, a subsequent decrease was identified. Consequently, we further investigated the interplay between mobilome and resistome in bvSP.

In bvSP, the predominant resistome is carried by the 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 serve as 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. Subsequently, we further investigated the HGT frequency of each ARG carried by bvSP isolated from China and explored the HGT frequency of resistome between various regions. Horizontal transfer of ARGs is defined as sequences with perfect identity (100% identity and 100% coverage) (Fig. 6a). Our findings indicate that horizontal transfer of the ARGs is widespread among Chinese isolates of bvSP. Specifically, 60% of the ARGs exhibited an HGT frequency of 100%, indicating that these ARGs had undergone horizontal transfer at least once (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 non-uniform resistome HGT frequency. Overall, bvSP from the southern region of China exhibited the highest HGT frequency (HGT frequency=97%). The HGT frequencies for bvSP within the eastern and northern regions of China are lower, at 91% and 90%, respectively (Fig. 6c). For specifically 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 bla_TEM-116 only undergo horizontal transfer within a single region of eastern China (Fig. 6d). Notably, as a localized transmission pathogen, resistome carried by bvSP exhibited a strong risk of cross-regional transmission, especially from northern region to southern region (HGT frequency=95%) (Fig. 6e).

Discussion

Over the past decades, the emergence of AMR has garnered significant attention from both public health agencies and the poultry industry (21). 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 (22). 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 L3c of bvSP being predominant global lineage types. The historical geographical transmission was verified using a spatiotemporal Bayesian framework. The result shows that both lineages were initially spread from China to Europe in the early 19th century, which correlates strongly with the European hen fever event in the mid-19th century (23). Subsequently, the practice of breeding Asian chickens gradually spread from Europe to South America, leading to the bvSP prevalence in the United States. Considering the economic losses caused by bvSP, the United States, Europe, and other major countries implemented eradication programs in the mid-19th century, thus eliminating the risk of pullorum disease. However, with the strengthening of China-EU relations in the 1980s, a transmission of L2b from Europe to China was observed. 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 (24, 25). As a manifestation of localization-related forces, the emergence of AMR allows pathogens to acquire additional fitness benefits (26). 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 (27). 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 (28–30), 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 L3b, 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 L3c 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 (31). The L3b, 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. 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 host. 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.

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. Secondly, in silico analysis relies heavily on sophisticated, continuously updated databases. Although we utilized the most up-to-date databases, the exact prevalence of 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 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 (Novogene, Beijing, China). Assembly of genome sequences was performed by SPAdes (32) 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 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 (33). Genomic data exceeding 500 contigs and having an N50 of less than 30,000 were excluded. Moreover, we considered GC content within the range of 51% to 53%. Lastly, the bacterial species associated with each genomic data set was confirmed by utilizing KmerFinder (34) 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, and Gubbins v.2.3 (35) was used to eliminate the recombination regions, 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 (36). The final phylogenetic tree was annotated and visualized using iTOL v.6.0 (37).

The global assembly of S. Gallinarum underwent a population structure analysis using the RhierBAPS v1.1.4 (38), 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 TreeTime (39) 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) (40) v.2.5. The BEAST analysis was performed on alignments of concatenated cgSNPs following previously described parameters (41). Maximum clade credibility trees were generated using TreeAnnotator v.2.6.7, and the tree was annotated and visualized using Figtree v.1.4.3.

To understand potential transmission events, we calculated the cgSNPs distance between distinct bacterial strains employing the SNP-dists, with cgSNPs distance less than 100 deemed indicatives of high similarity between two strains. Additionally, the geographic location and isolation date were considered to enhance the precision of speculation regarding potential transmission sites and the direction of transmission.

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 (42) database. Only plasmids with a similarity of more than 90% and a coverage over 95% were identified. BacAnt (43) 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). In a previous study (44), 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 (45) v.1.1 and Seqsero2 (46). Multilocus Sequence Typing was carried out using MLST v.2.22 with the senterica_achtman_2 scheme. ARGs were detected by ResFinder (47). Both similarity and coverage were set to a minimum value of 90. The invasiveness index was calculated using methods previously reported (48). Briefly, we utilized the pre-trained random forest model along with R v4.1.2 to compute the invasiveness index. This approach demonstrated the feasibility of calculating the invasiveness index across various Salmonella species. Specifically, we adopted the DeltaBS metric to detect mutations in protein-coding genes within the Whole Genome Sequence of the strains used for model training, resulting in the identification of 196 top predicted genes associated with invasiveness in S. enterica. The code used to calculate the invasiveness index is accessible at https://github.com/Gardner-BinfLab/invasive_salmonella.

Data Availability

The genome data of the newly isolated forty-five S. Gallinarum strains have been deposited into the China National GenBank Database (CNGBdb) under accession number CNP0005633. 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)

Additional information

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

The authors have declared no competing interests.

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), and the European Union’s Horizon 2020 Research and Innovation Programme under Grant Agreement No. 861917 – SAFFI.