A total of 441 NAS isolates (Table 1) were selected for WGS from 5507 isolates obtained from bovine staphylococcal IMI stored at the Mastitis Pathogen Collection of the Canadian Bovine Mastitis and Milk Quality Research Network (CBMQRN). Selection was randomized over the 25 NAS species previously identified (Mellmann et al., 2006; Condas et al., in review) which originated from 87 herds distributed over four geographic regions encompassing Atlantic Canada (Nova Scotia, Prince Edward Island, New Brunswick), Central Canada (Québec and Ontario), and Western Canada (Alberta) (Reyher et al., 2011). Inclusion criteria of isolates were: (1) inclusion of 68 NAS isolated from clinical mastitis samples; (2) inclusion of 26 NAS isolates with multi-drug resistance (MDR); (3) for uncommon species (defined as <20 unique isolations at cow level), inclusion of one arbitrarily selected isolate of that species per cow; (4) one isolate per cow for all other species, until achieving the final number. Our criteria allowed two isolates of the same NAS species isolated from the same cow to be included (i.e., S. chromogenes isolated from clinical mastitis and subclinical mastitis with MDR in the same cow), although the majority of isolates originated from different cows (Table 1).

All isolates obtained were grown on 5% defibrinated sheep blood agar plates (BD Diagnostics, Mississauga, ON, Canada) with subsequent incubation at 35°C for 24 h to yield single colonies. For each of these isolates, colonies were picked from blood agar plates and suspended in Bacto Brain Heart Infusion broth (BD Diagnostics), and incubated at 37°C overnight. Genomic DNA for all isolates was extracted with DNeasy Blood and Tissue Kit (Qiagen, Toronto, ON, Canada), according to the corresponding protocol for Gram-positive bacteria. Quality and quantity of DNA was checked using a NanoVue plus spectrophotometer (GE Healthcare Life Sciences, Mississauga, ON, Canada) and the Qubit 2.0 fluorometer (Invitrogen, Burlington, ON, Canada). Each DNA sample was diluted to a final concentration of 0.2 ng/μl. Sequencing of these samples was performed using the Illumina MiSeq platform (Illumina, San Diego, CA, USA); DNA libraries for sequencing were prepared using a Nextera XT DNA library preparation kit (Illumina, San Diego, CA, US). All sequencing steps, including cluster generation, paired-end sequencing (2 × 250 bp), and primary data analysis for quality control, were performed on the instrument.

Sequence reads obtained from the MiSeq platform were checked for poorly sequenced regions and Illumina adapters sequences were trimmed using cutadapt (Martin, 2011), implemented in Trim Galore! 0.4.0 (with default parameters). After trimming sequences, filtered reads were assembled into contigs using the de novo assembly program SPAdes version 3.6.0 (Nurk et al., 2013), employing built-in error correction and default parameters. To determine average coverage of sequencing, reads were mapped back to the assembled genome using BWA 0.7.12-r1039 (Li and Durbin, 2010) and depth of sequencing for each contig was plotted using BEDTools (Quinlan and Hall, 2010). Genome annotations and gene predictions for contigs larger than 200 bp were performed with Prokka 1.11 (Seemann, 2014), using the provided Staphylococcus database. Quality of the assembled genomes and assembly metrics was determined using Quast (Gurevich et al., 2013). The entire genome assembly process was automated using the Snakemake workflow engine (Köster and Rahmann, 2012). The data was submitted to NCBI under BioProject ID PRJNA342349.

Phylogenetic analyses were performed on 24 highly conserved housekeeping proteins (Supplementary Table 1). Sequences for each of these 24 proteins for all NAS species were retrieved from the genomes using BLAST+2.2.31 (Camacho et al., 2009). The sequence for Macrococcus caseolyticus, an outgroup species, was downloaded from NCBIs GenBank database. Multiple sequence alignments for these proteins were created using MUSCLE v3.8.31 (Edgar, 2004). The resulting alignments were used for phylogenetic analysis. Phylogenetic trees, based on 100 bootstrap replicates, were constructed by employing Maximum-Likelihood (ML), Maximum-Parsimony (MP), and Neighbor-Joining (NJ) methods using MEGA 6.0 (Tamura et al., 2013). Evolutionary distances for ML and NJ methods were computed using a JTT matrix-based model (Jones et al., 1992). Maximum-Parsimony trees were obtained using the Subtree-Pruning-Regrafting (SPR) algorithm (Nei and Kumar, 2000).

Phylogenetic trees were constructed based on full-length sequences of 16S rRNA, hsp60, rpoB, sodA, and tuf genes. Full-length sequences of these genes for NAS species were obtained using BLAST+2.2.31 (Camacho et al., 2009). Multiple sequence alignments for each of gene were created using MUSCLE v3.8.31 (Edgar, 2004). Maximum-Likelihood, MP and NJ trees based on these sequence alignments were created using 100 bootstrap replicates in MEGA 6.0 (Tamura et al., 2013).

Multilocus sequence analysis was performed on nucleotide sequences of the 16S rRNA, hsp60, rpoB, sodA, and tuf genes. Individual gene alignments were manually concatenated to create a combined dataset of these five genes. Poorly aligned regions from this concatenated alignment were removed using Gblocks 0.92 (Castresana, 2000) with default settings except the allowed gap position parameter which was changed to 0.5 (50%). Maximum-Likelihood, MP and NJ trees based on 100 bootstrap replicates of this dataset were constructed using MEGA 6.0 (Tamura et al., 2013). For all trees, evolutionary distances for ML, MP, and NJ methods were computed using the General Time Reversible model (Nei and Kumar, 2000), Subtree-Pruning-Regrafting (Nei and Kumar, 2000), and Kimura 2-parameter model (Kimura, 1980), respectively. Codon positions included were 1st+2nd+3rd+Noncoding. All positions containing gaps and missing data were eliminated from final alignments. All trees were rooted using Macrococcus caseolyticus, and gene sequences for this species were downloaded from NCBI GenBank database.

A genome-based phylogenetic tree of 441 NAS isolates was constructed using the published pipeline PhyloPhlAn (Segata et al., 2013). Briefly, the PhyloPhlAn approach is based on the use of 400 ubiquitous and phylogenetically informative proteins. Orthologs of these proteins in NAS genomes were detected using USEARCH v5.2.32 (Edgar, 2010). Multiple sequence alignments of these proteins were generated using MUSCLE v3.8.31 (Edgar, 2004). A final concatenated dataset containing 4231 aligned amino acid positions was generated. Final tree construction was performed using FastTree version 2.1 (Price et al., 2010). To determine effects of various evolutionary models on the resulting tree, Maximum-Likelihood, MP, and NJ trees based on 100 bootstrap replicates of the final alignment were also constructed using MEGA 6.0 (Tamura et al., 2013). For ML and NJ trees, evolutionary distances were computed using a JTT matrix-based model (Jones et al., 1992), whereas the Subtree-Pruning-Regrafting algorithm (Nei and Kumar, 2000) was used for MP tree construction.

A phylogenetic tree of all NAS isolates, rooted using Macrococcus caseolyticus, was created based on the core genome of the bovine NAS group. The core set of NAS proteins were identified using the UCLUST algorithm (Edgar, 2010). Protein families with at least 30% sequence identity and 50% sequence length in 441 NAS isolates were considered core. However, protein families present in <80% of the input genomes were excluded from further analysis. Proteins families which contained potential paralogous sequences (duplicated sequence in same genome) were also excluded from further analysis. Each protein family was individually aligned using MAFFT 7 (Katoh and Standley, 2013). Aligned amino acid positions which contained gaps in more than 50% of genomes, were excluded from further analysis. Remaining amino acid positions were concatenated to create a combined dataset consisting of 128,080 aligned amino acids. A maximum-likelihood tree based on this alignment was constructed using FastTree 2.1 (Price et al., 2010), using the Whelan and Goldman substitution model (Whelan and Goldman, 2001) and JTT matrix-based model (Jones et al., 1992).

Genome-based phylogenetic trees (PhyloPhlAn, ML, MP, NJ) and Core-Genome-Tree (CGT) were visually inspected and compared. Single gene and protein trees were compared with each other and with CGT. Topological differences among them were computed using Robinson-Foulds (RF) distance matrix (Robinson and Foulds, 1981; Makarenkov and Leclerc, 2000), implemented as webserver; T-REX (Tree and reticulogram REConstruction, www.trex.uqam.ca; Boc et al., 2012). To facilitate comparisons of our results, RF distance scores were normalized using the maximum possible distance between two trees, calculated according to 2(N-3) where N represents the total number of taxa on a tree (Kupczok et al., 2010; Bernard et al., 2016). We designated these as normalized-Robinson-Foulds (nRF) scores, with values ranging from 0 (0%) to 1 (100%), and nRF = 0 indicating that topologies of two trees under investigation are congruent. Consequently, higher nRF score indicate a low level of congruence (or high level of incongruence) between two tree topologies.

Five phylogenetic trees were constructed based on WGS data of NAS isolates. The first of these trees was constructed using a core genome dataset of NAS species (Figure 1), whereas the second tree (Supplementary Figure 1) was constructed using PhyloPhlAn (Segata et al., 2013). The same dataset from aligned sequences, obtained from PhyloPhlAn was also used to construct the ML (Supplementary Figure 2), MP, (Supplementary Figure 3), and NJ trees (Supplementary Figure 4). Within these trees, species of the NAS group branched into five distinct monophyletic and statistically well-supported clades (100% nodal support) in all phylogenetic reconstructions. We referred to them as Clade A, Clade B, Clade C, Clade D, and Clade E, which contained 3, 3, 1, 8, and 10 NAS species, respectively (Figure 1). The phylogeny of NAS species was also constructed by excluding 26 multidrug resistant isolates, as they could have introduced selection bias. However, no differences in phylogeny of the NAS species were seen (data not shown).

In all phylogenetic trees, Clade A branched as first or as the most ancient divergence of NAS group and was composed of three NAS species: S. sciuri, S. fleurettii, and S. vitulinus. Within this clade, S. sciuri formed a distant cluster basal to S. fleurettii and S. vitulinus, which appeared as sister taxa to each other (Figure 1 and Supplementary Figures 1–4). The second or Clade B of NAS group contained S. chromogenes, S. hyicus, and S. agnetis, with S. hyicus appearing as sister taxon to S. agnetis and S. chromogenes forming the basal group (Figure 1 and Supplementary Figures 1–3). However, there was a discrepancy in branching order among these species in the phylogenetic tree constructed using NJ method (Supplementary Figure 4), which placed S. chromogenes as sister taxon to S. hyicus and displayed S. agnetis branching separately. The next lineage of NAS to diverge was constituted by all isolates of a single species S. simulans (Clade C). All 42 sequenced isolates of this species diverged independently from other NAS members by a long branch, forming a strongly supported clade (Figure 1 and Supplementary Figures 1–4).

Clade D contained eight species comprised of S. hominis, S. devriesei, S. haemolyticus, S. pasteuri, S. warneri, S. epidermidis, S. caprae, and S. capitis (Figure 1 and Supplementary Figures 1–4). Staphylococcus aureus, a major pathogen of IMI, also shared common ancestor with species of Clade D. Phylogenetic placement of S. aureus is indicated on the tree (Figure 1). The eight NAS species of this clade were split into three clusters. The first cluster (D1) was composed of S. hominis, S. devriesei, and S. haemolyticus. Within this cluster, S. hominis appeared as a deeper branching species, whereas S. devriesei and S. haemolyticus branched together, representing divergence from a common ancestor (Figure 1 and Supplementary Figures 1–4). The second cluster (D2) was formed by two NAS species, S. pasteuri and S. warneri, grouping together with 100% support in all genomic phylogenies (Figure 1 and Supplementary Figures 1–4). Cluster D3 displayed a closer relationship of S. caprae to S. capitis, whereas S. epidermidis appeared as the basal lineage of this cluster. Monophyletic branching of these clusters within Clade D was supported by 100% bootstrap scores in all trees. The branching order and phylogenetic placement of these cluster groups was the same in all phylogenies (Figure 1 and Supplementary Figures 1–4).

The fifth clade or Clade E was the most recently diverged clade of NAS and contained 10 species, comprised of S. auricularis, S. kloosii, S. arlettae, S. gallinarum, S. succinus, S. nepalensis, S. cohnii, S. equorum, S. saprophyticus, and S. xylosus (Figure 1 and Supplementary Figures 1–4). Within this clade, S. auricularis was the deepest branching species, followed by a cluster group shared by S. kloosii and S. arlettae. Next, S. gallinarum and S. succinus diverged as two well-supported separate lineages (Figure 1 and Supplementary Figures 1–4). Within the remaining species, there was a closer relationship of S. nepalensis to S. cohnii and S. saprophyticus to S. xylosus (Figure 1 and Supplementary Figures 1–4). All isolates of S. equorum branched together as sister taxa to the S. xylosus—S. saprophyticus group.

Phylogenetic relationships among NAS species were also inferred using 16S rRNA, hsp60, rpoB, sodA, and tuf gene sequences separately and as a multilocus data set. The ML, MP, and NJ trees produced from these genes were compared with each other and with the CGT. The nRF distance-based comparisons among ML, NJ single trees and CGT (Tables 2A,B) and single genes MP trees and CGT (Table 2C) are presented in Table 2. Consistent with the CGT, multilocus ML (Figure 2) and NJ (Supplementary Figure 5) trees displayed branching of NAS species into five distinct clades (Clade A, Clade B, Clade C, Clade D, and Clade E). The MP-multilocus tree (Supplementary Figure 6), however, had a different pattern of clade divergence, and displayed branching of S. simulans with S. auricularis within Clade E. Hence, Clade C was absent from MP-based phylogeny. Except for the multilocus MP tree, species composition within each clade was identical to the CGT. However, branching order within Clade D and Clade E differed among each other and with CGT. Within Clade D, in the multilocus MP tree, the cluster group containing S. pasteuri and S. warneri suggested a deeper branching cluster compared to other species of this clade (Supplementary Figure 6). The notable difference within clade E was branching of S. succinus and S. equorum as sister taxa in all multilocus trees (Figure 2 and Supplementary Figures 5, 6).

Phylogenetic trees (ML, MP, and NJ) based upon 16S rRNA gene (Figure 3 and Supplementary Figures 7, 8), hsp60 (Supplementary Figures 9–11), rpoB (Supplementary Figures 12–14), sodA (Supplementary Figures 15–17), and tuf (Supplementary Figures 18–20) gene sequences were constructed. The nRF distance-based analyses for ML trees (Table 2A), NJ trees (Table 2B), and MP trees (Table 2C) based upon these five genes, indicated that these trees differed from each other and from the CGT (Table 2). Within the 16S rRNA gene based phylogenies, the only reliably identified relationship was the composition and branching order of S. sciuri, S. fleurettii, and S. vitulinus within Clade A (Figure 3 and Supplementary Figures 7, 8). All other species had varying patterns of evolutionary placement. Although S. simulans branched distinctly in ML and NJ trees, the order of its divergence was different compared to the CGT. In these cases, S. simulans branched after divergence of Clade A from the rest of the NAS species (Figure 3 and Supplementary Figure 8), as opposed to branching after separation of Clades A and B in the CGT (Figure 1 and Supplementary Figures 1–4). The branching of S. simulans after Clade A had 83 and 80% statistical support in ML (Figure 3) and NJ (Supplementary Figure 8) trees, respectively. In the 16S-rRNA-MP tree, S. simulans clustered with other NAS species (Supplementary Figure 7). Within Clade B, S. auricularis branched with S. chromogenes, S. agnetis, and S. hyicus (Figure 3). However, Clade D was not recovered as a monophyletic lineage. The cluster group containing S. hominis, S. haemolyticus, and S. devriesei appeared as separate taxa distinct from the rest of Clade D species (Figure 3). Except for S. auricularis, which clustered with clade B, species composition was invariable within Clade E. However, there was discordance among branching of these species. Additionally, phylogenetic positioning of S. kloosii and S. equorum was not resolved properly within Clade E (Figure 3).

Phylogenies constructed using hsp60 gene sequence differed significantly from each other, both in species composition and their hierarchy of branching pattern. For three hsp60 trees (ML, MP, and NJ), Clade A branches as monophyletic group with 100% bootstrap scores (Supplementary Figures 9–11). However, species of other clades were intermixed. The ML and NJ phylogenetic trees constructed using rpoB gene sequence revealed the same major clades. However, in both trees, S. simulans branched with S. chromogenes, S. agnetis, and S. hyicus, species of Clade B. Within Clade D, in both the ML and NJ trees, S. pasteuri and S. warneri formed a basal branching group compared to other Clade D species. Within Clade E, there were similar branching patterns among ML and NJ trees; the only difference was branching of S. succinus after the S. cohnii—S. nepalensis pair, as opposed to its divergence after S. gallinarum in CGT. Except for Clade A, the MP tree (Supplementary Figure 13) based on the rpoB gene, had a largely different pattern of evolutionary relationships among NAS species. Phylogenetic trees based on the sodA gene only revealed formation of Clades A and B, similar to CGT (Supplementary Figures 15–17), whereas remaining species of other clades (Clade C, D, and E) had a polyphyletic grouping pattern. For ML and NJ method-based phylogenies obtained from tuf gene, NAS species branched into the same five distinct clades as the CGT (Supplementary Figures 18, 20). However, in the MP (tuf) tree, S. simulans grouped with species of Clade D (Supplementary Figure 19). Hierarchy and clustering patterns among members of Clade D and Clade E differed among various phylogenies, based on the tuf gene.

Evolutionary relationships of NAS species were also investigated, based on 24 housekeeping protein sequences (Supplementary Table 1). Phylogenies obtained from each protein sequence were examined by various methods of evolutionary estimations (ML, MP, and NJ). Three trees (ML, MP, and NJ) produced from each single protein, were compared with each other and with CGT. The nRF distances between each of ML, NJ, and MP single protein trees and CGT are presented (Tables 3–5, respectively). Differences and similarities in species composition and branching order among species of various clades among trees are presented (Table 6).

Species composition and branching order for Clade A were retrieved in most phylogenetic reconstructions (Table 3) except MP trees [Elongation factor Tu (EF-Tu) and ribosomal protein S11 (RibS11)] and NJ tree [ribosomal protein S11 (RibS11)] which could not recover Clade A. Among protein trees that depicted branching of Clade A as monophyletic group, there were branching order differences among species of Clade A in MP trees, based upon cytidine diphosphate (CDP) diglyceride synthetase (CdsA), ribosomal protein S2 (RibS2) and RNA polymerase subunits α and β protein sequences (RpoA and RpoB). Discordance in branching of Clade A species was also apparent in NJ trees constructed from RecA, ribosomal protein S2 (RibS2) and ML tree reconstructed from CDP diglyceride synthetase (CdsA), elongation factor Tu (EF-Tu), cell division protein FtsY, RecA, ribosomal protein S2 (RibS2), and RNA polymerase α subunit (RpoA) protein sequences.

Similar to Clade A, most phylogenetic estimation methods produced Clade B as an independent and distinct monophyletic group, except in MP trees built from cell division protein FtsY, gyrase B (GyrB), ribosomal protein L2 (RibL2), RibS11, and RNA polymerase β′ (RpoC) subunit protein sequences. Within Clade B, dissimilarities in branching among its species were seen in MP [serine hydroxymethyl transferase (GlyA) and gyrase A (GyrA)], NJ [ribosomal protein L15 (RibL15)] and ML trees [cell division protein FtsY, serine hydroxymethyl transferase (GlyA), gyrase B (GyrB), and HSP60]. Clade C was considered accurate only when it branched after Clades A and B in phylogenetic trees. Several protein sequences, namely, CDP diglyceride synthetase (CdsA), gyrase A (GryA), HSP70 and ribosomal protein S8 (RibS8), produced precise branching of Clade C in ML, MP, and NJ trees. Accurate branching of Clade C was also recovered in MP trees constructed from ribosomal protein L1 (RibL1) and RibS2, preprotein translocase SecA, SecY, and NJ trees, based on cell division protein FtsY, ribosomal protein L2, L15, and S2 (RibL2, RibL15, and RibS2) and PcrA helicase. ML trees constructed using sequences from cell division protein FtsY and preprotein translocase SecY also produced correct branching of Clade C (identical to CGT).

Species composition for Clade D, matching to CGT, was obtained in ML, MP, and NJ phylogenies constructed using sequences from elongation factor G (EF-G), cell division protein FtsY, Gyrase A (GyrA), RecA, RNA polymerase subunit β (RpoB), preprotein translocase SecY and PcrA helicase. However, identical branching patterns among species of Clade D were only observed in MP (GyrA) and NJ trees constructed from serine hydroxymethyl transferase (GlyA), preprotein translocase SecA and PcrA helicase protein sequences. Accurate species composition within Clade E was recovered in all phylogenetic estimations (ML, MP, and NJ) obtained from CDP diglyceride synthetase (CdsA), cell division protein FtsY, serine hydroxymethyl transferase (GlyA), gyrase A (GyrA), HSP60, RNA polymerase subunit β (RpoB), and preprotein translocase SecY. However, the matching branching order among the species of Clade E to CGT was only present in one MP tree [Serine hydroxymethyl transferase (GlyA)] and three NJ trees [HSP60, RNA polymerase subunit β (RpoB), and preprotein translocase SecA].

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.