All strains used in this study except AB011P were described previously as Corallococcus members (Livingstone et al., 2017). AB011P was isolated from soil sampled from farmland at Penrhyncoch near Aberystwyth, using an Escherichia coli prey baiting method. Soil samples were placed on nutrient free medium next to a spot of prey organism. Swarming predators which consumed the prey were picked and isolated onto fresh growth media. 16S rRNA sequencing of AB011P and phylogenetic analysis, according to the method described by Livingstone et al. (2017) confirmed AB011P belonged to the genus Corallococcus. All bacterial strains (including AB011P) were cultivated as described previously (Livingstone et al., 2017) and their predatory activity was assayed against a panel of prey organisms. Myxobacteria were spotted onto lawns of each prey and the resulting zone of prey killing measured, exactly as described by Livingstone et al. (2017), except killing zones were measured after 7 days instead of 4 days in order to increase the sensitivity of the assay.

Strains of Corallococcus spp. were sequenced using 2 × 250 bp paired-end reads on an Illumina Hiseq 2500 platform by either the Centre for Genomic Research (CGR Liverpool, United Kingdom – strains AB011P and AB038B) or MicrobesNG (Birmingham, United Kingdom – another 21 strains). MicrobesNG were provided with cell masses, while DNA was extracted from AB011P and AB038B using the DNeasy Blood and Tissue kit (Qiagen) and quantified using a nanodrop spectrophotometer before provision to the CGR. Kraken 2 was used to identify the closest reference genomes for read mapping (Wood and Salzberg, 2014). BWA-MEM was used to check the quality of the reads (Li and Durbin, 2009) while de novo assembly was performed using SPAdes 3.7 (Bankevich et al., 2012). Annotation was undertaken using Prokka 1.1 (Seemann, 2014), and RAST 2.0 (Aziz et al., 2008) and contigs with coverage less than 7-fold were removed. The genome sequence of Corallococcus coralloides DSM2259^T was retrieved from the NCBI database. Genomes sequenced in this study are listed in Table 1 and available from the NCBI nucleotides database under accession numbers RAWS00000000 (AB004), RAVX00000000 (AB011P), RAWR00000000 (AB018), RAWQ00000000 (AB030), RAWP00000000 (AB032C), RAWO00000000 (AB038B), RAVW00000000 (AB043A), RAWN00000000 (AB045), RAWM00000000 (AB047A), RAWL00000000 (AB049A), RAWK00000000 (AB050A), RAWJ00000000 (AB050B), RAWI00000000 (CA031B), RAWH00000000 (CA031C), RAWG00000000 (CA040B), RAWF00000000 (CA041A), RAWE00000000 (CA043D), RAWD00000000 (CA047B), RAWC00000000 (CA049B), RAWB00000000 (CA051B), RAWA00000000 (CA053C), RAVZ00000000 (CA054A), and RAVY00000000 (CA054B).

Complete 16S rRNA gene sequences were obtained from sequenced genomes and used to query a BLASTn search of the NCBI database to identify their closest relatives. Query sequences along with those from relatives with >92% similarity and >82% nucleotide coverage were aligned and a neighbor-joining tree was constructed in MEGA-7.0 (Kumar et al., 2016), using the Kimura 2-parameter model with 500 bootstrap replicates. AMPHORA2 was used to predict phylogenetic relationships based on 31 concatenated and aligned marker genes (Kerepesi et al., 2014), and a maximum likelihood tree was generated using a Jones-Taylor-Thornton model on MEGA 7.0 (Kumar et al., 2016). We also used PhyloPhlAn to provide a more highly resolved phylogeny, based on 400 conserved genes (Segata et al., 2013). Newick format tree files were uploaded onto iTOL, a web based tool for annotating and editing trees (Letunic and Bork, 2016).

The average nucleotide identity (ANI) is a widely accepted genome-based method for species delineation and an ANI-based all-vs-all matrix and resulting clustering tree were constructed using the ANI-Matrix genome-based distance matrix calculator (Rodriguez and Konstantinidis, 2016). The digital-DNA/DNA hybridization (DDH) was calculated for pairs of genomes using the genome-to-genome distance calculator (GGDC 2.1) (Meier-Kolthoff et al., 2013). Results from GGDC Formula 2 are presented, as Formula 2 is the most robust Formula when applied to incomplete genomes (Auch et al., 2010).

Orthogroups (sets of orthologs) were compiled and trees based on orthogroup member presence/absence were constructed both using OrthoFinder 1.1.8 (Emms and Kelly, 2015). Pan-genome analysis was carried out using Roary (Page et al., 2015). Phandango (Hadfield et al., 2017) and recommended R scripts were used to view the resulting output graphs. Using RAxML (Stamatakis et al., 2004) core genome gene sequences were aligned and maximum likelihood trees produced for viewing and editing using iTOL (Letunic and Bork, 2016).

SPINE (Ozer et al., 2014) was used to separate core and accessory gene sequences. Initially genomes were aligned using the NUCmer function of the integrated MUMmer software package to identify core genes. Core gene sequences were then mapped to individual genomes using AGEnt (Ozer et al., 2014), allowing identification of the accessory (non-core) genes sequences in each genome. Accessory gene sequences were then clustered using ClustAGE (Ozer et al., 2014). Core and accessory gene sequences were then functionally annotated using eggNOG mapper and antiSMASH (Weber et al., 2015; Huerta-Cepas et al., 2017).

Only two “genomes” assigned to Corallococcus spp. were publically available at the time of this study, including the complete 10.08 Mbp genome of the type strain C. coralloides DSM 2259^T (Huntley et al., 2012). The second “genome” was a 1.56 Mbp contig assembled from metagenomic data, which we excluded from further comparisons. In a previous study, we isolated a large number of strains, identified as Corallococcus spp., and characterized their predatory activity (Livingstone et al., 2017). Twenty two of those Corallococcus strains, with diverse predatory profiles, and a novel Corallococcus spp. isolate (AB011P), with broad and potent antimicrobial activity, were chosen for genome sequencing. All strains grew well on VY/2 Agar and were sequenced variously by the Centre for Genomic Research (Liverpool, United Kingdom – strains AB011P and AB038B) and MicrobesNG (Birmingham, United Kingdom – 21 strains), using the Illumina Hi-Seq platform.

The number of contigs in each draft genome ranged from 104 to 1,546, with an average of 716 per genome. The L50 values (the largest X contigs which together constituted more than half of the genome sequence) had a median of 93, with an average N50 value (the size of the Xth contig) of 60,472 bp (Table 1). Mean sequence coverage was greater than 25-fold for all genomes (averaging 65-fold coverage), with a mean %GC content of 69.9. Assembled draft genomes were consistently large, having an average size of 10.22 Mbp, with the smallest genome (strain AB050B) being 9.40 Mbp and the largest (strain CA043D) being 10.79 Mbp (Table 1). Consequently, the number of coding sequences (CDS) also did not vary widely, with an average of 8,158, ranging between 7,638 and 8,724. Genome sequence data for this project have been deposited in the NCBI nucleotide database.

A phylogenetic tree (Supplementary Figure S1) was constructed using 16S rRNA gene sequences retrieved from the 23 genome sequenced strains, 37 Corallococcus spp. DSM strains and 8 non-Corallococcus myxobacterial strains from the NCBI database (using Polyangium cellulosum as the outgroup). All the non-Corallococcus sequences grouped together, separately from 58 of the 60 Corallococcus sequences. Two of the DSM Corallococcus sequences (C. coralloides DSM 52496 and C. coralloides DSM 52498) grouped into a clade with Myxococcus/Pyxidicoccus and querying them against the NCBI 16S database gave best hits to Myxococcus fulvus, suggesting they had been misclassified as Corallococcus members. The remaining Corallococcus sequences grouped together, with strong (86%) bootstrap support (Supplementary Figure S1). There was little difference between the placement of the type strains of C. coralloides (DSM 2259^T) and C. exiguus (DSM 14696^T), as expected given they are likely conspecifics (Garcia et al., 2010).

Trees were also constructed for sequenced organisms using AMPHORA2 on 31 conserved genes (Figure 1A), based on genomic nucleotide identity (Figure 1B), using PhyloPhlAn on 400 conserved genes (Supplementary Figure S2), using RAxML on core pan-genome genes, and based on the presence/absence of orthogroup members (Supplementary Figure S3). All trees had similar branching patterns, and Corallococcus genomes consistently grouped into two major clades. The first, and the larger clade, containing Corallococcus coralloides/exiguus, designated as the Group A Corallococci and the second smaller clade designated as the Group B Corallococci. From inspecting the 16S sequence tree, (Supplementary Figure S1), the following DSM strains are members of the Group B Corallococci: DSM 51298, DSM 51619, DSM 51620, DSM 51625, DSM 51639, DSM 51643, DSM 52500, and DSM 52501.

Average nucleotide identity and digital DNA-DNA hybridization (dDDH) values, as calculated using the genome-to-genome calculator (GGDC), are being increasingly used to benchmark genomospeciation. Currently an ANI of ≥95% is required to conclude that two organisms belong to the same species (Richter and Rosselló-Móra, 2009; Sangal et al., 2016). An ANI of ≤92% is taken as evidence that two genomes belong to different species (<75% for different genera), while pairs of genomes giving intermediate ANI values (93–94%), likely belong to separate sub-species (Sangal et al., 2016).

ANI values for all pair-wise comparisons of the 24 Corallococcus genomes are presented in Table 2. All pair-wise comparisons gave ANI values ≥82%, confirming Corallococcus is a single genus, although Group A and Group B may represent discrete sub-genera. Using a ≤92% ANI cut-off, nine “species” groups of genomes can be observed, with members from different groups belonging to separate species (Table 2). Applying an ANI value cut-off of ≥95% gives 11 groups of genomes, within which members definitely belong to the same species. However, pair-wise comparisons giving ANI values of 93–94% associated four of those 11 into two larger groups, suggesting that two of the species groups each contain two discrete sub-species.

Table 2 also provides the dDDH values for pair-wise comparisons of the Corallococcus genomes. The use of dDDH values for species delineation is not well benchmarked, however, cut-offs of >50 and >30% give species groups mirroring the results of ANI-based species definitions at ANI values of >95 and >93%, with 9 and 11 species groups, respectively. A dDDH cut-off of > 70% has been proposed by Richter and Rosselló-Móra (2009) to define membership of the same species. Applying a >70% cut-off gives 13 species groups, due to splitting the Group A subspecies containing C. coralloides DSM 2259^T into three separate species, which is not supported by ANI values (>95%).

From the ANI and dDDH values, we therefore propose conservatively that the genus Corallococcus comprises at least nine species, with five species in Group A and four in group B: Group A species 1 subspecies 1 (CA049B, CA054B, DSM2259^T, AB011P, and AB045) and subspecies 2 (AB004, AB018, AB030, AB032C, AB38B, and CA041A), Group A species 2 (AB047A, AB049A, and AB050B), Group A species 3 (AB043A), Group A species 4 (AB050A), Group A species 5 (CA043D), Group B species 1 subspecies 1 (CA051B) and subspecies 2 (CA053C), Group B species 2 (CA040B), Group B species 3 (CA054A), and Group B species 4 (CA031B, CA031C, and CA047B).

To analyze the Corallococcus pan-genome in detail, Roary was used (with a 90% BLASTp percentage identity cut-off) to cluster the genes encoding complete protein sequences into core (hard core and soft core) and accessory (shell and cloud) genomes (Page et al., 2015). Hard core genes are found in >99% genomes, soft core genes are found in 95–99% of genomes, shell genes are found in 15–95%, while cloud genes are present in less than 15% of genomes. Of 43,561 orthologous proteins, the Corallococcus pan-genome core comprised 2,486 genes (on average 30.5% of each genome), with the accessory genome containing 9,846 genes in the shell (22.6%) and 31,229 in the cloud (71.7%) (Table 3). Enough genomes were available in Group A and Group B to also undertake pan-genome analysis of the two groups separately (Table 3). As expected, due the smaller number of genomes in Groups A and B, they showed larger cores and smaller accessory pan-genomes than for the whole set of Corallococcus genomes. Surprisingly the Group A core pan-genome was substantially larger than that of Group B (4,195 and 3,282, respectively), despite including more than twice as many genomes (17 and 7, respectively), with each Group A member having a smaller complement of accessory genes than Group B.

Figure 2 displays the number of core genes (A) and the pan-genome size (B) as a function of the number of included genomes. The number of core genes plateaus and if the number of core genes is plotted against the reciprocal of the number of genomes included, the size of the core pan-genome at an infinite number of genomes would be expected to be 988, suggesting that 1498 (60.3%) of the current core genes will prove to actually be accessory genes. However, such conclusions should be taken with caution, as 100+ genome sequences should ideally be used for these types of comparisons (Whitworth, 2008a). Nevertheless, the core pan-genome at infinite genomes is much higher for Groups A and B (96.7 and 94.6% of the current core, respectively), presumably as a consequence of the greater sequence similarity of orthologs within the groups compared to between the groups.

The size of the pan-genome increases steadily with the addition of each further genome suggesting Corallococcus has a large, open pan-genome. By fitting the curve in Figure 2A to a power law (8127n^0.5481) it is possible to extrapolate to n genomes (Table 3). Even when 1,000 Corallococcus genomes have been sequenced, it is expected that sequencing another one would add a further 196 genes to the pan-genome. Only when ∼4,500 genomes have been sequenced would fewer than 100 novel genes be identified per further genome sequenced.

In comparison to the Corallococcus pan-genome, the Myxococcus pan-genome (which included genomes of five validly described species [fulvus HW1, macrosporus DSM 14697, stipitatus DSM 14675, virescens DSM 2260, and xanthus DK1622] and M. hansupus), (Supplementary File S1) had a substantially smaller core genome (609 genes) and a proportionately larger cloud/shell pan-genome (91.6% of each genome).

SPINE (Ozer et al., 2014) was used to separate core from accessory gene sequences, and the resulting protein sets annotated using eggNOG (Supplementary File S1), for all 24 Corallococcus spp., Group A Corallococci and Group B Corallococci. The percentage abundances of COG functional categories within a set were calculated and compared between protein sets to obtain a fold-difference in relative abundance. Inclusion criteria of log₂(relative abundance) > ±1, with an abundance > 2%, were applied to identify COG categories that were relatively enriched or impoverished in a protein set. There were no clear differences in the COG categorizations of the core pan-proteomes of Group A, Group B or all Corallococci. However, differences between the core and accessory pan-proteomes were observed in categories relating to conflict/predation (Table 3). The accessory pan-proteomes of Group A, Group B and all Corallococci were enriched for category Q (Secondary metabolites biosynthesis, transport, and catabolism) compared to their respective core pan-proteomes, while the Group A Corallococci and all Corallococci (but not Group B Corallococci) accessory pan-proteomes were also enriched for proteins of category V (Defense mechanisms). As expected the core pan-proteomes were variously enriched for housekeeping COG functions E (Amino acid transport and metabolism), F (Nucleotide transport and metabolism) and J (Translation, ribosomal structure and biogenesis).

Similarly to Corallococcus, the Myxococcus accessory pan-proteome (Supplementary File S1) was enriched in proteins assigned to COG category Q (secondary metabolism) compared to the core proteome. However, in contrast to Corallococcus, the Myxococcus core pan-proteome was relatively enriched in proteins of COGs J (Translation, ribosomal structure, and biogenesis) and U (Intracellular trafficking, secretion, and vesicular transport) and relatively impoverished in proteins of COGs G (Carbohydrate transport and metabolism) and P (Inorganic ion transport and metabolism).

AntiSMASH prediction of BGCs (Supplementary File S2) revealed that the average Corallococcus genome possesses 55 BGCs (1315 BGCs were found amongst the 24 genomes). CA051B had the highest number of BGCs (94), while the type strain C. coralloides DSM 2259^T had the lowest (35). Non-ribosomal peptide synthetases (NRPS) were the most prevalent BGCs (21 per genome on average), but there were also on average ten T1PKS-NRPS and five bacteriocin gene clusters per genome. Amongst these genomes, in four cases a single genome contained a BGC for a unique predicted compound class (AB043A is predicted to produce a lantipeptide-β-lactam, AB047A makes a NRPS-T1PKS-linaridin molecule, CA049B makes a bacteriocin-T1PKS-NRPS hybrid, while CA054A produces a nucleoside).

Many of the Corallococcus BGCs showed significant similarity to BGCs producing known “named” metabolites (Supplementary File S2). Amongst the 24 genomes, BGCs producing 50 named metabolites/classes were predicted. Five predicted “core” named metabolite BGCs were present in at least 23 of the 24 genomes (including alkylresorcinol, geosmin, and nostopeptolide), 25 “cloud” named BGCs were present in three or fewer genomes (including althiomycin, chivosazole, crocacin, cryptophycin, cyanopeptin, nannoystin, pellasoren, streptomycin, teixobactin, and thanamycin), while the other 20 “shell” named metabolite BGCs were present in between four and 22 genomes each. Group A and group B Corallococci also had a small core set of BGCs producing named metabolites (five and six, respectively). However, they were different in character between the two groups. Leupyrrin was predicted to be produced by each member of Group B, but was not predicted for any member of Group A, making the leupyrrin BGC a distinctive feature of Group B genomes.

The pattern of a small set of core BGCs and a large repertoire of accessory BGCs was also observed for BGCs predicted to produce unnamed, and thus potentially novel, metabolites. SPINE was used to separate core from accessory genes for the Corallococcus genomes, and antiSMASH predictions generated. There were twenty BGCs within the core pan-genome (including those producing four named metabolites) and 290 BGCs within the accessory pan-genome (including 42 predicted to produce 21 different named metabolites). Thus, named metabolites only represent a small proportion of both the core and accessory secondary metabolomes, and the core metabolome is only a small fraction of the metabolic diversity found across the pan-genome.

These trends were also observed for the Myxococcus secondary metabolome. AntiSMASH predictions of individual Myxococcus genomes revealed BGCs predicted to produce seven “core” named compounds (alkylresorcinol, carotenoid, dkxanthene, geosmin, myxoprincomide, puwainaphycin, and VEPE_/_AEPE_/_TG-1) in all six species. In addition, a further 37 BGCs were predicted to produce another 19 “accessory” named products. Using SPINE to split the core and accessory pan-proteomes gave just one BGC in the core (producing VEPE_/_AEPE_/_TG-1) and 184 BGCs in the accessory genome including 47 predicted to make 26 different named metabolites.

To investigate relationships between taxonomy, BGCs and predation, the predatory activity of all strains was assayed against 10 prey organisms, including yeast, Gram-negative bacteria and Gram-positive bacteria. Corallococcus strains were inoculated onto lawns of potential prey organisms and the diameter of the zone of killing was measured after incubating for 7 days (Supplementary File S3). The pattern of predatory activity and prey susceptibility was similar to previous descriptions with Corallococcus strains exhibiting broad but patchy activity against all prey (Livingstone et al., 2017).

Four trees (predation, ANI, 16S, and BGCs) were constructed based on the profiles of predatory activity (Supplementary File S3), genomic ANI (Table 2), 16S rRNA sequence similarity and ANI trees based on the DNA sequences of the BGCs from each strain, extracted from antiSMASH predictions. Differences in matrix structure for sequence similarity and predatory activity precluded generation of correlation coefficients. Therefore, similarities and differences between the trees were visualized as tanglegrams (Figure 3) where the branches of pairs of trees are rearranged to give maximum similarity in topology. Tanglegrams revealed the greatest congruence between taxonomy (both ANI and 16S) and BGC trees with the predation tree exhibiting much less congruence with either taxonomy or BGC trees. Our interpretation of this result is that predatory activity is not solely dependent on BGCs or non-BGC genes acquired by vertical descent (captured in the ANI). Therefore, genes acquired by horizontal gene transfer are likely to be important determinants of predatory activity against particular prey.

In studies isolating predatory myxobacteria from soil, strains belonging to the Corallococcus genus commonly predominate, often representing more than 70% of isolates (Livingstone et al., 2017; Mohr et al., 2017). Their colorless swarms are easily distinguished from the other most abundant myxobacterial genus, Myxococcus, which is brightly pigmented yellow/orange. Currently, two species of Corallococcus are recognized, C. exiguus and C. coralloides. Although diversity within the genus has been acknowledged (Spröer et al., 1999; Stackebrandt et al., 2007), it seems likely that the two type strains belong to a single group. Garcia et al. (2010), argue that C. exiguus and C. coralloides are the same species which is also supported from the 16S inferred phylogeny of the current work.

Traditionally, a polyphasic approach has been adopted for the identification and classification of myxobacteria, based on biochemical markers, cell/population physiology, and sequence data. However, with the recent explosion of sequence information, there is increasing acceptance of using sequence data as a basis for taxonomic assignment with whole genome sequence data seen as more reliable and with greater discriminatory power than 16S rRNA gene sequence analysis alone (Janda and Abbott, 2007; Richter and Rosselló-Móra, 2009). Increasingly, genome-based studies are revealing misidentifications and misclassifications of bacteria, highlighting problems arising from single-gene 16S rRNA phylogeny (Benevides et al., 2017; Tran et al., 2017).

Recent studies have proposed minimum standards for species identification based on an overall genome related index (OGRI) which incorporates dDDH and ANI values (Chun et al., 2018). When using the recommended thresholds for ANI and dDDH, the 23 Corallococcus genomes sequenced in the current work fall into at least nine genomospecies with potentially as many as 11 identified. Further biochemical and phenotypic analyses are thus required to propose members of these genomospecies as type strains for novel species of Corallococcus. Nevertheless, identifying 8 candidate novel species from amongst 23 novel isolates implies a huge diversity of Corallococcus in the environment.

Across the full set of 24 Corallococcus genomes, the average genome is around one-third core and two-thirds accessory genome, in line with other bacterial species and genera (Lapierre and Gogarten, 2009). For bacteria, the “average” organism has a core pan-genome of around 1650 genes when calculated within a species, reducing to ∼1,350 genes when analyzed across a genus (Vernikos et al., 2015). However, this core pan-genome will depend on the percentage identity cut-offs used, the numbers of genomes included and any sampling biases within the various taxa. Corallococcus genomes are large (around 10 Mbp) and the Corallococcus core pan-genome is therefore proportionately large (around 2,500 genes). For Corallococcus, three-quarters (76.0%) of the accessory genome was composed of cloud genes being identified in three or fewer of the 24 genomes.

It is apparent that Myxococcus has access to a larger pool of accessory genes than Corallococcus and that Corallococcus has a larger set of “core” functions. However, the relative diversity in the Myxococcus pan-genome could be a consequence of including single representatives from each Myxococcus species. Thus, it would be interesting to characterize the Myxococcus pan-genome if generated from 20 independently isolated wild-type strains.

Myxobacteria have been actively screened for natural products for several decades, with more than 100 core structures and over 500 derivatives now having been published (Herrmann et al., 2017). Secondary metabolites with potent antimicrobial properties have been obtained from Corallococcus cultures, including corallopyronin, corallorazine, and coralmycin (Schmitz et al., 2014; Schäberle et al., 2015; Kim et al., 2016). However, the analysis of Corallococcus genome sequences for BGCs with the potential to synthesize novel metabolites (Xiao et al., 2011) has not previously been undertaken.

Phylogenetic trees based on the sequences of Corallococcus BGCs are only partially similar to those based on Corallococcus genomic sequence similarity (Figure 3), implying gene gain or loss has a substantial impact on the complement of BGCs in a strain’s genome. This suggestion was supported by Corallococcus pan-genome analysis, which showed that only a small proportion of the BGCs in the pan-genome could be considered core. Of the 50 BGCs with similarity to those producing named compounds, fully half were found in three or fewer of the 24 genomes. This implies that the full pan-genome of Corallococcus will include a huge array of BGCs. Therefore, it will be interesting to determine to what extent the Corallococcus and Myxococcus pan-genomes share pools of BGCs and whether there are “Corallococcus-exclusive” BGCs awaiting discovery.

It also seems likely that non-BGC horizontally acquired genes will prove to be important determinants of predatory activity given that trees based on predatory activity were barely congruent with those based on genomic ANI (Figure 3). Individual strains of Corallococcus exhibit very different patterns of predatory activity against different prey species (Supplementary File S3), and only around 30% of any Corallococcus strain’s gene set is shared by all other Corallococcus strains (Table 3). Presumably, some recently acquired genes confer selective advantages during predation, or conversely, loss of predatory genes might have impaired a strain’s ability to predate. It is an important goal to identify such genes and to discover the mechanisms by which those genes affect predatory activity. To that end, making available 23 genome sequences, from at least nine species of Corallococcus, should prove to be an important resource for the continued genetic exploration and exploitation of these inherently antimicrobial bacteria.