The complete genomes of 345 Pasteurellaceae were retrieved from the National Center for Biotechnology Information (NCBI) non-redundant RefSeq database (last accessed in May 2020)¹ (Supplementary Table 1). Phylogenetic relationships between the genomes were determined using alignment of 16S rRNA sequences using ClustalW (Juraschek et al., 2019), and a maximum likelihood phylogenetic tree was constructed in MEGA X (Kumar et al., 2018), using the Generalized Time Reversible (GTR) and bootstrap confidence value of 1000. The phylogenetic tree was visualized with Interactive Tree of Life (iTOL) (Letunic and Bork, 2019), where the tree was edited and supplemented with genome information. The global distribution of genomes was determined using the ggmap R package version 3.0.0, a heatmap being plotted based on the number of genomes available for each region.

GenBank format (.gbk) of the genome dataset was used as input for insertion sequence (IS) prediction using ISsaga (Varani et al., 2011) with default parameters. In addition, the recommendations of the Everyman’s Guide to Bacterial Insertion Sequences (Siguier et al., 2015) were used to identify partial elements and provide IS family features.

Hierarchic organization of IS family distribution was visualized in the R environment using the ggplot2 circular package version 3.3.5. ISsaga was used to assess IS ORF genome context. Transposase and adjacent gene sequence were extracted and analyzed for conserved domains using CD-Search (Marchler-Bauer and Bryant, 2004) against the CDD v3.18 database with an expected value threshold of 0.01, and its gene product was inspected for gene ontology annotation through the QuickGO resource at EMBL-EBI (Binns et al., 2009). Next, we created a local database of genes flanking ISs (Supplementary Table 3) and divided them into four classes according to their function: stress response, AMR, adaptation, and virulence. Circular visualization of these classes was created in Circos Table Viewer (Krzywinski et al., 2009).

The genome sequences (.gbk format) of Mannheimia haemolytica (strains M42548, USDA-ARS-USMARC-184, and NCTC10643) and Aggregatibacter actinomycetemcomitans (strains HK_907, KK1651, and VT1169) were selected as examples of genomes carrying many or few IS copies, and a multiple genome alignment was performed using progressive Mauve (Darling et al., 2004).

We analyzed the prophage elements integrated into the genomes of Pasteurellaceae in three critical aspects: prophage-like elements, complete prophages, and novel putative prophage elements. First, we looked for prophage-like elements, i.e., candidate intact prophage that contained phage attachment sites, genes encoding structural phage proteins, genes coding for proteins involved in DNA regulation, insertion into the host genome, and lysis (Arndt et al., 2016; Czajkowski, 2019). We used PHASTER (Arndt et al., 2016) and Prophage Hunter (Song et al., 2019) to predict prophage-like elements. Raw data from the predictions were used for further analysis (Supplementary Tables 4, 5). Next, we screened for complete prophages and prophage-like elements in the classes of intact and active prophage-like elements from PHASTER and Prophage Hunter, respectively. Sequences were subject to BLASTN searches against reference viral genomes already described in the family (Supplementary Table 7) using MegaBLAST (Chen et al., 2015) alignment with a cut-off of 75% of cover and 85% identity, allowing identification of complete prophages found in the bacterial genomes (Supplementary Table 6). More than one complete prophage insertion into the genome was considered a poly-lysogenic event. Finally, we identified novel putative prophages considering the classes of incomplete, questionable, and ambiguous. We also manually inspected the results of mismatched prophage-like elements, according to our criteria, from previous analyses of complete prophages. We screened these sequences for phage-associated functions. Upper boundaries of the novel prophages were determined, wherever possible, by searching for phage integrases from the tyrosine recombinase family at the tRNAs. A sequence identity matrix was built using whole nucleotide sequences of putative novel prophages. These sequences were aligned by Clustal Omega (Sievers and Higgins, 2014) with the default parameters.

To show that putative novel prophages were different from those previously reported in the family (Supplementary Table 7), a synteny analysis using clinker and clustermap.js (Gilchrist and Chooi, 2021) among phages classified by genus was performed. Additionally, we classified a prophage to a family belonging to the order Caudovirales using the occurrence of head-neck-tail module genes detected by the Virfam (Lopes et al., 2014). Phylogenomic analysis of the novel prophages was performed by the ViPTree (Nishimura et al., 2017) webserver. A MegaBLAST analysis against the novel prophages to evaluate their dispersion among Pasteurellaceae genomes, using a cutoff of 90% of cover and identity, was carried out. A bipartite network was constructed using vConTACT (Bolduc et al., 2017) to estimate the relationship and clustering of prophages belonging to the Pasteurellaceae family.

A graphic representation of novel and previously reported prophages mapped onto the genomes of the family (grouped by genus), based on size, GC content, and number of ORFs was carried out. The annotation of protein sequences was done using GeneMarkS (Besemer et al., 2001) version 4.28 with the sequence type of the phage marked. BLASTP (McGinnis and Madden, 2004) was used to build a local database of the protein sequences, which provided a homologous protein cluster (HPC) with sequences > 30% amino acid identity, > 80% alignment coverage, and clustering E-value < 1E-5. Next, the functional annotation of HPC was done using HMMER69 v3.b2 searches with default parameters to the PFAM (Finn et al., 2014). We also analyzed the lifestyle of prophages using the phage Classification Tool Set (PHACTS) (McNair et al., 2012).

To identify integrative and conjugative elements (ICEs) in the genomes of the Pasteurellaceae, the genomes were inspected for MGE-encoding relaxases, type-IV coupling proteins (T4CP), and the type-IV secretion system (T4SS) gene cluster using the oriTfinder (Li X. et al., 2018) (Supplementary Table 9). ICE sequences experimentally validated for the family were retrieved from the ICEberg (Liu et al., 2018) database version 2.0. Comparative analyses with the ICEberg data allowed the identification of novel putative ICE elements within the Pasteurellaceae family (Supplementary Table 11). An element was considered as conjugative when it contained a relaxase, a T4CP, a T4SS gene cluster, and type-specific genes related to mating pair formation (Cury et al., 2017). We inspected attachment sites located between the tRNA and the integrase gene using Repeat Finder plugins on Geneious Prime^® version 2020. The opposite boundary of the element was delimited by BLASTN searching, and alignment of the ICE and attachment site sequences using minimum general parameters and filter low complexity region marking. The integrase family was classified by conserved domain searches using CD-Search (Marchler-Bauer and Bryant, 2004) tRNAscan-SE (Chan and Lowe, 2019) and MOBscan (Garcillán-Barcia et al., 2020) were used to identify tRNA genes and classify relaxase families, respectively.

We built a sequence identity matrix using whole nucleotide sequences of putative novel ICEs. These sequences were aligned by Clustal Omega (Sievers and Higgins, 2014) with the default parameters, and a heatmap was generated using the ggplot2 R package. To show that putative novel ICEs were different from those previously reported in the family (Supplementary Table 7), a synteny analysis using clinker and clustermap.js (Gilchrist and Chooi, 2021) among groups of ICEs classified by genus was performed. Next, we performed a global analysis of ICEs mapped onto the genomes, based on size, GC content, and number of ORFs. GeneMarkS (Besemer et al., 2001) was used to annotate predicted protein sequences. A BLASTP alignment to obtain HPC with sequences > 30% protein identity, > 80% alignment coverage, and clustering E-value < 1E-5 was done. Next, the functional annotation of HPC was done using HMMER (Wheeler and Eddy, 2013) v3.b2 searches with default parameters to the Pfam (Finn et al., 2014).

Key ICE genes were used for tracking the evolutionary history of these elements. The integrase (int), topoisomerase (parA), and coupling protein (traD) genes were found in almost all ICE sequences. Nucleotide sequences of these genes were aligned using ClustalW to construct maximum likelihood phylogenetic trees. The GTR model and a bootstrap confidence value of 1,000 were applied to each tree. The alignment and phylogenetic analysis were done using MEGA X (Kumar et al., 2018) and edited using iTOL (Letunic and Bork, 2019).

Finally, we constructed a network to evaluate the dispersion of the ICEs among the Pasteurellaceae host genomes. A local database with nucleotide sequences of the ICEs was created and MegaBLAST alignments were done, with a cut-off of 75% of cover and identity to consider the interaction between individual ICE and the host genomes, and a tabular representation of this interaction was constructed. The network was visualized using Cytoscape (Shannon et al., 2003).

A total of 162 plasmid sequences from the Pasteurellaceae were retrieved from the NCBI RefSeq database (last accessed on Nov 2020)² (Supplementary Table 12). We manually curated the plasmid database, eliminating partial plasmid DNA sequences, redundant nomenclature, and unassignable hosts. Relaxases were classified using MOBscan (Garcillán-Barcia et al., 2020), and the nucleotide sequence of mob genes was aligned using ClustalW to construct a maximum likelihood phylogenetic tree. The Tamura-Nei (TN) model and a bootstrap confidence value of 1,000 were applied. Alignment and phylogenetic analyses were done using MEGA X and edited using iTOL. To analyze potentially non-mobilizable plasmids, amino acid sequences were aligned using Muscle (Madeira et al., 2019) with the output format adjusted to Phylip sequential. A distant matrix was constructed using the output file through EMBOSS (Madeira et al., 2019). A multi-dimensional graph was created in R. Lastly, to correlate the redundant nomenclature plasmids mentioned above with the One-health concept, an arc diagram was created using the common plasmids found in different species/genomes.

AMR genes were identified by ResFinder 4.0 (Bortolaia et al., 2020) and the Comprehensive Antibiotic Resistance Database (CARD 2020) (Alcock et al., 2020). The latter was also used to identify any synonymous AMR genes indicated in different published elements using different names. Firstly, we used the context genes associated with AMR provided by ISsaga analysis (Supplementary Table 3) to investigate the association of transposable elements (TEs) with AMR genes. Subsequently, the impact of IS elements in three classes (upstream, downstream, and interruption of AMR gene) according to the transposase position in relation to the AMR genes, was classified. Sequences were extracted from transposons associated with ICEs to analyze the possible transfer of AMR gene sequences. Individual transposases were aligned using the ISfinder database (Siguier et al., 2006), using default parameters (E-value ≤ 10-5), and a minimum alignment coverage of 50% and with at least 70% identity was considered. The direct repeat and terminal inverted repeats were manually identified and annotated using Geneious Prime^® based BLASTn searches against ISfinder to identify known IS elements.

Fasta sequences of all genomes and MGEs of the Pasteurellaceae were used as input for the AMR gene predictions. To show the contribution of the MGEs in carrying AMR genes, we compared the content of AMR genes of these elements to the total number of AMR genes in their respective genomes. For comparison, we used localization of the ICE in the genomes and verified whether the genes were located within the element, and also looked for AMR genes in plasmids belonging to isolates from our genome dataset. Also, individual plasmid and prophage sequences were inspected for AMR genes. The distribution of AMR genes and classes among ICEs and plasmids was represented using the Sankey diagram.³ Finally, sequences of the most prevalent AMR genes detected in MGEs were evaluated for selective pressure. AMR gene sequences were aligned in MEGAX (Kumar et al., 2018) and exported in meg format. The program DnaSP v6 (Rozas et al., 2017) was used to carry out Tajima (1989) and Fu and Li (1993) tests.

We surveyed 345 publicly available Pasteurellaceae complete genomes, spanning 34 species belonging to 14 genera (Figure 1A). Our genome dataset contains representatives from diverse veterinary species (212), humans (119), and unknown sources (14) (Figure 1A, Supplementary Figure 1, and Supplementary Table 1). Genome sizes ranged from 1.5 Mb in Haemophilus spp. to 2.8 Mb in Mannheimia spp. (Figure 1A and Supplementary Table 1). We inferred an evolutionary tree using 16S rRNA (Figure 1A). Similar to previous studies (Naushad et al., 2015), our phylogeny indicates that many Pasteurellaceae genera have different paraphyletic clades and tend to form different separated clusters, which may reflect misclassifications awaiting resolution. Four different species were assigned to the genus Mannheimia (i.e., Mannheimia granulomatis, M. haemolytica, Mannheimia varigena, and Mannheimia succiniciproducens), and eight different species have names indicating they are members of the genus Actinobacillus [i.e., (Actinobacillus) delphinicola, Actinobacillus equuli, (Actinobacillus) indolicus, Actinobacillus lignieresii, Actinobacillus pleuropneumoniae, (Actinobacillus) porcitonsillarum, (Actinobacillus) succinogenes, and Actinobacillus suis], though square brackets around the genus names indicate species that require reclassification, as they are not Actinobacillus sensu stricto (Blackall and Turni, 2020). Our dataset shows a global distribution of sources, though America, Europe and Asia are more highly represented (Figure 1B). H. influenzae, Pasteurella multocida, A. pleuropneumoniae, and A. actinomycetemcomitans had the most genomes available per species, presumably reflecting their importance for human and veterinary health.

We used ISsaga to predict, map and annotate a total of 10,820 ISs, belonging to 19 different families and divided into 12 subgroups (ssgr), in our genome dataset (Supplementary Table 2 and Supplementary Figure 2A). Of the known elements, 37.8% (4,016) were intact and 62.2% (6,611) were partial. The most frequent IS families found were: IS481, ISL3 and IS1595 ssgr IS1016 with 2,802 (25.9%), 1,588 (14.67%), 1,211 (11.2%) occurrences, respectively. These elements were broadly disseminated among the Pasteurellaceae genomes, with a higher diversity of IS families and an average number of 53, 23, and 94 ISs per genome in Mannheimia, Haemophilus, and Glaesserella genera, respectively (Supplementary Table 2 and Supplementary Figure 2B). Despite the majority of ISs belonging to IS481, this family was only found in eight genera, whereas ISL3, IS1595, and IS3 were the most widespread within the Pasteurellaceae (Figure 2A). Mapping of ISs with regards to genomic locations (Supplementary Table 3) identified 14 IS families inserted within/adjacent to 337 genes associated with stress responses (n = 145), AMR (n = 89), adaptation (n = 42), and virulence (n = 61), with IS481 and ISL3 being the most represented families (Figure 2C).

The abundance of IS elements was positively correlated with genome size for Mannheimia spp., Haemophilus spp., and Histophilus spp., but not for Glaesserella spp., Actinobacillus spp., Aggregatibacter spp., and Pasteurella spp. (Figure 2B). Alignment of M. haemolytica genomes (species with many ISs) compared to alignment of A. actinomycetemcomitans genomes (species with few ISs) revealed, based on synteny analysis, numerous internal rearrangements in genomes possessing a higher number of IS copies compared with those with a lower number of ISs (Supplementary Figure 3), which might have an impact on the genetic organization of these species.

Using PHASTER and Prophage Hunter, we mapped 2,939 prophage-like elements within our genome dataset (Supplementary Table 4, 5). Of these, 1,398 were classified as intact, questionable, or incomplete by PHASTER, and 1,541 were classified as active or ambiguous using Prophage Hunter (Supplementary Figure 4A). For prophage-like elements classified as intact and active from PHASTER and Prophage Hunter, respectively, we delimited complete prophages comprising eleven different viral species, identified in 193 insertions in 91 Pasteurellaceae genomes (Supplementary Table 6 and Supplementary Figure 4B).

Manual inspection of sequences identified by PHASTER and Prophage Hunter as incomplete/questionable/ambiguous, as well as complete prophages not matching any reported phage for the family, further identified possible novel prophage species, beyond those previously reported in the Pasteurellaceae (Supplementary Table 7). This approach identified 23 putative novel prophages in ten different genera (Supplementary Table 8), which exhibit the potential to encode proteins (Open Reading Frames- ORFs) necessary for their assembly and replication (Figure 3A). Whole-sequence alignment and synteny analysis using known Pasteurellaceae phages demonstrated that the predicted novel phages are new for this family (Supplementary Figures 4C, 5). These novel phages, predicted to belong to the Myoviridae (69%) and Syphoviridae (31%) families (Supplementary Table 8), are disseminated in different genomes in our dataset (Supplementary Figure 4D), with the most common being those designated here as Haemophilus phage GHA9 and Mannheimia phage 38599, and several representing the first report of phage in some of the Pasteurellaceae species (Supplementary Figure 4E). In some M. haemolytica and A. actinomycetemcomitans genomes, poly- insertions were noted (Supplementary Figure 6).

Combining our predicted novel prophage sequences with those previously reported for the Pasteurellaceae, overview analysis of the complete dataset revealed that the average size ranged from 31 kb (in Bibersteinia spp.) to 47 kb (in Haemophilus spp.), with between 30 and 136 ORFs and GC contents of 38–44% (Figure 3B). The average GC contents were similar to those of the genomes (40% ± 2) and the encoded proteins were predicted to be primarily involved in structure/assembly, transcriptional regulation, and recombination/replication, as well as some unknown function (Figure 3B). All prophages belonged to either the Myoviridae (73.8%) or Syphoviridae (26.2%) family (Supplementary Figure 4F). By comparison to closely related prophages, 83.3% are predicted to be temperate and 16.7% lytic (Figure 3C). Phylogenomic analysis revealed two distant clades, both containing monophyletic groups. None of the putative novel prophages grouped with previously reported Pasteurellaceae phages (cluster with 100% identity), confirming the novelty of our findings (Figure 3D). Results of VConTACT analysis further support the viral relationships among the published Pasteurellaceae phages, indicating four different clusters (according to the viral family), with homogeneous groups, found for Aggregatibacter and Mannheimia, and heterogeneous groups for Haemophilus, Mannheimia, and Pasteurella (Supplementary Figure 6). These results agree with the phylogenomic analysis, suggesting a strong evolutionary relationship amongst Pasteurellaceae phages.

OriTFinder analysis identified possible ICEs in our dataset. Results indicated several regions with high potential for self-transferability in 126 genomes comprising ten genera (Supplementary Table 9 and Supplementary Figure 8A). Using ICEfinder and MegaBLAST against known Pasteurellaceae ICEs (Supplementary Table 10), we identified nine previously reported ICEs in the genomes. We additionally found evidence of 20 putative novel ICEs (Supplementary Table 11), which exhibit typical modular structures, with blocs of core genes related to ICE replication and dissemination interspersed with cargo genes encoding functions that may benefit the host bacterium (such as resistance or virulence-related genes) (Figure 4A). Most ICEs are preferentially integrated into tRNA sites and encode an integrase belonging to the Xer family and a MOB_H relaxase (Supplementary Table 11 and Supplementary Figure 8B). Whole sequence alignment and synteny analysis with previously reported Pasteurellaceae ICEs demonstrated the novelty of the 20 putative ICEs (Supplementary Figures 8C, 9).

Collective analysis of the previously reported and novel ICEs indicated that they ranged from 16 to 78 kb, encoding 23–107 ORFs with GC contents from 37 to 41% (averaging 40%, similar to the genomes). Most of the encoded proteins have predicted functions related to ICE replication/dissemination, as well as antimicrobial and/or heavy-metal resistance, however, few virulence genes were found (Figure 4B). Using the key genes int, parA, and traD (where present), we generated phylogenies to determine evolutionary relationships between the various ICE and identified conserved clusters within some genera which encompasses species that may cohabit in the host, such as Pasteurella spp., Mannheimia spp. and Histophilus spp., and other groups in Haemophilus spp. and Aggregatibacter spp. (Figure 4C). The parA gene was the most informative for tracing the evolutionary history of these ICEs. Results also indicate the presence of putative novel ICEs in seven of the Pasteurellaceae species (Supplementary Figure 8D).

Network analysis revealed that some related ICEs found in the Pasteurellaceae are heterogeneous and widespread among species, e.g., the cluster of ICEs (not yet assigned an ICE family designation) found in M. haemolytica, P. multocida, and H. somni, and the ICEHin1056 family of ICEs found in H. influenzae, Haemophilus parainfluenzae, A. pleuropneumoniae, and Haemophilus ducreyi. By contrast, some ICEs were found exclusively in A. actinomycetemcomitans genomes and do not connect with any other group (Figure 4D). In several H. ducreyi isolates (GHA3, GHA5, GHA8, and GHA9), two coexisting ICEs (ICEHdu2 and ICEHin1056) were identified (data not shown).

Plasmids are extrachromosomal MGEs that contribute significantly to AMR gene dissemination by HGT (Rodríguez-Beltrán et al., 2021). For our analysis, we retrieved 162 complete Pasteurellaceae plasmid sequences from the NCBI RefSeq database (Supplementary Table 12), comprising plasmids from worldwide representatives of 22 species, with Glaesserella parasuis, P. multocida, and A. pleuropneumoniae having the most entries (Supplementary Figure 10). Although all available plasmids are listed in Supplementary Table 12, where different accession numbers exist for the same plasmid, only one representative sequence of each was included in our analysis to avoid bias. Of the 151 different plasmids analyzed, the smallest was 1 kb and the largest 325 kb. The number of ORFs ranged from 1 to 353, with GC contents from 31 to 61% (average of 41%). Using MOBscan, plasmids were classified according to the type of relaxase encoded, with most assigned to the MOB_P, MOB_Q, MOB_V, and MOB_F families (Figure 5A), whereas 52 were designated as potentially non-mobilizable. Results of phylogenetic analysis of the different MOB genes revealed that, even within families, there are polyphyletic clusters, indicating the wide diversity of these sequences (Smillie et al., 2010; Figure 5B and Supplementary Figure 11). Whole-sequence alignments of the potentially non-mobilizable plasmids revealed four clusters of related, but not identical, plasmids present in multiple species (Figure 5C), possibly indicating host-specific divergence following the loss of mobilization function. Notably, four particular plasmids were present in multiple Pasteurellaceae species known to infect different animal hosts, i.e., pIG1 (Wright et al., 1997), pB1002 (San Millan et al., 2010), pB1001/p780 (San Millan et al., 2010; Li Y. et al., 2018), and pB1000 (San Millan et al., 2010; Gangaiah et al., 2016; Bossé et al., 2017), the last of which is found in five species including the human pathogens H. influenzae and H. ducreyi, indicating a One Health concern (Figure 5D).

Comprehensive analysis of all Pasteurellaceae chromosome and plasmid sequences in our datasets revealed a total of 33 different AMR genes were mapped to 478 locations onto 131 genomes from 10 genera (Supplementary Figure 12A). Genes encoding resistance to aminoglycosides, tetracyclines, and sulfonamides were the most represented classes identified (Supplementary Figure 12B). Multiple AMR genes were found in many genomes, and 371 of the 478 AMR genes were associated with MGEs (Supplementary Figure 12C). Also, a wide diversity of AMR genes has been found in species from animal reservoirs such as P. multocida, M. haemolytica and Bibersteinia trehalosi. However, the same was not observed for human pathogens.

Determination of TE genetic context with regards position upstream, downstream, or interrupting ORFs mediating AMR (using the dataset shown in Figure 2C and Supplementary Table 3), we identified ten IS families (ISL3, IS481 and IS30 being most prevalent) mapping to such sites (Figure 6A. In some M. haemolytica genomes, IS481 interrupts and potentially impacts the macrolide-resistance gene macB and the fluoroquinolone (norfloxacin and enoxacin) resistance gene mdtH. IS30 in B. trehalose interrupts marB, a repressor of the marRAB operon involved in activation of AMR and oxidative stress genes (Ariza et al., 1994; Figure 6A). Members of the Tn3, Tn5, and Tn10 transposon families were mapped to ICEs carrying sulfonamide, phenicol, and aminoglycoside resistance genes (Figure 6B).

Although we found prophage-like elements carrying AMR genes in the genomes of H. influenzae and A. indolicus, no AMR genes were found in our dataset of complete prophages (Supplementary Table 13). In contrast, a total of 33 different AMR genes were mapped to 195 locations onto 103 plasmids (68% of the plasmids), and a total of 26 different AMR genes were mapped to 126 locations onto 28 ICEs (65% of the ICEs), including ICEs disseminated among different species of the family, with genes encoding resistance to 10 and 8 classes of antimicrobials, respectively (Supplementary Figures 13A,B and Supplementary Table 13). In plasmids, genes for resistance to aminoglycosides and sulfonamides were most prevalent, whereas genes for resistance to aminoglycosides and tetracyclines were the most common in ICEs (Supplementary Table 13). We assessed the contribution of the ICEs and plasmids carrying AMR genes compared to their respective chromosome, and we found that the majority of AMR genes are exclusively associated with ICEs (Figure 6C).

Some classes of AMR genes were exclusively associated with specific elements, e.g., genes for resistance to rifamycin, lincosamide, and N-glycosides were only found in plasmids, and fluoroquinolone resistance genes only in ICEs (Figure 6D). But most classes of AMR genes were distributed among both ICEs and plasmids (Figure 6D) and several plasmids and ICEs encoded multiple AMR genes (Supplementary Figure 13C).

Copies of the chromosomally encoded kpnH gene, mediating resistance to various antimicrobials such as azithromycin, ertapenem, imipenem, norfloxacin, and polymyxin B (colistin) (Srinivasan et al., 2014), found in human and animal pathogens share 75% nucleotide identity. Whereas, some plasmid and ICE-encoded resistance genes, such as tetB, catII, and bla_ROB–1, show > 99% nucleotide identity between human and animal pathogens. However, only the pB1000 plasmid (carrying bla_ROB–1) shows clear evidence of transfer, with identical sequences found in pathogens from five different host species, including humans (Figure 5D). In contrast, many AMR genes were found exclusively in human pathogens (e.g., hmrM, catS, lpsA, and bla_TEM–1) or animal pathogens (e.g., floR, strA, strB dfrA, ermF, and bla_OXA), despite being MGE-associated. In total, 87.24% (n = 629) of AMR genes were found in animal and 12.76% (n = 92) in human pathogens.

Some AMR genes associated with MGEs were more globally distributed than others in our dataset, e.g., the bla_ROB–1, floR, strA, sul2, and tetB genes were more abundant in isolates from Europe and Asia (Figure 6E). In contrast, 28 AMR genes were associated with specific localities (Supplementary Figure 13D). Neutrality analysis of the most prevalent MGE-associated AMR genes per geographic location indicated possible recent selection for resistance to sulfonamides (sul2) in Europe (Tajima’s D –2.54153, p < 0,001; Fu and Li D*: –4.3437, p < 0.02; Fu and Li’s F* test statistic: –4.43062, p < 0.02) and tetracyclines (tetB) in Asia (Tajima’s D: –1.27745; Fu and Li’s D*: –1.27336, P < 0.02; Fu and Li’s F* test statistic: –1.38944, P < 0.02).