# Phylogenomics picks out the par excellence markers for species phylogeny in the genus *Staphylococcus*.

Lucia Graña-Miraglia<fnref n="1" />, César Arreguín-Pérez<fnref n="2" />, Gamaliel Lopez-Leal<fnref n="3" />, Alan Muñoz<fnref n="1" />, Angeles Pérez-Oseguera

¹, Estefan Miranda-Miranda<fnref n="2" />, Raquel Cossío-Bayúgar<fnref n="2" />, Santiago Castillo-Ramírez Corresp. <fnref n="1" />

Corresponding Author: Santiago Castillo-Ramírez

Email address: iago@ccg.unam.mx

Although genome sequencing has become a very promising approach to conduct microbial taxonomy, few labs have the resources to afford this especially when dealing with data sets of hundreds to thousands of isolates. The goal of this study was to identify the most adequate loci for inferring the phylogeny of the species within the genus *Staphylococcus*; with the idea that those who cannot afford whole genome sequencing can use these loci to carry out species assignation confidently. We retrieved 177 orthologous groups (OGs) by using a genome-based phylogeny and an average nucleotide identity analysis. The top 26 OGs showed topologies similar to the species tree and the concatenation of them yielded a topology almost identical to that of the species tree. Furthermore, a phylogeny of just the top 7 OGs could be used for species assignment. We sequenced four staphylococcus isolates to test the 26 OGs and found that these OGs were far superior to commonly used markers for this genus. On the whole, our procedure allowed identification of the most adequate markers for inferring the phylogeny within the genus *Staphylococcus*. We anticipate that this approach will be employed for the identification of the most suitable markers for other bacterial genera and can be very helpful to sort out poorly classified genera.

<footnote n="1">¹ Programa de Genómica Evolutiva, Centro de Ciencias Genómicas, Universidad Nacional Autónoma de México, Cuernavaca, Morelos, Mexico</footnote>

<footnote n="2">² Centro Nacional de Investigación Disciplinaria en Parasitología Veterinaria del Instituto Nacional de Investigaciones Forestales, Agrícolas y Pecuarias, Jiutepec, Morelos, Mexico</footnote>

<footnote n="3">³ Departamento de Microbiología Molecular, Instituto de Biotecnología, Universidad Nacional Autónoma de México, Cuernavaca, Morelos, Mexico</footnote>
1 Phylogenomics picks out the par excellence markers for species phylogeny in the genus

2 *Staphylococcus*.

3

4

5 Authors

6 Lucía Graña-Miraglia¹, César Arreguín-Pérez², Gamaliel López-Leal³, Alan Muñoz¹, Ángeles Pérez-

7 Oseguera¹, Estefan Miranda-Miranda², Raquel Cossío-Bayúgar² and Santiago Castillo-Ramírez¹*.

8

9 ¹Programa de Genómica Evolutiva, Centro de Ciencias Génomicas, Universidad Nacional Autónoma de México, Cuernavaca, Morelos 62210, México.

11 ²Centro Nacional de Investigación Disciplinaria en Parasitología Veterinaria del Instituto Nacional de Investigaciones Forestales, Agrícolas y Pecuarias, Carr. Fed. Cuernavaca-Cuautla No. 8534, Jiutepec, Morelos, 62550 México.

14 ³Departamento de Microbiología Molecular, Instituto de Biotecnología, Universidad Nacional Autónoma de México, Cuernavaca, México.

16

17

18 *Correspondence: Santiago Castillo-Ramírez, iago@ccg.unam.mx

19

20

21

22

23

24

25
26 ABSTRACT

27 Although genome sequencing has become a very promising approach to conduct microbial taxonomy,
28 few labs have the resources to afford this especially when dealing with data sets of hundreds to thousands
29 of isolates. The goal of this study was to identify the most adequate loci for inferring the phylogeny of
30 the species within the genus *Staphylococcus*; with the idea that those who cannot afford whole genome
31 sequencing can use these loci to carry out species assignation confidently. We retrieved 177 orthologous
32 groups (OGs) by using a genome-based phylogeny and an average nucleotide identity analysis. The top
33 26 OGs showed topologies similar to the species tree and the concatenation of them yielded a topology
34 almost identical to that of the species tree. Furthermore, a phylogeny of just the top 7 OGs could be used
35 for species assignment. We sequenced four staphylococcus isolates to test the 26 OGs and found that
36 these OGs were far superior to commonly used markers for this genus. On the whole, our procedure
37 allowed identification of the most adequate markers for inferring the phylogeny within the genus
38 *Staphylococcus*. We anticipate that this approach will be employed for the identification of the most
39 suitable markers for other bacterial genera and can be very helpful to sort out poorly classified genera.
40

40

41

42 INTRODUCTION

43 The study of ecology and evolution has been substantially transformed by the omics technologies.

44 Remarkably, the use of these technologies has allowed essential questions in evolutionary biology to be
45 addressed and has advanced our knowledge in many biological processes (López-Leal, Tabche et al.

46 2014, Joseph, Marti et al. 2015, Joseph, Cox et al. 2016). The precise identification of the different
47 species within any given genus is highly valuable for many branches of microbiology. For instance, as far
48 as clinical microbiology is concerned, this is instrumental in establishing which species are human/animal
49 pathogens and which are just regular commensal organisms; whereas in microbial ecology species
50 assignation is very helpful in defining the niche-range of the different species (Becker, Margos et al.

51 2016).
52

Although whole-genome sequencing is the best method for genotyping bacterial isolates and, therefore, to conduct species assignation/identification, this approach is not yet affordable for routine use in many laboratories all over the world. This is especially true for low and middle-income countries in which the amount of money invested in science is not as much as in developed countries. For this particular purpose, and from a theoretical point of view, orthologous genes are the perfect candidates for inferring the phylogeny of species, as they should reflect the species tree of the taxa considered. As per definition, orthologous genes are the ideal candidates to track the sequence of past of speciation events within a given lineage (Fitch 2000). Therefore, in order to infer the phylogeny of the species, a clear-cut ascertainment of orthologous genes is of paramount importance and phylogenomic pipelines can be implemented to try to identify the potential orthologous genes. Orthologous relationships are context-sensitive and a gene family could be mainly composed of orthologous genes at one taxonomic level but if one considers a higher taxonomic level many more non-orthologous genes will appear. Furthermore, orthologous relationships could involve one to many relationships (Gabaldón and Koonin 2013). The word co-orthologue was coined to describe that exact situation in which a genome has more than one orthologous gene (Gabaldón and Koonin 2013). Hence, from a practical point of view (i.e. operational implementation), orthologous genes with just one gene per genome among the taxa considered should be the best markers to delineate the history of the species.

70

The genus *Staphylococcus* has a few dozens of species of Gram-positive bacteria, which are commensals colonizing the skin and mucous membranes of some mammals and birds. However, some of these species have a clear clinical and economic relevance, as they are a frequent cause of infection in humans, livestock, and domestic animals. Although the genus is very well known for the human opportunist pathogen *S. aureus*, which is famous worldwide as major source of nosocomial infections (Challagundla, Reyes et al. 2018, Frisch, Castillo-Ramirez et al. 2018), there are some other species such as *S. epidermidis*, *S. lugdunesis*, *S. saprophyticus* and *S. schleiferi* that have been associated with human
78 infections. Several molecular-based methods have been introduced to carry out species identification
79 within the genus *Staphylococcus* (Kwok and Chow 2003, Ghebremedhin, Layer et al. 2008, Sasaki,
80 Tsubakishita et al. 2010, Lamers, Muthukrishnan et al. 2012). PCR along with sequence analysis of the
81 genes *dnaJ*, *tuf*, *sodA*, *rpoB*, *hsp60* and *nuc* have been used to differentiate *Staphylococcus* species (Kwok
82 and Chow 2003, Ghebremedhin, Layer et al. 2008, Sasaki, Tsubakishita et al. 2010, Lamers,
83 Muthukrishnan et al. 2012). Although they have been very useful and demonstrate a better resolution
84 than the 16S rRNA gene (Ghebremedhin, Layer et al. 2008), these genes show different amounts of
85 genetic diversity and, therefore, varying levels of discriminatory power for the different species. Thus,
86 depending on the gene and the species, low-level resolution or even mis-identification can occur. Clearly,
87 an ideal solution would be to conduct whole-genome sequencing of all the isolates considered to tell apart
88 the different species. In terms of genome sequences, *Staphylococcus* is one of the few bacterial genera for
89 which many genomes are available.

90

91 In this study, we address the question of which genes are the best inter-species markers (*bona fide*
92 orthologous genes) for the genus *Staphylococcus* and, by applying a phylogenomic approach, we provide
93 a list of the top candidates for species assignation within this genus. Underfunded groups could use these
94 top candidates to accurately conduct species assignation without the necessity of using whole-genome
95 sequencing. Furthermore, this approach could be used for many other species/genera to identify the most
96 suitable markers for inferring the phylogeny of the taxa under investigation.

97

98 METHODS

99 Genomes and homologous groups

100 We downloaded 265 publically available complete genomes (see Supplementary Table 1). These cover 46
101 different species from the genus *Staphylococcus* and are a good representation of the host range within
102 this genus. We ran CheckM (Parks, Imelfort et al. 2015) on these genomes to discard poorly sequenced
103 genomes (i.e. with contamination or incomplete) and only one of them (marked in red in Supplementary
Table 1) did not pass the criteria (≥95% complete and with ≤5% contamination) and thus were not included in the rest of the analyses. We also sequenced the genomes of 4 bacterial isolates (see Supplementary Table 2) obtained from hemolymph or hypostome exudates from adult ticks showing signs of bacterial infection and reared over experimentally infested bovines in the Centro Nacional de Investigación Disciplinaria en Parasitología Veterinaria (CENID-PAVET-INIFAP), Jiutepec, Morelos, México. Of note the 4 isolates sequenced were confirmed to be Staphylococcus spp by biochemical tests and 16S rRNA gene sequence analysis. The isolate INIFAP 005-08 was initially identified as *S. saprophyticus* by its physical and biochemical characteristics (Gram-positive coccus bacteria, catalase positive; novobiocin resistant, absence of coagulase, gelatinase and caseinase activity). This isolate was also able to produce acids by fermenting glycerol, lactose, D(+) mannose, sucrose and turanose; however it was unable to use L(+) arabinose or L(+) lactose. This isolates was classified as *S. xylosus* as per 16S rRNA gene analysis. The isolates INIFAP 002-16 and INIFAP 004-15 were identified as Gram-positive coccus bacteria, catalase positive and ribotyped by the 16S rRNA gene, which exhibited positive identity for *S. xylosus*. INIFAP 009-16 was identified as Gram-positive coccus bacteria, catalase-positive and identified as *S. succinus* using the 16S rRNA gene. This same isolate was identified as *S. carnosus* using the API20E (bioMérieux 2010) system. The Nextera XT DNA Library Prep Kit was used for the sequencing libraries and Agilent High Sensitivity DNA Kit was employed for quality control. The isolates were sequenced using an Illumina MiSeq platform, with a 2X250 bp configuration; the genome sequencing was conducted at Instituto Nacional de Medicina Genómica (http://www.inmegen.gob.mx) in Mexico City. Prior to assembling the genomes, we employed the SolexQA v.3.7.1 program (Cox, Peterson et al. 2010) to trim the reads, employing the default parameters. We used Velvet version 1.2.09 (Zerbino and Birney 2008) and Spades v3.11.0 (Bankevich, Nurk et al. 2012) to carry out *de novo* assembly setting the option careful and we tested the following k-mer sizes: 51,61,71,81,91,101,111,121,127. We did not consider contigs smaller than 300 bp in the final assemblies. In all cases, Spades outperformed Velvet in terms of the number of contigs and the N50 statistics. Again, we used CheckM (Parks, Imelfort et al. 2015) to evaluate the quality of the genomes
130 sequenced by us; notably all the 4 isolates were shown to be ≥95% complete and with ≤5%

131 contamination and therefore reliable for downstream analyses. The details of the genomes assemblies are
132 provided in Supplementary Table 2. These seem to be good genome assemblies as judged by the high
133 median coverage and the low number of contigs for each one of the newly sequenced isolates. These 4
134 genomes have been submitted to the GenBank and have the following accession numbers:

135 PIZQ01000000, PIZN00000000, PIZP00000000, PIZO00000000 (BioProject PRJNA421192). Then, we
136 employed PROKKA v1.11 (we set the genus option to *Staphylococcus*) to annotate all the genomes
137 sequences (newly sequenced ones and the publically available) to have a consistent annotation for all of
138 them. To get the orthologous groups (OGs) with one-to-one relationships among the species, we only
139 considered single gene families (SGF); these were constructed running BLASTP searches with an e-value
140 of 1.0 e⁻³⁰ between the genome of *S. aureus* TW20 and the rest of the genomes. We kept all the cases
141 where there was only one hit per genome and requiring that the seed from TW20 and the hit aligned ≥ 60
142 % of their lengths and were ≥ 45 % identical - the rationale behind these two criteria was to make sure
143 that we had whole genes and not just domains. Then, for the all the SGF we constructed DNA alignments
144 in frame via the program Fast statistical alignment version 1.5.9 (Bradley, Roberts et al. 2009), setting the
145 option --nucprot that align nucleotide sequences taking into account the protein space. We conducted
146 recombination analysis on each of the SGFs using PhiTest (Bruen, Philippe et al. 2006) that was
147 implemented via the PhiPack program, using a window size of 50 bp. Additionally, we also determined
148 the nucleotide diversity for each one of the SGFs using R pegas library function nuc.div() with the default
149 parameters. We carried out a function enrichment analysis as follows: first, we conducted a Gene
150 Onthology (GO) annotation via InterProScan version 5 using the genes from the *S. aureus* TW20 genome
151 as a reference. Then, we conducted a GO enrichment analysis via Blast2GO PRO using a Fisher's Exact
152 Test (https://www.blast2go.com/), the analysis was conducted at three different GO levels: Molecular
153 Function (MF), Biological process (BP) and Cellular Component (CC). Then to evaluate the results, we
154 used a False Discovery Rate of 0.05 and the Benjamini-Hochberg correction was used to account for
155 multiple testing.
156

157 Phylogenetic reconstructions, average nucleotide identity analysis and neighbour nets.

158 We performed phylogenetic reconstructions for every SGF that did not show recombination signals as per
159 PhiTest (Bruen, Philippe et al. 2006) and for the super alignment (read below). All the gene trees were
160 constructed with RaxML version 8.2.11 (Stamatakis 2014) executing 10 inferences on the alignment,
161 using 10 distinct randomized Maximum Parsimony (MP) trees and with the GTR+G+I model. Both the
162 Shimodaira-Hasegawa topology test and the Robinson-Foulds distance were also conducted via RAxML
163 with the default settings. It is known that species tree estimation is not a trivial matter and we employed a
164 previous approach that gave trustable results (Castillo-Ramirez and Gonzalez 2008). We created a super
165 alignment concatenating all the SGFs that did not have signals for recombination and on this alignment a
166 ML phylogeny was constructed also through RAxML, this time executing 20 independent inferences
167 starting from 20 different MP trees and with GTR+G+I model. For this ML phylogeny bootstrap
168 replicates were generated again employing RaxML, using the -x option that implements a fast algorithm
169 for bootstrapping. The average nucleotide identity analysis (ANI) was run via the python module pyani
170 (http://widdowquinn.github.io/pyani/), specifically we employed the ANIm method (Richter and
171 Rosselló-Móra 2009). To determine the GC content and the proportion of variable sites for each SGF, we
172 used the function summary from the program AMAS (Borowiec 2016). In order to visualize the
173 conflicting phylogenetic signals in the genes *rpoB* and *tuf* we used SplitTree4 (Huson and Bryant 2005) to
174 construct Neighbour nets with uncorrected P distances.

175

176 RESULTS

177 Defining the species tree, confirming genome affiliations and the set of orthologous genes

178 We focused on the genus *Staphylococcus* as it has clear clinical and veterinary relevance and, due to that,
179 it has been extensively covered in terms of genome sequences. The data set employed for this study
180 represents 46 species, incorporating 40 type strains, and included a total of 269 genomes for the analyses
181 (see Supplementary Tables 1 and 2). Notably, this data set has species with a wide host-range covering
182 many of the niches described for this genus. First, we determined the single gene families (SGFs), as
183 these are potential candidates to be orthologous genes with one-to-one relationships, and found 208 SGFs.

184 However recombination could have affected them and, therefore, we conducted recombination tests on
185 them and around 15% of them (31 SGFs) showed signals of recombination and were discarded, which left
186 177 SGFs as good candidates to be orthologous groups (OG); the list of these 177 SGFs is provided in
187 Supplementary Table 3. We employed a previously used strategy (Castillo-Ramirez and Gonzalez 2008)

188 to approximate the species tree. This is the total evidence approach, in which all the 177 SGFs without
189 signals of recombination are concatenated and treated as if they were a single marker on which a
190 Maximum Likelihood (ML) phylogeny was constructed (see Figure 1) – this ML phylogeny was our
191 Proxy for the Species Tree Topology (PSTT). From Figure 1 one can see that type strains are scattered
192 throughout the tree and that most of the isolates from individual species tend to form monophyletic
193 groups and these groups are very well-supported as most of them have bootstrap values higher than 80
194 (see orange dots in the phylogeny). However, we also noted a region on the tree (shaded areas), where
195 different species intermingle together. This is caused by two species, namely *S. warneri* and *S*.

196 *saccharolyticus*, which do not form monophyletic groups and clearly some strains from these species
197 have been mislabelled. Additionally, as an independent strategy, we also carried out an ANI analysis to
198 calculate the relatedness of the strains to the type strains included in this study (see Figure 2). This
199 analysis shows that in terms of their taxonomy most of the genomes have been properly labelled. For
200 instance, all the strains designated as *S. lugdunensis* clustered with the type strain, *S. lugdunensis* NCTC
201 12217, with identity percentages well above 95% and thus could be assigned to this species (Chun, Oren
202 et al. 2018). The same applies to 78.26 % of the other species; see for example *S. gallinarum* or *S*.

203 *simulans* where genomes labelled as belonging to each of these species clustered tightly (again above the
204 95%) with the type strains. On the other hand, the mis-classification previously noted in the phylogeny is
205 also evident in this analysis, as the strains of both *S. warneri* and *S. saccharolyticus* do not cluster all
206 together in each case and their identity percentages are well below 95% (see green labels in Figure 2). We
207 then carried out a GO enrichment analysis to have an idea of the overrepresented functions of the 177
208 SGFs (see Supplementary Table 4), as expected many of the enriched biological processes and molecular
209 functions have to do with housekeeping functions; the top 2 GO molecular functions were ATP binding
210 and GTP binding, whereas the top 2 Biological processes were DNA repair and fatty acid biosynthesis
211 (see Supplementary Table 4 for more details). Taken together, these results demonstrate that the 177
212 SGFs seem to be real OGs and that most of the genomes have a proper affiliation regarding their
213 taxonomy.

214

215 Ranking the orthologous groups.

216 Then we analysed which of the 177 SGFs, from now on OGs, could be the best markers to infer the
217 evolutionary relationships among the species. For that end, we established how similar each one of the
218 single gene trees (from the 177 OGs) was compared to the PSTT; this was done employing the Robinson
219 and Foulds (RF) distance between the single gene trees and the PSTT – this distance gives the number of
220 bipartitions that are not shared by the two trees under consideration. We also used π, a common measure
221 of genetic diversity, to establish the amount of genetic variation for each OG. Figure 3 gives the
222 percentage of similarity of the 177 OG trees to the PSTT and the nucleotide diversity for each of the 177
223 OGs; this figure shows that no single OG yielded the same topology as the PSTT. Furthermore, we also
224 noted that every single OG has its own topology – not shared by any other OG - as none of all the
225 pairwise comparisons of the OGs trees gave a RF equal to 0 (a RF distance of 0 implies that the two
226 topologies in comparison are the same, see Supplementary Figure 2). Figure 3 also shows that 167 OGs
227 (94%) have nucleotide diversity values higher than 0.15, which make them good candidates for
228 phylogenetic markers. Of note, for two species (S. aureus and S. epidermidis) we also computed the intra-
229 species diversity and it seems that these OGs even at this level have a good amount of genetic diversity
230 (see Supplementary Figure 3). We found that the top 26 OGs were similar to the PSTT, as all of them
231 have RF distances below 214, indicating that they are more than 60 % identical to the PSTT. Notably
232 these 26 OGs present good values of nucleotide diversity (all but two showing values higher than 0.2
233 nucleotide changes per nucleotide site and the average for the 26 being 0.26 nucleotide changes per
nucleotide site), which is very convenient for their use as to phylogenetic markers. Table 1 provides details about these OGs such as RF distance PSTT, nucleotide diversity, function and the ID in the *S. aureus* TW20 genome. Remarkably, when we concatenated these 26 OGs and constructed a phylogeny, this tree showed a topology pretty similar to the PSTT (see Figure 4, panel A), the RF distance was 74 being 86.1% identical to the PSTT. Furthermore, using only the best 7 OGs (the bold candidates from Table 1), we got a similar result – that is the concatenated ML phylogeny of these 7 OGs is almost 80% identical to the PSTT (see panel B, Figure 4). Here it is worth mentioning that although these two phylogenies (Figure 4) are not 100% identical to the PSTT, the two phylogenies recovered almost all the species as monophyletic groups with very good bootstrap values (> 80%), see blue dots in both phylogenies - the only exceptions being the miss-classifications noted also in the PSTT and ANI analysis. In considering this part, although no single gene tree of the 177 OGs yielded the PSTT, just using the best 26 (and even just the top 7) OGs we were able to recover a topology pretty similar to the PSTT.

246

247 **The inadequacy of usual markers and the soundness of our strategy**

Then, we wanted to know how a set of commonly used markers for inferring the phylogenetic relationships within this genus compared to our list. We chose the genes *dnaJ*, *tuf*, *sodA*, *rpoB* and *hsp60* as these are the set of genes most used in previous studies (Kwok and Chow 2003, Ghebremedhin, Layer et al. 2008). Most of them had clear issues to be considered suitable markers given the criteria we employed to define our genuine OGs. Notably two of them, *tuf* and *rpoB*, did have signals for recombination as shown by their Neighbour nets and the PhiTest (see Supplementary Figure 1) and their nucleotide diversity was very low (below 0.15, see Figure 3 orange dots). On the other hand, *sodA* and *hsp60* were not single gene families, as they did have more than one gene per genome in some species. Among these usual markers, only *dnaJ* seems to be a good candidate in as much as it is a SGF, did not have signals of recombination and has a nucleotide diversity value above 0.2 (see Figure 3). The shortcoming of most of these genes is not totally unexpected, as these genes were not selected based on
phylogenetic criteria; however, it is clear that most of these genes do not seem to be a good option to infer the history of the species.

Finally, to prove the validity of our phylogenomic approach, in the data set here employed we included the genome sequences of four bacterial isolates collected by us from cattle ticks and phenotypically classified as members of the genus *Staphylococcus* (see Supplementary Table 2) - these isolates were sequenced using an Illumina MiSeq platform (see methods). Initially these isolates were classified using the 16S rRNA gene sequence (see Supplementary Table 2); however, in 2 of these cases, namely INIFAP 009-16 and INIFAP 002-15, the initial classification was incorrect. For instance, INIFAP 009-16 was classified as *S. succinus* by its 16S rRNA gene sequence but the PSTT placed this isolate with *S. xylosus*. In the case of INIFAP 002-15, whereas PSTT assigned it to *S. succinus*, the starting classification via the 16S rRNA gene was *S. xylosus*. Importantly, in the concatenate alignments using our 26 OGs (Table 1) the isolates were properly classified (see panel A, Figure 4). Furthermore, the ML phylogeny based on the concatenated alignment of the top 7 OGs also recovers the true affiliation of these isolates (see panel B, Figure 4). To sum up, our genuine OGs are a much better option than the usual markers to establish the phylogeny of the species and even just a few of these (the top 7) are able to adequately carry out taxonomy classification of the newly sequenced bacterial isolates.

## DISCUSSION

The main goal of our study was to establish the best markers for species identification in the genus *Staphylococcus*. To try to be as comprehensive as possible, we used a data set with 46 species and a total of 269 genomes from this genus that broadly represents most if not all the niches covered by this genus. Here we identified 177 OGs and ranked them according to their potential as phylogenetic markers; clearly, this set of markers should be useful not only for ecological and evolutionary studies but also for clinical and biotechnological purposes. In addition, our phylogenomic approach allowed us to identify two species (*S. warneri* and *S. saccharolyticus*) in which mislabelling of deposited genomes has occurred.
This issue of misclassification is not that rare, as genome sequences submitted in public databases many times do not undergo a proper inspection in terms of microbial taxonomy. However, the combination of these phylogenomic approaches (genome-based phylogeny and ANI analysis) can be very useful to pinpoint those cases of misclassification.

We want to emphasise that our study has two major contributions: one is the list of adequate markers for inferring the phylogeny of the species within the genus *Staphylococcus* and the other is the phylogenomic strategy that we used to identify such markers. Regarding the first major contribution, given that the average gene content of the species here analysed is 2413, our 177 OGs represent barely the 7 % of the genome of these species. However, we need to emphasise here that, for merely practical reasons, we focused on orthologous genes with one-to-one relationships and, very likely, there are many more orthologous genes within these species that did not pass our criteria given that they are not SGF. Despite the fact that all these OGs are potential good markers, in as much as their histories do not deviate significantly from the history of the species, they do differ among them in the amount of phylogenetic signal they each present. Importantly, no single OG was able to exactly represent the PSTT and no two OGs yielded the same topology; clearly, we found a cloud of topologies. This is in agreement with a previous study that found that only one of the several hundreds of OGs reflected the species tree (Castillo-Ramirez and Gonzalez 2008), although this focused on different evolutionary scales and different bacteria. It is worth mentioning, that the top 26 OGs have good values of nucleotide diversity (all but two have values higher than 0.21), and are localized in different parts of the chromosome and present different functions. We think these 26 OGs represent a really good set of markers for inferring the PSTT; actually concatenating the best 7 it was possible to get a very good approximation of the PSTT. Here, we want to highlight that using the phylogeny of the top 7 OGs all but the 2 species with issues of mis-classification were recovered as monophyletic groups and thus it seems that just using these 7 markers species assignation can be carried out. Clearly, this set of markers could be very useful for underfunded laboratories working with *Staphylococcus* isolates. Although from a genotyping perspective, ideally one
311 would want to use whole-genome sequencing for species identification, many laboratories in the world
312 (specially in developing countries) still cannot afford the sequencing of tens to hundreds of isolates. Thus,
313 this short list of markers should be extremely useful for those who only can sequence a few loci, which is
314 often the case for clinical, environmental and evolutionary biology microbiologists in developing
315 countries. Clearly, the markers highlighted in this study are a much better option that the commonly used
316 markers, as the latter seem to exhibit obvious flaws (more than one copy per genome, signals of
317 recombination, low nucleotide diversity values) for inferring the history of the species. Our strategy, and
318 in turn the list of genes found, proved to be factually sound, even when using just a few markers for the
319 taxonomic assignment of newly sequenced bacterial isolates.

320

321 Maybe more important than the first major contribution is the fact that our phylogenomic approach will
322 allow the identification of adequate markers in many different genera – not only from bacteria but also
323 from archaea. This is very important, as gene families showing only orthologous relationships very likely
324 do not extend all over the tree of life – or just for a few gene families. Along these lines, it is very likely
325 that many of the markers found here will not work for other genera (i.e. they will not show one-to-one
326 orthologous relationships) as the homologous genes within those genera might have been affected by
327 horizontal gene transfer or duplication and differential loss, or some other molecular event that prevent
328 the history of the gene to reflect the speciation events. We want to highlight that the biology of the
329 species under consideration could have a very important effect on the number of orthologous genes found.
330 For instance, highly recombinogenic species, such as *Neisseria gonorrhoeae* (Ezewudo, Joseph et al.
331 2015) and *Acinetobacter baumannii* (Grana-Miraglia, Lozano et al. 2017), would have fewer orthologous
332 genes than more clonal species. Nonetheless, the strategy employed by us should work even in these
333 species, although the number of potential orthologous genes should be much less.

334

335 CONCLUSIONS
336 In summary, here we devised and applied a phylogenomic approach that allowed us to define the most
337 suitable markers for inferring the species phylogeny within the genus *Staphylococcus*. We acknowledge
338 that the effectiveness of these markers for other bacterial genera remains open to investigation. However,
339 we are confident that the phylogenomic approach here implemented can be employed to identify the most
340 suitable markers for other genera. On a broader level, this study has very practical implications, as it
341 provides a general framework for mapping out the best markers to construct the phylogeny of the taxa
342 considered not only at the genus level but also at other taxonomic levels.

343

344

345

346

347 Acknowledgments

348 We are extremely grateful to the editor and the reviewers, as their suggestions have drastically improved
349 our manuscript.

350

351 REFERENCES

352 Bankevich, A., S. Nurk, D. Antipov, A. A. Gurevich, M. Dvorkin, A. S. Kulikov, V. M. Lesin, S. I.
353 Nikolenko, S. Pham and A. D. Prjibelski (2012). "SPAdes: a new genome assembly algorithm
354 and its applications to single-cell sequencing." *Journal of computational biology* **19**(5): 455-477.

356 Becker, N. S., G. Margos, H. Blum, S. Krebs, A. Graf, R. S. Lane, S. Castillo-Ramírez, A. Sing and
357 V. Fingerle (2016). "Recurrent evolution of host and vector association in bacteria of the
358 Borrelia burgdorferi sensu lato species complex." *BMC genomics* **17**(1): 734.

359 Borowiec, M. L. (2016). "AMAS: a fast tool for alignment manipulation and computing of
360 summary statistics." *PeerJ* **4**: e1660.

361 Bradley, R. K., A. Roberts, M. Smoot, S. Juvekar, J. Do, C. Dewey, I. Holmes and L. Pachter
362 (2009). "Fast statistical alignment." *PLoS computational biology* **5**(5): e1000392.

363 Bruen, T. C., H. Philippe and D. Bryant (2006). "A simple and robust statistical test for
detecting the presence of recombination." *Genetics* **172**(4): 2665-2681.

365 Castillo-Ramirez, S. and V. Gonzalez (2008). "Factors affecting the concordance between
366 orthologous gene trees and species tree in bacteria." *BMC Evol Biol* **8**: 300.
367 Challagundla, L., J. Reyes, I. Rafiqullah, D. O. Sordelli, G. Echaniz-Aviles, M. E. Velazquez-

368 Meza, S. Castillo-Ramírez, N. Fittipaldi, M. Feldgarden and S. Chapman (2018).

369 "Phylogenomic Classification and the Evolution of Clonal Complex 5 Methicillin-Resistant

370 Staphylococcus aureus in the Western Hemisphere." *Frontiers in Microbiology* **9**: 1901.

371 Chun, J., A. Oren, A. Ventosa, H. Christensen, D. R. Arahal, M. S. da Costa, A. P. Rooney, H. Yi, X. W. Xu, S. De Meyer and M. E. Trujillo (2018). "Proposed minimal standards for the use of genome data for the taxonomy of prokaryotes." *Int J Syst Evol Microbiol* **68**(1): 461-466.

374 Cox, M. P., D. A. Peterson and P. J. Biggs (2010). "SolexaQA: At-a-glance quality assessment of Illumina second-generation sequencing data." *BMC bioinformatics* **11**(1): 485.

376 Ezewudo, M. N., S. J. Joseph, S. Castillo-Ramirez, D. Dean, C. Del Rio, X. Didelot, J.-A. Dillon, R.

377 F. Selden, W. M. Shafer and R. S. Turingan (2015). "Population structure of Neisseria

378 gonorrhoeae based on whole genome data and its relationship with antibiotic resistance."

379 *PeerJ* **3**: e806.

380 Fitch, W. M. (2000). "Homology a personal view on some of the problems." *Trends Genet*

381 **16**(5): 227-231.

382 Frisch, M., S. Castillo-Ramirez, R. Petit, M. Farley, S. Ray, V. Albrecht, B. Limbago, J.

383 Hernandez, I. See and S. Satola (2018). "Invasive methicillin-resistant Staphylococcus

384 aureus USA500 strains from the US Emerging Infections Program constitute three

385 geographically distinct lineages." *mSphere* **3**(3): e00571-00517.

386 Gabaldón, T. and E. V. Koonin (2013). "Functional and evolutionary implications of gene

387 orthology." *Nature Reviews Genetics* **14**(5): 360.

388 Ghebremedhin, B., F. Layer, W. König and B. König (2008). "Genetic classification and

389 distinguishing of Staphylococcus species based on different partial gap, 16S rRNA, hsp60,

390 rpoB, sodA, and tuf gene sequences." *Journal of clinical microbiology* **46**(3): 1019-1025.

391 Grana-Miraglia, L., L. F. Lozano, C. Velazquez, P. Volkow-Fernandez, A. Perez-Oseguera, M.

392 A. Cevallos and S. Castillo-Ramirez (2017). "Rapid Gene Turnover as a Significant Source of

393 Genetic Variation in a Recently Seeded Population of a Healthcare-Associated Pathogen."

394 *Front Microbiol* **8**: 1817.

395 Huson, D. H. and D. Bryant (2005). "Application of phylogenetic networks in evolutionary

396 studies." *Molecular biology and evolution* **23**(2): 254-267.

397 Joseph, S. J., D. Cox, B. Wolff, S. S. Morrison, N. A. Kozak-Muiznieks, M. Frace, X. Didelot, S.

398 Castillo-Ramirez, J. Winchell and T. D. Read (2016). "Dynamics of genome change among

399 Legionella species." *Scientific reports* **6**: 33442.

400 Joseph, S. J., H. Marti, X. Didelot, S. Castillo-Ramirez, T. D. Read and D. Dean (2015).

401 "Chlamydiaceae genomics reveals interspecies admixture and the recent evolution of Chlamydia abortus infecting lower mammalian species and humans." *Genome biology and evolution* **7**(11): 3070-3084.

404 Kwok, A. Y. and A. W. Chow (2003). "Phylogenetic study of Staphylococcus and

405 Macrococcus species based on partial hsp60 gene sequences." *International journal of*

406 *systematic and evolutionary microbiology* **53**(1): 87-92.

407 Lamers, R. P., G. Muthukrishnan, T. A. Castoe, S. Tafur, A. M. Cole and C. L. Parkinson (2012).

408 "Phylogenetic relationships among Staphylococcus species and refinement of cluster

409 groups based on multilocus data." *BMC evolutionary biology* **12**(1): 171.

410 López-Leal, G., M. L. Tabche, S. Castillo-Ramírez, A. Mendoza-Vargas, M. A. Ramírez-Romero

411 and G. Dávila (2014). "RNA-Seq analysis of the multipartite genome of Rhizobium etli CE3
412 shows different replicon contributions under heat and saline shock." BMC genomics 15(1):

413 770.

414 Parks, D. H., M. Imelfort, C. T. Skennerton, P. Hugenholtz and G. W. Tyson (2015). "CheckM:

415 assessing the quality of microbial genomes recovered from isolates, single cells, and

416 metagenomes." Genome research: gr. 186072.186114.

417 Richter, M. and R. Rosselló-Móra (2009). "Shifting the genomic gold standard for the

418 prokaryotic species definition." Proceedings of the National Academy of Sciences 106(45):

419 19126-19131.

420 Sasaki, T., S. Tsubakishita, Y. Tanaka, A. Sakusabe, M. Ohtsuka, S. Hirotaki, T. Kawakami, T.

421 Fukata and K. Hiramatsu (2010). "Multiplex-PCR method for species identification of

422 coagulase-positive staphylococci." Journal of clinical microbiology 48(3): 765-769.

423 Stamatakis, A. (2014). "RAxML version 8: a tool for phylogenetic analysis and post-analysis

424 of large phylogenies." Bioinformatics 30(9): 1312-1313.

425 Zerbino, D. and E. Birney (2008). "Velvet: algorithms for de novo short read assembly using

426 de Bruijn graphs." Genome research: gr. 074492.074107.

427
**Figure 1**(on next page)

Maximum likelihood phylogenetic tree based on the super alignment of the orthologous groups.

Maximum likelihood phylogenetic tree based on the super alignment of the 177 non-recombinant SGFs. Strains of the same species are labelled with the same colour and black stars indicate type strain for some species. Green and violet rectangles highlight strains with conflicting clade assignment, *S. warneri* and *S. saccharolyticus* respectively. Bootstrap values above 80 are shown in orange dots. The scale bar shows the number substitutions per site. The code for the strain identifiers is in Supplementary Table 1.
Tree scale: 0.1

★ Type Strains

<figure><img src="image_1.png" /></figure>
**Figure 2 (on next page)**

Heat map of the Average Nucleotide Identity (ANI) analysis of the 269 strains from
*Staphylococcus spp.*

The reddish cells show identity percentages above 95% implying that the strains belong to the same species, whereas non-reddish colors denote identity percentages below 95% (see ANI key). The top and left side rows on the heat map show the species assignation (see species key). The dendograms across the side and top show the clustering of the strains; strain identifiers as in Figure 1. Color of the strain identifiers (bottom row) is as follows: red gives the type strains, whereas blue shows the strains with issues of mis-classification. This analysis indicates 44 clearly discernible *Staphylococcus spp.*
<figure><img src="image_2.png" /></figure>

**Species**
* magnetis
* argenteus
* arlettae
* aureus
* auricularis
* capitis
* caprae
* carnosus
* chromogenes
* cohnii
* condimenti
* delphini
* devriesei
* epidermidis
* equorum
* felis
* fleurettii
* gallinarum
* haemolyticus
* hominis
* hyicus
* intermedius
* kloosii
* lentus
* lugdunensis
* lutrae
* massilliensis
* microti
* muscae
* nepalensis
* pasteurii
* petrasii
* pettenkoferi
* pseudintermedius
* saccharolyticus
* saprophyticus
* schleiferi
* schweitzeri
* sciuri
* simiae
* simulans
* stepanovici
* succinus
* vitulinus
* warneri
* xylosus
* Unknown

**ANI**
* 1
* 0.98
* 0.95
* 0.92
* 0.89
* 0.86
* 0.83
* 0.8
Figure 3

Similarity to the species tree and nucleotide diversity

Percentage of similarity between gene trees and the PSTT and the nucleotide diversity for each of the 177 SGF (green dots). The most similar gene tree topologies correspond to genes with high levels of nucleotide diversity. The commonly used marker genes *dnaJ*, *rpoB* and *tuf* are highlighted (orange dots), they all show low to moderate nucleotide diversity values and none of them show a similarity percentage to the PSTT above 60 %.

<figure><img src="image_3.png" /></figure>

Nucleotide Diversity
Figure 4

Phylogenies of the top orthologous groups

Phylogenetic trees based on the concatenated alignments of the 26 orthologous groups listed in Table 1 (A) and on the concatenated alignments of the top seven orthologous groups (B). The tree in panel A has a percentage of similarity with the PSTT of 86% and the tree in panel B of 79%. Bootstrap values above 80 are shown in blue dots. Scale bar shows the number of substitutions per site.

<figure><img src="image_4.png" /></figure>
**Table 1**(on next page)

Top 26 genuine orthologous genes

These are the most adequate markers for inferring the species phylogeny. * Gene names were taken from the UniProt (www.uniprot.org) database when possible, otherwise the proteins were blasted against NCBI protein database (www.ncbi.nlm.nih.gov) and a well-annotated *Staphylococcus* species was used for annotation.
<table><thead><tr><th>Description*</th><th>NCBI Reference Sequence</th><th>GC content</th><th>Nucleotide Diversity</th><th>Proportion of Variable sites</th><th>Similarity with Species Tree (Robinson-Foulds Distance)</th></tr></thead><tbody><tr><td>DNA mismatch repair protein (MutS)</td><td>WP_000073352.1</td><td>0.342</td><td>0.2370335</td><td>0.412</td><td>65.4135 (184)</td></tr><tr><td>ATP-dependent RecD-like DNA helicase (recD2)</td><td>WP_001283311.1</td><td>0.345</td><td>0.2638535</td><td>0.588</td><td>65.4135 (184)</td></tr><tr><td>Cytosol aminopeptidase (pepA)</td><td>WP_001009697.1</td><td>0.361</td><td>0.3404309</td><td>0.760</td><td>64.6617 (188)</td></tr><tr><td>putative ABC transporter ATP-binding protein (YheS)</td><td>WP_000602071.1</td><td>0.339</td><td>0.2515735</td><td>0.688</td><td>64.6617 (188)</td></tr><tr><td>Phenylalanine-tRNA ligase beta subunit (pheT)</td><td>WP_000908982.1</td><td>0.366</td><td>0.2683997</td><td>0.696</td><td>63.9098 (192)</td></tr><tr><td>Penicillin-binding protein H (pbp3)</td><td>WP_000919772</td><td>0.340</td><td>0.2553319</td><td>0.623</td><td>63.9098 (192)</td></tr><tr><td>DNA-directed RNA polymerase subunit beta' (rpoC)</td><td>CBI48492.1</td><td>0.381</td><td>0.1488365</td><td>0.442</td><td>63.5338 (194)</td></tr><tr><td>Pyruvate kinase (pyk)</td><td>WP_001232648.1</td><td>0.361</td><td>0.2139191</td><td>0.584</td><td>63.5338 (194)</td></tr><tr><td>DNA polymerase III subunit alpha (DnaE)</td><td>WP_000226911.1</td><td>0.336</td><td>0.2895834</td><td>0.737</td><td>63.1579 (196)</td></tr><tr><td>UDP-N-acetylmuramate-L-alanine ligase (murC)</td><td>WP_000150163.1</td><td>0.335</td><td>0.2352931</td><td>0.622</td><td>63.1579 (196)</td></tr><tr><td>fibronectin-binding domain-containing protein (YloA)</td><td>WP_000312763.1</td><td>0.338</td><td>0.2644582</td><td>0.667</td><td>62.782 (198)</td></tr></tbody></table>
<table>
  <tr>
    <td>UDP-N-acetylmuramoyl-L-alanyl-D-glutamate-L-lysine ligase (murE)</td>
    <td>WP_000340119.1</td>
    <td>0.374</td>
    <td>0.2366963</td>
    <td>0.638</td>
    <td>62.0301 (202)</td>
  </tr>
  <tr>
    <td>Homoserine dehydrogenase (dhoM)</td>
    <td>WP_000735864.1</td>
    <td>0.344</td>
    <td>0.2720552</td>
    <td>0.579</td>
    <td>62.0301 (202)</td>
  </tr>
  <tr>
    <td>Acetyl-coenzyme A carboxylase carboxyl transferase subunit alpha (accA)</td>
    <td>WP_000883645.1</td>
    <td>0.357</td>
    <td>0.2383029</td>
    <td>0.584</td>
    <td>62.0301 (202)</td>
  </tr>
  <tr>
    <td>Molybdopterin molybdenumtransferase (MoeA)</td>
    <td>WP_000259718.1</td>
    <td>0.386</td>
    <td>0.2844009</td>
    <td>0.679</td>
    <td>62.0301 (202)</td>
  </tr>
  <tr>
    <td>S-adenosylmethionine tRNA ribosyltransferase-isomerase (queA)</td>
    <td>WP_001019171.1</td>
    <td>0.357</td>
    <td>0.2468871</td>
    <td>0.620</td>
    <td>61.6541 (204)</td>
  </tr>
  <tr>
    <td>ATP-dependent DNA helicase (RecQ)</td>
    <td>WP_000983677.1</td>
    <td>0.345</td>
    <td>0.2654353</td>
    <td>0.659</td>
    <td>61.2782 (206)</td>
  </tr>
  <tr>
    <td>hypothetical protein (YibE/F-like)</td>
    <td>WP_001794550.1</td>
    <td>0.363</td>
    <td>0.3498756</td>
    <td>0.723</td>
    <td>60.9023 (208)</td>
  </tr>
  <tr>
    <td>DNA polymerase III subunit tau (dnaX)</td>
    <td>WP_001109047.1</td>
    <td>0.368</td>
    <td>0.2727606</td>
    <td>0.645</td>
    <td>60.9023 (208)</td>
  </tr>
  <tr>
    <td>Ribonuclease R (rnr)</td>
    <td>WP_001050064.1</td>
    <td>0.367</td>
    <td>0.2284088</td>
    <td>0.626</td>
    <td>60.5263 (210)</td>
  </tr>
  <tr>
    <td>Ktr system potassium uptake protein B (ktrB)</td>
    <td>WP_000021864.1</td>
    <td>0.342</td>
    <td>0.2674265</td>
    <td>0.646</td>
    <td>60.5263 (210)</td>
  </tr>
  <tr>
    <td>Dihydrolipoyl dehydrogenase (pdhD)</td>
    <td>WP_001291535.1</td>
    <td>0.353</td>
    <td>0.3373007</td>
    <td>0.759</td>
    <td>60.5263 (210)</td>
  </tr>
  <tr>
    <td>3-oxoacyl-[acyl-3-oxoacyl]-</td>
    <td>WP_001100524.1</td>
    <td>0.377</td>
    <td>0.2367141</td>
    <td>0.601</td>
    <td>60.1504 (212)</td>
  </tr>
</table>
<table><tr><td>carrier-protein]<br/>synthase 3 (fabH)</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Signal recognition<br/>particle receptor<br/>(FtsY)</td><td>WP_000007682.1</td><td>0.366</td><td>0.1979450</td><td>0.506</td><td>60.1504 (212)</td></tr><tr><td>6-<br/>phosphogluconolac<br/>tonase<br/>(hypothetical)</td><td>WP_000181322.1</td><td>0.358</td><td>0.3401093</td><td>0.737</td><td>60.1504 (212)</td></tr><tr><td>Aspartate<br/>aminotransferase<br/>(AAT-like)</td><td>WP_001068542.1</td><td>0.342</td><td>0.2853785</td><td>0.678</td><td>60.1504 (212)</td></tr></table>