The Photobacterium genomes ranged in size from 4.2 to 6.4 Mb (Figure 3). The genome size range was in accordance with the ones reported for other genera from the Vibrionaceae family (Thompson et al., 2009). The largest genomes found belonged to P. profundum, which is an extremely versatile species. P. profundum SS9 is able to grow at cold temperatures and in pressurized environments (Eloe et al., 2008; Lauro et al., 2008, 2014). It seems therefore logical that such environmental versatility comes associated with a larger genome size (Konstantinidis and Tiedje, 2004).

The GC content of the genomes varied between 38.7 and 50.9%, and the strains clustered in two groups based on GC%: one of approximately 40% and another of approximately 50%. The only exceptions were P. swingsii CAIM1393 and Photobacterium sp. AK15, which had a GC content of 43.4 and 46.2%, respectively. The GC content has been linked to a division of strains that also is reflected in their environmental conditions as well as to amino acid usage (Lightfield et al., 2011). Independently of the method used, the species with higher GC content clustered together, suggesting an evolutionary association of GC content and phylogenetic proximity within the genus (Figure 2).

Several theories for genome evolution in bacteria explaining smaller cell size and genomes, the different GC content, the genomic reductions and expansions have been presented and discussed (Moran et al., 2007; Hunt et al., 2008a; Newton et al., 2010; Morris et al., 2012; Fernández-Gómez et al., 2013; Giovannoni et al., 2014; Luo and Moran, 2015). It has, for example, been shown that genome size and GC content are related with the ecological strategies of the different marine bacteria, with free-living bacteria having lower GC content and smaller genomes, as compared to patch-associated bacteria (Luo and Moran, 2015). Also, it is expected that symbionts, parasites and commensals would experience genome reduction due to specialization (Morris et al., 2012; Giovannoni et al., 2014). Here, we observe that the smallest genomes and lower GC content are indeed found in the known symbiotic organisms (P. iliopiscarium, P. damselae, P. phosphoreum), but also in P. galatheae and P. halotolerans. This could be explained by genetic drift and streamlining, respectively (Giovannoni et al., 2014; Luo and Moran, 2015).

The only two closed genomes (P. profundum SS9 and P. gaetbulicola Gung47) have a larger and a smaller chromosome of approximately 4 and 2 Mb, plus megaplasmids of 80 and 35 Kb, respectively. The presence of 2 chromosomes is a trend of the Vibrionaceae family and it is assumed that the draft genomes used also include two chromosomes and large plasmids, but fully closed genomes would be required to assert this. The number of genes per genome was estimated using prodigalrunner and ranged from 4,041 to 7,027.

To visualize the protein-coding gene content conservation in the Photobacterium genus, we constructed a BLAST atlas using P. profundum SS9 and P. gaetbulicola Gung47 as reference genomes (Figure S6). As previously described for other members of the Vibrionaceae family (Thompson et al., 2009), the Photobacterium large chromosome seems to be more conserved between species than the smaller chromosome or plasmids (Figure S6). The smaller chromosomes are highly variable and harbor the main genomic differences between strains of the same species. These secondary chromosomes and plasmids therefore seem to be the source of the genetic plasticity, mirrored in the different phenotypes observed within members of the Photobacterium genus and Vibrionaceae family in general (Vesth et al., 2010; Lukjancenko and Ussery, 2014). A well-studied example is V. splendidus, where the genomes of strains from the same species vary and result in distinct phenotypic capabilities (Thompson J. R. et al., 2005; Hunt et al., 2008a). This diversity is also evident for the Photobacterium species in the several genetic features analyzed, such as virulence, bioluminescence, histamine production, secondary metabolism and CRISPR-Cas operons. All these genomic features presented random distributions rarely correlated with phylogenetic relatedness (Figure 4), raising the issue of species-phenotype association within the Photobacterium genus. Nevertheless, genomic features such as genome size and GC content seem to be associated with different lifestyles adopted by Photobacterium species.

Bioluminescence for example was initially thought to be a widespread feature of the genus (Urbanczyk et al., 2011). We identified the genetic basis for bioluminescence only in three species P. kishitanii, P. phosphoreum, and P. leiognathi (Figure S5), isolated from different marine animals (Table S1). This also compared with the phylogenetic analysis (Figure 2A), with the exception of P. iliopiscarium, that may have lost this trait due to niche adaptation (Figure 4). Analysis of the lux-rib operon further showed that strains ATCC 25521, ATCC 33979 and Irivu.4.1 are most likely P. leiognathi subsp. leiognathi, according to their lux-rib gene organization (Figure S5) (Ast and Dunlap, 2004). Similarly, histamine production seems to be specific for P. kishitanii and P. damselae, however genes responsible for this feature could be identified also in one P. angustum strain (ATCC 33977) (Figure S7). Scombrotoxin fish poisoning (high levels of histamine) is the most frequent cause of fish poisoning incidents within the United States (Pennotti et al., 2013). Recently, high-histamine producing Photobacterium strains have been isolated from freshly caught fish (Bjornsdottir-Butler et al., 2016), raising awareness for the need to revise food safety rules regarding sea food. There are also several reports describing P. phosphoreum as histamine producing species (Kanki et al., 2004, 2007), but we did not identify the genes responsible for histamine production in P. phosphoreum ATCC 11040, although it may be present in other P. phosphoreum strains.

The pan-genome refers to the total number of genes in all the 35 strains, while the core genome represents the number of orthologous genes shared between them. Using 35 genomes we identified a pan-genome of 28,951 genes and a core-genome of 1,232 genes (Figure 5). Taking into consideration the average gene number of 4,750 for the Photobacterium strains, 1,232 genes, represents approximately 25% of the total genome, meaning that approximately 1/4 of the genome is conserved in all the strains. The number of core-genome genes is in agreement with what has been found in other marine Gammaproteobacteria (Qin et al., 2013), although this is two times higher than what has been found for the Vibrio genus, where the core-genome comprised approximately 500 genes (Thompson et al., 2009). This difference can be explained by the high number of Vibrio species (more than 120), which reflects the genomic and ecological diversity of the genus.

The pan-genome for the Photobacterium genus alone is greater than what has been reported in a study of the Vibrionaceae family, where 43 genomes of a total of 13 species from three different genera were used (Thompson et al., 2009). This suggests Photobacterium harbors high genomic diversity, reflected in their ability to colonize different environmental niches (Konstantinidis and Tiedje, 2004; Konstantinidis et al., 2009) and supports the theory that high gene content variation exists in environmental marine strains (Tettelin et al., 2008; Konstantinidis et al., 2009). Using a power-law regression, it is possible to evaluate the openness of a pan-genome (Tettelin et al., 2008). The Photobacterium pan-genome is open, with a γ parameter of 0.62 in a power-law regression fitting relatively well the data analyzed (R² = 0.89) (Figure S8). This genetic variation and uniqueness of each strain is also evident in the number of singletons per strain, which for some strains represent almost 20% of the genes identified in the genome (Table S3).

Most Photobacterium species have been described as symbiotic or associated with other marine organisms (Ast et al., 2007; Dunlap, 2009; Gomez-Gil et al., 2011; Urbanczyk et al., 2011, 2013). Furthermore, some strains of P. damselae are pathogens of marine organisms, especially fish (Hundenborn et al., 2013; Andreoni and Magnani, 2014). The key virulence genes of P. damselae are a phospholipase-D damselysin gene (dly) and a pore-forming toxin gene (hlyA), which act in a synergistic manner (Rivas et al., 2011, 2013; Le Roux et al., 2015). We conducted a homology search to evaluate the possible virulence of other Photobacterium species. The genes could only be identified in P. damselae subsp. damselae CIP 102761. In this strain, two copies of the hlyA gene were identified, one in contig_1, close to an IS4 transposase and a phage integrase and the other in contig_4, next to the dly gene.

Often draft-whole genome sequences contain plasmid sequences. Plasmids are important mediators of several physiological traits of Photobacterium, such as virulence, drug resistance and biosynthetic capabilities (Kim et al., 2008; Rivas et al., 2011; Nonaka et al., 2012; Osorio et al., 2015). We compared known Photobacterium plasmids to the here studied genomes. Plasmid pPHDD1 showed high similarity to contig_4 of P. damselae subsp. damselae CIP 102761 genome. Some of the contigs of P. damselae subsp. piscicida DI21 also had high similarity to plasmids pPHDP60, pPHDP10 and pPHDP70, which have been previously isolated from this strain (Osorio et al., 2008, 2015).

Using Island Viewer (Dhillon et al., 2015), we searched for genomic islands in the fully sequenced genomes of P. profundum SS9 and P. gaetbulicola Gung47, and compared these to genomes of other strains from the same species (Figure 6). Some of the major genetic differences between the strains of the same species seem indeed to be related to the presence or absence of specific genomic islands. Genomic regions only present in the reference strain are placed close to identified genomic islands.

Additionally, strains varied in their number of transposase genes (Table S3). While P. profundum SS9 had 219 transposase genes, P. leioghnathi subsp. mandapamensis svers.1.1. had none. No correlation between number of prophages and number of transposase genes was observed. A high number of transposase genes could indicate transposon-mediated exchange of genetic material from any source (both plasmid-borne and random environmental DNA). P. profundum SS9 and Photobacterium sp. SKA34 are cases where transposase genes represent almost 4% of the total genes in the genome (Table S3). This is in accordance with metagenomics data that showed that transposase coding genes are the most abundant genes in nature, most likely accelerating biological diversification and evolution (Aziz et al., 2010) and that a higher number of transposase genes is present in higher oceanic depths (Konstantinidis et al., 2009).

Using the PHAge Search Tool (PHAST) (Zhou et al., 2011) we identified 33 intact prophage sequences and 72 incomplete ones (Table S3). From the 33 intact ones, 17 were unique, while the other 16 were re-occurrences of three prophages. One was present only in P. profundum strains; another in all P. kishitanii strains, plus one P. damselae and one P. angustum; and the third was randomly distributed (Figure 4).

The bacterial and archaeal adaptive immune systems entail CRISPR-Cas modules (Makarova et al., 2015) and we queried the genomes for the architecture of the CRISPR-Cas systems including the cas gene organization, the direct repeats (DR) and the protospacers in the CRISPR locus. We divided the architecture of the CRISPR-Cas systems into seven clusters (Figure 7). Most of the architectures identified were similar to the ones described in Yersinia pestis, Escherichia coli and Desulfovibrio vulgaris (Figures 7A–C; Haft et al., 2005). Also, we identified two clusters in P. profundum SS9 similar to the ones of Y. pestis and E. coli, encoded in the chromosome and plasmid, respectively (Figures 7D,E). These clusters had different gene arrangement and/or included genes coding for unknown proteins in the operon. Two other clusters containing CRISPR-associated genes were identified (Figures 7F,G).

The estimation of DRs and protospacers in the CRISPR arrays is very difficult in draft genomes due to the short sequencing reads and the repetitive nature of the sequences, which make them difficult to assemble correctly. In most of the cases, both direct repeats and protospacers could be identified both upstream and downstream of the cas-operon. These were also identified elsewhere in the genome or in distinct contigs, usually consisting of one of the ends of the contig. For example, further upstream from the Tn7 in the P. galatheae S2753 (Figure 7G), an array of 16 protospacers could be identified beside a gene coding for a DNA nicking enzyme.

The DRs that flanked the protospacers were extracted and compared. The similarity of the DRs correlated with the type of clusters the strain had. DRs from clusters (a) and (d) were similar, as were the ones within (b) and (f). The protospacers are short sequences that are derived from bacteriophages or other foreign DNA, such as conjugative plasmids (Attar, 2015; Makarova et al., 2015). Therefore, these sequences may provide a history of encounters of a specific bacterium with phages and/or plasmids. The number of spacers varied considerable between different strains, from 1 to 64 protospacers in the same CRISPR array. Arrays with extensive number of protospacers in some strains indicate numerous bacteriophage infections, which is the case of e.g., P. profundum SS9, P. angustum S14 and P. aquae CGMCC 1.12159 with 64, 45, and 33 protospacers, respectively (Figure 7). The 4 spacers in P. leiognathi lrivu.4.1 were 100% identical to the first 4 spacers in P. angustum S14, which had an array of 45 spacers. Also high similarity (>93% Identity) could be identified in the first spacers of P. angustum ATCC 25915, ATCC 33975 and P. damselae subsp. damselae ATCC 33539 (which should be classified as P. angustum), although for the P. angustum strains the protospacers have not just been identified downstream of the cas operon, but elsewhere.

We have also compared the protospacers to the prophage sequences identified in these genomes. The bacterial immune system (CRISPR-Cas) allows the protection against re-infection by the same bacteriophage (Attar, 2015), nevertheless, for P. angustum ATCC 33975 an array of 9 sequential spacers were 100% identical to an intact prophage sequence identified within the same genome, suggesting multiple re-infection events.

Phage infection seems to be a frequent event in some strains, supported by the extensive number of protospacers in the CRISPR arrays and the number of prophages within the genomes (Figure 7). Remarkably, the cas operons previously associated with E. coli, D. vulgaris and Y. pestis were identified in Photobacterium strains (Haft et al., 2005). Yet, presence/absence and type of cluster seems randomly spread across the different species (Figure 4).

Another indication of high genomic exchange is the number and distribution of secondary metabolism biosynthetic clusters, many of which are believed to be acquired by HGT (Khaldi et al., 2008). Secondary metabolites such as non-ribosomal peptides and polyketides are known to have antagonistic properties, which may constitute an advantage in several ecological niches. Although Photobacterium strains have mostly been studied due to their association with marine animals, their potential in drug discovery and other applications has recently been reported (Wietz et al., 2010b; Machado et al., 2014, 2015a,b; Nielsen et al., 2014; Mathew et al., 2015a,b).

Using antiSMASH 3.0 (Weber et al., 2015), different biosynthetic gene clusters were identified (Table S4) with some clusters being present across the genus and others being species specific (Figure 4). The same terpene biosynthetic cluster could be identified in all the P. angustum strains (including the misidentified P. damselae subsp. damselae ATCC 33539). The function of this terpene cluster is not known, but it may be related to the planktonic lifestyle of this species, since all the P. angustum strains were isolated from seawater. Another species-specific cluster was the polyunsaturated fatty-acid (PUFA) cluster present in P. profundum strains. This cluster could be involved in the high pressure and cold temperature adaptation of this species, since these PUFAs are known to modify membrane fluidity in response to hydrostatic pressure and temperature (Campanaro et al., 2005).

Other biosynthetic clusters such as siderophore, aryl-polyene and ectoine were widely distributed across the genus (Figure 4). A siderophore cluster was present in 11 strains, although two distinct siderophore biosynthetic clusters were identified. The ectoine cluster was present in 18 strains (two had only two out of three genes needed for its biosynthesis). The most widely distributed cluster was the aryl-polyene biosynthetic cluster, identified in 23 out of the 35 strains. Interestingly strains associated with marine animals did not have a siderophore cluster nor were they prolific in other secondary metabolite clusters (Figure 4, Figure S4). Strains isolated from seawater or sediments seem to be enriched in secondary metabolite clusters, when compared to marine animal associated species. These strains have therefore the genetic capability of possibly antagonizing microbial competitors in their environment.

The distribution of secondary metabolite clusters suggests events of gain and loss of these clusters throughout the evolution of Photobacterium species. The random presence of specific traits, i.e., not phylogenetically related, has been reported for other genus. For example, in Pseudovibrio the random distribution of type IV secretion systems was attributed to the frequency with which those genes are horizontally acquired (Cascales and Christie, 2003; Romano et al., 2016). In fact, several of the interesting phenotypic traits identified in Photobacterium strains, such as bioluminescence, virulence, histamine production and piezophilia seem to be acquired by HGT (Ast and Dunlap, 2004; Campanaro et al., 2005; Dunlap, 2009; Urbanczyk et al., 2011; Hundenborn et al., 2013; Andreoni and Magnani, 2014; Osorio et al., 2015; Bjornsdottir-Butler et al., 2016).