Piscirickettsia salmonis, a Gram-negative and facultative intracellular bacterium, is responsible for salmonid rickettsial septicemia, or piscirickettsiosis (Fryer et al., 1990; Yañez et al., 2012). This disease causes high mortality rates in the three most important farmed salmonid species for the Chilean aquaculture industry, resulting in significant economic losses (Almendras and Fuentealba, 1997). P. salmonis was originally isolated and identified from diseased Oncorhynchus kisutch (Coho salmon) in southern Chile in 1989 (Bravo and Campos, 1989; Fryer et al., 1990), but this pathogen is not restricted to Chile. Indeed, piscirickettsiosis has subsequently been reported in Canada (Brocklebank et al., 1993), Ireland (Rodger and Drinan, 1993), Scotland (Grant et al., 1996), Norway (Olsen et al., 1997), southern USA (Olsen et al., 1997; Arkush et al., 2005), and Turkey (Öztürk and Altinok, 2014). P. salmonis causes systemic, chronic septicemia in fish, which exhibit a variety of internal symptoms, such as a discolored kidney, enlarged spleen, and pale liver with numerous nodules measuring 5–6 mm in diameter (Branson and Nieto Diaz-Munoz, 1991; Cvitanich et al., 1991). Despite research in several in vitro and in vivo models (Rozas and Enríquez, 2014), the mechanisms involved in P. salmonis pathogenicity are not entirely known. As such, current vaccination strategies, antibiotic treatments, and other biotechnological tools are not effective in controlling piscirickettsiosis (SERNAPESCA, 2016b). Furthermore, the continued high incidence of piscirickettsiosis outbreaks reported in Chile (SERNAPESCA, 2016a) is likely due to antibiotic resistance acquired by the pathogen (Cartes et al., 2016; Henríquez et al., 2016; Sandoval et al., 2016).

P. salmonis was originally grouped in the Rickettsiaceae family and reported as a rickettsia-like organism (Fryer et al., 1990, 1992). Further analyses of the nucleotide sequence and secondary structure of the 16S ribosomal RNA (rRNA) resulted in the relocation of P. salmonis into the Gammaproteobacteria in cursive class and Piscirickettsiaceae in cursive family (Fryer et al., 1992; Fryer and Hedrick, 2003). More recent research by Reid et al. (2004) using specific primers for 16S ribosomal DNA and the intergenic transcribed spacer (ITS) region of Irish and Scottish isolates revealed the potential existence of two geographically and phylogenetically separate P. salmonis groups. Similar approaches evidenced the potential existence of two representative genogroups isolated in Chile (Mauel et al., 1999; Bohle et al., 2014; Mandakovic et al., 2016; Otterlei et al., 2016), as specifically associated with the EM-90 and LF-89 strains. Differences within these genogroups have been reported in relation to geographic distribution, antibiotic susceptibility, and host specificity (Saavedra et al., 2017) Furthermore, these genogroups have been computationally compared using a limited number of fully sequenced genomes. For example, Bohle et al. (2014) used two complete genomes (strains PM15972A1, EM-90-like; and PM32597B1, LF-89-like) to genomically characterize the potential existence of genogroups. Additionally, Bravo and Martinez (2016) used 11 genomes (five complete and six draft genomes) to obtain a first approximation of the P. salmonis pan-genome. Nevertheless, these studies are limited in that they were performed with a restricted quantity of fully sequenced genomes. For instance, only one representative of the EM-90 strain was used by Bohle et al. (2014) and Bravo and Martinez (2016).

As with phylogenetic relationships and functional annotations, studies evaluating genomic structural divergences in virulence factors and metabolism within genogroups are scarce. Comparative genomic analyses of fully sequenced genomes are fundamental for defining the entire core- and pan-genomes of different isolates from the same species. The core-genome is defined as the entire repertoire of translated genes conserved among all isolates, with conservation suggesting that these genes are essential for basic bacterial survival. In turn, the pan-genome is the sum of the core genes and those within the “accessory genome,” i.e., the set of unique genes shared by or exclusive to certain strains. The presence of a pan-genome is associated with increased intra-species diversity, environment/host adaptations, and differentiations in pathogenic mechanisms (Tettelin et al., 2008). Several pan-genomic studies exist for different fish microorganisms, including marine pathogens such as Aeromonas salmonicida (Vincent and Charette, 2017), the Alteromonas ssp genus (López-Pérez and Rodriguez-Valera, 2016), and Vibrio harveyi (Espinoza-Valles et al., 2015). To define the entire core- and pan-genomes of an organism, a large quantity of fully sequenced genomes is essential. Advances in recent years mean that enough P. salmonis genome sequences are now available for pan-genome analysis.

The aim of this study was to differentiate and classify the number of P. salmonis genogroups by using phylogenetic and phylogenomic analyses, as well as by identifying structural and functional genomic divergences within the 19 fully sequenced genomes available for this pathogen. The obtained phylogenies and genomic structural divergences among the assessed P. salmonis strains strongly support a clear existence of two genogroups, termed herein as LF and EM. Moreover, pan-genome characterization identified genes involved in the P. salmonis core-genome that play roles in the life cycle, invasion, and pathogenic processes. Differences in virulence factors associated with established genogroups were elucidated by identifying and categorizing genogroup-exclusive genes. The detection of different genes potentially associated with pathogenicity, infection mechanisms, and host-pathogen interactions in the LF and EM genogroups is an important step toward understand the nature of variants of P. salmonis that needs to be explored in more detail in in vivo challenge studies.

Different studies using 16S rRNA-based phylogenetic approaches suggest the existence of two clusters among P. salmonis strains, as associated with the representative LF-89 and EM-90 strains (Mauel et al., 1999; Bohle et al., 2014; Mandakovic et al., 2016; Otterlei et al., 2016). Differences between these genogroups include geographic distributions, antibiotic susceptibilities, and host specificities (Saavedra et al., 2017). Recently, Bravo and Martinez (2016) reported comparative genomics results supporting the existence of these two genogroups. More specifically, seven draft genomes and five complete genome sequences (LF genogroup: 4 and EM genogroup: 1) were used to describe the P. salmonis core-genome, with a total of 1,683 proteins. In the present study, comparative genomics analyses were performed using 19 complete genome sequences (LF genogroup: 13 and EM genogroup: 6). Phylogenetic analysis was conducted using seven conserved single-copy genes. In turn, whole core-genome phylogenomics analysis used the core proteins from all 19 genomes, resulting in strong evidence for the existence of two genogroups. This analysis of all 19 isolates also provided a comprehensive representation of the P. salmonis core- (1,732 proteins) and pan-genomes (3,463 proteins), as well as of the genogroup-specific core- (LF genogroup: 2,170 proteins and EM genogroup: 2,228 proteins) and pan-genomes (LF genogroup: 2,924 proteins and EM genogroup: 2,778 proteins).

The “open” or “close” state of pan-genome it will depend of the capacity of acquire exogenous DNA (Medini et al., 2005; Diene et al., 2013), This phenomenon is commonly seen in species living in bacterial communities, tend to have large genomes and an open pan-genome, a high horizontal gene transfer range and several ribosomal operons (Georgiades and Raoult, 2011; Diene et al., 2013). Exist evidence that the pathogen P. salmonis generates biofilm in stages of stress (Marshall et al., 2012; Albornoz et al., 2017), which are essential for the organization of bacterial communities (Filloux and Vallet, 2003). The obtained results indicate that the P. salmonis pan-genome has an open state. Any efforts to “close” the pan-genome should focus on obtaining new complete genomes from various isolates. The higher number of core proteins found in the present study might be due to the use of more complete genomes than have been analyzed in prior reports. Indeed, comparative genomic assays using incomplete genomes run the risks of only detecting partial gene sequences, which may erroneously be classified as novel gene families, or not detecting certain genes at all. For instance, publicly available genomic P. salmonis annotations present potentially erroneous protein-coding gene predictions. Several incomplete genes, i.e., containing only partial fragments of protein-coding genes, are spread across all available incomplete genome sequences (Denton et al., 2014). The presence of these small gene fragments could be the result of inherent genomic traits (e.g., chromosome rearrangements) (Darling et al., 2008; Jackson et al., 2011), gene duplications (Jackson et al., 2011), or bioinformatics-based biases in errors (e.g., assembling complications, false positive gene predictions) arising from the incompleteness of these genome sequences (Schnoes et al., 2009; Denton et al., 2014).

The annotation process of all predicted genes revealed functional similarity and divergences within strains and genogroups according to KEGG metabolic pathways and gene ontology categories. Previous studies confirm that metabolic pathways of nucleotides, followed by carbohydrates, cofactors and vitamins, amino acids and lipids are highly conserved in gammaproteobacteria (Poot-Hernandez et al., 2015). These same metabolic pathways are highly conserved in both P. salmonis genogroups, as can be seen in the “Core-genome” item in Figure 5.

On the other hand, a difference might exist at the colony-phenotype level for both LF-like and EM-like P. salmonis strains (Bohle et al., 2014). For instance, LF-89-like isolate colonies present a stick phenotype (Supplementary Figure 7A) when left for several days on culture plates, whereas EM-90-like strain colonies present a mucoid phenotype (Supplementary Figure 7B). This phenotype difference between EM and LF may be due to the varied distribution of genes encoding glycosyltransferase in association with lipopolysaccharides formation. More specifically, distinct protein sets appear to be preferably distributed within genogroups. Moreover, studies by Saavedra et al. (2017) demonstrate that this differential phenotype is confirmed by differential growing kinetics at 18 and 22°C. Interestingly, only the EM genogroups is able to grow at 22°C whereas the LF genogroups did not grow (Saavedra et al., 2017). In that work, they use different field isolated that belong to both genogroups. Similar analysis of the fully sequence genomes strains used in this in silico analysis shown the same results indicating that both genogroups have divergences that change the ability to grow at different temperature (Supplementary Figure 6).

This is the first phylogenomic and structural genomics study suggesting the existence of two P. salmonis genogroups (denominated LF and EM), as well as at least two subgroups within the EM genogroup (determined via phylogenomic and structural localization of the ribosomal operon in the chromosome). These genomic analyses are interesting in light of previous studies and particularly the recent work by Rozas-Serri et al. (2017), a study in which controlled challenge assays resulted in different clinical signs for the EM and LF strains. These differences demonstrate the need for genetically understanding the differences between the genogroups characterized in the present study. The data-supported existence of these genogroups represents material and basic knowledge that might impact on vaccine development where it will be important to show if there are immunogenic differences between EM and LF genogroups, particularly since current vaccines are mainly based on LF variants. Additionally, these findings should also encourage comparison of EM and LF genogroup with regard to pathogenic differences both in vitro and in vivo. Moreover, the study of multiple genomes of pathogen species combined with reverse vaccinology can let to identify common antigens that could potentially protect against more than one subgroup of P. salmonis and other bacterial species. Further studies need to be done in order to stablish the potential of this vaccine technology (Gourlay et al., 2017).

Common virulence factors were found in the P. salmonis core-genome, and the identified P. salmonis virulence factors were grouped into the following six functional categories: (i) endotoxins, including lipopolysaccharides and lipooligosaccharides; (ii) iron uptake, including Fur, siderophores, vibrioferrins, and Heme biosynthesis; (iii) capsule and other surface components, including capsule biosynthesis, exopolysaccharides, the type IVB secretion system, and the type VI secretion system; (iv) adherence, colonization, and invasion factors, including type IV pili, flagella, and polar flagella; (v) enzymes, including Pantothenate biosynthesis and GAPDH; and (vi) stress response category, including Hsp40, Hsp60, Hsp70, and ClpP. Interestingly, some virulence factors were genogroup-specific. For example, LF isolates evidenced an exclusive presence of the Patatin-like (Pat) protein, which exhibits phospholipase activity, which are associated with host entry, phagosomal escape, and host lysis for spreading cell-to-cell (Rahman et al., 2013); transposases, which are responsible for transposons movement (Jackson et al., 2011); and a glycosyltransferase group 1, which catalyzes the transfer of glucose residues involved in synthesizing the outer core region of lipopolysaccharides (Ünligil and Rini, 2000; Leipold et al., 2007). In turn, EM strains uniquely evidenced the Listeria Adhesion Protein, that interacts with the Hsp60 host-cell receptor to promote bacterial adhesion in Listeria ssp. (Jagadeesan et al., 2010); the fic/DOC family protein, a toxin that mediates post-translational modifications of host cell proteins, thus specifically interfering with signaling, cytoskeletal processes, or translation (Roy and Cherfils, 2015); and a putative glycosyltransferase, which is essential for O-polysaccharide synthesis (Dabral et al., 2015). These genogroup-specific virulence factors could grant special traits. In a similar context, specific Rickettsia typhi groups possessing genes encoding for the Patatin-like protein (Rahman et al., 2013) have advantages against the host during the infection process.

Gene duplication is an important factor for the formation and evolution of bacterial genomes (Romero and Palacios, 1997). Six complete rRNA operons were found in all presently analyzed P. salmonis genomes. All six rRNA operons shared the same basic structure (16S–ITS−23S−5S), and two carried a 276-nucleotide insert on the respective ITS regions. Differences in the genomic organization of rRNA operons were able to associate both types of ribosomal operons with the two defined P. salmonis genogroups (LF or EM). In a recently published study, Saavedra et al. (2017) described marked differences between genogroups, with the EM strains being able to grow in a wider temperature range and having lower nutritional demands. Another structural divergence within genogroups was in the organization of ribosomal operons. The rRNA operon is composed by the 16S, 23S, and 5S rRNA genes, as well as by an ITS between the 16S and 23S genes (Lecompte, 2002). The rRNA operon is essential for protein biosynthesis and, since the 1980s, has been the gold standard for bacterial identification and for establishing evolutionary relationships between bacterial species (Woese, 1987; Hashimoto et al., 2003; Pei et al., 2010). rRNA operon analyses have also helped elucidate how different bacteria grow and thrive. Indeed, the rRNA operon copy number is considered a measure of growth rate and/or adaptability, as organisms with a higher number of rRNA operons can achieve faster growth under certain conditions (Klappenbach et al., 2000; Yano et al., 2013). Despite this, not all rRNA operons are transcribed at the same rate, and in Escherichia coli, the loss of certain operons does not affect growth (Condon et al., 1995). In turn, diverse promoters of rRNA operons are differentially regulated under stringent conditions in Bacillus subtilis (Samarrai et al., 2011), while in Vibrio cholerae, rearrangements of constitutive genes, including ribosomal operons, increase the bacterial transcription rate if able to move to a location closer to the replication origin (Couturier and Rocha, 2006). Organisms with more chromosomal rRNA operons might have a higher potential for “fine-tuning” expressions to adapt to different environmental conditions (Condon et al., 1995). This would reinforce the postulation that while some operons are important for growth rate or physiology, others have more adaptive purposes. This mechanism also accounts for the quantity of type IVB secretion system clusters spread over the P. salmonis genomes, as well as for the degree of observed synteny. Bacterial genomes may carry multiple divergent copies of T4BSS, which is reflective of functional diversity (Gillespie et al., 2015). Intriguingly, this was not the case for P. salmonis T4BSSs, which was found to present copy-number divergences mainly within genogroups. For instance, Dot/Icm T4BSS in LF isolates was predominantly distributed in three copies, each composed by 16 genes. Exceptions were the PSCGR01 isolate, distributed into two copies. In EM strains, the Dot/Icm secretion system was distributed in two copies, each composed by 16 genes. The expression of Dot/Icm T4BSS genes has been detected in P. salmonis, supporting participation in transcriptional activity (Gómez et al., 2013). The data suggest that the chromosomal divergences in genomic localization and quantity of genetic cassettes found in the P. salmonis genome, indicated that is a common driving force for the differential evolution of LF and EM genogroups.

Comparative genomic analysis of 19 fully sequenced P. salmonis genomes provided evidence for the existence of two genogroups, as linked to the EM and LF type strains. This hypothesis was further supported by phylogenetic and phylogenomics analysis; and chromosome distribution assessments for the ribosomal and type IVB secretion systems. Each genogroup presented differences in DNA rearrangements, with more divergences observed within isolates from the EM genogroup. Divergences included variations in ribosomal operons, suggesting a potential evolution into two novel subgroups. On the other hand, the P. salmonis pan-genome was determined in an open state, which might account for the observed genomic variability. Pan-genome analysis also allowed for identifying auxiliary and unique genes specific to each strain and genogroup. These genes might provide diverse advantages to the respective genogroups.

To the best of our knowledge, the present study is the most comprehensive comparative genomics characterization of complete genome sequences for P. salmonis. The characterized divergences and similarities represent a notable contribution toward understanding the biology and evolution of P. salmonis. Additionally, the obtained data serve as a starting point in developing innovative methodologies for discriminating between P. salmonis strains and genogroups, which will be useful for more specifically diagnosing infection and for selecting potential candidate genes that could be applied in the control of piscirickettsiosis.

VM and AY supervised the study. GN, CM, VM, and AY participated in the design and discussion of the research. GN, CM, and DE performed bioinformatics analyses. GN, VM, AY, PR, LV, JC, JF, and MM contributed to the acquisition, analysis, and interpretation of data for the work. GN, PS, CO, VM, and AY drafted the manuscript. All authors read and approved the final manuscript.