The classical Bordetella species include Bordetella bronchiseptica, Bordetella parapertussis, and Bordetella pertussis. B. pertussis is a severe human pathogen and the causative agent of pertussis, a respiratory disease that can be particularly severe in infants and the elderly, up to a fatal outcome [1–3]. Bordetella parapertussis comprises two host-associated lineages: a human-adapted lineage that causes pertussis and a sheep-adapted lineage that causes respiratory disease in sheep [4]. B. bronchiseptica has a wider range of hosts, causing acute and chronic respiratory infections in various mammals [3, 5]. Human infections are relatively uncommon, but their clinical manifestations range from asymptomatic carriage to bronchitis, pneumonia, and sepsis [5].

Classical Bordetella species possess broadly similar sets of core virulence factors, including surface adhesins, fimbriae, toxins, and secretion systems [3, 4]. The composition and functionality of these factors vary among species and across phylogenetic lineages. In B. bronchiseptica, the ptx–ptl region, which encodes pertussis toxin and its secretion apparatus, is absent from some phylogenetic lineages and may be lost or highly divergent in individual genomes [4]. In strains retaining this locus, pertussis toxin production may be disrupted by changes in the promoter region. However, it has been experimentally shown that the introduction of four nucleotide substitutions into the ptx promoter of B. bronchiseptica was sufficient to activate the expression of pertussis toxin to the level of B. pertussis [6].

The historical differentiation of classical Bordetella into separate species was based on a set of phenotypic, clinical, and environmental characteristics, including the range of hosts, the nature of the disease caused, and cultural and biochemical properties [7]. Later, comparative genomic studies revealed that the differences between these species resulted from gene loss and inactivation, the accumulation of insertion sequences, and genomic rearrangements [8, 9]. In particular, the changes affected the genes of surface structures, metabolic pathways, and regulatory systems, and could have contributed to adaptation to different hosts [8].

Despite the preservation of three independent species names, classical Bordetella are characterized by high genomic similarity. The average nucleotide identity (ANI) values between the genomes ranged from 98.3% to 99.4%, which is above the commonly used 95–96% threshold for species delineation in bacteria [2]. That is why Bridel et al. assigned B. bronchiseptica, B. parapertussis, and B. pertussis to the same group, the B. bronchiseptica genomic species [2]. This grouping does not cancel out the traditional names of the species but rather emphasizes their kinship according to modern genomic criteria.

Phylogenomic studies have also shown that B. pertussis and the human lineage B. parapertussis independently originated from different B. bronchiseptica-like lineages [8–10]. Their adaptation to a restricted host range was accompanied by genome reduction, pseudogene accumulation, insertion sequence proliferation, and the loss of certain metabolic and regulatory functions [4, 8–10].

The collection of the G. N. Gabrichevsky Research Institute of Epidemiology and Microbiology in Moscow contains 31 previously uncharacterized isolates of B. bronchiseptica, including 25 from human lineages and 6 from animal hosts. Comparing this collection with the reference genomes of classical Bordetella species makes it possible to determine the position of the studied isolates relative to the known phylogenetic lineages of B. bronchiseptica and to assess whether the observed phylogenetic differences are associated with variation in functionally significant genomic traits. This study investigated the phylogenetic placement, genomic features, and diversity of these isolates relative to known classical Bordetella lineages.

METHODS

We studied a collection of 31 B. bronchiseptica isolates. Twentyfive isolates were obtained from human patients by the Pertussis Monitoring Reference Center of the G. N. Gabrichevsky Research Institute of Epidemiology and Microbiology. The remaining six isolates were collected from animals by specialists at the K. I. Scriabin Moscow State Medical University: three from pigs, two from rabbits, and one from a domestic cat.

The array used for comparative analysis also included 30 reference genomes from NCBI: 10 genomes each of B. bronchiseptica, B. pertussis, and B. parapertussis from human-associated lineages. The strains were selected using the BacDive database (DSMZ); for each species, we included a typical strain. The genomic assemblies of the selected strains were downloaded from NCBI using the identifiers provided by BacDive. Table S2 (supplement) lists reference strains and their identifiers.

Culturing of isolates

The isolates were cultured in accordance with methodological guidelines MUK 4.2.3701-21 "Laboratory diagnostics of pertussis and other diseases caused by other Bordetella." Lyophilized cultures were stored at 2-8 °C. Before the study, the cultures were rehydrated in 0.5 mL of sterile 0.9% NaCl solution, after which 0.1 mL of each suspension was plated onto Bordetelagar (State Research Center for Applied Microbiology and Biotechnology, Obolensk) supplemented with 10% defibrinated sheep blood, and incubated at 37 °C for 24 h.

Morphological and biochemical identification

Colony morphology was assessed using a SteREO Discovery V12 stereomicroscope (Carl Zeiss, Germany). Staining characteristics were determined by Gram staining followed by light microscopy. Biochemical identification included assessment of catalase, oxidase, tyrosinase, and urease activities, as well as nitrate reduction and citrate assimilation. Motility was determined based on the growth pattern in semisolid nutrient agar. All studies were performed in accordance with MUK 4.2.3701-21.

Genome-wide sequencing, assembly, and annotation of genomes

Genomic DNA was isolated using the ExtractDNA Blood & Cells kit (Eurogen, Russia) according to the manufacturer's instructions. DNA integrity was assessed by electrophoresis in 1.5% agarose gel, and concentration was assessed with a Spectra Q HS Plus kit (Raissol, Russia) on a Qubit 2 fluorimeter (Invitrogen, USA).

Libraries were prepared from 220 ng of DNA using the MGIEasy Fast PCR-Free FS DNA Library Prep Set V2.0 (MGI, China), pooled at equimolar ratios, and circularized using the MGIEasy Dual Barcode Circularization Kit (MGI, China). Sequencing was performed on the DNBSEQ-G50 platform (MGI, China) in 2 × 150-bp paired-end mode using an FCL PE150 flow cell.

Quality control and read filtering were performed using fastp v0.24.1 [11], de novo assembly using Unicycler v0.5.1 [12], and assembly quality appraisal using QUAST v5.2.0 [13] run with default settings. The completeness and contamination of assemblies were assessed using CheckM2 v1.1.0 [14]. Assemblies meeting the quality thresholds of > 95% completeness and < 5% contamination were selected for pangenomic and metabolic analyses. Table S1 presents assembly quality indicators, including the number of contigs, the N50 value, completeness, and contamination. Genome annotation was performed in Bakta v1.12 with default settings [15].

ANI and phylogenetic analysis

To assess the genetic similarity, we calculated ANI using OrthoANI v.1.40 [16]. Phylogenetic analysis based on 120 BAC120 marker genes was performed using GTDB-Tk v2.6.1 [17] to reconstruct the phylogeny. After the 120 marker proteins had been aligned, the standard canonical GTDB-Tk mask was applied, reducing the length of the concatenated amino acid alignment from 41,084 to 5,036 positions. The resulting tree was rooted by Achromobacter xylosoxidans DSM 2402 (GCA_000508285.1) as a closely related outgroup selected based on a BLAST search for the Bordetella bronchiseptica rpoB gene homolog with exclusion of representatives of the genus Bordetella [18]. This choice is further supported by published phylogenetic data [19, 20].

To assess the stability of the inferred positions of the major phylogenetic groups with respect to marker set selection, we additionally reconstructed phylogenies based on 42 ribosomal proteins. Using UBCG2 v2.0 [21], we identified 81 conserved genes for phylogenetic analysis. Forty-two genes encoding ribosomal proteins belonging to the rpl, rpm, and rps groups were selected from this set for subsequent analysis. Next, we aligned the amino acid sequences using MAFFT v.7.48 (L-INS-i) [22] and concatenated the resulting alignments using AMAS v.1.0 [23]. The maximum likelihood tree was inferred using IQ-TREE v2.4.0 [24], with automatic model selection and 1000 bootstrap replicates, and visualized using iTOL [25].

Search for antibiotic resistance genes and virulence factors

The genes of antibiotic resistance and virulence factors were identified in ABRicate v.1.0.1 [26], using CARD and VFDB databases. The matrices of the number of detected hits were built in R using tidyverse, and visualized with ComplexHeatmap.

The analysis included matches with a nucleotide identity of ≥ 60% and coverage of the reference sequence of ≥ 80%.

Search for antiphage systems

Antiphage systems were identified in 61 Bordetella genomes using PADLOC v.1.1.0 run under default settings [27].

Search for prophage regions

Prophage regions in 61 Bordetella genomes were identified using PHASTEST [28] and Phigaro v.1.1.0 [29] run under default settings. PHASTEST classified prophages as intact (score > 90), doubtful (70–90), or incomplete (< 70), and Phigaro refined their coordinates and characteristics.

Analysis of the pangenome and metabolism

We analyzed the pangenome and metabolic potential of 20 genomes, including 13 laboratory isolates of B. bronchiseptica selected based on assembly quality and reference genomes from NCBI representing B. bronchiseptica, B. pertussis, B. parapertussis, and A. xylosoxidans DSM 2402. The analyses were performed using Anvi'o v9 [30]. The pangenome was constructed with anvi-pan-genome (--minbit 0.5, --minpercent-identity 30, and --mcl-inflation 2, NCBI BLAST). To compile functional annotation, we used COG24 and KEGG/KOfam, and anvi-estimate-metabolism to assess the completeness of metabolic modules.

RESULTS

Phenotypic features

On Bordetelagar, colonies of B. bronchiseptica were round, convex, moist, smooth, shiny, grayish-cream in color, with smooth edges. Their consistency was soft and oily; the colonies could easily be removed with a loop (fig. 1A). When Gramstained, B. bronchiseptica appeared as short Gram-negative rods, occurring singly, in pairs, or in groups in the smears (fig. 1B). In biochemical tests, B. bronchiseptica strains were urease-positive, catalase-positive, and oxidase-positive but tyrosinase-negative; they grew on Simmons citrate agar and reduced nitrate to nitrite. In addition, B. bronchiseptica grew on 10% blood agar and MPA and exhibited motility, as indicated by diffuse growth along and around the inoculation stab in semisolid agar.

Average nucleotide identity (ANI) analysis

Comparative analysis of the ANI of 61 Bordetella genomes, including 31 genomes of the studied B. bronchiseptica isolates and 30 reference genomes, revealed several highly similar groups (fig. 2). Within these groups, the ANI values were 99.2–100%, and in less similar groups they were 97.2–99.9%. All the obtained values exceeded the commonly used threshold of 95–96% for distinguishing bacterial species, confirming the high genomic similarity among classical Bordetella species. The overlap of the ranges reflects the absence of a single threshold value of ANI separating the selected groups; their differentiation was based on the structure of the dendrogram of pairwise ANI values, rather than on a fixed threshold.

GTDB-Tk Analysis

A phylogenetic analysis based on 120 marker proteins from 61 Bordetella genomes, with the A. xylosoxidans DSM 2402 genome used as an outgroup, showed that the B. pertussis genomes and the human-associated B. parapertussis genomes formed separate, compact clades, whereas the B. bronchiseptica genomes did not form a single monophyletic group (fig. 3).

Analysis of ribosomal proteins

Phylogenetic analysis based on a concatenated amino acid sequence alignment of 42 ribosomal proteins from the same 61 Bordetella genomes and the A. xylosoxidans DSM 2402 genome has also shown that the B. bronchiseptica genomes did not form a single compact monophyletic group (fig. 4). The studied isolates formed several phylogenetic groups, some of which clustered near the B. pertussis and B. parapertussis clades in the phylogenetic tree. The B. pertussis and B. parapertussis clades received bootstrap support values of 99–100%. Individual groups of B. bronchiseptica were also characterized by high support, while bootstrap support values for a number of internal branches were below 50%.

Antibiotic resistance genes and virulence factors

Isolates 882, 883, 884, and 937-1 had the highest number of resistance determinants. Isolates 882 and 883 had the same profile of six genes: BOR-1, APH(3')-Ia, aadA, QACEID1, sul1, and tet(A). Fragmented sequences of APH(3')-Ia and sul1, as well as aadA and qacEΔ1, were found in isolate 884. BOR-1, TEM-112, APH(3')-Ia, catA1, and the evgA regulatory gene were detected in isolate 937-1.

Of the determinants analyzed, only BOR-1 was detected in all 10 B. parapertussis genomes, whereas acquired antibioticresistance genes were not detected in the B. pertussis genomes at the specified thresholds. The prn gene, which encodes pertactin, was found in all B. pertussis and B. parapertussis genomes, as well as in 39 of the 41 B. bronchiseptica genomes; prn was not detected in the assemblies of isolates 884 and 885.

The ptx–ptl locus encoding the pertussis toxin subunits and its secretion system was found in all 10 genomes of B. pertussis and B. parapertussis, and in 31 of the 41 genomes of B. bronchiseptica. Its distribution corresponded to the phylogenetic structure of the sample: ptx–ptl-positive B. bronchiseptica genomes were distributed across several branches, whereas the 10 ptx–ptl-negative genomes (14–23, 19–24, 21–23, 25–24, 69–23, 69–24, 7–13, 899t, F709, and FDAARGOS_634) formed a separate branch.

In addition to the ptx–ptl locus, a wide range of other virulence factors were identified in the studied genomes. The bvgA and bvgS genes encoding a two-component virulence regulation system were found in all the genomes of B. pertussis and B. parapertussis, and in 40 and 38 of the 41 genomes of B. bronchiseptica, respectively. The components of the type III secretion system (bsc, bopB, bopD, bopN, bsp22, bteA/ bopC) involved in interaction with host cells and delivery of effector proteins were widely represented in all three species; most of the relevant genes were identified in 39-41 of the 41 genomes of B. bronchiseptica. The adenylate cyclase toxin locus was also predominantly conserved: CyaA was detected in 38 of the 41 genomes of B. bronchiseptica, as were cyaC and cyaBDE, the locus genes associated with toxin activation and secretion. The dnt gene encoding the dermonecrotic toxin was found in 32 of the 41 genomes of B. bronchiseptica and in all the studied genomes of B. pertussis and B. parapertussis. Greater variability was observed among the adhesion factors: prn was identified in 39 of 41 B. bronchiseptica genomes, fim3 in 39 of 41, and fim2 in 28 of 41. The fim2 gene was present in all B. pertussis genomes but absent from all studied human B. parapertussis genomes. The greatest coverage variability was observed for fhaB, which encodes filamentous hemagglutinin. Although the initial ABRicate analysis identified fhaB matches in 17 of the 41 B. bronchiseptica genomes, only eight met the predefined coverage threshold of ≥ 80%.

Prophage regions

The analysis of the prophage regions revealed their presence exclusively in the B. bronchiseptica and B. pertussis genomes. B. parapertussis genomes had no prophage regions.

Among the 47 detected prophage regions, 42 were classified as intact. The analysis revealed that the prophages in the studied isolates belong mainly to phages infecting Pseudomonas, Burkholderia, and Bordetella. Most prophages were identified as siphoviruses. In B. bronchiseptica 20-24, 69-24, 7-13, and S-55 isolates we found coexisting myoviruses, podoviruses, and siphoviruses (phages of three different morphological types).

Antiphage protection systems

The range of predicted protective loci was the widest (from 4 to 29 per genome) in the B. bronchiseptica isolates. In particular, we detected RM systems of types I and II, AbiE, dXTPase, Lamassu, PDC protective cassettes (M06, M40, S02, S34), and Zorya systems of types I–III. All 10 genomes of the human B. parapertussis lineage had the same set of 13 protective loci, including RM systems of types I and II, AbiE, dXTPase, PDC_S19, and Zorya types I–III. Seven protective loci were identified in the genomes of B. pertussis; the protective systems identified there were dXTPase, PDC_M06, and SoFic, while the RM systems were absent.

Pangenomic analysis and metabolism analysis

We identified 9348 gene clusters in the pangenome constructed from 20 selected genomes. The functions of 6736 clusters were annotated, whereas the remaining 2612 could not be functionally characterized (fig. 6). Within the studied subsample, B. bronchiseptica genomes differed in the presence or absence of several gene clusters, whereas some clusters were absent from the B. pertussis and B. parapertussis genomes.

In all 20 genomes included in the metabolic analysis, 58 main modules were ≥ 75% complete (fig. 7; see Figure S1 in the supplement for a full-size version). Glycolysis, gluconeogenesis, the pentose phosphate pathway, the tricarboxylic acid cycle and the glyoxylate cycle were preserved, as well as the main biosynthetic pathways of amino acids, nucleotides, fatty acids, phospholipids, cofactors and vitamins, and the main components of the respiratory chain.

The most pronounced differences among the signature modules involved the pertussis pathogenicity signature, specifically the pertussis toxin module, which was detected in B. pertussis and in genomic regions of B. bronchiseptica and B. parapertussis. In the clade that included B. bronchiseptica isolates 69-23, 69-24, 25-24, 21-23 and 14-23, this module was not detected. A number of modules present in B. bronchiseptica and B. parapertussis, including those related to nitrogen metabolism, degradation of aromatic compounds, and biosynthesis of nucleotide sugars, were also not found in the three B. pertussis genomes included in the analysis.

DISCUSSION

A genome-wide comparison confirmed the close phylogenetic relationship of classical Bordetella. The ANI values between the studied and reference genomes exceeded the species-level threshold, supporting the classification of B. bronchiseptica, B. parapertussis, and B. pertussis as related taxa within the B. bronchiseptica genomic species [2]. The analysis of the ANI heat map and phylogenetic trees showed heterogeneity within this group. Phylogenetic reconstructions based on BAC120 markers and ribosomal proteins produced a similar picture. B. pertussis and B. parapertussis formed compact clades, whereas B. bronchiseptica isolates were distributed among several lineages. This is consistent with the published data on the origin of B. pertussis and the human lineage of B. parapertussis from different B. bronchiseptica-like ancestors [8–10]. Therefore, in this sample, B. bronchiseptica is represented not as a single monophyletic group, but as a heterogeneous complex of close lineages.

In addition, bootstrap support below 50% for several internal branches in B. bronchiseptica indicates limited resolution of the relationships among the most closely related genomes based on the conserved ribosomal markers. Therefore, the topology of these branches should be interpreted with caution.

The distribution of virulence factors was also uneven. The ptx–ptl locus was absent from some B. bronchiseptica isolates, and genomes lacking this locus formed a distinct phylogenetic branch. Thus, the distribution of the ptx–ptl locus was associated with the phylogenetic positions of the studied B. bronchiseptica genomes. However, the presence of the locus alone does not demonstrate the production of functional pertussis toxin; assessment of the relevant regulatory regions and experimental confirmation of gene expression and toxin secretion are required to evaluate pertussis toxin production [6].

The variability also affected other toxin-associated factors. The CyaA gene encoding adenylate cyclase toxin was detected in most, but not all, B. bronchiseptica genomes, and dnt encoding dermonecrotic toxin was absent in some of the studied genomes of this species. At the same time, most components of the type III secretion system (bsc, bopB, bopD, bopN, bsp22, bteA/bopC) were widely represented in the studied sample. The bvgA and bvgS genes encoding a twocomponent system of global virulence regulation were also preserved in most genomes.

Differences were also found among the surface factors involved in adhesion and colonization. The genes prn, encoding pertactin, and fim3, encoding one of the main subunits of fimbriae, were found in most of the genomes of B. bronchiseptica, whereas fim2 was not detected in the studied human genomes of the B. parapertussis lineage. The greatest variability in coverage was observed for the fhaB gene. Incomplete coverage of this gene in some isolates may reflect both real interlinear variability and fragmentation of individual genomic assemblies, so the result requires careful interpretation. The variable distribution of fim2 also indicated differences in the set of adhesive factors between individual groups of classical Bordetella species.

In the studied sample, genomes of B. bronchiseptica carried a more diverse set of potential antibiotic-resistance determinants than did those of B. pertussis and B. parapertussis. The largest number of such determinants was found in isolates 882, 883, 884 and 937-1. These results should be interpreted with caution: the study's main purpose was to investigate the presence of relevant genes rather than phenotypic sensitivity to drugs.

The greatest variability in prophage regions and antiphage protection systems was also observed among the B. bronchiseptica isolates. This may reflect the broader ecological niche of the species and more frequent contacts with different phages. In B. pertussis and B. parapertussis, such profiles were more of the same type, which is consistent with the specialization of these species to a specific host and the reduction of their genome. However, the functionality of the predicted prophages and protective systems requires a separate experimental verification.

Pangenomic and metabolic analyses supplemented the phylogenomic data. The studied subsample of "classical" Bordetellae retained a common genomic core and the main pathways of central metabolism. At the same time, B. pertussis and B. parapertussis lacked some of the gene clusters and metabolic modules in comparison with B. bronchiseptica, which aligns with the previously described reduction of the genome during adaptation to a narrower host range [8, 9].

Thus, the studied B. bronchiseptica collection isolates represent several phylogenetic variants within the complex of "classical" Bordetellae. The revealed differences between the lineages affect not only their phylogenetic position, but also the composition of genes associated with virulence, antibiotic resistance, prophages, antiphage protection systems, and individual metabolic pathways. This suggests that these differences are evolutionary in nature and may reflect adaptation to various ecological niches, including the host.

CONCLUSIONS

The analyzed B. bronchiseptica isolates were distributed across several phylogenetic lineages, whereas B. pertussis and human-lineage B. parapertussis formed compact clades. The B. bronchiseptica lineages differed in the composition of their virulence factors, putative antibiotic-resistance determinants, prophages, and antiphage defense systems. The distribution of the ptx–ptl locus was consistent with the phylogenetic structure of the isolates. Despite the pronounced phylogenetic and functional-genomic heterogeneity of B. bronchiseptica, the classical Bordetellae exhibited high genomic similarity and conserved core metabolic functions.