With time more Brucella genomes are constantly becoming available. In this study, Brucella genomes with “complete” and “chromosome” assembly levels were collected from the Genome database of NCBI (accessed on 11 November 2022). After quality control was performed using CheckM 1.2.1 (Parks et al., 2015), 255 genomes with completeness greater than 90% were retained for subsequent analysis. The sample information of the 255 genomes was collected (Supplementary Figure S1; Supplementary Table S1). These genomes were obtained from 244 strains of 14 Brucella species and 11 strains that have not yet been identified. Moreover, these strains were collected from 2001 to 2022 in multiple countries, including Italy, China, Norway, and others. Additionally, these Brucella strains were isolated from a wide range of animal hosts and environments, such as humans, water buffaloes, cattle, pigs, dogs, sheep, and so on. The genomes of B. sp. 09RB8471 (GCA_001971625.1) and B. melitensis 6144 (GCA_023796775.1) were separately compared for collinearity with the genome of B. melitensis 16M (GCA_000007125.1) using Mauve 20150226 (Darling et al., 2004), and chromosome numbering errors in these two genomes were corrected (Supplementary Figure S2).

To ensure comparability, 255 Brucella genomes were annotated using Prokka 1.14.6 (Seemann, 2014). EggNOG-mapper 2.1.5 was utilized for functional annotation and for providing COG (clusters of orthologous groups) information for all the genes (Cantalapiedra et al., 2021). Virulence factor genes (VFGs) and antibiotic resistance genes (ARGs) in the Brucella genome were identified using the Virulence Factors of Pathogenic Bacteria Database (VFDB) (accessed on 27 March 2023) (Liu et al., 2022) and the Comprehensive Antibiotic Resistance Database (CARD) (version: card-data.3.2.9) (Alcock et al., 2023), respectively.

Gene families in the pangenome of 255 Brucella strains were identified using OrthoFinder 2.5.4 (-S diamond-M msa-T fasttree) (Emms and Kelly, 2015; Emms and Kelly, 2019). Core gene families (distributed in all 255 genomes), accessory gene families (distributed in at least 2 genomes, at most 254 genomes), and strain-specific genes (distributed only in one genome) were identified based on the OrthoFinder results. The pangenome growth curve was modeled using Heaps’ law (y = Ax^B + C) (Heaps, 1978; Tettelin et al., 2008), where y represents the number of gene families and x represents the number of genomes. Additionally, a 0–1 matrix file representing the presence/absence pattern of gene families was generated based on the results of OrthoFinder. To obtain a pan tree file, clustering analysis was conducted on the 0–1 matrix file using the R package “pheatmap” (clustering_method = “ward. D2,” clustering_distance = “manhattan”). Moreover, to compare the core tree and pan tree topology, the normalized matching-cluster (nMC) and normalized Robinson–Foulds (nRF) scores were calculated using TreeCmp 2.0-b76 (Bogdanowicz et al., 2012). These two scores range from 0 to 1, with higher values indicating greater differences in the topology between the two trees. Finally, a tanglegram of the two trees was generated using the R package “ape.”

The average nucleotide identity (ANI) values of the pairwise comparison of 255 Brucella genomes were calculated using Pyani 0.2.12 (−m ANIm) (Pritchard et al., 2016). The digital DNA–DNA hybridization (dDDH) values between each of the 255 genomes and B. melitensis 16M were calculated using GGDC 3.0 (Meier-Kolthoff et al., 2013, 2022). The dDDH values of pairwise comparison of the 25 reference genomes from different Brucella species were also calculated using GGDC 3.0. A phylogenetic tree of the 255 Brucella strains was constructed based on the core single-copy orthologous genes. Falsochrobactrum ovis DSM 26720 (GCA_003259955.1) was selected as an outgroup. Specifically, all orthologous groups within the 255 Brucella genomes and one outgroup genome were identified using OrthoFinder 2.5.4. Then, all the core single-copy orthologous gene sequences were aligned using Muscle 3.8.1551 (Edgar, 2021). The alignment products were subsequently trimmed using the Trimal 1.4.rev15 (Capella-Gutierrez et al., 2009). Subsequently, AMAS.py (Borowiec, 2016) was used to concatenate all the alignment results. Finally, using the concatenated multiple sequence alignment file, a maximum likelihood (ML) tree was built using Iqtree 2.0.3 (−m MFP-B 1000 --bnni-T AUTO --seed 123) (Minh et al., 2020). The phylogenetic tree was visualized using ITOL (Letunic and Bork, 2021).

HGTector 2.0b3 (Zhu et al., 2014) was used to identify horizontal transfer gene families (HTGFs) in 255 Brucella genomes that have transferred from various organisms outside the family Brucellaceae to the genus Brucella. Brucella was utilized as the self-group, and Brucellaceae was utilized as the close group.

Insertion sequences (ISs) in 255 Brucella genomes were identified using Isfinder (Siguier et al., 2006). Prophages were detected using PHASTER (Arndt et al., 2016). Integrative conjugative elements (ICEs) and integrative mobilizable elements (IMEs) were identified using ICEfinder (Liu et al., 2019).

Xiaowen Yang et al. and Jagadesan Sankarasubramanian et al. separately conducted studies on the pangenomes of eight (Yang et al., 2016) and 11 (Sankarasubramanian et al., 2017) intracellular Brucella species, respectively. The 255 strains in our study included not only 10 species of intracellular Brucella but also four species of extracellular Brucella. Therefore, compared to those in the aforementioned two studies, the pangenome in our study more comprehensively reflected the species diversity and genetic polymorphisms of the Brucella genus. A total of 9,091 gene families were identified in the pangenome analysis, including 1,404 core gene families, 4,928 accessory gene families, and 2,759 strain-specific genes (Figure 1A; Supplementary Table S2). Several Brucella strains harbored more accessory gene families or strain-specific genes than did the other strains, resulting in a high proportion of accessory gene families and strain-specific genes in the pangenome (Figure 1B). In particular, B. intermedia, B. pseudintermedia, B. anthropi, and B. pseudogrignonensis harbored more accessory gene families and strain-specific genes. With the inclusion of additional genomes, the pangenome expanded continuously (Figure 1C), following Heaps’ law (y = Ax^B + C). Notably, the parameter B was greater than 0; thus, the pangenome in this study was open, suggesting that a large number of genomes need to be analyzed to fully characterize the complete gene pool of the Brucella genus. The open pangenome reflected the genetic diversity of the Brucella genus, indicating the possibility of genetic exchange among the gene pool of Brucella and the external environment. In contrast, the core genome exhibited the opposite trend (Figure 1C), as depicted by a fitted model (y = Ae^Bx + C), in which the size gradually decreased with the addition of new genomes.

To elucidate the impact of core gene families, accessory gene families, and strain-specific genes on the phylogeny of the Brucella genus, a core tree and a pan tree were constructed and compared (Figure 1D). The nMC and nRF scores were computed to assess the differences in tree topology between the core tree and the pan tree. In this study, the nMC and nRF scores were 0.7345 and 0.7984, respectively, indicating substantial disparities between the core tree and pan tree of the Brucella genus. The presence of a large number of accessory gene families and strain-specific genes in the open pangenome of the Brucella genus led to significant differences between the core tree and the pan tree. Furthermore, on the core tree and pan tree, extracellular Brucella and intracellular Brucella were significantly separated. Classical intracellular Brucella and nonclassical intracellular Brucella were separated from each other on the core tree, but there was a crossover phenomenon on the pan tree. This indicates a certain difference between extracellular Brucella and intracellular Brucella in core-pangenome assemblies.

The genomes of Brucella species are composed of two chromosomes. Pan-chromosome analyses were conducted on different chromosomes of Brucella. Based on Heaps’ law (y = Ax^B + C), the openness of chromosome II (Chr II) was greater than that of chromosome I (Chr I) (B_{Chr II} > B_{Chr I}) (Figure 1E), suggesting that Chr II may be more dynamic than Chr I. The two chromosomes of Brucella are considered to have different evolutionary origins, and the smaller Chr II could be derived from a mega-plasmid possessing many dispensable genes (Bavishi et al., 2010). The functional and distribution characteristics of the core and accessory gene families on Chr I and Chr II were further explored, respectively (Figure 1F). Many core and accessory gene families were distributed on Chr I, with the core gene families being dominant, indicating that Chr I of Brucella was relatively conserved and stable. Except for the core gene families associated with “N: cell motility,” “U: intracellular trafficking, secretion, and vesicular transport,” and “G: carbohydrate transport and metabolism,” the abundance of core gene families related to other functions was significantly higher than that of accessory gene families on Chr I (t-test, p < 0.05). In contrast, the accessory gene families were dominant on Chr II, indicating that Chr II of Brucella was more dynamic and plastic and served as the primary reservoir for genetic diversity in the Brucella genus. Except for the accessory gene families associated with “J: translation, ribosomal structure and biogenesis” and “F: nucleotide transport and metabolism,” the abundance of accessory gene families related to other functions was significantly higher than that of core gene families on Chr II (t-test, p < 0.05).

To further explore the genetic diversity of the Brucella genus, a comparison was conducted on the genome sizes and GC contents of different Brucella species. It was evident that, with an average size of 4.69 Mb and an average of 4,327 protein-coding genes per genome, the genomes of extracellular Brucella species, such as B. pseudogrignonensis, B. anthropi, B. pseudintermedia, and B. intermedia, were significantly larger than those of facultative intracellular Brucella species (Figure 2A; Supplementary Table S1). The average genome size of facultative intracellular Brucella was 3.36 Mb, with an average of 2,971 protein-coding genes per genome. Furthermore, the genomic GC content of extracellular Brucella was significantly different from that of facultative intracellular Brucella (Figure 2A, Supplementary Table S1). The genomic GC content of facultative intracellular Brucella was highly conserved, ranging from 57.1 to 57.3%. However, different extracellular Brucella species had significant variations in genomic GC content, with the highest being 57.96% (B. pseudintermedia strain ASAG-D25) and the lowest being 53.71% (B. pseudogrignonensis strain K8). In summary, there were significant differences in genome size and GC content between facultative intracellular Brucella and extracellular Brucella. The genomes of different facultative intracellular Brucella species were conservative, while greater heterogeneity was observed among the genomes of different extracellular Brucella species.

To further assess the phylogenetic relationships among the different Brucella species, the core genome phylogenetic tree was constructed for 255 Brucella strains (Figure 2B). Within the extracellular Brucella, a separate branch was formed by B. pseudogrignonensis. B. intermedia, B. pseudintermedia, and B. anthropi clustered together. All the intracellular Brucella species clustered together. These results are consistent with the findings of Leclercq et al. They found that the entire Brucella genus was comprised of two major branches. One branch encompassed extracellular Brucella not involved in opportunistic infections, such as B. pseudogrignonensis. The other branch contained opportunistic pathogenic extracellular Brucella species (including B. intermedia, B. pseudintermedia, and B. anthropi) alongside intracellular Brucella capable of causing brucellosis (Leclercq et al., 2019). By calculating the ANI and dDDH values between each Brucella strain and the reference strain B. melitensis 16M, significant differences were identified among extracellular Brucella, nonclassical intracellular Brucella, and classical intracellular Brucella (Figure 2B). The ANI values between the extracellular Brucella strains and the reference strain B. melitensis 16M ranged from 83.9 to 85.7%. The ANI values between nonclassical intracellular Brucella strains and the reference strain B. melitensis 16M ranged from 97.7 to 98%. The ANI values between classical intracellular Brucella strains and the reference strain B. melitensis 16M were greater than 99%. The dDDH values exhibited the same characteristics as the ANI values (Figure 2B; Supplementary Figure S3). The dDDH values between extracellular Brucella strains and the reference strain B. melitensis 16M ranged from 22.2 to 25.3%. The dDDH values between nonclassical intracellular Brucella strains and the reference strain B. melitensis 16M ranged from 80.4 to 82.4%. The dDDH values between classical intracellular Brucella strains and the reference strain B. melitensis 16M were greater than 95%. ANI and dDDH values are widely used for prokaryotic taxonomic identification (Sentausa and Fournier, 2013). An ANI value greater than 96% or a DDH value greater than 70% was considered to indicate a single species (Konstantinidis and Tiedje, 2005; Goris et al., 2007). However, these two thresholds were not applicable for species identification within the Brucella genus. Nonetheless, the ANI and dDDH values still revealed significant genomic differences between extracellular Brucella, nonclassical intracellular Brucella, and classical intracellular Brucella (Supplementary Figure S4; Supplementary Figure S5).

Subsequently, the phylogenetic positions of all the undefined Brucella strains were examined. B. sp. 2002734562 clustered together with other B. abortus strains (Figure 2B), suggesting that B. sp. 2002734562 may be a member of B. abortus. The branch closest to B. sp. 2002734562 was occupied by B. abortus RB51-AHVLA (Supplementary Figure S3). B. abortus 2308, B. abortus S19, and B. abortus A19 were also close to the branch of B. sp. 2002734562. B. abortus RB51-AHVLA, B. abortus S19, and B. abortus A19 are vaccine strains of B. abortus (Refai, 2002; Moriyón et al., 2004; Wang et al., 2022). B. abortus 2308 is a virulent strain of B. abortus (Kellerman et al., 1970). Moreover, the other undefined strains of nonclassical Brucella exhibited a close phylogenetic relationship with B. inopinata (Figure 2B; Supplementary Figure S3). B. sp. 1315, B. sp. 2594, B. sp. 09RB8910, B. sp. 2716, and B. sp. 10RB9215 clustered together. Interestingly, these five strains were all isolated from frogs or toads. B. sp. 09RB8471 was also isolated from a frog and clustered together with B. inopinata 141012304 and 4 undefined strains isolated from humans. Although these 10 undefined strains and B. inopinata were separated into two clades, considering that the ANI and dDDH values of the members from the two clades were generally similar, we inferred that these undefined strains may come from the same species, such as B. inopinata.

Horizontal gene transfer (HGT) refers to the transfer of genetic material between different organisms that are not related to the parent or offspring. HGT is an important mechanism in bacterial evolution that can promote the diversity and adaptability of bacterial populations (Jain et al., 1999; Ochman et al., 2000). In our study, the open pangenome of the Brucella genus contained a large number of accessory gene families and strain-specific genes, suggesting that HGT may have played an important role in the evolution of Brucella. Using HGTector, a total of 261 horizontal transfer gene families (HTGFs) in 255 Brucella genomes were identified, which had been transferred from various organisms outside the family Brucellaceae to the genus Brucella (Supplementary Table S3). All genomes of different Brucella species were found to contain HTGFs (Figure 3A). Variations were found in the number of HGTFs among different Brucella species (Figure 3B). The highest number of HTGFs was found in the genome of B. pseudogrignonensis, followed by B. intermedia, with 73 and 58 HTGFs, respectively. The genomes of B. ovis, B. vulpis, and B. inopinata contained fewer HTGFs, averaging 21–25 HTGFs each. The genomes of other Brucella species each had an average of 33–44 HTGFs.

Subsequently, the distribution of the 261 HTGFs in the genome and their relationships with the pangenome of the Brucella genus were further analyzed (Figure 3C). Among the 261 HTGFs, 52 (19.9%) were exclusively distributed on Chr I, 183 (70.1%) were exclusively distributed on Chr II, and the remaining 26 (10%) were distributed on both chromosomes. Furthermore, 19 (7.3%) were core gene families, 203 (77.8%) were accessory gene families, and the remaining 39 (14.9%) were strain-specific genes. The size of Chr II of the Brucella genus was significantly smaller than that of Chr I (Rajendhran, 2021). However, Chr II contained more HTGFs than Chr I, suggesting that Chr II of the Brucella genus might have experienced more HGT events than Chr I. HGT is an important mechanism for bacterial evolution and adaptation (De la Cruz and Davies, 2000). Using Alien_hunter, Anish Bavishi et al. found that the level of HGT on Chr II of B. melitensis 16M and B. anthropi ATCC 49188 was greater than that on Chr I (Bavishi et al., 2010). Our study examined HGT in the whole Brucella genus, aligning with Anish Bavishi et al.’s findings, and showing a higher HGT level in Chr II compared to Chr I in Brucella. Anish Bavishi et al. proposed a plasmid-mediated origin for Chr II as the preferred model explaining this bias of HGT (Bavishi et al., 2010). In this model, a resident plasmid would acquire essential gene clusters through intragenomic gene transfer from Chr I and HGT from other species, leading to the formation of Chr II. This accessory chromosome (Chr II) may then acquire new genes via various mechanisms like HGT. With Chr I housing crucial housekeeping genes, Chr II likely evolved as an adaptive reservoir under selection pressure to accommodate different ecological niches and nutritional challenges and played a vital role in Brucella’s evolution and adaptation.

To further elucidate the functional characteristics of these HTGFs, the COG information of these gene families was analyzed. These HTGFs were involved in diverse biological functions (Figure 3D). Nearly half (48%) of the HTGFs were associated with metabolism. Approximately one-fifth (21%) of the HTGFs were related to information storage and processing. Twelve percent of the HTGFs were related to cellular processes and signaling. However, the functions of the remaining 19% of the HTGFs are still unclear. In particular, the proportion of HTGFs associated with “E: amino acid transport and metabolism” was the highest (19%). The proportion of HTGFs associated with “L: replication, recombination, and repair” was the second highest (14%). Brucella are amino acid prototrophs that can grow with ammonium and contain many amino acid/peptide/polyamine uptake genes (Ronneau et al., 2016). Our research revealed that HGT played a significant role in the acquisition of essential biological functions such as amino acid transport and metabolism in the Brucella genus.

The functional characteristics of the HTGFs were further compared between different Brucella species. These 261 HTGFs were found to encode a total of 134 different proteins by a BLAST search using the NCBI nr database (Figure 3E; Supplementary Figure S6). The HTGFs associated with gene products such as ABC transporter, type IV secretion system proteins, APC family permease, and 15 other types were widely distributed throughout the Brucella genus (Figure 3E). The fact that different Brucella species shared these HTGFs suggested that relevant HGT events may have occurred in the ancient ancestor of the Brucella genus. The HTGFs associated with the transposase family proteins were distributed in all Brucella species except for B. vulpis. These HTGFs may represent historical traces of previous HGT in Brucella, and they may also provide a certain material basis for Brucella genomes to undergo additional HGT.

In classical intracellular Brucella, the HTGFs of B. abortus and B. suis encoded the most gene products, with 41 types (Figure 3E). The HTGFs of B. ovis encoded the fewest gene products, with 23 types. Different classical intracellular Brucella species commonly contained HTGFs associated with 22 gene products, such as ABC transporter permease, ABC transporter ATP-binding proteins, type IV secretion system proteins, and so on. Compared with other classical intracellular Brucella species, B. ovis lacked HTGFs associated with aldehyde dehydrogenase family proteins, fumarylacetoacetate hydrolase family proteins, and MDR/zinc-dependent alcohol dehydrogenase-like family proteins, and so on. Additionally, the abundance of HTGFs associated with ABC transporter permease, ABC transporter ATP-binding proteins in B. ovis was relatively lower than in other classical intracellular Brucella, and these ABC transport systems were primarily associated with amino acid transport (Supplementary Figure S6). ABC transporters play a significant role in the antibiotic resistance and virulence of bacteria (van Veen and Konings, 1998; Paulley et al., 2007). Jenner et al. reported that B. abortus, B. melitensis, B. suis, B. canis, and B. ovis all contain a large number of ABC transport systems and that B. ovis has relatively fewer ABC transport systems than the other four species, primarily lacking polyamine transporters. Considering that B. ovis does not cause zoonotic infections, Jenner et al. speculated that the absence of ABC transport systems in B. ovis may be related to the virulence of Brucella (Jenner et al., 2009). In our study, HTGFs associated with ABC transport systems were relatively less abundant in B. ovis compared to other classical intracellular Brucella species. The differential distribution of HTGFs associated with ABC transport systems and other functional proteins between B. ovis and other classical intracellular Brucella species may be linked to the significant differences in zoonotic potential and host preferences between them. Additionally, B. suis and B. canis have a close phylogenetic relationship, with B. canis often nested within the branch of B. suis on the phylogenetic tree (Fretin et al., 2008). However, significant differences were found in their HTGFs. The HTGFs associated with ACGS family proteins, arginase family proteins, and lysozyme family proteins were unique to B. suis (Figure 3E). These HTGFs were not found in B. canis.

In nonclassical intracellular Brucella, B. vulpis harbored HTGFs related to 20 gene products, among which those associated with DEAD/DEAH box helicase were unique to B. vulpis (Supplementary Figure S6). B. inopinata harbored HTGFs related to 24 gene products but did not contain any unique HTGFs. The unclassified Brucella strains harbored HTGFs related to 59 gene products, among which those associated with 14 gene products such as FAD-dependent oxidoreductase, GNAT family N-acetyltransferase, and helicase RepA family proteins were specific to these unclassified Brucella strains.

Different extracellular Brucella species had HTGFs that correspond to 34–54 gene products (Supplementary Figure S6). The HTGFs associated with 14 gene products such as the AraC family transcriptional regulator, and M20 aminoacylase family proteins were unique to B. intermedia. The HTGFs related to 10 gene products including ACR3 family arsenite efflux transporter, and LysE family translocator were specific to B. pseudintermedia. The HTGFs associated with 11 gene products such as arginine deiminase, and histidinol dehydrogenase were unique to B. anthropi. Additionally, HTGFs related to 19 gene products like amidohydrolase, and ferrochelatase were specific to B. pseudogrignonensis.

In summary, significant differences in the distribution of HTGFs were found among classical intracellular Brucella, nonclassical intracellular Brucella, and extracellular Brucella. In particular, extracellular Brucella was found to contain a greater number of unique HTGFs than intracellular Brucella. Intracellular Brucella has a host preference and an intracellular lifestyle, which to some extent hinders the occurrence of HGT (Rouli et al., 2015). On the contrary, the free distribution of extracellular Brucella in the environment was more conducive to genetic exchange with other organisms. These findings indicate that HGT enhanced the genetic diversity of Brucella and facilitated the acquisition of various important biological functions.

Virulence factors are properties (i.e., gene products) that enable a microorganism to establish itself on or within a host of a particular species and enhance its potential to cause disease (Liu et al., 2022). In our study, a total of 69 VFGs were identified in 255 Brucella genomes (Figure 7; Supplementary Table S8). Among these VFGs, 50 were widely distributed on Chr I of Brucella, while the remaining 19 were primarily distributed on Chr II of Brucella. Several VFGs, including acpXL, bmaB/omaA, and others, were distributed on MGEs in some Brucella genomes. Specifically, wbdA and wboA in B. abortus, B. suis, B. canis, B. ceti, B. pinnipedialis, and B. vulpis, as well as ricA in B. melitensis, B. pinnipedialis, and B. vulpis, were widely distributed on prophages of Chr I. Furthermore, manBcore and manCcore in B. suis, B. canis, and B. microti were mainly distributed on the ICEs of Chr II. The lpxK and waaA genes of B. vulpis were distributed on the IMEs of Chr II.

The classical intracellular Brucella possessed a relatively rich repertoire of VFGs, with each species containing 67–69 different VFGs (Figure 7). B. ovis lacked wbdA and wboA. B. canis did not possess btpA. Among the strains of nonclassical intracellular Brucella, B. vulpis strains, B. sp. 1315, and B. sp. 2594 possessed a considerable number of VFGs, with each strain harboring 65–68 different VFGs (Figure 7). Other nonclassical intracellular Brucella strains possessed fewer VFGs, with each strain containing 46–50 VFGs (Figure 7). In comparison to classical intracellular Brucella strains, these nonclassical intracellular Brucella strains lacked some VFGs associated with adherence and immune modulation (LPS) (Figure 7). The extracellular Brucella strains, which included B. intermedia, B. pseudintermedia, B. anthropi, and B. pseudogrignonensis, harbored the fewest VFGs, with each species possessing only 11–13 VFGs (Figure 7). These VFGs were predominantly located on Chr I of extracellular Brucella, with very few VFGs located on Chr II.

Intracellular Brucella is an important human and veterinary pathogen that is capable of causing brucellosis (Roop et al., 2021). In contrast, extracellular Brucella species represent a diverse group of free-living Brucella, among which only a few species are opportunistic pathogens that occasionally infect humans (Leclercq et al., 2019). In our study, intracellular Brucella possessed significantly more VFGs than did the extracellular Brucella. These VFGs were mainly involved in adherence, host immune modulation, and effector delivery systems and may be key factors leading to the difference in the pathogenicity of intracellular and extracellular Brucella. Additionally, some strains of nonclassical intracellular Brucella lacked a considerable number of VFGs associated with adherence and host immune modulation. These VFGs may cause differences in the pathogenicity between nonclassical intracellular Brucella strains and classical intracellular Brucella strains. Furthermore, although the distribution of VFGs is relatively conserved in classical intracellular Brucella species, B. ovis lacked wbdA and wboA, and B. canis lacked btpA in this study. Both wbdA and wboA are related to the biosynthesis of the lipopolysaccharide O-chain (Prasolova et al., 2023). btpA is related to the regulation of the host immune response (Salcedo et al., 2013). Vizcaino et al. reported the absence of wboA in B. ovis (Vizcaino et al., 2004). To our knowledge, there have been no other reports on the absence of wbdA in B. ovis or the absence of btpA in B. canis, and our research is the first to identify this phenomenon. In summary, different Brucella species exhibit diverse distribution patterns of VFGs, and these polymorphisms are likely closely related to variations in virulence, pathogenicity, and lifestyle among different Brucella species.

The inappropriate use of antibiotics has promoted a widespread increase in antibiotic resistance in bacteria. The HGT of ARGs encoded on MGEs is one of the key driving factors in this process (Botelho and Schulenburg, 2021). In our study, a total of 14 ARGs were identified in the 255 Brucella genomes (Figure 7; Supplementary Table S9). These ARGs confer antibiotic resistance through mechanisms such as antibiotic efflux, antibiotic inactivation, and antibiotic target alteration. mprF confers resistance to peptide antibiotics. Ant(9)-Ic and aph(3′)-Ia confer resistance to aminoglycoside antibiotics. catB10 confers resistance to phenicol antibiotics. It should be noted that catB10 only exhibits a query coverage of 73% in this study, and this incomplete sequence of catB10 may result in non-functional resistance of Brucella. ceoB confers resistance to aminoglycoside and fluoroquinolone antibiotics. bla_OXA-919 confers resistance to penam antibiotics. The remaining 8 ARGs confer resistance to penem, penam, cephamycin, cephalosporin, and monobactam antibiotics. There are nine ARGs that are beta-lactamase genes, including bla_OCH-1, bla_OCH-2, bla_OCH-3, bla_OCH-4, bla_OCH-5, bla_OCH-6, bla_OCH-7, bla_OCH-8, and bla_OXA-919. Moreover, mprF and ceoB were located primarily on Chr II of Brucella, while the remaining 12 ARGs were distributed only on Chr I (Figure 7). Furthermore, except for aph(3′)-Ia and mprF, the remaining 12 ARGs were distributed only among extracellular Brucella. The genome of B. abortus clpP contained both aph(3′)-Ia and mprF. The other intracellular Brucella strains contained only mprF.

In this study, various MGEs in Brucella genomes, including prophages, ICEs, and IMEs, were devoid of any ARGs, suggesting that the role of MGEs in the dissemination of ARGs might be limited in the Brucella genus. Furthermore, all the intracellular Brucella strains commonly harbored only one ARG, mprF. The remaining 12 ARGs were exclusively present in extracellular Brucella, suggesting that these two types of Brucella may have significantly different antibiotic resistance. Intracellular Brucella and extracellular Brucella exhibit distinct survival habits. The parasitism of intracellular Brucella may have limited their acquisition of ARGs through HGT. In contrast, extracellular Brucella species exist mainly in extracellular environments, which facilitates frequent genetic exchange and the acquisition of ARGs. In addition, recent studies have revealed that the resistance of Brucella spp. to various antibiotics, such as rifampicin, azithromycin, tetracycline, and so on, has been increasing (Wareth et al., 2021, 2022; Celik et al., 2023; Rezaei Shahrabi et al., 2023; Beig et al., 2024). These studies indicate an increase in antibiotic resistance in intracellular Brucella, which urgently requires attention. However, the distribution of ARGs within the Brucella genus was found to be very limited in our study. This may be attributed to the limitations of the CARD database that was used in our research, which may lack sufficient information on Brucella ARGs. In addition, more complex mechanisms may also contribute to drug resistance in Brucella.