Air samples were collected at multiple indoor and outdoor locations in Singapore in 2015 (including residential areas, offices, parks, and beaches) using Andersen single-stage impactors (SKC, USA). The air was impacted onto multiple agar types (M9 minimal salts, R2A, potato dextrose, malt extract, brain heart infusion, and marine agar) which were then incubated overnight at 30 °C. The resulting colonies were sub-cultured on Tryptic Soy Agar (Sigma-Aldrich, USA) and the strains were individually inoculated in lysogeny broth (LB, Becton–Dickinson, USA) at 30 °C overnight to obtain axenic cultures.

The sequencing and assembly of the P. aryabhattai whole genomes were performed using the same pipeline as described in Kalsi et al (Kalsi et al., 2019). Briefly, genomic DNA libraries for each genome were sequenced on a PacBio RS II sequencer, whilst whole-genome shotgun libraries were sequenced at 300bp length using the Illumina MiSeq. De novo draft assemblies were assembled using the PacBio reads and then improved using the MiSeq reads in order to create high-fidelity genome assemblies. The completed assemblies consist of 4 to 14 contigs each (Table 1), with GC content ranging from 37.18 to 38.14%.

To compare the nucleotide similarity and synteny between the recombinant clades and the major megaterium and aryabhattai clades, one genome from each of the four clades (Figure 1) was aligned against the reference genome assemblies, P. aryabhattai K13 and P. megaterium ATCC 14581. The genomes from each clade were selected for having few, well assembled contigs and a high N50. The whole genome alignments were performed using progressive Mauve 2.4.0 (Darling et al., 2010) and the results were visualized with the R package genoPlotR (Guy et al., 2010).

Isolates were freshly revived from –80 °C storage and cultured on trypticase soy agar at 30 °C overnight. After 24 hr, a colony was sub-cultured in trypticase soy broth at 30 °C for 48 hr. Biochemical phenotype tests were then performed using the API 20 E and 50 CH kits (bioMerieux, USA) following the manufacturer’s standard protocol. After incubating for 48 hr, the results were manually observed according to the manufacturer’s standard instruction. API 50 CHB (version 4.1) in APIWEB was used as the reference for species interpretation by matching the isolates’ biochemical profiles.

From the 190 genomes used for the phylogenetic trees, we selected the 18 which were listed on GenBank as assembled at the ‘Complete genome’ level, plus the 12 new genome assemblies that we generated. These 30 genomes included 13 from the megaterium clade, 14 from the aryabhattai clade, and 3 from recombinant clade 1. The genomes were annotated using Prokka (Seemann, 2014) in order to extract their protein-coding sequences. SwiftOrtho (Hu and Friedberg, 2019) was then used to compare the genes from the 30 genomes in an all-vs-all manner, grouping them into sets of orthologs using a reciprocal best-hits method.

The orthologous gene clustering revealed that many of the species-specific orthologs are collocated in species-specific blocks of adjacent genes. Eight such blocks were identified in all P. megaterium genomes, and three were found in all P. aryabhattai genomes (Supplementary Table 9). These species-specific gene blocks are ideal targets for a PCR test to discriminate between the two very similar clades (Figure 2). By using primers that straddle three of the adjacent genes, it can be determined that the PCR reaction correctly amplifies the adjacent gene block that is exclusive to the relevant Priestia clade, rather than any single gene that may have other homologs elsewhere in the genomes of the other clade. PCR primers were designed with Primer3Plus (Untergasser et al., 2012) using sequences from the largest orthologous gene blocks of each species. The target sequences were searched for in all available P. megaterium and P. aryabhattai assemblies in the GenBank database using BLASTN (Altschul et al., 1990) to confirm their specificity. See supplementary information for details of PCR experiments.

In order to identify orthologs that have a high degree of divergence between the two species, we examined orthologs that were present in at least 95% of genomes from both species. dN and dS were calculated for pairs of sequences from the same species and from pairs of sequences from opposite species. Orthologs were excluded if the average dN within one of the species was 0.15 or greater, or if the average dS within one of the species was 0.175 or greater. The gene trees made from orthologs with dN or dS values above these thresholds were unable to recreate the two species clades (Figure 1) due to lack of sequence conservation within species.

The remaining 4,197 orthologs were ranked using the following equation:

Where, for an ortholog with sequences present in both species:

dN(betw) is the average dN of pairs of sequences from opposite species.

dN(mega) is the average dN of pairs of sequences both from P. megaterium.

dN(arya) is the average dN of pairs of sequences both from P. aryabhattai.

We compared the global distributions of P. aryabhattai and P. megaterium in a dataset of 1,171 metagenomics samples from 33 countries (Gusareva et al., no date), with a latitude range of -45.033° to 65.685°. The air sample collection and sequencing methods are described in Gusareva et al. (Gusareva et al., no date). A total of 3.9 billion reads were produced for analysis.

An initial taxonomic classification of the reads was done by DNA-protein alignment to the NCBI NR protein database (downloaded 30/03/2021) using Kaiju v1.7.2, yielding 21–416 reads assigned to P. aryabhattai and 94,663 reads assigned to P. megaterium. Because of the high degree of nucleotide sequence similarity between the two species (Figure 3), we performed a second pass of DNA-DNA alignment to classify reads of the two species. Reads that had been assigned to either of the two species by Kaiju were aligned to both reference genomes (P. aryabhattai K13, NCBI accession GCF_002688605.1, and P. megaterium NBRC 15308 = ATCC 14581, NCBI accession GCF_006094495.1) using BWA v0.7.17. Reads were then assigned to the species that produced the longest alignment, or else the best alignment by percent nucleotide identity if the alignment length was equal between species.

Out of 116,079 reads that were classified as P. aryabhattai or P. megaterium by Kaiju, 43,424 (37%) could not be aligned to either genome by BWA and were discarded from the analysis. A further 5,866 (5%) aligned equally well to both genomes in alignment length and nucleotide identity, and were also removed.

12 new high-fidelity whole-genome assemblies were generated, of isolates from outdoor air samples collected in various locations in Singapore (Table 1). The assemblies have 4 to 14 contigs, with an average total size of 5,762,982 bp and an average N50 of 5,095,237 bp, indicating thorough assembly of the main chromosomes. The genomes were determined to have low contamination and heterogeneity by CheckM 1.0.7 (Parks et al., 2015) (Supplementary Table 2), and completeness was found to be over 98% by BUSCO 5.0.0 (Manni et al., 2021).

Species identification of the 12 genomes was performed using FastANI version 1.32 (Jain et al., 2018) to search for the closest matching genome in the NCBI RefSeq database. Nine out of the 12 assemblies have ANI hits above 98% to database genomes named P. megaterium, and one to P. aryabhattai, whilst two genomes showed ANI matches of just over 96% to P. megaterium. However, the 16S rRNA genes of all the genomes match closest to P. aryabhattai (Table 1). These mismatches between the two methods of ANI and 16S RNA indicate an issue with species identity in the Priestia genus. The 12 genome assemblies in Table 1 were all named Priestia aryabhattai due to the findings of the phylogenetic study below. The two genomes which were the closest GenBank matches to the 12 new genomes — P. megaterium Q3 and YC4-R4 — are part of a separate P. aryabhattai clade and are recommended to be renamed as such (Supplementary Table 7).

To determine the phylogenetic relationship between the two species P. aryabhattai and P. megaterium, we conducted a thorough phylogenetic study of all 190 available genomes of the two species, comprising the 12 genomes we generated in this study (Table 1) in addition to 178 whole-genome assemblies that were downloaded from GenBank. We constructed phylogenetic trees based on protein-coding gene sequences (Figure 1a), whole-genome ANI (Figures 1b, c), and 16S rRNA (Figure 1d).

The trees (Figures 1a-c) separate most of the genomes (165 genomes) into two consistent clades, which we refer to as the megaterium clade (130 genomes) and the aryabhattai clade (35 genomes). The classification of the two clades is supported by the protein-coding gene trees and the ANI trees, with the members of each clade being identical by each method. A phylogenetic tree constructed using the 16S rRNA sequence, a commonly used bacterial marker gene, was unable to distinguish the two clades due to high sequence similarity within the gene (Figure 1d).

Out of the 190 genomes, 24 were not placed in the megaterium nor aryabhattai clade but instead formed two smaller clades which showed inconsistent positions between trees. For example, one of these clades was grouped with the megaterium clade in the protein-coding gene tree but with the aryabhattai clade in the whole-genome ANI tree (Figures 1a, b). The network tree (Figure 1c) placed the two small clades in an intermediate position between the megaterium and aryabhattai clades. We defined these two smaller groups as recombinant clades 1 and 2 (14 and 10 genomes each). Although recombinant clade 1 is not monophyletic in the protein-coding gene tree (Figure 1a), the clade is monophyletic in the ANI and network trees (Figures 1b, c). This clade includes P. aryabhattai B8W22, the first described isolate of P. aryabhattai (Shivaji et al., 2009).

Most of the genome assemblies retrieved from GenBank have already been assigned the species name that corresponds to their clade in our trees. We found 17 genomes whose assigned species names do not match with our classifications based on this phylogenetic analysis, and the suggested species identification for these genomes are provided in Supplementary Table 7.

We calculated the average distances within and between the four clades by whole-genome ANI and by nucleotide p-distance of the 74 genes (Supplementary Table 6). The megaterium clade had greater diversity (2.67% +/- 1.23 ANI distance) than the aryabhattai clade (1.37% +/- 0.39 ANI distance), suggesting a more recent common ancestor for aryabhattai than for megaterium. The whole-genome ANI is 1.28 times larger than the p-distance of protein-coding genes (0.78% in the megaterium clade), suggesting that ANI is more effective in showing the extent of genetic diversity at the strain level than the p-distance of the protein-coding genes, which are house-keeping genes with conserved sequences between genomes.

Additionally, TYGS (Meier-Kolthoff and Göker, 2019) was used to further investigate the taxonomy of 3 strains from the two recombinant clades. No closer matches than P. aryabhattai and P. megaterium were identified (Supplementary Table 5).

We aimed to determine whether the whole-genome genetic distance between the megaterium and aryabhattai clades is great enough to consider them as separate species. Using the dDDH (Auch et al., 2010) and ANI (Jain et al., 2018) methods, genomes are conventionally defined as belonging to the same species if they are at least 70% or 95% similar, respectively (Klappenbach et al., 2007).

The two metrics of dDDH and ANI showed inconsistent results for species identification (Figure 3). The pairwise dDDH values for all comparisons between the megaterium clade and the aryabhattai clade (orange & blue points, lower left) were below 70% similarity, with an average of 64%, indicating that they are separate species. However, most of the pairwise ANI values were within the 95–96% threshold, with an average of 95.2%. The majority of genome pairs (174 out of 182) were above 95% ANI, suggesting that P. megaterium and P. aryabhattai are the same species. Thus, the two distance measurements cannot consistently separate the two clades as different species because the two clades are closely related, and because the intersection of the dDDH and ANI thresholds is not placed along the lined formed by the genome pairs.

The pairwise values of both dDDH and ANI (Figure 3) between the main megaterium/aryabhattai clades and the recombinant clades were higher than the values between the megaterium and aryabhattai clades, but lower than the values within the megaterium clade or within the aryabhattai clade. The dDDH values between the recombinant clades and the other clades were lower than the species identification threshold, while the ANI values were higher than the threshold. The majority of the recombinant genomes (19 out of 20) are isolates from outdoor environments; the existence of naturally occurring recombinants with intermediate genomic identity between P. megaterium and P. aryabhattai suggests a recent split or ongoing speciation between the two clades.

To assess the level of synteny between clades and the similarity of the recombinant clades to the aryabhattai and megaterium clades, a genome from each of the four clades was aligned against the reference genomes of P. aryabhattai and P. megaterium. The results (Figure 4) showed high identity between genomes within the megaterium clade and within the aryabhattai clade, as expected, and also a high level of synteny between these two clades. The recombinant clades both showed a loss of this synteny, therefore their intermediate similarity to the two larger clades cannot be explained as them being the common ancestor of megaterium and aryabhattai. These clades also showed a mosaic pattern of genomic segments with higher identity to both the megaterium and aryabhattai reference genomes, suggesting a history of extensive horizontal gene transfer between these groups.

The results of 61 biochemical phenotype tests on 13 Priestia isolates (1 megaterium, 10 aryabhattai, 2 in recombinant clade 1) are shown in Supplementary Table 4. No single test or combination of tests were found that could distinguish between the aryabhattai and megaterium clades that were shown in Figure 1. Out of 61 tests, 53 tests gave identical results for all 13 isolates, while 7 tests were inconsistent within the 10 aryabhattai isolates. The remaining one test showed a different result only in one of the recombinant isolates. APIWEB successfully identified all 13 strains as belonging to the larger Priestia megaterium group and was not designed to differentiate between the two species P. megaterium and P. aryabhattai.

To facilitate the identification and research of strains from these species, we aimed to design PCR primers that can unambiguously place any new strain into one of the two species, without the need for whole-genome sequencing. To create such primers, we compared the protein coding gene sequences of annotated genome assemblies and identified the gene regions which are specific to each of the aryabhattai and megaterium clades.

We used SwiftOrtho (Hu and Friedberg, 2019) to find clade-specific orthologs in 30 high-quality assemblies (listed as ‘Complete genome’ in GenBank). The analysis identified 46 orthologs which were present in every genome in the megaterium clade but in none of the genomes in the aryabhattai clade. Conversely, 21 orthologs were found in every genome in the aryabhattai clade but in none of the genomes in the megaterium clade. About half of these species-specific genes are located adjacent to each other in the genome in several blocks of genes up to 10kb in length (Supplementary Table 9). The consistent presence and absence of these gene blocks in the megaterium and aryabhattai clades was confirmed by BLAST searching for each gene block in all 190 Priestia genomes that were used in the phylogenetic study. Genomes in the recombinant clades often contained a combination of gene blocks from both the megaterium clade and the aryabhattai clade, with no clear pattern.

We identified 3 blocks of adjacent genes which were specific to the aryabhattai clade, and 8 blocks of adjacent genes that were only found in the megaterium clade. By placing forward and reverse PCR primers across different genes in these genomic regions, it can be ensured that these consecutive clade-specific genes are amplified, rather than single homologs elsewhere in the genome. Here we provide two sets of PCR primers (Supplementary Table 10); one for identifying P. megaterium strains, and another for identifying P. aryabhattai strains. Researchers can identify which clade their strain belongs to by running two PCRs, once with each pair of primers, since strains of one clade do not react to the other clade’s primers. Figure 2 shows the species-specific orthologs which are targeted by the PCR primers. The gel electrophoresis image shows the successful identification of 4 megaterium and 10 aryabhattai strains (3 strains from culture collections and 11 of the new isolates) using these primers.

In order to identify the orthologs with the greatest sequence divergence between the two clades, we calculated the average pairwise dN between the two clades minus the average pairwise dN within each of the two clades. This method identifies orthologs that have highly conserved sequences within a clade but have also diverged in sequence between clades. The top 100 identified orthologs by this method are given in Supplementary Table 12 with the gene names identified by Prokka and UniProt. A large number of the top 100 are unidentified proteins with less than 90% identity to their closest database match. The distribution of the sequence divergence for 4,197 orthologs is shown in Supplementary Figure 3.

The top 100 orthologs by sequence divergence between the species includes several which are related to the synthesis and use of cobalamin a.k.a. vitamin B12 (Supplementary Figure 3). These include the cobS, cobD, and cobU genes from the cobalamin synthesis pathway and a vitamin B12-dependent ribonucleotide reductase. Other orthologs that have diverged between the species include seven GNAT family acetyltransferases, three components of the PTS sugar transport system, three genes related to spore formation and germination, two prophages, two Shikimate kinases, and three genes related to flagellum formation. These functional groups each had orthologs showing sequence divergence between the two clades of more than two standard deviations above the mean of the 4,197 total orthologs.

To search for broad patterns of differentiation in geographical or climatological adaptation between the two species, we identified metagenomic reads belonging to P. aryabhattai and P. megaterium in global air samples and compared the relative abundances of the two species over different latitudes. The results in Figure 5 show a significant difference in the most abundant species at latitudes greater than or less than 16° from the equator (Chi-squared test, p = 1.48e-51). In samples taken at locations further than 16° from the equator, P. megaterium was usually the more abundant species, whereas at locations less than 16°from the equator, P. aryabhattai was usually the more abundant species.