Bacterial strains were grown at 28 °C using a Tryptone Yeast Extract (TY) medium containing 5 g Tryptone and 3 g yeast extract per liter of water (Beringer, 1974). Pikovskaya (PVK) medium was used to visualize phosphate solubilization. This contained (per liter): 10 g glucose, 5 g Ca₃(PO₄)_2–5 g (NH₄)₂SO₄, 0.5 g NaCl, 0.2 g MgSO₄, 0.2 g KCl, 0.2 g yeast extract, 0.5 g MnSO₄, and 0.002 g Fe SO₄ (Nautiyal, 1999). To isolate phosphate-solubilizing bacteria from soil samples, 5 g of soil was mixed with 50 ml of sterile water, diluted, and plated onto PVK agar (Nautiyal, 1999). Colonies forming larger zones of clearing were picked and single colony purified three times.

For long-term storage, purified cultures were grown overnight in TY to late log phase and mixed with a freezing solution (TY containing 24% vol/vol DMSO) such that the final concentration of DMSO was 8%. These were stored at -80 °C.

Soil samples were collected from the University of Manitoba Ian N. Morrison Research Farm near Carman, Manitoba (49.492261 N, 98.042497 W) in 2018. The plots sampled were part of a rotation study focused on incorporating soybean into a rotation with canola, corn, and wheat. In 2018, all four crops were grown. Soils were collected as 2 x 9 cm cores taken immediately adjacent to plant roots at crop maturity. Samples were stored at 4 °C, with subsamples frozen at -80 °C the same day.

To extract bacterial genomic DNA from isolated strains, colonies were picked with a sterile onoculating stick from an agar plate and resuspended in 0.1 ml ddH₂O. From this, 50 μL was extracted using Qiagen’s DNeasy PowerSoil kit using the manufacturer’s protocol. Initial identification of bacterial isolates was done by sequencing the V4 region of the 16S rRNA. Briefly, the V4 region was amplified using the 515F; 5’-GTGCCAGCMGCCGCGG-3’ (Parada et al., 2016) – 806R; 5’-GGACTACNVGGGTWTCTAAT-3’ (Apprill et al., 2015) primer pair as previously described (Yudistira et al., 2021). The isolated PCR products sent to Centre d’expertise et de services Génome Québec for sequencing. Initial taxonomy was assigned to sequences using BLAST (2.11.0) against the NCBI nucleotide database.

To quantify phosphate solubilization by bacterial cultures, a QuantiChrom™ Phosphate Assay Kit was utilized. Briefly, three independent colonies were grown overnight in 5 mL TY broth. The optical densities (OD600) of the cultures were normalized to 1.0 and 1 mL was added to 5 mL of PVK broth and incubated at 28 °C for three days. The cultures were then pelleted for one minute at top speed in a microfuge, and the supernatants were assayed. The phosphate solubilization assays were carried out in a 150 µL final volume in a 96-well clear bottom plate that included 50 µL of sample or standard and 100 µL of the provided reagent which contained Maclchite Green and molybdate for color development. The contents were mixed and incubated in the dark for 30 minutes at room temperature. Color development was assayed at 620 nm using a plate reader. The phosphate concentration was calculated as follows:

ODBlank, ODStandard, and ODSample are the OD620nm values of the 0 Pi Blank, the 0.28 mg dL-1 Pi Standard, and the sample, respectively. Conversions: 1 mg/dL Pi equals 105.3 µM, 0.001%, or 10 ppm.

Plant growth-promoting activity was determined using a vigor index (Abdul‐Baki and Anderson, 1973). Briefly, this index measures growth promotion using the mean root length, mean shoot length, and seed germination.

To carry out these assays, seeds were first surface sterilized for 20 minutes using a 1% hypochlorite solution. Subsequently, they were washed with approximately 10 volumes of sterile ddH₂O over a period of approximately 30 minutes. Bacteria being tested were first grown in 5 mL broth culture of TY and grown overnight. The overnight culture was then used to inoculate a 100 mL broth culture of TY. The OD600 of these cultures was adjusted to an optical density of 1. A volume of 20 ml of the adjusted cultures was used to soak the surface sterilized seeds. Seeds were treated with 20 mL of the bacterial culture for 30 minutes. Control seeds were soaked in 20 mL of sterilized water ddH₂O. Finally, the seeds were individually planted using sterile forceps into seed-starting trays that contained a sterile sand/vermiculite (1:1) mixture. These were subsequently covered with plastic cling film to reduce the drying of the seedlings. A total of 24 seeds were planted for each treatment. Seedlings were harvested 15 days after planting, and root and shoot lengths were measured. Significance was determined using a two-way ANOVA using Dunnett’s test for post hoc analysis.

Whole genome sequencing was carried out essentially as previously described using a MinION sequencer (Hawkins et al., 2022). To extract genomic DNA, bacteria were grown overnight in TY broth, and genomic DNA was extracted using the PureLink Genomic DNA mini kit (Invitrogen). Library preparation was performed using SQK-LSK-109 and EXP-NBD-104 kits following the manufacturer’s instructions (Oxford Nanopore Technologies, Oxford, UK). Sequencing was carried out using a Nanopore MinION Mk1B system with R10.3 flow cells. Sequencing was stopped once sufficient data were obtained to achieve approximately 100x coverage across all genomes. Sequencing reads were then base-called using Guppy-GPU (Wick et al., 2019). The generated FASTQ files were subjected to quality control and adaptor trimming. Adaptor trimming was performed using BBduk (Bushnell et al., 2017). de novo genome assembly was achieved using Flye, followed by three rounds of polishing using minimap2 (Kolmogorov et al., 2019; Li, 2018). Genome completion was estimated using CheckM (Parks et al., 2015) and found to be above 99%, with an estimated 1% contamination in all assemblies. Default parameters were used for all software in the analysis.

Phylogenetic analysis was carried out using programs and packages available on the Type Strain Genome Server (TYGS), found at https://tygs.dsmz.de/background/show. Assembled genomes were submitted to the Type Strain Genome Server (TYGS) to obtain DNA-DNA Hybridization (dDDH) values as well as G+C content (Meier-Kolthoff et al., 2022). Phylogenetic trees were also generated on the TYGS server using the Genome BLAST Distance Phylogeny (GBDP) approach, with the ‘coverage’ algorithm and distance formula d5.n (Meier-Kolthoff et al., 2013). Average Nucleotide Identity (ANI) values were calculated using fastANI (Jain et al., 2018). ANIclustermap (Shimoyama, 2022), found at https://github.com/moshi4/ANIclustermap) was used to visualize ANI data.

The WK1.1f assembly was scaffolded using RagTag (v2.1.0). The scaffold module was utilized with SI02 designated as the reference genome. The default minimap2 module was used to align the contigs from WK1.1f. Unplaced contigs were combined into a single scaffold labeled “ch0” with the -C option. All other parameters were left at their default settings. To compare the genomes, BRIG (v0.95) was employed. The concatenated WK1.1f fasta sequence formed the backbone for the alignment. GC content and GC skew rings were added, followed by a BLASTn alignment of the SI02 and SI08 genomes. The BLASTn thresholds were set at 50% and 70% identity. Symbiosis genes were annotated separately, following the layout of the WK1.1f backbone.

Genome synteny was assessed using SYNY (Julian and Pombert, 2024). Assembled genomes were annotated with Prokka (v1.14.6) using default settings, and the resulting GenBank files were used as input for the analysis. All isolates were subjected to pairwise comparisons using default settings for modification parameters. Synteny plots were generated based on collinearity inference using DIAMOND (Buchfink et al., 2021).

Pangenome analysis was conducted using Roary (v3.13). Genome assemblies were annotated with Prokka (v1.14.6). The GFF3 output files served as input for Roary. The pangenome was calculated with default settings, using the -e and -n options. The resulting presence-absence file was used to visualize the pangenome with an UpSet plot generated by the R package Complexupset version 1.33.

The genomes of SI01, SI02, SI04, SI08, and SI09 were aligned using Mauve within Geneious. Each aligned block was manually inspected for significant gaps in homology. Protein-coding regions unique to a cluster and containing more than 10 genes were extracted and concatenated. Each replicon was then subjected to cluster of orthologous groups (COGs) classification using eggNOG-mapper (v2.1.12), with DIAMOND selected for annotation. The total counts for each COG category per replicon were used to generate a heatmap in R.

Pathway completeness was assessed using eggNOG-mapper (v2.1.12). Protein-coding FASTA files generated by Prokka (v1.14.6) were used as input. The DIAMOND BLASTx option was selected to annotate coding sequences. To group and visualize metabolic pathways, KEGGaNOG (v0.5.6) was used. The annotation files produced by eggnog-mapper served as input to KEGGaNOG, using all default options.

Bacteria that have phosphate-solubilizing capabilities often exhibit plant growth-promoting traits (Cheng et al., 2023; Shen et al., 2021). Many of the phosphate-solubilizing bacteria previously isolated have been isolated from the rhizosphere of several crop plants. These tended to be more metabolically active than those isolated from non-rhizosphere sources (Gyaneshwar et al., 2002). We collected soil samples adjacent to the roots from a rotation study focused on incorporating soybeans into rotations with canola, corn, and wheat.

After soils from each crop were initially plated onto PVK agar, 20 putative phosphate-solubilizing bacteria were found. Following single-colony purification and retesting for phosphate-solubilizing activity, nine strains consistently produced larger zones of clearing and were selected for further analysis. These were labeled as SI01-SI09. Based on the sequencing results of the V4 regions of the isolates, as well as the sample/plate from which each of the strains was initially isolated, and the growth characteristics, it was deemed likely that SI02 and SI03 were likely siblings, so SI03 was not further characterized. This left 8 isolates capable of solubilizing phosphate on PVK agar (Figures 1A, B). SI01 and SI09 were isolated from the corn plot, SI02-SI05 and SI07 from the soybean plot, and SI06 and SI08 from the canola and wheat plots, respectively. Based on the V4 sequences, except for SI05 and SI06, the strains were tentatively identified as Paraburkholderia strydomiana. SI05 and SI06 were tentatively identified as Burkholderia ambifaria and Paraburkholderia graminis, respectively.

To quantify their ability, these strains were grown overnight in TY, subcultured into PVK broth for three days, and the supernatants were assayed for soluble P_i (Figure 1C). The results show that in all cases, the isolates could solubilize between 95–144 mg dL^-1 phosphate, in contrast to the 10 mg dL^-1 normally found in PVK medium. The quantification also corroborated the visual observation that each of the strains looked similar on PVK agar.

Seed germination and seedling emergence are crucial in plant development and linked to better later health (Gardarin et al., 2016). To assess if the strains promoted growth, they were tested on soybean seeds.

When soybean seeds were inoculated with the phosphate-solubilizing strains, it appeared that seeds treated with SI02, SI04, SI07, and SI09 looked larger and greener than the control plants that were treated with only water (Figure 2). Seeds treated with SI05 and SI06 looked smaller and were negatively impacted compared to the control plants (Figure 2). When these observations were used to calculate a vigor index, it was found that treating soybean seeds with these strains resulted in significant differences (Figure 3). Seeds treated with SI02, SI04, SI07, and SI09 increased the vigor index, whereas SI05 and SI06 decreased the vigor index (Figure 3C). Looking through the individual measures that were used to determine the vigor index, it was noted that changes that were seen in the soybean seedlings were related to the percentage of seeds that germinated as well as the differences in the root and shoot lengths (Figures 2, 3). We note that all the strains that affected the germination rate of soybean were P. strydomiana (Figure 3A). Due to inconsistent growth, strains SI01 and SI08 were excluded from these assays.

P. strydomiana strains have previously been isolated from nitrogen-fixing nodules taken from Hypocalyptus sophoroides in South Africa (Beukes et al., 2019). Additionally, Paraburkholderia is strongly associated with plant growth-promoting traits. Since the isolates in this study were not derived from nodules but were selected based on their ability to solubilize phosphate, sequencing these isolates was of interest to see if their genome sequences might provide insights into their capabilities.

Whole genome sequencing using the MinION yielded greater than 35x sequencing depth per sample (Table 1). As expected, based on preliminary sequencing of the V4 region, six of the eight isolates belonged to the species P. strydomiana (SI01, SI02, SI04, SI07, SI08, and SI09), and the other two were P. graminis (SI06) and B. ambifaria (SI05).

The genome assembly of B. ambifaria showed 5 contigs with a total size of 8.076 Mbp, an N50 of 3.02 Mbp, a GC content of 66.3%, and coverage of 75x. Using the Rapid Annotation Subsystem Technology (RAST), this genome was predicted to contain 7130 protein-coding genes.

The genome of P. graminis is 7.243 Mbp with 3 contigs, an N50 of 4.256 Mbp, 60.4% GC content, and 50x coverage, comprising 6385 protein-coding and 69 tRNA genes. P. strydomiana assemblies mostly had 3 contigs, except SI07 with 9. (Table 1). If other assemblies were used as scaffolds, SI07 could also be resolved to have 3 contigs. Their sizes ranged 7.2–7.791 Mbp, N50 3.9–4.33 Mbp, GC 58.2–61.8%, coverage 35x–118x, with 6915–7226 protein-coding and 65–68 tRNA genes.

Using assembled genomes, we calculated pairwise ANI and dDDH values for each isolate and its closest relatives with the Type Strain Genome Server (TYGS) (Table 2). The nearest related species for each query is shown in Table 2. Each strain’s ANI compared to its match was over 95%, and its dDDH exceeded 70%, which are the currently accepted cut-offs for species identity, when compared to the type strain.

Using the genome sequences from the three identified subject strains (Table 2), an all-against-all ANI analysis was conducted. As expected, isolates SI05 and SI06, which were identified as P. ambifaria and B. graminis, respectively, had no relation to the any of the other sequences in the comparison (Figure 4). The isolates identified as P. strydomiana were divided into three groups with 100% identity (Figure 4). These groups consisted of SI01, SI108, and SI109, SI02 and SI04, as well as SI07. We note that while all of these had greater than 95% identity with Wk1.1fT, these groups appeared to be distinct based on the analysis (Figure 4).

To corroborate the ANI analysis, a phylogenetic analysis was performed using the Genome Blast Distance Phylogeny method on the TYGS server. The results show that the resultant tree supported the ANI analysis (Figure 5). Collectively, the data support the hypothesis that these isolates appear to represent three distinct strains, which are in turn distinct from the type strain Wk1.1f (Figure 5).

The genome of strain WK1.1fT, with over 120X coverage, was assembled into 259 contigs (Beukes et al., 2019). Other P. strydomiana genomes in NCBI (https://www.ncbi.nlm.nih.gov/datasets/genome/?taxon=1245417) also vary from 38 to 259 contigs, showing an incomplete understanding of their genome architecture. We aimed to see if WK1.1fT shared the same genome structure as our isolates.

Since WK1.1fT shared the greatest similarity with SI02, the WK1.1fT contigs were aligned to SI02 using the minimapp2 option in RagTAG (Alonge et al., 2022), and the existing assembly was scaffolded. Contigs that were not placed were concatenated into a single scaffold. Utilizing BRIG (Alikhan et al., 2011), the resultant contigs were then subsequently compared with SI02, SI07, and SI08 (Figure 6). From the results, it is clear that the genomes show high similarity (Figure 6). Aside from some minor differences, the WK1.f1 genome was similar to our sequenced isolates across the two largest contigs.

Additionally, 127 sequences that did not align well from WK1.f1 were placed in the concatenated scaffold showed identity with the smallest replicon from SI02 and SI08. This region exhibits multiple regions of low or no similarity to the other isolates and contains numerous insertion elements and recombinases, suggesting it may be highly dynamic. Notably, SI07, which we also could not close, aligned well with WK1.1f and exhibited similar gaps to the other isolates. Generally, based on our analysis P. strydomiana genomes appear to consist of 3 replicons, Contig 1, which is between 3.8 - 4.3 Mb, Contig 2 which is between 3.2 - 3.3 Mb, and Contig 3 which shows the greatest variability with strains S102 and S104 having a replicon of approximately 0.24 Mb, whereas S101, S108, and S109 had a replicon of approximately 0.5 Mb.

Using SYNY to visualize genome collinearity (Julian and Pombert, 2024), the analysis revealed high concordance within P. strydomiana groups (SI02 vs. SI04, and SI08 vs SI09) across all three replicons (Supplementary Figures 1A, B, respectively). When comparing SI02 to SI08, there is substantial synteny between the two largest replicons, and very limited synteny with the smallest replicon (Supplementary Figure 1C). When SI02 and SI08 are compared to SI07, there is significant collinearity across the two largest replicons; however, qualitatively, each has unique regions (Supplementary Figures 1D, E). Similar results were seen when comparing these strains with WK1.1fT (Supplementary Figures 1F–H). Collectively, these analyses corroborate the average nucleotide identity as well as the phylogenetic analysis (Figures 4, 5) and suggest the major differences in the genome make-up are due to the variability of the smallest replicon.

Although the two largest replicons were largely syntenic, it was of interest to determine if these regions of variability affected the organism’s overall physiology. To determine the functionality in these variable regions, the genomes of SI01, SI02, SI04, SI08, and SI09 were first aligned with Mauve. Since SI07 was not a closed genome, it was excluded from the analysis. Each aligned block from the larger replicons was manually checked for gaps in homology. Protein-coding regions unique to each cluster (SI01, SI08, SI09, or SI02 and SI04), with more than 10 genes, were extracted and concatenated. Due to the dissimilarity of the smallest replicons, the entire replicon was considered. The concatenated variable regions were then classified using eggNOG-mapper (Cantalapiedra et al., 2021).

Comparison of the gaps in the largest replicon revealed that both strain groups mainly contain genes in COG groups S, P, M, and V. Group M covers cell wall/membrane/envelope biogenesis, Group P involves inorganic ion transport and metabolism, Group V includes defense mechanisms, and Group S has genes for unknown proteins (Figure 7).

Comparing the variable regions of replicon two shows that both strain groups contain COG groups C, G, K, and L. Group L includes genes for replication, recombination, and repair, while Group K involves proteins related to transcription. Group G encodes carbohydrate transport and metabolism, and group C includes genes for energy production and conversion (Figure 7).

Comparing the plasmids, the highest proportion of shared genes are from Group S, which encodes proteins of unknown function. Plasmids from SI01, SI08, and SI09 have more recognized COG groups than SI02 and SI04. Notably, these plasmids are nearly twice as large as those in other strains (Figure 7).

Genome structure analysis indicated some variability among strains. To evaluate this, a pangenome analysis was performed (Figure 8). The pangenome includes 12,123 genes, with the core genome consisting of 5,195 genes. The type strain, WK1.1f, has 2551 additional genes not present in the other isolates. The number of unique shared genes aligns with hierarchical clustering from ANI analysis and phylogeny (Figures 4 and 5). Strains SI01, SI08, and SI09, which are closely related phylogenetically, share 920 unique genes; SI02 and SI04 share 723; SI07 has an additional 598 unique genes. Aside from SI01, each strain also contains unique genes: SI02 and SI04 have 14 and 12, respectively, while SI08 and SI09 have 69 and 73, respectively. Notably, the isolates share 267 unique genes not found in the type strain.

Using annotated genomes, P. strydomiana was analyzed with KEGGaNOG to visualize pathway completeness (Popov et al., 2025). Core metabolic pathways, including glycolysis, pentose phosphate pathway, Entner-Doudoroff pathway, TCA cycle, glyoxylate shunt, and gluconeogenesis, are present in all strains (Figure 9). All strains contain RuBisCO and genes for the Calvin-Benson-Bassham cycle (Figure 9). In addition, all amino acid biosynthetic pathways are present (Figure 9). The presence of complete cytochrome c, cytochrome bd, and electron transport oxidoreductases suggests an aerobic lifestyle. All strains also have many genes for Type II, III, and VI secretion systems (Figure 9). Regarding phosphate use, all strains have phosphate transport and bidirectional polyphosphate pathways (Figure 9).

There are notable differences between the type strain and the strains isolated in this study. The type strain WK1.1f contains a phosphonate transporter and genes for breaking carbon-phosphate bonds. It also possesses genes necessary for nitrogen fixation and has been demonstrated to be capable of symbiotic nitrogen fixation (Beukes et al., 2019). We note that the majority of the genes associated with nitrogen fixation are located on the third replicon (Figure 6). Additionally, the annotated genome predicts that strains SI01, SI02, SI04, SI07, SI08, and SI09 are capable of performing mixed acid fermentation to produce lactate.

Although many plant growth-promoting bacteria have been reported (Abou Jaoudé et al., 2024; Khoso et al., 2024), very little is known about the molecular basis for their ability to promote plant growth. Several genes have been correlated with plant growth promotion, including pyrroloquinoline quinone encoding genes (pqq), nifHDK, which encodes nitrogenase, acdS, which encodes 1-aminocyclopropane-1-carboxylic acid (ACC) deaminase, hcnABC, which encode enzymes that produce hydrogen cyanide; phlABCD, which encode genes necessary for the biosynthesis of 2,4-diacetylphloroglucinol, and budABC, which has been shown to affect stomatal activity (Bruto et al., 2014). To determine if any of these might be present, the genomes of each strain were queried for their presence. The results show that all the strains contained pqqBCDE, hcnABC, and acdS. The strains appeared to contain two genes necessary for the synthesis of acetoin/2,3-butanediol but were missing the initial step of the pathway. The strains did not contain genes indicative of nitrogenase or 2,4-diacetylphloroglucinol (Table 3).