Comparative Genomics of Staphylococcus Reveals Determinants of Speciation and Diversification of Antimicrobial Defense.

The bacterial genus Staphylococcus comprises diverse species with most being described as colonizers of human and animal skin. A relational analysis of features that discriminate its species and contribute to niche adaptation and survival remains to be fully described. In this study, an interspecies, whole-genome comparative analysis of 21 Staphylococcus species was performed based on their orthologues. Three well-defined multi-species groups were identified: group A (including aureus/epidermidis); group B (including saprophyticus/xylosus) and group C (including pseudintermedius/delphini). The machine learning algorithm Random Forest was applied to prioritize orthologs that drive formation of the Staphylococcus species groups A-C. Orthologues driving staphylococcal intrageneric diversity comprised regulatory, metabolic and antimicrobial resistance proteins. Notably, the BraSR (NsaRS) two-component system (TCS) and its associated BraDE transporters that regulate antimicrobial resistance showed limited distribution in the genus and their presence was most closely associated with a subset of Staphylococcus species dominated by those that colonize human skin. Divergence of BraSR and GraSR antimicrobial peptide survival TCS and their associated transporters was observed across the staphylococci, likely reflecting niche specific evolution of these TCS/transporters and their specificities for AMPs. Experimental evolution, with selection for resistance to the lantibiotic nisin, revealed multiple routes to resistance and differences in the selection outcomes of the BraSR-positive species S. hominis and S. aureus. Selection supported a role for GraSR in nisin survival responses of the BraSR-negative species S. saprophyticus. Our study reveals diversification of antimicrobial-sensing TCS across the staphylococci and hints at differential relationships between GraSR and BraSR in those species positive for both TCS.

Background

Staphylococcus Species and Genomics

The existence of taxonomically distinct species groups was first proposed for Staphylococcus based on differential DNA–DNA hybridization methods (Kloos et al., 1991). These groups were supported by 16S rDNA sequence analysis of 38 taxa (Takahashi et al., 1999) and multilocus sequence data of around 60 species and subspecies (Lamers et al., 2012).

A comparative analysis that utilized next generation genome sequencing data of staphylococci to probe phylogenetic relationships with 491 shared orthologues across 12 Staphylococcus species (Suzuki et al., 2012) proposed S. pseudintermedius and S. carnosus as the most basal lineages. Moreover, with 10 species in their analysis being residents of human skin, the authors proposed that evolution selected for human adaptation after branching from S. carnosus. The relationships between the strains generated from shared orthologues were maintained using total gene content (Suzuki et al., 2012). However, in contrast to the conclusions of 16S rDNA and multilocus data (Takahashi et al., 1999; Lamers et al., 2012) their analysis revealed discrete clustering of Staphylococcus species. In contrast with this analysis, no distinct clustering of S. hominis with S. haemolyticus was observed, and S. saprophyticus was assigned to the S. epidermidis group of species (Suzuki et al., 2012). Currently, there is a knowledge gap in Staphylococcus species comparisons with a need to determine if this clustering of staphylococcal species is supported using whole genome data. Our findings here begin to close this gap.

Two Component Systems

Prokaryotes are receptive to environmental stimuli through diverse sensory and transducing two component systems (TCS). These TCS archetypically comprise a sensor histidine kinase (HK) that spans the cell membrane to interact with the external environment. Stimulus perception causes conditional autophosphorylation that is relayed to an interacting response regulator (RR) to enable DNA-binding directed transcription modulation (Stock et al., 2000).

While TCS are widespread and diverse across prokaryotes, the intramembrane-sensing histidine kinases (IM-HK) are specific to the Firmicutes. This family of small HKs has a short, 25 amino acid linker region between each 400 amino acid transmembrane helix. S. aureus GraSR uses a IM-HK to regulate a global network responsible for resistance to antimicrobial peptides (AMPs). GraSR modulates the expression of DltABCD and MprF that in concert alter the S. aureus surface charge to evade electrostatic interaction-mediated targeting of cationic AMPs (Falord et al., 2011).

An orthologous TCS to GraSR described in S. aureus was concurrently designated BraSR and NsaSR by two different groups (Blake et al., 2011; Hiron et al., 2011). Serial passage in sub-MIC concentrations of the lantibiotic nisin was shown to select increased nisin MIC due to a SNP in nsaS gene encoding sensor histidine kinase of NsaRS (nisin susceptibility-associated sensor regulator) (Blake et al., 2011). The TCS was separately designated BraSR (bacitracin resistance-associated sensor regulator) from the reduced MIC of bacitracin and nisin determined for the TCS gene mutant (Hiron et al., 2011). BraR binding sites were revealed upstream of the ABC transporter genes braDE and vraDEH (Popella et al., 2016) that were not transcribed in the mutant but induced in the presence of bacitracin. The transporter BraDE contributes to the detection of nisin and bacitracin and subsequent signal transduction via BraSR, whereas VraDE is more directly involved in detoxification by efflux (Hiron et al., 2011). Transcription of braSR is increased following exposure to multiple antibiotics, including ampicillin, phosphomycin and nisin. Inactivation of braS (nsaS) revealed differential transcription of 245 genes (Kolar et al., 2011), revealing the TCS might report cell envelope stress to directly regulate biofilm formation, cellular transport and responses to anoxia.

In this study, a comparative genome analysis of 21 Staphylococcus species was performed based upon their orthologous gene content. Species groups were revealed and then interrogated using the Random Forest algorithm to identify group-contributing genes. The operon encoding the BraSR TCS was found to differentiate the S. aureus/S. epidermidis species group from other species groups determined in the study and the TCS was found to have restricted distribution across 49 species of Staphylococcus. Experimental evolution of representative braSR-positive and -negative species with nisin selection identified differential selection of BraSR and GraSR to produce resistance to this AMP.

Results and Discussion

Analysis of Orthologous Gene Content Across the Staphylococci

The orthologous gene content of 21 sequenced staphylococcal species’ genomes (Table 1) was determined using OrthoMCL to group orthologous genes (homologues separated by speciation) into clusters across the different species. The number of shared orthologous clusters between the different species’ genomes was then represented as a heatmap (Figure 1). The output from this analysis revealed the assembly of three major groups of species, each with high numbers of shared orthologous clusters. An associated cladogram supported three groups (groups A, B, and C) when defined as containing three or more species (Figure 1). This supported previous reported groupings from 16S rDNA and multilocus analyses (Takahashi et al., 1999; Lamers et al., 2012). Additionally, three species pairs showed a high degree of shared orthologous clusters of genes and branched together in the cladogram: S. aureus/S. simiae, S. simulans/S. carnosus, and S. lentus/S. vitulinus. S. aureus and S. simiae were proposed as members of the S. aureus group of staphylococci from gene content (Takahashi et al., 1999).

Table 1

Staphylococcus species and strains included in OrthoMCL analysis.

Staphylococcus	Strain	Genome	Sequence
species		accession	status
			(Reference)
S. arlettae	CVD059	ALWK01000000	Draft
		(Uid175126)	(Dinakaran et al., 2012)
S. aureus	Newman	AP009351	Complete
		(Uid58839)	(Baba et al., 2008)
S. capitis	SK14	ACFR01000000	Draft
		(Uid55415)
S. carnosus	TM300	NC012121	Complete
		(Uid59401)	(Rosenstein et al., 2009)
S. cohnii	ATCC29974	LT963440	Draft
			(This study)
S. delphini	8086	CAIA00000000	Draft
		(Uid199664)
S. epidermidis	ATCC_12228	NC005008	Complete
		(Uid57861)	(Zhang et al., 2003)
S. equorum	Mu2	CAJL01000000	Draft
		(Uid169178)
S. haemolyticus	K8	LT963441	Draft
			(This study)
S. hominis	J6	LT963442	Draft
			(This study)
S. intermedius	NCTC_11048	CAIB01000000	Draft
		(Uid199665)
S. lentus	F1142	AJXO01000000	Draft
		(Uid200144)	(Nam et al., 2012b)
S. pettenkoferi	VCU012	AGUA00000000	Draft
		(Uid180074)
S. lugdunensis	HKU09	CP001837	Complete
		(Uid46233)	(Tse et al., 2010)
S. pseudintermedius	HKU10	Uid62125	Complete
			(Tse et al., 2011)
S. saprophyticus	ATCC_15305	AP008934	Complete
		(Uid58411)	(Kuroda et al., 2005)
S. simiae	CCM_7213	AEUN0000000	Draft
		(Uid77893)	(Suzuki et al., 2012)
S. simulans	ATCC 27848	LT963435	Draft
			(This study)
S. vitulinus	F1028	AJTR0000000	Draft
		(Uid200114)	(Nam et al., 2012a)
S. warneri	SG1	CP003668	Complete
		(Uid187059)	(Cheng et al., 2013)
S. xylosus	ATCC29971	LT963439	Draft
			(This study)

FIGURE 1

Heat map representation of shared orthologous proteins across Staphylococcus species. Presence is indicated using a color scale from red (highest number of shared clusters of orthologous proteins) to white (lowest number). Major groups of species observed in the analysis are highlighted as groups A–C.

The largest and least well-defined, species group comprises S. epidermidis, S. capitis, S. warneri, S. haemolyticus, S. hominis, S. lugdunensis, S. pettenkoferi, S. aureus, and S. simiae (Figure 1). Designated group A, it is dominated by species that colonize human skin (Kloos, 1980; Schleifer and Bell, 2015). The likelihood of a strain-dependent effect structuring group A was investigated by substituting S. epidermidis, S. hominis, and S. aureus strains based on multiple available genomes (Table 1 and Supplementary File S1). Substituting these individual species with alternative strains and repeating the OrthoMCL analysis did not alter species groupings. Groups B and C were similarly unaffected by switching strains of S. saprophyticus and S. pseudintermedius, respectively.

The smaller species group B comprises S. equorum, S. arlettae, S. cohnii, S. saprophyticus, and S. xylosus (Figure 1). Though not universal, a frequent lifestyle identified in the group B species is human or animal host colonization; several species are associated with meat products and novobiocin resistance (Devriese et al., 1985; Bannerman et al., 1991) with commonalities in their cell wall composition (Schleifer et al., 1984).

Species group C comprises S. pseudintermedius, S. delphini, and S. intermedius and this collective was previously designated the S. intermedius group (SIG); the species cause opportunistic infection of companion animals and equids (Devriese et al., 1985). Emerging antibiotic resistance in the SIG species group is a clinical veterinary concern (Stull et al., 2014) and their routine speciation is complicated by their high degree of 16S rRNA locus sequence identity (Slettemeås et al., 2010).

While preserving known species groupings, the whole genome analysis identified discrete species groups of staphylococci (A-C) and an explanation for their formation was sought. Genetic determinants directing the formation of species group A were tested in R using machine learning with the Random Forests algorithm for classification (Breiman, 2001). This algorithm was used to identify variables, in this case OrthoMCL clusters (data not shown), that contributed to formation of the groups, based on a forest of trees generated from these variables. A gene from each cluster was then determined and mapped back to a representative genome and the PROKKA annotation of each protein coding sequence was verified using BLAST, for group A this representative was S. epidermidis. Contributing variables were investigated for group A, based on the strain set described in Table 1, where permutations were used to verify the existence and reproducibility of species groups (Supplementary File S1).

Clusters Driving Formation of Group a Species

The presence of 7 and absence of 6 OrthoMCL clusters collectively contribute to defining group A, with differing levels of support (Mean Decrease in Accuracy [MDA] values) (Table 2 and Supplementary File S2). Four orthologs that are sequentially encoded in the genome as an operon (epi_02134 - epi_02137; MDA 3.2, 3.0, 2.6, 2.2, respectively) were also the most strongly supported in this analysis (Table 2). The latter cluster pair epi_02136/epi_02137 was annotated by PROKKA as a TCS sensor/regulator (Table 2 and Supplementary File S2) and shares ∼100% similarity with BraSR (SA2417/SA2418 of S. aureus N315), a TCS associated with resistance to AMPs nisin and bacitracin (Hiron et al., 2011). The adjacent clusters encoded in the same operon (epi_02134, epi_02135) comprise the BraD/BraE ABC transporter subunits with 98 and 99% similarity with SA2415/SA2416 of S. aureus N315, respectively (Hiron et al., 2011). We demonstrate as a key finding of our analysis that BraSR and BraDE are associated with genomes of group A Staphylococcus species.

Table 2

Proteins driving formation of species group A.

Group A staphylococci Random Forest Output
PROKKA annotation	BLAST homology	MDA
Presence
Putative cell wall associated hydrolase (epi_00542)	Hypothetical protein	2.2
Hypothetical protein (epi_02098)	Cell wall surface anchor protein	2.7
Hypothetical protein (epi_02108)	Hypothetical protein	1.9
FtsX like permease family protein (epi_02134)	ABC transporter permease	3.2
Macrolide export ATP binding/permease protein MacB (epi_02135)	Bacteriocin ABC transporter ATP-binding protein	3.0
Sensor histidine kinase GraS (epi_02136)	TCS histidine kinase	2.6
Glycopeptide resistance associated protein R (epi_02137)	TCS transcriptional regulator	2.2
Absence
Succinate semialdehyde dehydrogenase NADP+ (sap_00201)	Succinate-semialdehyde dehydrogenase	2.9
Putative membrane protein putative toxin regulator (sap_00203)	PTS sugar transporter subunit IIC	2.7
Putative multidrug resistance ABC transporter ATP binding/permease protein YheI (sap_00398)	Multidrug ABC transporter ATP-binding protein	3.3
Putative multidrug resistance ABC transporter ATP binding/permease protein YheH (sap_00399)	Multidrug ABC transporter ATP-binding protein	1.4
L-lactate utilization operon repressor (sap_00760)	Transcriptional regulator	2.7
Glutamate aspartate carrier protein (sap_01003)	Sodium:dicarboxylate symporter	2.7

The presence of orthologue epi_00542 (MDA 2.2; Table 2 and Supplementary File S2) contributes to species group A, with support that the protein functions as a putative cell wall hydrolase from the Nlp-P60 family hydrolase domain that is associated with hydrolysis of peptidoglycan. Also, contributing to defining group A are the absences of two orthologue clusters (sap_00398; MDA 3.3 and sap_00399; MDA 1.4; Table 2 and Supplementary File S2) that are annotated as multidrug ABC transporters. A range of cytotoxic molecules are mobilized across the cell membrane by multidrug ABC transporters where certain families of these can also act as sensors (Stock et al., 2000; Rietkötter et al., 2008). Across staphylococcal groups, differential repertoires of ABC transporters associated with antimicrobial survival are consistent with the importance of community competition in species evolution.

Sequence variation of the NADP-dependent succinate semialdehyde dehydrogenase (SSADH) between group A staphylococci versus groups B and C was identified by the association of cluster sap_00201 (MDA 2.9, Table 2 and Supplementary File S2) with group A species; this variation might be allied to differences in glutamate metabolism across the genus. Glutamate is involved in multiple metabolic processes and bacterial glutamate dehydrogenase catabolizes glutamate, which contributes to acid tolerance. NADP-SSADH catalyzes catabolism of γ-aminobutyrate, a product of glutamate dehydrogenase activity (Feehily and Karatzas, 2013); this pathway is oxidative stress sensitive owing to the catalytic cysteine residue of SSADH.

With respect to clusters driving formation of group B and C species, the size of species input groups B and C (Figure 1) limit use of the random forest algorithm. Consequently, a similar species-defined analysis of groups B and C is not included and a broader species comparison of staphylococci could be considered in future.

Diversity of Cationic AMP Survival Loci Across the Staphylococci

The described comparative genomic analysis revealed that while BraSR TCS is associated with group A species of staphylococci, the GraSR TCS is distributed across all species groups. Supporting predictions from the Random Forest analysis, low sequence identity of BraR/BraS with GraR/GraS was confirmed. BraR mean sequence identity with GraR of group A (44%) and group B/C species (40%) was greater than that of BraS compared with GraS of group A and groups B/C (mean ∼30% and ∼26%, respectively) (Table 3).

Table 3

Comparative sequence identity of the BraRS and GraRS TCS across species groups A–C.

TCS protein	Mean identity within group A	Mean identity to groups B and C GraR or GraS	Mean identity to group A GraR or GraS
BraR	77.1	39.6	44.3
BraS	62.9	26.4	29.8
TCS protein	Mean identity within group A	Mean identity within group B	Mean identity within group C	Mean identity across groups
GraR	87.8	84	97.9	66.7
GraS	69.4	66	88.2	48.2

High mean sequence identity (84–98%) of GraR regulator protein occurs within each of the three species groups (Table 3) with divergence of GraR between species groups identified by lower mean sequence identity (67%). GraS sensor histidine kinase was less conserved within species groups A (mean 69%) and B (mean 66%), compared with GraS of species group C that shared greatest mean sequence identity (88%), albeit that group C defined here is a small, related species set. Both BraS and GraS sensor proteins have lower sequence conservation across staphylococci than BraR and GraR (Table 3). The reduced divergence of these response regulators might reflect their relative isolation from selection by the external environment and differential stimuli.

Responses to cationic AMPs in the staphylococci are complex (Peschel and Collins, 2001; Herbert et al., 2007) and ligand specificity could account for species divergence of GraSR and BraSR TCS. This evolutionary outcome could be explained with strong selection pressure driven by ubiquity and diversity of cAMPs in staphylococcal niches. One intrigue in our analysis is the absence of GraSR and presence of only BraSR TCS in the group A species, S. pettenkoferi, with the sole related sensor protein having a mean sequence identity of 27% with group A GraS but 58% with group A BraS. S. pettenkoferi BraR has a mean sequence identity of 47% with GraR and 73% with BraR from group A. These values support the S. pettenkoferi TCS is a BraSR ortholog. GraSR was also absent and BraSR present in the four additional publicly available S. pettenkoferi genome sequences (strains 1286_SHAE, 589_SHAE, UMB0834 and CCUG 51270). The absence of GraSR in S. pettenkoferi raises questions about the evolution of BraSR in group A staphylococci. Gene duplication of GraSR in a group A species, with subsequent sequence divergence over time to BraSR and spread throughout group A species by horizontal gene transfer, is tempting to suggest. S. pettenkoferi having BraSR but not GraSR presents a challenge to this paralog hypothesis. We propose two possibilities; S. pettenkoferi may have suffered deletion of graSR following acquisition of braSR, or S. pettenkoferi never acquired braSR, but rather its TCS evolved from ancestral genes. Such a scenario would enable group A organisms to acquire braSR from S. pettenkoferi as an additional and sufficiently divergent TCS locus.

Staphylococcus species genomes sequenced recently were investigated for their encoded GraS and BraS protein homologues, which supported the limited distribution of BraS in staphylococci as identified in the Random Forest analysis (Table 4 and Supplementary File S3). Furthermore, it revealed additional species encoding BraSR but not GraSR (S. agnetis, S. auricularis, S. chromogenes, S. hyicus, S. massiliensis). Regardless of the origins of both TCSs, the divergence between and within GraSR and BraSR likely reflect specificities for their ligands and selection driven by the niches to which the staphylococci are specialized.

Table 4

Presence and absence of GraS and BraS across 49 Staphylococcus species.

GraS-encoding	GraS and BraS-	BraS-encoding
only	encoding	only
S. arlettae	S. argenteus	S. agnetis
S. carnosus	S. aureus ∗	S. auricularis
S. cohnii	S. capitis ∗	S. chromogenes
S. condimenti	S. caprae ∗	S. hyicus
S. delphini	S. devriesei	S. massiliensis
S. edaphicus	S. epidermidis ∗	S. pettenkoferi ∗
S. equorum	S. haemolyticus ∗
S. felis	S. hominis ∗
S. fleurettii	S. lugdunensis ∗
S. gallinarum	S. pasteuri
S. intermedius	S. petrasii
S. kloosii	S. saccharolyticus
S. lentus	S. schweitzeri
S. lutrae	S. simiae ∗
S. microti	S. warneri ∗
S. muscae
S. nepalensis
S. piscifermentans
S. pseudintermedius
S. rostri
S. saprophyticus
S. schleiferi
S. sciuri
S. simulans
S. stepanovicii
S. succinus
S. vitulinus
S. xylosus

GraSR and BraSR-Associated ABC Transporters

Both GraSR and BraSR, as members of the BceS-like IM-HK family of TCS, are activated by AMP ligand bound to an associated ABC transporter (Mascher, 2014). Given the important function of these TCS, the conservation of their associated transporter protein sequences was compared across the staphylococci.

VraFG is the GraSR-associated ABC transporter (Meehl et al., 2007) and in the genomes encoding VraFG (absent from group B species and S. pettenkoferi) there is a high degree of shared protein sequence conservation. VraF has a mean sequence identity of 68% across the staphylococci examined (Table 1), with greatest conservation within species groups (group A, 79% identity; group B, 85.3% identity; group C, 96.8% identity). Shared sequence identity among the VraG proteins was 47.5%, with 88%, 65.2 and 61.9% identity within groups A, B, and C, respectively. The BraDE ABC transporter associated with BraSR was identified in group A species and, similar to VraFG, revealed greater identity (68.4%) across BraD sequences compared with BraE (38.9%) protein sequences. Divergence within BraSR and GraSR-associated transporters has likely arisen from concurrent evolution of the ABC transporter specificities for AMPs.

Experimental Evolution of Nisin Resistance in S. aureus, S. hominis and S. saprophyticus

Previous studies demonstrated that selection by experimental evolution identified mutations conferring antimicrobial resistance in overarching regulators, notably SNPs in braS revealed roles for BraSR in nisin sensing and survival (Hiron et al., 2011). Following our identified species association of BraSR to group A staphylococci, we adopted an experimental evolution strategy to interrogate the contributions of GraSR and BraSR TCS under selection for nisin resistance.

Strains of group A species, S. aureus and S. hominis plus group B S. saprophyticus were each serially passaged in triplicate cultures with increasing concentrations of nisin using a microtiter plate method, with an equivalent sodium citrate buffer control passaged in parallel. Stepwise increases in nisin MIC were observed for all strains tested with no obvious pattern in the rate of resistance acquisition between the species. After selection, both S. aureus 171 and S. aureus SH1000 strains exhibited ∼100-fold increases in nisin MIC, a greater fold increase in resistance than that observed by Blake et al (2011; Hiron et al., 2011), which may be due to experimental design differences. Selection of both S. hominis strains increased nisin MIC ∼25-fold, and S. saprophyticus strains CCM_883 and CCM_349 showed 80-fold and 5-fold increases, respectively. Multiple clones of S. aureus 171, S. hominis J31 and S. saprophyticus CCM883 were genome sequenced to identify sequence variants that potentially contributed to increased nisin MIC. T0 genomes were assembled and annotated, then reads from three pools (each comprising 5 independent clones) and one individual clone of each experimentally evolved species were aligned to their respective assembled genomes to identify sequence variants (SNPs, insertions/deletions) specific to nisin selection (Tables 5–7).

Table 5

Non-synonymous, homozygous SNPs from independent clone pools of staphylococci after nisin selection.

Gene ID (Prokka)	Protein ID	Pool	Position	Base change	Amino acid change	Allele frequency
S. aureus 171
walk	WalK Sensor kinase	1	17119	A - > G	H364R	1
gltB_1	Glutamate synthase	1	437361	A - > T	Q797L	1
rpoB	DNA-directed RNA polymerase subunit beta	1	520176	C - > T	H506Y	1
mraY	Phospho-N-acetylmuramoyl-pentapeptide transferase	1	1103071	G - > A	V266I	1
yhcF	Transcriptional regulator	1	1998279	G - > A	G73R	1
Hypothetical	Membrane protein	2	615003	C - > T	Q57∗	0.35
rpoC	DNA-directed RNA polymerase subunit beta	3	523690	G - > T	A448S	0.99
femB	FemB	3	1323450	G - > A	R215H	0.63
Phage terminase	Terminase	3	1491745	G - > C	G240A	1
greA	GreA	3	1625376	A - > T	L76∗	0.37
bras	BraS sensor histidine kinase	3	2627088	C - > T	T175I	1
S. hominis J31
rpoB	DNA-directed RNA polymerase subunit beta	1	41317	G - > T	D1046Y	0.79
graS	GraS sensor histidine kinase	2	147503	C - > T	S120L	0.42
ftsH	FtsH Zinc metalloprotease	2	2172470	C - > A	D171E	1
gmk	Guanylate Kinase	3	552630	G - > A	R135H	1
yhcF	Transcriptional regulator	3	1212777	C - > T	Q47∗	1
S. saprophyticus 883
graS	GraS sensor histidine kinase	1	2068987	G - > T	G209C	1
codY	CodY regulator	2,3	1537681	C - > A	L79F	1,1
pitA	Phosphate transporter	2,3	2066724	C - > T	A195V	0.99,1
graS	GraS sensor histidine kinase	2,3	2069134	C - > A	R160S	1,1
marR family	MarR regulator	2,3	2136907	C - > T	T62I	1,1

Table 7

INDELs from nisin selection pools and single clones of S. aureus, S. hominis, and S. saprophyticus.

Source	Gene ID (Prokka)	Protein ID	Location	Base change	Effect
S. aureus 171
Single, Pool 3	lpl2_2	Lipoprotein	403168	195_196 insGG	I66fs §
Single, Pool 1,3	sdrE_1	MSCRAMM family adhesin	554844	2672_2673 insC	K891fs §
Single			1990375	-1_-1 insCC	uv
Single, Pool 2,3			2052262	∗1530_∗1530 delTG	fs
Single, Pool 3	deoC2	Deoxyribose-phosphate aldolase 2	2126552	450_451 insAG	K151fs §
Single, Pool 1	Hypothetical	Hypothetical protein	2471437	116_117 insT	E40fs
Single, Pool 3	fnbA_2	Fibronectin binding protein	2490720	1763_1764 delCG	S588fs §
Pool 1	Hypothetical	Hypothetical protein	1556011	201 delT	S67fs §
Pool 2			1337112	-1_-1 insATG
Pool 2,3	Hypothetical	Transposase	1337623	209_211 delAAG	E70del
Pool 2	Hypothetical	Transposase	1819467	1047_1048 insC	∗350fs §
Pool 2	leuA_2		2052046	∗1530 delC	dv
Pool 2	Hypothetical	Hypothetical	2471438	115_116 insATA	P39del, insHT
Pool 3	lpl2_1	Hypothetical	402338	212_213 insCT	Q71fs §
Pool 3	sdrD_1	MSCRAMM family adhesin	550856	3237_3238 insC	M1080fs §
Pool 3	Hypothetical		555082	-1_-1 insG
Pool 3	Hypothetical	LPXTG surface protein	2480756	216 delA	T72fs §
S. hominis J31
Single, Pool 2	ftsH	Zinc metalloprotease	2172279	325 delA	S109fs §
Pool 1	relA	GTP Pyrophosphokinase	951655	764 delA	Q255fs §
Pool 2	ssaA2_2	CHAP domain containing protein	1437614	575_578 delGTTA	G192fs §
S. saprophyticus 883
Single, Pool 1,2,3	ptsG	PTS alpha-glucoside transporter subunit IIBC	604054	-1_-1 insAA	uv

Nisin-Selected SNPs in Staphylococci

Experimental evolution of S. saprophyticus identified a SNP in graS (GraS: A160S; Table 5) that was present in two clone pools, and SNP graS G209C in a third pool. A single clone sequenced from the latter pool identified only one SNP in graS (GraS: G209C) and an upstream variant associated with ptsG (Table 5). These data provide support for GraSR contributing to nisin resistance in S. saprophyticus given the absence of the BraSR TCS in this group B Staphylococcus species. Aside from TCS, other regulators may contribute to the nisin response in S. saprophyticus as evidenced by an identical SNP identified in two separate nisin resistance selections (pools 2 and 3) corresponding to a T62I change in an uncharacterized MarR transcriptional repressor.

In both S. aureus and S. hominis there are multiple pathways to high-level nisin resistance. Each species revealed SNPs in TCS systems, but these differed across the parallel selection experiments (Tables 5–7). In S. aureus, a non-synonymous SNP in braS (BraS: T175I) was present in 100% of reads from one sequenced pool, differing from previous work that identified a discrete braS SNP (BraS: A208E) (Hiron et al., 2011). Evidence for a second TCS contributing to nisin resistance arose from a walK non-synonymous SNP (WalK: H364R) within the diverse and flexible signal sensing PAS domain of WalK in S. aureus (Hefti et al., 2004). WalKR is essential and functions to maintain cell wall metabolism (Delaune et al., 2012) and SNPs in this TCS contribute to vancomycin and daptomycin resistance due to cell wall-thickening (Howden et al., 2011). Should this cell wall phenotype be associated with the H364R WalK variant it could similarly limit nisin interaction with its lipid II target to abrogate pore formation. A large overlap was reported between the WalKR and GraSR regulatory networks in S. aureus (Falord et al., 2011).

In S. hominis, a graS SNP (GraS: S120L) was present in 2 clones of sequence pool 2 and no SNPs or other sequence variants were identified in braSR (Tables 5–7). S. hominis has both braSR and graSR loci and therefore it is intriguing nisin resistance selection resulted in SNPs in a different TCS to S. aureus despite encoding both, potentially reflecting differences in their contribution across group A staphylococci. A further transcriptional regulator might contribute to nisin resistance in both S. aureus and S. hominis, where the uncharacterized yhcF revealed SNPs producing G73R and N47∗, respectively; the presence of SNPs in yhcF of both species supports a role for this regulator. The YhcF transcriptional regulator proteins of S. aureus and S. hominis have 75% similarity and their cognate genes are adjacent to an ABC transporter locus with potential specificity for GlcNAc, which might catalyze recycling of cell wall substrates from nisin damage. The role of this operon is currently being investigated.

In summary, we have identified differential encoding and diversity of antimicrobial resistance regulators and their associated transporters across the staphylococci. Our previous studies of the nasal microbiome correlated cumulative antimicrobial production with community structure, limitation of invasion and S. aureus exclusion (Coates et al., 2014; Libberton et al., 2014, 2015). Further dissection of antimicrobial sensing and discrimination via the TCS systems BraSR and GraSR combined with analysis of their associated transport specificities will provide information that can be layered with niche-relevant antimicrobial activities from competing species. Such analyses are now emerging and will provide a more holistic determination of Staphylococcus ecology.

Materials and Methods

Staphylococcus Orthologous Gene Content

Representative genomes of 21 different Staphylococcus species available at the time of analysis (Table 1) were either sequenced (see later section) or retrieved from the NCBI FTP repository[1]. Complete genomes were used where possible. Draft genomes available as NCBI scaffolds were reordered against an appropriate reference using a bespoke perl script. Genomes were annotated using PROKKA (version 1.5.2) (Seemann, 2014) to ensure consistent gene calling and annotation. OrthoMCL (version 1.4) was used to cluster orthologous proteins (Li et al., 2003), with input parameters, e-value cut-off: 1e-5, percentage identity cut-off: 30, percentage match cut off: 20. Briefly, initial BLAST steps of orthoMCL used the latter two low stringency cut-off values; these values were used to retain more proteins for clustering from these BLAST stages. Inparalog, ortholog and co-ortholog pairwise relationships were generated through reciprocal best and better hits in subsequent stages that used the p-value cut-off of 1e-5. Finally, the MCL (Markov clustering) aspect of the tool was applied to these pairwise relationships to allow clustering into orthologous groups (Li et al., 2003; Fischer et al., 2011). A bespoke python script was used to create a table describing the presence or absence of each OrthoMCL cluster within every genome. These data were converted to a matrix for analysis in the statistical package R and a heatmap was generated from the matrix. To control for gross strain-specific effects on the heat map (and thus OrthoMCL clusters), this step was repeated by substituting with alternative strains (Supplementary Table S1) and all permutations were analyzed in subsequent steps of the analysis.

Drivers of OrthoMCL Group Formation

The R library, Random Forest (version 4.6-7) (Liaw and Wiener, 2003) was used to investigate the genetic inputs directing classification of the species into their OrthoMCL groups. A presence/absence table of each of the orthologous groups obtained from the United States300 permutation of the OrthoMCL analysis was generated using a bespoke python script and used as the input data for the Random Forest algorithm.

The data was split into a test and training data set with both sets including equal proportions of group A species. The optimum value for mtry was found to be 66 using the tuneRF function (ntree = 1001, stepFactor = 1.5, improve = 0.001). These mtry and ntree parameters resulted in a model with an out of bag (OOB) error rate of 9.09% and area under ROC curve (AUC) of 0.96.

Data output was summarized using the variable importance plot function and the numeric mean decrease in accuracy (MDA) resulting from the permutation of each variable was obtained through the importance function; these data were used as the measure of the importance of each variable. The maximum MDA in this analysis was 3.3. Clusters were mapped back to the genome and the annotation of protein sequence for a species representative of each cluster was retrieved. Protein sequences of clusters identified as important were retrieved and their annotations curated and verified against published annotations. In addition, outputs were generated by substituting strains of species in the analysis to compare conservation of identified clusters between the variable importance plots. Sequences of protein clusters from the single species representative in Table 2 and identified by Random Forest output are listed in Supplementary Files S2, S3. Protein sequences were retrieved from their respective genomes and alignments were performed using ClustalW2 (version 2.1).

Minimum Inhibitory Concentration Assay

Nisin (Sigma-Aldrich Company Ltd, United Kingdom) was prepared as a 20 mg mL-1 solution in 10 mM sodium citrate (Sigma–Aldrich Company Ltd, United Kingdom) at pH 3 and stored at 4°C. MIC assay used microtiter plates with doubling dilutions of nisin in BHI (Thermo Scientific) inoculated 1 in 2 with 100 μL bacterial suspension adjusted to OD600 0.2 ± 0.005. The lowest concentration with an optical density ≤ to that of the initial optical density was taken as the minimum inhibitory concentration (MIC).

Selection for Nisin Resistance

Experimental evolution was performed by serial passage in broth containing doubling dilutions of nisin in triplicate wells of a microtiter plate. For selection of S. aureus and S. saprophyticus, the maximal assay concentration of nisin was 5 mg mL-1 and for S. hominis 50 μg mL-1. Control selection experiments with equivalent sodium citrate concentrations were performed in parallel. Experiments were initiated with inoculation of bacteria to OD600 = 0.2 for the first passage and plates were incubated static at 37°C. Bacteria growing at the highest concentration of nisin after 24–48 h were passaged forward to the next plate; subsequent passages were inoculated with a 1:1000 dilution of culture. Serial passage was continued until growth occurred at the maximal nisin concentration (for strains S. saprophyticus = 10 mg ml-1, S. aureus = 10 mg ml-1, and S. hominis = 250 μg ml-1) or for a period of 12 days. All passaged cultures were collected and stored at -80°C in 20% (v/v) glycerol (Fisher Scientific) after each passage and the T0 time point served as comparator strain.

Colonies were randomly selected for sequencing after plating from independent biological replicate cultures that had reached an equivalent maximum level of nisin resistance. Clones from each repeat were selected and cultured in 10 mL of BHI at 37°C with shaking at 200 rpm overnight. Increased MICs were confirmed by using the MIC assay described above at the highest nisin concentrations. Selection was performed for a corresponding citrate control time point for each of the three species.

DNA Extraction, Library Preparation and Sequencing

Cells were harvested from overnight culture and lysed in buffer containing 12.5 μg ml-1 lysostaphin (Sigma–Aldrich) and 10 U mutanolysin (Sigma–Aldrich). DNA was purified using a DNeasy Blood and Tissue Kit (Qiagen). DNA (30 ng) from each of five selected clones was pooled to make Illumina Truseq DNA libraries with an insert size of 350 bp. In addition to three separate clone pools, a single clone was selected for sequencing from the clones used to constitute the pools. Single clones were selected on the basis of the highest DNA quality. The single clones and the T0 isolates were also sequenced using Illumina Truseq nano DNA libraries with 350 bp inserts.

Identification of SNPs and INDELS

T0 comparator strains were assembled using VelvetOptimiser (version 2.2.5; Victoria Bioinformatics Consortium) with Kmer sizes from 19 to 99 and Velvet version 1.2.06 (Zerbino and Birney, 2008). Annotation was carried out using PROKKA version 1.5.2 (Seemann, 2014). The PacBio assembly of S. hominis strain J31 (Accession FBVO01000000) (Coates-Brown and Horsburgh, 2017) was used as the comparator assembly for this strain. Good quality filtered reads from experimentally evolved pools and single clones were aligned to respective comparator strains using the BWA (version 0.5.9-r16) (Li and Durbin, 2009) packages aln and sampe, and also using BWA (version 0.7.5a-r405) mem package. SAM files were converted to bcf (binary variant call) files with samtools for SNP calling using the mpileup package. The bcf output file from mpileup was then converted to vcf (variant call format) files and quality filtered. For SNPs, only this quality filtered vcf file from the pooled clones, along with mpileup output without base data, were used to further filter the SNPs to include only those present in 33.33% of reads, which equates to the SNP being present in more than one clone. To reduce falsely called SNPs, SNPs not called from both alignments (from either BWA aln and sampe or BWA mem) were removed from the data set, as recommended by Li (Li, 2014). SNPs called in the control data and evolved isolates were filtered from the data.

Availability of Data and Materials

Genomes resulting from this work can be retrieved from the ENA database at EMBL-EBI (https://www.ebi.ac.uk/ena/data/view) under the bioproject accession PRJEB22856, including data from experimental evolution of S. aureus 171; Parental S. aureus 171 data accession: LT963437. Individual genome assembly accessions used in Figure 1 are listed in Table 1 and Supplementary File S1. Strains not already publicly archived are available on request. This manuscript was submitted to bioRxiv ahead of review (Coates-Brown et al., 2018).

Author Contributions

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

Conflict of Interest Statement

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.

Table 6

Non-synonymous, homozygous SNPs from single clones of S. aureus, S. hominis, and S. saprophyticus after nisin selection.

S. aureus strain 171 (single clone from pool 1)
Gene ID (Prokka)	Protein ID	Position	Base change	Amino acid change
walK	WalK Sensor histidine kinase	17119	A - > G	H364R
gltB_1	Glutamate synthase	437361	A - > T	Q797L
rpoB	DNA-directed RNA polymerase subunit beta	520176	C - > T	H506Y
mraY	Phospho-N-acetylmuramoyl-pentapeptide transferase	1103071	G - > A	V266I
msrR	Regulatory protein MsrR	1309845	G - > T	E181∗
yurK	Transcriptional regulator	1998279	G - > A	G73R
S. hominis strain J31 (single clone from pool 2)
ftsH	Zinc metalloprotease FtsH	2172470	C - > A	D171E
S. saprophyticus strain 883 (single clone from pool 1)
graS	GraS sensor histidine kinase	2068987	G - > T	G209C