Genomic insights into the virulence and salt tolerance of Staphylococcus equorum

To shed light on the genetic background behind the virulence and salt tolerance of Staphylococcus equorum, we performed comparative genome analysis of six S. equorum strains. Data on four previously published genome sequences were obtained from the NCBI database, while those on strain KM1031 displaying resistance to multiple antibiotics and strain C2014 causing haemolysis were determined in this study. Examination of the pan-genome of five of the six S. equorum strains showed that the conserved core genome retained the genes for general physiological processes and survival of the species. In this comparative genomic analysis, the factors that distinguish the strains from each other, including acquired genomic factors in mobile elements, were identified. Additionally, the high salt tolerance of strains enabling growth at a NaCl concentration of 25% (w/v) was attributed to the genes encoding potassium voltage-gated channels. Among the six strains, KS1039 does not possess any of the functional virulence determinants expressed in the other strains.

Haemolysis activity test. TSA supplemented with 5% sheep blood (v/v) (BBL Microbiology Systems, Sparks, MD, USA) was used for βand δ-haemolytic activity tests. β-Haemolytic activity was determined by cold shock at 4 °C for 24 h after incubation at 30 °C for 24 h; δ-haemolytic activity was determined by cross-streaking the test strains perpendicularly to S. aureus RN4220 at 30 °C for 24 h. The same method was implemented in the selection of S. equorum strains having βand δ-haemolytic activities 13 . Experiments were conducted three times, on separate days.
Growth monitoring in the presence of antibiotics and NaCl. S. equorum strains cultured in TSB were normalized to 0.5 turbidity at OD 600 and then diluted 1:100 in TSB supplemented with antibiotics to confirm the function of annotated antibiotic resistance genes. The antibiotics chloramphenicol, ciprofloxacin, erythromycin, lincomycin, methicillin, penicillin G and tetracycline were purchased from Sigma (St. Louis, MO, USA) and employed at concentrations of 30,5,15,30,8,6 and 30 μg/ml, respectively, based on previous research results 11,13 and the guidelines of methicillin for S. aureus set out by the Clinical and Laboratory Standards Institute 28 . The salt tolerance of S. equorum strains was determined by examining their growth in TSB supplemented with NaCl at concentrations of 15% (w/v), 20% and 25%. To determine the salt tolerance of E. coli transformants, NaCl was employed at a final concentration of 3% or 6% in LB broth. Cell growth was monitored by measuring OD 600 using a Varioskan Flash (Thermo Scientific, Waltham, MA, USA).
Cloning of the potassium voltage-gated channel genes. The potassium voltage-gated channel genes of strains KM1031 and C2014 were amplified with specific primer sets (Supplementary Table S1) designed from their genome sequences. PCR amplifications were performed using a T-3000 thermocycler (Biometra, Göttingen, Germany) and the PCR mixture consisted of the template DNA, 0.5 μM of each primer, 1.25 units of Inclone Taq polymerase (Inclone Biotech, Daejeon, South Korea), 10 mM dNTPs and 2 mM MgCl 2 . Samples were preheated for 5 min at 95 °C and then amplified using 30 cycles of 1 min at 95 °C, 30 s at 55 °C and 1 min at 72 °C. The amplified PCR products were cloned into pGEM-T Easy vector (Promega, Madison, WI, USA) under the control of the T7 promoter. The successful integration of the cloned fragments was confirmed by sequencing using GenoTech (Daejeon, South Korea) with the T7 primer set.
Nucleotide sequence accession numbers. The complete genome sequences and annotation data of S. equorum strains KM1031 and C2014 have been deposited in GenBank under the accession numbers CP013980-CP013983 and CP013714-CP013719, respectively. The average genome sequence length of the six strains is 2,828,805 bp. S. equorum UMC-CNS-924 exhibits the smallest genome (2,700,865 bp), while strain C2014 possesses the largest one (2,930,519 bp). All S. equorum strains display average G + C content of 33%.

Results and Discussion
To facilitate a coherent comparative analysis, we performed consistent ORF prediction for the six S. equorum (complete and incomplete) genome sequences. In this way, comparable numbers of genes were obtained for each genome, with an average of 2,672 ORFs per genome (Table 1). Notably, (BLAST-based) functional in silico prediction could be performed for 89.2% of the identified ORFs, while the remaining 10.8% not assigned a COG functional classification were predicted to encode hypothetical proteins.
Analysis using the SEED subsystem categorization and COG functional categorization predicted the existence of an average of 1,998 CDSs and 2,383 CDSs per genome, respectively (Table 1). Based on the SEED subsystem, over 321 CDSs accounting for 15.8-16.8% of the S. equorum genomes were allocated to the genes for amino acid biosynthesis and utilization (Fig. 1). The next most abundant subsystem category is related to carbohydrate utilization (15.5-16.7%), followed by protein metabolism. Major COG subsystems are related to amino acid transport and metabolism, as well as carbohydrate transport and metabolism. Both analyses enabled coherent conclusions to be drawn regarding the subsystem category of S. equorum.
Comparative analysis of S. equorum genomes. Whole-genome comparison of the six S. equorum strains showed that the genomes are highly homologous in terms of functional category (Fig. 1). MAUVE alignment of the six genomes allowed the identification of approximately 8-10 locally collinear blocks (LCBs), regions without rearrangement of the homologous backbone sequence ( Supplementary Fig. S1). However, the LCBs are interspaced by specific DNA stretches of various lengths. MAUVE analysis showed an overall collinear relationship across S. equorum strains KM1031, C2014 and KS1039, which were isolated from different jeotgal samples 9,13 . When the genome of strain KS1039 was established as the standard, a large-scale chromosomal reorganization by a single recombination event was found to have occurred in the genome of strain G8HB1, which resulted in the inversion of a genomic region; in addition, the results indicated that the genomes of strains Mu2 and UMC-CNS-924 might have been generated by complex rearrangements. The collinear relationship and complex rearrangement found in the six S. equorum genome structures can be explained by their differences in isolation source and geographic location. However, the contigs in genomic data for strains Mu2, UMC-CNS-924 and G8HB1 might have distorted the MAUVE analysis results. The gene pools shared by the genomes of the five S. equorum strains KM1031, C2014, KS1039, Mu2 and UMC-CNS-924 are depicted in a Venn diagram (Fig. 2). These five strains share 2,166 CDSs in their core genome, corresponding to approximately 76.1-86.7% of their ORFs. Many of the CDSs in the core genome are assigned via COG annotation to functions relating to metabolism and the transport of amino acids and carbohydrates. The genome of strain UMC-CNS-924 has the smallest proportion (2.0%) of unique CDSs that are absent from the four other S. equorum genomes. In contrast, the proportions of unique CDSs in the genomes of strains KM1031, C2014, KS1039 and Mu2 are 4.5%, 10.9%, 6.0% and 8.1%, respectively. The majority of singleton-specific genes  are associated with hypothetical proteins (Supplementary Table S2). Meanwhile, functional singletons in the genomes of strains KM1031, C2014, KS1039, UMC-CNS-924 and Mu2 are allocated to transporter, transposase, CRISPR-associated protein, tetracycline resistance and phage-related genes, respectively.

Insights into virulence.
The well-known food pathogen S. aureus produces several virulence factors: a clumping factor to protect against phagocytosis, an extracellular adhesion protein for adhesion, an enterotoxin, a toxic shock syndrome toxin, and cytotoxins such as αand β-haemolysins for tissue invasion 29 . However, genomic analysis revealed that the six S. equorum strains do not possess any of the virulence determinants for adhesion, enterotoxins and pathogenicity islands that are found in S. aureus. Nonetheless, strain C2014 induced haemolysis on sheep blood-supplemented agar and strain KM1031 displayed resistance to chloramphenicol, erythromycin, lincomycin and penicillin G 13 ; therefore, we focused on genetic analysis to explain the differences in phenotype among the S. equorum strains.
Haemolysis. The αand β-haemolysins are prevailing toxins of S. aureus, but their homologs were not identified in any of the six S. equorum genomes. However, three CDSs annotated as haemolysin, haemolysin III and haemolysin activation protein genes were identified in all of the genomes of five S. equorum strains, excluding strain KS1039 (Table 2). Strain KS1039 does not possess a predicted haemolysin gene. These three CDSs also exist in all of 15 S. equorum strains from cheeses (Supplementary Table S3).
The deduced amino acid sequence of annotated haemolysin gene (AVJ22-RS03205) from strain C2014 has 98% identity with the TlyC amino acid sequence of Streptococcus equi. TlyC was described as a haemolysin in Brachyspira hyodysenteriae and was considered to be a virulence factor contributing to the disease caused by this spirochete because tlyC-expressing E. coli displayed haemolytic activity in vitro 30 . Meanwhile, Carvalho et al. 31 concluded that TlyC of Leptospira does not exert a direct haemolytic effect but may contribute to Leptospira binding to the extracellular matrix during host infection. In addition, Turner and Helmann 29 reported that YhdP (TlyC) in Bacillus subtilis is a multidrug efflux protein and its amino acid sequence has 67% similarity with that of a membrane protein in Bacillus isronensis. These results suggest that the haemolysin homolog may not contribute to haemolysin activity directly but may instead be a membrane protein or enhance haemolytic activity.
The integral membrane protein gene hlyIII was identified from Bacillus cereus and Vibrio vulnificus; its involvement in haemolysis was verified by its expression in a nonhaemolytic E. coli strain 30,32 . However, its homolog found in Bacteroides fragilis was not linked to haemolytic activity 33 . Among the S. equorum strains tested for  Table 2. Potential virulence determinants identified in six S. equorum genomes. *Genes located in a plasmid.
haemolytic activity, only strain C2014 exhibited β-haemolytic activity, despite all strains having the haemolysin III gene, which means that this gene is not an independent determinant of β-haemolysis.
All of the six S. equorum strains possess a gene encoding a putative haemolysin activation protein composed of 43 amino acids. However, only strain UMC-CNS-924 has three genes encoding three distinct haemolysin activation proteins under the control of a promoter ( Table 2) and exhibited δ-haemolysis ( Supplementary Fig. S2B). The three putative haemolysin activation proteins presented ≤45.5% amino acid sequence identities with each other and ≤37.2% sequence identity with that of strain C2014 (Supplementary Table S4). Staphylococcus lugdunensis was reported to secrete three 43-amino-acid peptides with synergistic haemolytic activity, phenotypically similar to the δ-haemolysin of S. aureus, and their genes are located in an operon 30 . The identification of δ-haemolysin activity only in strain UMC-CNS-924 may be attributable to these three genes, while their homologs found in the other five strains may contribute to other functions.
Among the strains tested for haemolytic activity, only strain C2014 exhibited β-haemolytic activity ( Supplementary Fig. S2A). Comparative genomic analysis revealed an additional CDS only found in strain C2014, a putative haemolysin family calcium-binding region gene (AVJ22_RS14095) ( Table 2) and a unique gene harboured in the plasmid named pC2014-5 (12.0 kb). The product of this gene was reported to induce stabilization of the β-sheet structure of haemolysin and be enhanced by Ca 2+ binding 34 . Therefore, we assumed that the expression of the gene AVJ22_RS14095 may enhance the haemolytic activity of other gene products that cannot induce haemolysis on their own. Based on the SEED subsystem, the genes AVJ22_RS03205 and AVJ22_RS09420 in strain C2014 were classified as encoding magnesium/cobalt efflux protein and the membrane protein haemolysin III, respectively. Thus, we cautiously suggest that the gene AVJ22_RS09420 might be a determinant of the haemolysis activity, with the support of AVJ22_RS14095, in strain C2014.
Acquired antibiotic resistance. Based on the functional categories as determined by the SEED subsystem, an average of 65.5% of the annotated CDSs in the virulence, disease and defence category for the six S. equorum strains are predicted to be genes for resistance to antibiotics and toxic compounds (Fig. 1B). Meanwhile, the genes for resistance to antibiotics and toxic compounds in strains KM1031 and G8HB1 constitute 71.8% and 73.7% of the CDSs in this category, respectively.
Putative efflux pump genes for chloramphenicol, lincomycin, quinolone and multiple drugs were identified across the six S. equorum chromosomes as well as those of 15 S. equorum strains from cheeses ( Table 2,  Supplementary Table S3). KM1039 is the only strain that does not harbour the putative lincomycin-resistance gene lmrB. Conversely, the phenotype of lincomycin resistance was exhibited in strains KM1031 and UMC-CNS-924 (Fig. 3). The lincomycin-resistant S. equorum KM1031 harbours a plasmid encoding the lnuA gene, which it can transfer to Gram-positive recipients 11 . Strain UMC-CNS-924 also harbours the lnuA-encoding plasmid. The lincomycin resistance of both strains might thus have been acquired via horizontal transfer of the resistant plasmid. The six commonly identified putative efflux pump genes including the lmrB homolog may thus not function in lincomycin resistance.
Resistance to ciprofloxacin, a kind of quinolone, was not identified in any of the test strains. In our previous antibiotic resistance test of 126 S. equorum strains, two strains showed independent ciprofloxacin and ofloxacin resistance and both strains also exhibited multidrug resistance 14 . Their quinolone resistance might be an acquired trait and the chromosomal quinolone resistance homologs may be weakly related to the phenotypic quinolone resistance of S. equorum.
Among the test strains, only the strain KM1031 exhibited resistance to chloramphenicol as well as erythromycin. Comparative genomic analysis has highlighted the putative antibiotic ABC transporter ATP-binding protein gene (AWC34_RS11115) and the antibiotic biosynthesis monooxygenase gene (AWC34_RS01805) identified in KM1031 as the possible determinants of chloramphenicol and erythromycin resistance. Nguyen and Nguyen 35 reported that an E. coli transformant containing antibiotic ABC transporter ATP-binding protein homolog showed resistance to cefalotin, kanamycin, ampicillin, erythromycin and chloramphenicol. In addition, the antibiotic ABC transporter ATP-binding protein was reported to confer resistance to several antibiotics via a ribosomal protection mechanism 36 . Therefore, we assumed that the antibiotic ABC transporter ATP-binding protein gene may confer the chloramphenicol and erythromycin resistance to strain KM1031. However, the existence of an antibiotic ABC transporter ATP-binding protein gene (SEQU_RS24115) does not confer resistance to chloramphenicol and erythromycin to the strain UMC-CNS-924. The putative antibiotic ABC transporter ATP-binding protein of strain KM1031 showed 100% sequence identity to that of UMC-CNS-924 (Supplementary  Table S4), which means that the homologs do not contribute to the chloramphenicol and erythromycin resistance if their transcriptional regulators function properly. Uniquely, an antibiotic biosynthesis monooxygenase gene (AWC34_RS01805) was identified only in strain KM1031; this gene was reported to alter polyketide antibiotics such as macrolide antibiotics through oxidation 37,38 . Meanwhile, strain-specific possession of ABC transporter ATP-binding protein or antibiotic biosynthesis monooxygenase homologs in the 15 strains from cheeses hampered clarification of which gene determines the resistance (Supplementary Table S3). Interestingly, two IS6 family transposase genes (COG AWC34_RS01810 and AWC34_RS01800) were identified in the flanking regions of the antibiotic biosynthesis monooxygenase gene for strain KM1031. Antibiotic biosynthesis monooxygenase gene was also found in the contigs of strains RE2.35 and 908_10 and an IS6 family transposase gene (A4A32_12915) was identified upstream of the homolog in strain RE2. 35. The contig of strain RE2.35 composed of a pseudogene exists between the transposase and antibiotic biosynthesis monooxygenase genes. The contig of strain 908_10 composed of a pseudogene and antibiotic biosynthesis monooxygenase gene suggests the existence of transposase in the flanking region. Therefore, we cautiously assumed that AWC34_RS01805 in KM1031 inserted by the action of transposase conferred the chloramphenicol and erythromycin resistance in this strain.
The penicillin G-resistant strain KM1031 harbours a plasmid encoding a putative β-lactamase gene ( Table 2, Fig. 3). β-Lactamase homologs are found in the chromosomes of the other four strains, except strain KS1039. The putative amino acid sequence of the plasmid-encoded β-lactamase gene (AWC34_RS13020) shows 79.7-95.3% identity with those from homologs encoded chromosomally (Supplementary Table S4). Meanwhile, the plasmid-encoded β-lactamase has 97% identity with the Zn-dependent hydrolase of S. aureus (Supplementary  Table S5). This suggests that the chromosomal β-lactamase homologs contribute not to penicillin G resistance but to other functions. The absence of a β-lactamase homolog in strain KS1039 implies that this gene is not critical to the survival of this strain.
Strain UMC-CNS-924 was resistant to tetracycline and was shown to contain a plasmid-encoded gene whose product has 99% sequence identity with tetracycline resistance MFS (major facilitator superfamily) efflux pumps of Staphylococcus species (Table 2, Fig. 3, Supplementary Table S5). This gene may thus encode an efflux pump that directly promotes resistance to tetracycline.
It is already well known that the mecA gene, which encodes the low-affinity penicillin-binding protein PBP 2A, confers methicillin resistance 39 . The mecA gene was not identified in the six strains, and methicillin resistance was also not identified. Meanwhile, three genes annotated to encode methicillin resistance proteins were commonly identified in the six strains. Two putative methicillin resistance proteins show ≥88% sequence identity with the aminoacyltransferase of Staphylococcus species (Supplementary Table S5). For these two genes, there is the possibility of mis-annotation. The third gene has ≥66% sequence identity with FemA of the femAB operon. FemA and FemB are known as additional components for methicillin resistance, which enhance the methicillin resistance of MecA. Therefore, these three genes may not confer methicillin resistance to S. equorum strains.
In the current study, all plasmid-mediated antibiotic resistance genes were linked to antibiotic resistance. However, most of the putative antibiotic resistance genes encoded chromosomally were not linked to their expected resistance. The potential antibiotic resistance genes encoded chromosomally may not contribute to resistance owing to their low activity 40,41 . It is already well known that plasmid-mediated antibiotic resistance genes confer higher resistance than genes on chromosomes 42,43 . However, antibiotic-resistant strains with chromosome-generated adaptation have been reported to survive at a higher rate than strains without antibiotic resistance genes upon exposure to antibiotics 41,44 . Therefore, antibiotic resistance genes on chromosomes may allow organisms to survive under various conditions. Two-component systems. The success of S. aureus as a pathogen is in part due to the precise regulation of genes for survival in various environments. Two-component systems (TCSs) serve as a basic stimulus-response coupling mechanism to allow organisms to sense and respond to changes in the environment. S. aureus has been reported to possess 16 TCSs within its relatively small genome, two of which, saeRS and agrCA, are known to regulate virulence [45][46][47] . In the genomes of S. equorum strains, three types of putative TCS involved in the regulation of ion acquisition, cell-wall synthesis and nitrate reduction have been identified (Table 3). arlSR (yhcSR) was reported to be essential for survival under various environmental conditions 48 but was not reported to be involved in the regulation of virulence factors. arlSR TCS was also reported to regulate the biofilm formation of Staphylococcus epidermidis in an ica-dependent manner 49 . However, the ica gene involved in biofilm formation was not identified in the six S. equorum genomes. These results imply that S. equorum arlSR is involved not in the regulation of virulence genes but in adaptation to the environment. As another example, the vraSR TCS system was reported to be a positive modulator of cell wall biosynthesis 50 . Other TCSs in S. equorum have not been reported to be related to virulence factors. In this context, S. equorum may have little possibility of expressing the pathogenicity seen in S. aureus, but it can exhibit virulence simply by acquiring virulence determinants including those for haemolysis and antibiotic resistance.
Scientific RepoRts | 7: 5383 | DOI:10.1038/s41598-017-05918-5 Salt tolerance of C2014 and KS1039. Bacteria respond to hyperosmotic stress either by controlling the flux of ions across their cellular membrane or by accumulating osmolytes called compatible solutes 51 . In ion homeostasis, potassium plays a pivotal role and is the most abundant ion in the cytoplasm of bacteria 52 . Although CNS are frequently identified in foods with a high salt concentration, the mechanism behind their salt tolerance is not well understood. Conversely, a number of characteristics allowing S. aureus to survive osmotic stress have been reported. For example, osmoprotectants such as choline, glycine, betaine and proline accumulate in S. aureus in response to osmotic stress 53,54 . In addition, multiple genes, including the branched-chain amino acid transporter gene brnQ 55 and the arsenic operon regulatory gene arsR 56 , have been reported to cooperate in conferring salt tolerance to S. aureus. Furthermore, the involvement of a very large cell-wall protein, Ebh, in the tolerance to transient hyperosmotic pressure was reported 57 . The phospholipid cardiolipin, an important component of the cell membrane, was also reported to be necessary for the prolonged survival of S. aureus in high-salt conditions 58,59 . The six S. equorum strains also possess two cardiolipin synthetase genes (Supplementary  Table S6), as well as a sodium/potassium transport system (Supplementary Table S7), and all strains exhibited growth on TSA plates supplemented with 15% (w/v) NaCl. Meanwhile, two strains, KS1039 and C2014, exhibited growth at a NaCl concentration of 25%; the growth rate of KS1039 was slightly higher than that of C2014 under these conditions (Supplementary Fig. S3).
Strain C2014 possesses an ortholog (AVJ22_RS01775) of an ion transporter that encodes a protein having a ball domain and a potassium voltage-gated channel in its ORF ( Supplementary Fig. S4, Supplementary Table S7). The ball domain was shown to be responsible for the inactivation of voltage-gated ion channels 60,61 . DNA sequences encoding ball domains have been widely identified in Staphylococcus species as well as in a Bacillus species and its relatives; the potassium voltage-gated channel-encoding sequences are located in the downstream parts of these ball domains ( Supplementary Fig. S4). Strain KS1039 possesses two unique orthologs, SE1039_RS01900 and SE1039_RS01905, which encode a potassium voltage-gated channel and a protein with a ball domain, respectively; we supposed that these genes contribute to the high salt tolerance of strain KS1039 16 . The ion selectivity of potassium voltage-gated channels was reported to be associated with a conserved sequence motif TVGYG located in a re-entrant loop present in-between two predicted transmembrane regions 62 . This conserved sequence was identified in both AVJ22_RS01775 and SE1039-RS01900 genes.
To investigate the effect of the potassium voltage-gated channel and ball domain on salt tolerance, the SE1039_ RS01900/RS01905 and AVJ22_RS01775 genes were amplified and then cloned into the pGEM-T easy vector under the control of the T7 promoter. The resulting plasmids were designated as p1039P for the gene SE1039_RS01900, p1039B for the gene SE1039_RS01905, p1039BP for the genes SE1039_RS01905/RS01900 and p2014BP for the gene AVJ22_RS01775 ( Supplementary Fig. S4). Under IPTG induction, the effect of NaCl on the growth was pronounced when NaCl was applied at a concentration of 6% (Fig. 4). The transformant harbouring p1039P showed the highest growth, followed by the transformant harbouring p1039BP. E. coli harbouring p1039B showed much slower growth than the control containing pGEM-T easy vector. These results suggest that the ball domain downregulates the activity of voltage-gated ion channels. Higher salt tolerances were exhibited in strain KS1039 and the transformant harbouring p1039BP than in strain C2014 and the recombinant harbouring p2014BP, implying that salt tolerance can be increased when potassium voltage-gated channels and ball domains exist in separate ORFs.
The flanking regions of the potassium voltage-gated channel genes for strains C2014 and KS1039 do not provide any clues about if or where any insertion might have occurred. However, transposase genes are found at distant loci on both sides of the potassium voltage-gated channel genes (COG AVJ22_RS01765 and AVJ22_RS03035 in strain C2014; SE1039_RS03265 and SE1039_RS01160 in strain KS1039). We thus assume that strains C2014 and KS1039 acquired the potassium voltage-gated channel gene by a random insertion event, enabling them to survive in the high-salt conditions of jeotgal. The potassium voltage-gated channel gene allows S. equorum strains to survive in high-salt fermented foods with a NaCl concentration of over 15%.

Conclusions
The previously reported analyses of the three genomes of S. equorum strains Mu2, UMC-CNS-924 and KS1039 revealed that these three strains do not possess any clear virulence factors found in the well-known pathogen S. aureus 10,16,17 . Genomic analysis should help to confirm the lack of known safety hazard determinants in strains, but cannot guarantee the safety of candidate industrial strains. This comparative genomic analysis identified the strain-specific determinants of S. equorum strains, which are linked to phenotypic characteristics including virulence and salt tolerance. The identification of factors related to phenotypic safety hazards should help  Table 3. Putative two-component systems identified in six S. equorum genomes.
in guaranteeing the safety of industrial strains. In this context, this study strongly supports the safety of strain KS1039, which was selected as a starter candidate for jeotgal fermentation 13 because it does not possess any of the functional virulence determinants identified in the other strains. Additionally, this study identified the gene that contributes to the strain-specific high-salt tolerance of S. equorum strains. This study supports the usefulness of comparative genomics in the safety assessment and functional analysis of industrial strains including starter candidates for food fermentation.