Cases of Burkholderia infections detected in patients with CF attending The Prince Charles Hospital and the Mater Hospital, Queensland, Australia, between 2017 and 2022 were included in the study. Clinical specimens were collected from the Prince Charles Hospital as per the Research Collaboration Agreement HREC/17/QPCH/277, and the Mater Hospital as per HREC/14/QPAH/275. Bacterial isolation was performed from expectorated sputum samples. Briefly, a direct inoculum from sputum was plated out onto Burkholderia cepacia selective agar (BCSA) and incubated at 37°C under aerobic conditions for 48 h. The medium, supplemented with crystal violet, vancomycin and gentamicin, inhibits growth of gram-positive bacteria and gram-negative bacilli Pseudomonas, for increased selectivity of Burkholderia. Bacteria were grown under aerobic conditions in lysogeny broth (LB) (1% (w/v) tryptone, 1% (w/v) yeast extract, 0.5% (w/v) sodium chloride) at 37°C, with shaking at 220 rpm overnight.

For genomic DNA extraction, bacterial cultures were grown in 5 mL of LB, followed by overnight incubation at 37°C with shaking at 220 rpm. Genomic DNA was extracted using a QIAamp DNA Mini Kit (Qiagen) as per the manufacturer’s instructions. DNA concentration and purity were assessed using a Nanodrop 1,000 spectrophotometer. Whole genome sequencing (WGS) was performed by MicrobesNG¹ using the Illumina sequencing platform with 2 × 250 bp pair-end reads. Genome assembly was conducted by MicrobesNG using their standard analysis pipeline. Briefly, Kraken was used to identify the closet available reference genome (Wood and Salzberg, 2014). The reads were mapped to the reference genome using BWA mem to assess the quality of the data (Li and Durbin, 2009). De novo assembly of the reads was conducted using SPAdes (Bankevich et al., 2012). Genome annotation was completed using Prokka (Seemann, 2014). Reference genomes were downloaded from the NCBI database² (Supplementary Table 1).

Genome quality assessment and assembly metrics were assessed using CheckM (v1.1.2) (Parks et al., 2015). Genomes with an estimated completeness of ≥90% and contamination of ≤5% were determined high-quality and included for comparative analysis. Multi-locus sequence typing (MLST) was determined by querying genomes against the Burkholderia cepacia complex typing database³ available on PubMLST (Jolley et al., 2018). Under the MLST scheme seven independent loci comprising of genes atpD, gltB, gyrB, recA, lepA, phaC, and trpB are assessed to identify differences in the allelic profile to define the sequence type (ST). Taxonomic classification was assigned using the Genome Taxonomy Database Toolkit (GTDB-Tk; v2.2.4) with reference database r207. A concatenated set of 975 core genes were produced using BPGA and the resultant output used to predict phylogeny (Parks et al., 2018; Supplementary Data 1). Phylogenetic analysis was conducted using RAxML (v8.2.13) (Stamatakis, 2014) and visualised in iToL (v6.8.1).⁴ Average nucleotide identity (ANI) was determined using FastANI (v1.1) and visualised using the ‘pheatmap’ (v1.0.12) and ‘ggplot’ (v3.1.3) packages in Rstudio (v4.1.3) and R (v4.2.3). Genomes were annotated through EnrichM (v0.4.15⁵) using the ‘annotate’ function, with the input parameter ‘--orthologs’ to annotate orthologous genes.

A total number of 34 B. anthina genomes were included in the comparative genomic analyses. Statistical analyses for comparing distribution of genome size, number of CDSs, and coding density were determined using an unpaired t-test. The Bacterial Pan Genome Analysis (BPGA) pipeline using the USEARCH tool was used to cluster protein sequences into gene families (Chaudhari et al., 2016). A sequence identity cut-off value of 80% was used for orthologous clustering. The pan-genome profile analysis tool was used to determine the pan and core genome size, calculated with 100 iterations. Data are calculated as a function of genomes introduced into the analysis. The pan-genome functional analysis tool was used to assign KEGG Ortholog pathways to protein sequences, and execute comparative functional analysis for core, accessory, and unique genes. Principle component analysis (PCA) of KEGG orthologs was performed using the ‘factoextra’ (v1.0.7), ‘ggplot2’ (v3.4.4), ‘purrr’ (v1.0.2), ‘scatterplot3d’ (v0.3-44), ‘ade4’ (v1.7-22), and ‘gdata’ (v3.0.0) packages in Rstudio (v4.1.3) and R (v4.2.3). Comparative genomics based on KEGG Orthology was completed using EnrichM (v0.4.15; see footnote 5). Genomes were annotated using the ‘annotate’ function, with the input parameter ‘--ko’ used to annotate KEGG Orthology. Statistical analyses between clades were performed in EnrichM by Fisher’s exact test. Corrected p-values of <0.05 were considered significant. Secondary metabolite gene clusters were predicted using anti-SMASH (v7.0.0⁶) using the parameter ‘relaxed’ to identify hits (Blin et al., 2023). In silico characterisation of the ethylenediaminesuccinic acid hydroxyarginine (EDHA) biosynthetic gene cluster was conducted through BLASTn searches. Streptomyces sp. MA5143a was used as the reference sequence and encodes for an EDHA operon containing four genes denoted AesA, AesB, AesC, and AesD. The genetic organisation of the EDHA gene cluster was visualised using the ‘gggenes’ (v0.5.1) and ‘ggplot2’ (v3.4.4) packages in Rstudio (v4.1.3) and R (v4.2.3). Identification of iron genes and iron gene operons was predicted using the bioinformatic tool FeGenie (Garber et al., 2020). The RAST server (v2.0) was used to annotate genomes and identify genes associated with iron metabolism, specifically siderophore biosynthesis (Aziz et al., 2008; Overbeek et al., 2014; Brettin et al., 2015).

Taxonomic analysis of Burkholderia previously relied on recA gene sequence analysis, with the molecular approach proving to be a valuable tool in the early distinction between species. Although widely implemented and used in epidemiological investigations, its utility to resolve in depth phylogenetic relationships among closely related species within the BCC is limited. To further understand the evolutionary relationships of Burkholderia, we performed genomic analyses among and within species of the BCC. Here, we genome sequenced 18 clinical isolates cultured from patients with CF in Queensland, Australia. The quality of each genome assembly was assessed using CheckM, with genomes more than 95% complete, and < 5% contaminated retained for analysis (Table 1). A total of 177 Burkholderia genomes were included as part of this study, including 159 publicly available genomes, with B. gladioli selected as the outgroup species. The genome size of Burkholderia strains isolated in this study ranged from 6.3 to 8.3 Mbp, with a G + C content of 67%, and predicted number of coding genes ranging from 5,600–7,600 genes (Table 1).

Phylogenetic analysis was performed using a set of 975 core genes produced using BPGA and inferred using RAxML. MLST analysis was performed to inform strain relatedness among isolates. As part of this study, we obtained four isolates belonging to B. cenocepacia III-A, and one isolate belonging to B. cenocepacia III-B. All B. cenocepacia isolates cultured in this study clustered separate from each other and displayed a unique ST (Figure 1). While five of the seven cultured B. multivorans isolates (BCCIQ02A, BCCIQ03A, BCCIQ04C/D and BCCIQ07B) displayed the same ST-622 and clustered tightly together. Phylogenetic placement of Burkholderia genomes identified several distinct phylogenetic clades within B. anthina. Notably, all five B. anthina isolates clustered within Clade II, displayed the same ST-2133, and represent the first clinical isolates reported within the group. Based on a well-supported phylogeny, it is evident that B. anthina Clade II represents a truly distinct group to B. anthina Clade I. As shown in Figure 1, B. anthina is polyphyletic with Clade I, as well as B. vietnamiensis, and B. ambifaria, diverging earlier than B. anthina Clade II. Phylogenetic placement of B. anthina Clade II suggest that this lineage shares a more recent common ancestor with B. cenocepacia, B. cepacia, B. stabilis, and B. pyrrocinia.

The phylogenetic separation in B. anthina is further reflected in the genomic architecture between the two clades. B. anthina Clade I genomes have a significantly larger genome than B. anthina Clade II (p < 0.0001), with a mean size of 7.7 Mbp as compared to 7.4 Mbp (Figure 2A). Accordingly, B. anthina Clade I have a significantly greater number of predicted coding sequences than B. anthina Clade II (p = 0.0029), with an average of 472 more coding sequences (Figure 2B). The coding density of both clades was comparable with an average of 87% (Figure 2C). The pan-genome of B. anthina demonstrated characteristics of an “open” pan genome, composed of a core genome of 3,371 orthologous genes (12.2%) shared by all strains, and a pan genome of 24,274 genes (87.8%) (Figure 2D). The gene accumulation rarefaction curve increased with the addition of each new genome, reflecting an increase in the number of gene families, and expanding the gene pool diversity. When analysed within their respective clades, Clade I (Figure 2E) and Clade II (Figure 2F) genomes also revealed an “open” profile. At an analysis of 16 genomes, Clade I demonstrated the greatest diversity among strains with only 4,237 core genes (19.4%) and 17,624 (80.6%) pan genes. This contrasts with Clade II which had 5,464 core genes (38%) and 9,020 (62%) pan genes. The pan and core genome profile of B. anthina Clade I highlight the genetic variability within the group and reflects the differences among genomes.

Whole genome average nucleotide identity (ANI) was calculated to determine strain relatedness and assess the species boundaries of B. anthina and closely related species as indicated by phylogenetic placement. A threshold of 98% ANI is considered the cut-off for strain classification within species, whereas a threshold of 95% ANI is considered the cut-off for species delineation (Goris et al., 2007; Jain et al., 2018). Consistent with the relationships inferred by phylogenetic analyses, ANI evaluation confirmed the identification of two distinct lineages within B. anthina (Figure 3). Pairwise ANI of B. anthina Clade I and Clade II demonstrated 90% pairwise identity, indicating that the two groups represent distinct species according to the operational 95% ANI threshold commonly used for species designation and appropriate for Burkholderia. Genomes clustering with B. anthina Clade I appeared to be genetically variable with further clusters designated I-A, I-B, and I-C identified within the group. Among these clusters, genomes displayed a pairwise ANI between 94.7–95.8%, bordering the threshold for species demarcation. In comparison B. anthina Clade II harboured highly similar strains with a pairwise ANI ≥ 98%. Assessment of species boundaries between B. anthina Clade II and B. cenocepacia confirmed that these groups are distinct to each other, with less than 95% pairwise identity (Supplementary Table 2). Taken together, these analyses corroborate our findings that there are two distinct lineages within B. anthina.

To explore the functional differences of the different B. anthina clades, the pangenome was analysed according to their KEGG ortholog pathways. Both groups demonstrated a similar KEGG ortholog distribution profile regarding the relative proportion of core, accessory, and unique genes assigned to each KEGG pathway. However, B. anthina Clade II, displayed a higher relative proportion of unique genes associated with cellular processes and motility. KEGG categories associated with metabolism, including amino acid and carbohydrate metabolism were the most abundant categories across both B. anthina clades for core, accessory, and unique genes (Figure 4A). Interestingly, B. anthina genomes displayed substantial variance in KEGG orthologs between and within the identified clade groups (Figure 4B). To highlight clade specific adaptations and identify distinct biological interactions, we performed a KEGG ortholog enrichment analysis using a Fisher’s exact test. An examination of enriched KEGG orthologs revealed 722 genes differentially enriched among the two groups, with 451 (62.5%) genes enriched in Clade I and 271 (37.5%) genes enriched in Clade II (Figure 4B; Supplementary Table 3). KEGG orthologs associated with metabolism represented the most prominent KEGG category enriched in both groups. KEGG categories associated with genetic information processing, environmental information processing, and cellular processes were also enriched, albeit at a lower proportion. Genes associated with global overview and map, carbohydrate metabolism, amino acid metabolism, and xenobiotics biodegradation and metabolism showed frequent enrichment in both groups. In terms of genes encoding functions for metabolism, Clade I showed greater functional enrichments linked to metabolic pathways, microbial metabolism in diverse environments, and degradation of aromatic compounds. Of ecological importance, xenobiotics biodegradation and metabolism were largely enriched in Clade I, with gene functions specific to benzoate degradation, xylene degradation, chlorocyclohexane and chlorobenzene degradation, dioxin degradation, toluene degradation, steroid degradation, and styrene degradation.

Analysis of secondary metabolites identified 18 different biosynthetic gene clusters (BGC) in B. anthina with clear differences between Clades I and II (Figure 5). Genes clusters inferred to encode the biosynthetic pathway for terpene production (186), non-ribosomal peptide synthases (NRPS) (72), type I polyketide synthase (T1PKs) (65), and NRP-metallophone (56) were most abundant. BGC counts associated with terpene production, NRPS, and T1PKs were higher in Clade II than that of Clade I. While Clade I demonstrated a broader distribution of specialised metabolite cluster types, with ranthipeptide, betalactone, and butyrolactone exclusively present in Clade I, albeit not detected across all genomes. Aminopolycarboxylic-acid and other were exclusively present in all genomes within Clade II. The aminopolycarboxylic-acid gene cluster is predicted to encode the biosynthesis of ethylenediaminesuccinic acid hydroxyarginine (EDHA), an iron responsive regulated zwitterionic siderophore. To validate the presence of the predicted EDHA gene cluster we bioinformatically identified the operon within the genome, using Streptomyces sp. MA5143a as the reference sequence. As predicted, all Clade II B. anthina genomes harboured homologues for all four genes, allowing for the biosynthesis of EDHA a newly described siderophore (Supplementary Figure 1). To gain insight into the overall iron metabolism potential of B. anthina, the bioinformatic tool FeGenie was used identify genes related to iron acquisition, storage, and redox cycling. Genes associated with iron gene regulation and iron storage were comparable across the two clades, and gene families related with heme oxygenase, iron reduction, or magnetosome formation were absent across all genomes (Figure 6A). Generally, a greater number of gene families associated with iron transport were identified in Clade I (Figure 6B), while gene families associated with siderophore synthesis were higher in B. anthina Clade II (Figure 6C). Genome annotations by RAST predicted on average 31 more genes associated with siderophore synthesis in Clade II compared to Clade I, highlighting the extended siderophore assembly capabilities of the group (Supplementary Table 4). All B. anthina genomes encoded genes required for ornibactin synthesis, while only Clade II genomes encoded genes for pyochelin synthesis, and at least 20 genes putatively involved in siderophore biosynthesis. Interestingly, four phylogenetically closely related strains (B. anthina 1CH1, B. anthina AU37523, B. anthina DSM 16086, and B. anthina LMG 20980) identified in Clade I harboured the genes FhuB, FhuC, and FhuD required for uptake of ferric aerobactin. The synthesis and transport of the aerobactin siderophore is frequently produced by several Enterobacteriaceae, including Escherichia, Vibrio, Salmonella, and Shigella (Waters and Crosa, 1988; Sheldon et al., 2016). Given that aerobactin biosynthesis is not an established feature of Burkholderia it is likely these genes were horizontally acquired. Overall, these results highlight clear genomic and metabolic differences between the two clades identified in B. anthina and provides strong evidence that they should be considered two different species.

Burkholderia anthinoferrum sp. nov. (an.thin.no.ferrum: anthino, meaning ‘from the garden’ or ‘from flowers’, ferrum; iron). Cells are gram negative, non-motile, non-sporulating, coccobacilli bacteria. Growth is observed in LB, MacConkey and BCSA agar at 30, 37, and 42°C. Colonies in LB and BCSA are circular, convex, and cream coloured, while colonies in MacConkey are light pink coloured. Growth in BCSA produces an alkaline medium due the fermentation of lactose and sucrose and is strain dependent. The G + C content for the members of the species is 67.0%. Strains have been isolated predominately from environmental samples, but recently from clinical samples as described in this study. We propose that Burkholderia anthinoferrum sp. nov. XXVI be used as the type strain, a rhizosphere bacterium isolated from a mango orchid in Mexico. The strain holds a central position in the phylogeny of its clade and has been characterised in the literature (de Los Santos-Villalobos et al., 2012, 2018).