Comprehensive Virulence Gene Profiling of Bovine Non-aureus Staphylococci Based on Whole-Genome Sequencing Data

Abstract

ABSTRACT

Non-aureus staphylococci (NAS) are the most frequently isolated pathogens from intramammary infection (IMI) in dairy cattle. Virulence factors (VFs) and mechanisms by which NAS cause IMI are not fully known. Herein, we analyzed the distribution of 191 VFs in 441 genomes of 25 NAS species, after classifying VFs into functional categories: adherence (n = 28), exoenzymes (n = 21), immune evasion (n = 20), iron metabolism (n = 29), and toxins (n = 93). In addition to establishing VF gene profiles, associations of VF genes between and among functional categories were computed, revealing distinctive patterns of association among VFs for various NAS species. Associations were also computed for low, medium, and high somatic cell count (SCC) and clinical mastitis (CM) isolates, demonstrating distinctive patterns of associations for low SCC and CM isolates, but no differences between high SCC and CM isolates. To determine whether VF distributions had any association with SCC or CM, various clustering approaches, including complete linkages, Ward clustering, and t-distributed stochastic neighbor embedding, were applied. However, no clustering of isolates representing low SCC, medium SCC, or high SCC or CM was identified. Regression analysis to test for associations with individual VF functional categories demonstrated that each additional toxin and host immune evasion gene increased the odds of having high SCC or CM, although an overall increase in the number of VFs was not associated with increased SCC or occurrence of CM. In conclusion, we established comprehensive VF gene profiling, determined VF gene distributions and associations, calculated pathogenic potentials of all NAS species, and detected no clear link between VF genes and mastitis.

IMPORTANCE Non-aureus staphylococci (NAS) are the most frequently isolated pathogens from milk in dairy cattle worldwide. The virulence factors (VFs) and mechanisms by which these bacteria cause udder infection are not fully known. We determined the distribution and associations of 191 VFs in 25 NAS species and investigated the relationship between VFs and disease. Although the overall number of VFs was not associated with disease severity, increasing numbers of toxin and host immune evasion genes specifically were associated with more severe disease outcomes. These findings suggest that the development of disease and the interactions of VFs with the host are complex and determined by the interplay of genes rather than just the presence of virulence genes. Together, our results provide foundational genetic knowledge to other researchers to design and conduct further experiments, focusing on understanding the synergy between VFs and roles of individual NAS species in IMI and characterizing species-specific effects on udder health.

INTRODUCTION

Non-aureus staphylococci (NAS), most of which are coagulase-negative staphylococci (CNS), are the most frequently isolated bacteria from bovine milk (1–3). Although NAS are often considered minor mastitis pathogens (3, 4), they are increasingly recognized as dominant pathogens of bovine mastitis worldwide (1, 3, 5). The genus Staphylococcus (as of October 2018) includes 53 species and 28 subspecies (http://www.bacterio.net/staphylococcus.html), of which 25 NAS species are commonly isolated from milk from dairy cows in Canada and other countries. Interspecies relationships and prevalence of these species were recently reported by our group (1, 6). However, pathogenesis of these bacteria is not fully understood. Therefore, it is not clear whether NAS should be considered commensal bacteria or opportunistic pathogens. Additionally, the effects of individual NAS species on udder health are not well characterized (7–9). Mechanisms that allow these organisms to colonize and cause mastitis are not well-known (2, 3, 10). Generally, obligate or opportunistic pathogenic bacteria sense host signals and adapt gene expression to match environmental conditions (11–13). A subset of these genes has a key role in the ability of the bacterium to cause disease (14, 15). Products of such genes that facilitate successful colonization and survival of the bacterium in a host environment and damage that host are considered pathogenicity determinants or virulence factors (VFs) (13, 15, 16). The VFs are either coded within bacterial genomes, mostly on specific genomic loci, referred to as pathogenicity islands, or coded within transmissible genetic elements that can spread among bacteria (13, 17–20). Numerous VFs are known for Staphylococcus aureus (SAU), an important pathogen of various animals, including dairy cattle, where it is recognized as a major udder pathogen, commonly associated with (sub)clinical mastitis (21, 22). Staphylococcus aureus has an array of genes involved in adhesion, invasion, and host defense evasion (23). The VFs and mechanisms by which SAU causes intramammary infections (IMI) in dairy cattle are well characterized, compared to NAS-related IMI (24–27). Few studies have focused on VFs of NAS isolated from bovine mastitis (28–31). All VFs and mechanisms by which various bovine NAS species survive, multiply, and cause disease in the host are yet to be fully identified. Thus, to determine pathogenic potential (number of VFs and their associations) of individual or closely related NAS species and to investigate whether the presence and absence of VFs have associations with occurrence of mastitis, comprehensive identification of putative VFs from the 25 most common NAS species is essential. Genome-based phylogeny of 25 bovine NAS species was recently established by our group (32), which divided NAS species into five distinct clades (A to E). Clade A contained SVI, SFL, and SSC (for NAS species and NAS species abbreviations, see “Classification and distribution of virulence factors” in Materials and Methods). Clade B included SAG, SHY, and SSC, clade C was represented by SSI. Clade D was divided into D1 (SHO, SDE, and SHA), D2 (SPA and SWA), and D3 (SEP, SCR and SCI), and clade E contained SAC, SAR, SKL, SSU, SGA, SCO, SNE, SEQ, SSA, and SXY.

In this study, genomic data for 441 isolates from 25 NAS species from 87 herds were used. Our objectives follow: (i) to identify and determine the distribution of putative VFs (pVFs) among 25 NAS species; (ii) to investigate relationships among pVFs; (iii) to identify distinct pVF associations for low, medium, and high somatic cell count (SCC) and clinical mastitis (CM) isolates; and finally (iv) to investigate the association between the presence of pVFs and CM.

RESULTS

Distribution and associations of VF genes involved in adherence.

The 28 VFs of the adherence category tested in this study included accumulation-associated protein (aap), biofilm-associated surface protein Bap (bap), autolysin (atl), clumping factors (clfA and clfB), collagen adhesion (cna), elastin binding protein (ebp), fibronectin binding proteins (ebh, efb, uafA, fnbA, and fnbB), extracellular adherence/MHC analogous protein (eap/map), cell wall surface anchor family proteins (sasC, sasG, and sasP), intercellular adhesins (icaA, icaB, icaC, icaD, and icaR), and Ser-Asp-rich fibrinogen binding proteins (sdrC, sdrD, sdrE, sdrF, sdrG, sdrH, and sdrI) (Fig. 1A). Among these genes, atl was the most frequently detected gene in 20 out of 25 NAS species in at least one isolate. The atl gene was present in SCH, SVI, SSU, SCO, SGA, SHA, SHY, SHO, SDE, SPA, SWA, SEP, SCR, SCI, SNE, SEQ, SSA, and SXY but was not detected in isolates of SFL, SSC, SSI, SAC, SKL, and SAR. The icaC gene of the ica operon, believed to be involved in biofilm formation (33) was the second most frequent gene and was present in 17 of 25 NAS species (Fig. 1A). Other genes of the ica operon, namely, icaA, icaB, icaD, and icaR were detected in eight, seven, eight, and seven NAS species, respectively. Another biofilm-related gene, bap (encodes biofilm-associated surface protein Bap) was present in only six NAS species (Fig. 1A).

FIG 1

Distributions and pairwise associations of 28 adherence-related virulence factors. (A) Distribution of adherence-related virulence factors of 25 NAS species, arranged into five clades (A, B, C, D, and E) according to their placement in phylogenetic trees (32). Values and different colors within cells represent the percentage of isolates containing a given gene. The color gradient from red to dark green indicates increasing percentage of isolates containing the particular gene within species (light yellow, 1 to 25% of isolates; light blue, 26 to 50% of isolates; light green, 51 to 75% of isolates; dark green, 76 to 100% of isolates). (B) Pairwise associations of genes, computed using phi coefficients (143), with a color gradient representing the type of association (blue for positive; red for negative), while the intensity of the color and size of the circle show the strength of the association.

Among the clade A species, SVI contained only 2 (alt and icaC) of 28 adhesion genes, whereas SFL and SSC contained 4 and 5 genes from the ica operon (Fig. 1A). Among species of clade B (SAG, SHY, and SCH), atl, clfB, cna, efb, uafA, fnbA, and fnbB were the most widespread genes. SSI, the sole representative species of clade C (in our data), contained only icaA and icaR. Many genes from the adherence category, including aap, clfA, eap/map, sasG, sraP, sdrD, sdrE, and sdrI, were not detected in any of the 441 NAS isolates. One or more genes from the intracellular adhesin family (ica operon) were present in many NAS species, except species from clades B and D1. Two genes, fnbA and fnbB, responsible for encoding fibronectin binding proteins, were identified in only three species (SAG, SHY, and SCH) from clade B. The cna gene, which contributes to collagen adhesion, was detected in two (of three) species from clade B (SAG and SHY) only. The distributions and frequencies of other clade-specific VF genes are shown in Fig. 1A. In this study, a VF was considered present in the NAS species, even if that VF was detected only in one isolate of that particular species.

Associations between 28 adherence genes in NAS species are shown (Fig. 1B). A strong positive association was detected between fnbB and cna, with both fnbB and cna present in 100% isolates of SAG and SHY (Fig. 1A). There were also positive associations for a second fibronectin binding fnbA and cna, however, with lower intensity compared to the previous example shown for fnbB, corresponding to frequencies of 62 and 67% in isolates of SAG and SHY (Fig. 1A), although cna was present in 100% of isolates of these species. Similarly, strong positive associations were present among icaA, icaB, and icaD of the ica operon (Fig. 1B). Generally, within the adherence category, most associations were either neutral or positive; no strong negative association was detected within this category, and all negative associations were very weak.

Distribution and associations of exoenzymes.

The second category of VFs (21 genes), exoenzymes (Fig. 2A), consisted of adenosine synthase A (adsA), aureolysin (aur), cysteine proteases (sspA, sspB, sspC, sspD, sspE, and sspF), hyaluronate lyase (hysA), lipases (lip and geh), serine proteases (splA, splB, splC, splD, splE, and splF), staphylocoagulase (coa), staphylokinase (sak), thermonuclease (nuc) and von Willebrand factor binding protein (vWbp). Thermonuclease (nuc) was present in almost 100% of NAS isolates, except SVI, where it was identified in 17% of isolates. The second most frequent exoenzyme genes were aur and lip, which were detected in 19 and 18 NAS species, respectively. The gene for another lipase, geh, was detected in six NAS species. The enzyme adsA was identified in all species from clades B and C and in all species from clade D, except SPA and SEP, but in only one species from clade E (SKL) (represented by a single isolate). Clade-specific distributions of enzymes were observed for sspB, which was identified in all isolates of SPA and SWA (species from clade D2) and vWbp, detected in SAG (100%), SHY (100%), and SCH (94%) species from clade B. Hyaluronate lyase (hysA) was present in only two clade B species (SAG and SHY). Seven of the listed exoenzymes (sspD, sspE, and sspF and splA, splB, splC, and coa) were not detected in any NAS species. Distributions and frequencies of each exoenzyme described above are shown (Fig. 2A). Within exoenzymes, there was a strong positive association between splF and geh, whereas there was a strong negative association between sspA and adsA (Fig. 2B).

FIG 2

Distributions and pairwise associations of 21 exoenzyme genes in NAS species. (A) Distribution of exoenzymes genes of 25 NAS species, arranged into five clades (A, B, C, D, and E) according to their placement in phylogenetic trees (32). Values and different colors within cells represent the percentage of isolates containing a given gene. The color gradient from red to dark green indicates increasing percentage of isolates containing the particular gene within species(light yellow, 1 to 25% of isolates; light blue, 26 to 50% of isolates; light green, 51 to 75% of isolates; dark green, 76 to 100% of isolates). (B) Pairwise associations of genes, computed using phi coefficients (143), with a color gradient representing the type of association (blue for positive; red for negative), while the intensity of the color and the size of the circle show the strength of the association.

Distribution and associations of VF genes involved in host immune evasion.

The host immune evasion category consisted of 16 genes involved in capsular synthesis (capA, capB, capC, capD, capE, capF, capG, capH, capI, capJ, capK, capL, capM, capN, capO, and capP), chemotaxis inhibitory protein (chp), staphylococcal complement inhibitor (scn), staphylococcal protein A (spa), staphylococcal binder of immunoglobulin (sbi) gene. The distribution and frequencies of these VFs are shown (Fig. 3A). The chp gene was not detected in any NAS species. The scn and sbi genes were exclusively identified in species from clade B (SAG, SHY, and SCH) and were absent in other NAS species. The spa gene was present only in SHY (100%) and SSI (62%) species. Capsular genes were the most commonly detected and frequently distributed genes in this category. In this study, cap genes were considered present if hits were detected for either isoform cap5 or cap8 (34, 35). Among capsular genes, capM was detected in all 441 NAS isolates, and most NAS species contained ≥8 capsular genes, except SAC, SSA, and SNE which contained capM and capP only. SCH contained most of the capsular genes, except capN. However, except for capM, capO, and capP, which were detected in all isolates of SCH, frequencies of other genes were within 7 to 14% of isolates. Most capsular genes had positive associations, except capP which had a negative association with capD and capL. The strongest associations were between capA and capB and among capE, capF, and capG (Fig. 3B).

FIG 3

Distributions and pairwise associations of 16 host immune evasion genes. (A) Distribution of host immune evasion genes of 25 NAS species, arranged into five clades (A, B, C, D, and E) according to their placement in phylogenetic trees (32). Values and different colors within cells represent the percentage of isolates containing a given gene. The color gradient from red to dark green indicates increasing percentage of isolates containing the particular gene within species (light yellow, 1 to 25% of isolates; light blue, 26 to 50% of isolates; light green, 51 to 75% of isolates; dark green, 76 to 100% of isolates). (B) Pairwise associations of genes, computed using phi coefficients (143), with a color gradient representing the type of association (blue for positive; red for negative), while the intensity of the color and the size of the circle show the strength of the association.

Distribution and associations of VF genes involved in iron uptake and metabolism.

Twenty-nine genes involved in iron uptake and metabolism were detected, including 9 iron-regulated surface determinant genes (isdA, isdB, isdC, isdD, isdE, isdF, isdG, isdH, and isdI), 7 ABC transporter (also known as siderophore receptors) genes (htsA, htsB, htsC, sfaA, sfaB, sfaC, and sfaD), 12 staphyloferrin A and B synthesis-related genes (sirA, sirB, sirC, sbnA, sbnB, sbnC, sbnD, sbnE, sbnF, sbnG, sbnH, and sbnI) and 1 sortase B gene (srtB). Among 9 iron-regulated surface determinant genes, isdI was detected in 22 NAS species but was not present in SEP, SAR, and SSU. The isdC, isdE, and srtB genes were detected in SSC, SSI, SPA, SCR, SCI, and SAC only. The isdH gene was identified only in SPA (100%). Two genes from iron-regulated surface determinants (isdB and isdD) were not detected in any NAS species. Seven genes of ABC transporters (siderophore receptors) were uniformly distributed, with few exceptions, such as the absence of htsA from SVI, SFL, and SAR and the absence of htsB and htsC from SAR and SFL, respectively. Staphyloferrin A synthesis-related genes (sirA, sirB, and sirC) were detected more frequently than staphyloferrin B synthesis-related genes. sbnA, involved in staphyloferrin B synthesis, was the only gene from the iron uptake and metabolism category identified in all 441 NAS isolates. Distributions and frequencies of all 29 VFs involved in iron uptake and metabolism are shown (Fig. 4A). Both negative and positive associations were detected among VFs of iron uptake and metabolism. For instance, siderophore receptor genes (ABC transporters) were mostly positively associated (Fig. 4B). Similarly, staphyloferrin B synthesis genes also had positive associations. Negative associations were present among the ABC transporter htsA and htsB (Fig. 4B). Iron-regulated surface determinant isdC and isdE had strong positive associations, whereas isdG and isdI had a negative association (Fig. 4B). The srtB gene had strong associations with isdC and isdE (Fig. 4B).

FIG 4

Distributions and pairwise associations of 29 virulence factors related to iron uptake and metabolism. (A) Distribution of iron uptake- and metabolism-related virulence factors of 25 NAS species, arranged into five clades (A, B, C, D, and E) according to their placement in phylogenetic trees (32). Values and different colors within cells represent the percentage of isolates containing a given gene. The color gradient from red to dark green indicates increasing percentage of isolates containing the particular gene within species (light yellow, 1 to 25% of isolates; light blue, 26 to 50% of isolates; light green, 51 to 75% of isolates; dark green, 76 to 100% of isolates). (B) Pairwise associations of genes, computed using phi coefficients (143), with a color gradient representing the type of association (blue for positive; red for negative), while the intensity of the color and the size of the circle show the strength of the association.

Distribution and association of toxin, type IV secretion, and phenol-soluble modulin genes.

The identification and distributions of 93 toxin genes (Fig. 5 and 6) were determined in NAS isolates. Figure 5 includes 36 toxin genes from various categories, including 6 genes for alpha, beta, delta, and gamma hemolysins (hly/hla, hlb, hld, hlgA, hlgB, and hlgC), 4 genes for leukocidins, including leukocidin M (lukM and lukF-like) and Panton-Valentine leukocidins (lukS-PV and lukF-PV), 2 leukotoxins (lukD and lukE), toxic shock syndrome toxin (tsst), 4 exfoliative toxins (eta, etb, etc, and etd), 8 genes of type VII secretion system (esaA, esaB, esaC, essA, essB, essC, esxA, and esxB), and 11 genes for phenol-soluble modulins, including the 5 alpha modulins (PSMα1, PSMα2, PSMα3, PSMα4, and PSMmec) and 6 beta modulins (PSMβ1, PSMβ2, PSMβ3, PSMβ4, PSMβ5, and PSMβ6). The majority of these toxin genes were not identified in NAS species. Among hemolysin genes, beta-hemolysin (hlb) was identified from species of clade B (SAG, SHY, and SCH) and clade D3 (SEP, SCR, and SCI) and in 4% of isolates of SXY. The delta-hemolysin gene (hld) was detected in 16% of isolates of SWA only. Genes for other hemolysins (alpha and gamma), leukocidins, and leukotoxins were not detected in any NAS isolates. Similarly, tsst, etc, and etd were also not detected in NAS. Among the eight known genes of the type VII secretion system, esaC and esxB were not identified in NAS. The other six genes of secretion system (esaA, esaB, essA, essB, essC, and esxA) were detected in SVI, SAG, SHY, SCH, SEP, SCR, SAR, SGA, SCO, SSA, and SXY in frequencies ranging from 4% to 100%. Among PSMs, none of the PSMα genes were identified in NAS. In contrast, PSMβ genes were present in various NAS species. Detection and distribution patterns of 21 enterotoxins and 36 staphylococcal exotoxins (SETs) genes are shown (Fig. 6). These genes were not identified in the majority of the NAS isolates except for some species in clade B, especially SAG and SHY. For instance, enterotoxins, sed, seg, seh, sei, sej, seln, selo, and selp were identified in 15% of SAG isolates, whereas 33% of isolates of SHY contain seb, seg, seh, sej, sell, and selm. Similarly, from 36 SETs, only 5 SET genes were detected in 1 or more species of clade B. Of 93 toxin genes, 32 were detected in 1 or more NAS isolates (Fig. 5 and 6). There were strong positive associations among type VII secretion system genes (esaA, esaB, essA, essB, essC, and esxA) (Fig. 7). In addition, there were positive associations for some enterotoxins and staphylococcal exotoxins (Fig. 7). However, these genes were mostly identified only in SAG and SHY species (Fig. 6).

FIG 5

Distributions of toxin genes of different categories. Distribution of 36 toxin genes of 25 NAS species, arranged into five clades (A, B, C, D, and E) according to their placement in phylogenetic trees (32). Values and different colors within cells represent the percentage of isolates containing a given gene. The color gradient from red to dark green indicates increasing percentage of isolates containing the particular gene within species (light yellow, 1 to 25% of isolates; light blue, 26 to 50% of isolates; light green, 51 to 75% of isolates; dark green, 76 to 100% of isolates).

FIG 6

Distributions of enterotoxin and staphylococcal exotoxin genes. Distribution of 21 enterotoxin (A) and 36 staphylococcal exotoxin (B) genes of 25 NAS species, arranged into five clades (A, B, C, D, and E) according to their placement in phylogenetic trees (32). Values and different colors within cells represent the percentage of isolates containing a given gene. The color gradient from red to dark green indicates increasing percentage of isolates containing the particular gene within species (light yellow, 1 to 25% of isolates; light blue, 26 to 50% of isolates; light green, 51 to 75% of isolates; dark green, 76 to 100% of isolates).

FIG 7

Pairwise associations of toxin genes. Pairwise associations of toxin genes, computed using the phi coefficient (143). Colors represent the type of association (blue for positive; red for negative), while the intensity of the color and the size of the circle show the strength of the correlation.

Associations between VF genes of different categories.

Relationships between the VFs from five categories were investigated using an association plot in which VFs from the VF functional categories adherence, exoenzymes, host immune evasion, iron uptake and metabolism, and toxins were analyzed concurrently (Fig. 8A). Most associations among the VF functional categories were neutral or very weak (Fig. 8B), with some moderate to strong associations. For instance, adhesin genes (icaA, icaB, icaC, and icaD) were positively associated with iron regulatory genes (Fig. 8A). Similarly, capsular genes (capE, capF, capG, capH, and capI) were positively associated with staphyloferrin B synthesis-related genes (Fig. 8A). To assess differences in associations among four SSC classes, association plots were generated separately for isolates from low, medium, and high SCC and from CM cases. Differences in association patterns between the isolates from low SCC and isolates from CM are shown (Fig. 9). There were many distinct positive and negative association patterns in CM isolates. For example, PSMβ5 and PSMβ6 had very strong positive associations with the icaA, icaB, and icaD genes of the ica operon. Similarly, there were strong positive associations for staphyloferrin B synthesis-related genes and cna, hysA, and sbi, whereas these positive associations were absent in low SCC isolates. Also, there were many strong negative associations in CM isolates, but absent in low SCC isolates. For instance, sirB and sirC from the staphyloferrin A synthesis family had strong negative associations with collagen adhesion genes (icaC and icaR) with isdE, isdF, and isdF. Other important differences in gene associations between low SCC and CM isolates are shown (Fig. 9). Medium and high SCC-specific associations are shown (see Fig. S1 and S2 in the supplemental material). Finally, to discover species-specific differences in VF associations between NAS species, graphs were generated for each NAS species (n = 25). An example showing the species-specific difference in association between SCH and SSI is presented (Fig. 10). In this case, no strong positive and negative associations were observed in SSI, except for capL, which was strongly negatively correlated with capH, whereas in SCH, there were strong positive associations among capsular genes and type VII secretion system genes (Fig. 10). However, there were no strong negative associations in SCH VF genes.

FIG 8

Pairwise associations and histogram of virulence factors from different functional categories. (A) Pairwise associations between virulence genes of five functional categories: adherence, exoenzymes, host immune evasion, iron uptake, and toxins. The associations were computed using the phi coefficient (143). (B) Histogram shows the total numbers of associations that fall within a particular association type. Colors in both panels represent the type of association (blue for positive; red for negative). Red boxes in panel A represent examples of clusters of positive associations between different functional categories.

FIG 9

Pairwise associations and comparison of low somatic cell count isolates with clinical mastitis isolates. Mirror image showing difference in VF associations between low somatic cell count and clinical mastitis isolates. The top right triangle shows the pairwise associations of genes from clinical mastitis isolates. The bottom left triangle represents associations of virulence genes detected in low somatic cell count isolates. The associations were computed using phi coefficient (143). Colors represent the type of association (blue for positive; red for negative), while the intensity of the color and the size of the circle show the strength of the correlation.

FIG 10

Comparison of pairwise associations between Staphylococcus chromogenes and Staphylococcus simulans. Mirror image showing differences in VF associations between two NAS species. The top right triangle shows the pairwise associations of virulence genes identified in Staphylococcus chromogenes. The bottom left triangle represents associations of virulence genes detected in Staphylococcus simulans. The distinctive associations between two species are marked with red boxes. The associations were computed using the phi coefficients (143). Colors represent the type of association (blue for positive; red for negative), while the intensity of the color and the size of the circle show the strength of the correlation.

Association between the presence of virulence factors and mastitis.

Overall, an increase in the number of putative VFs was not associated with an increase in log SCC, although upon stratification by type of VFs, the presence of each additional toxin gene was associated with a 0.024 increase in log SCC (P = 0.006). None of the other VF types were associated with changes in log SCC (Table 1). With inclusion of CM samples, association of numbers of VF genes were tested, with the probability of greater inflammation (i.e., higher SCC), with CM having the strongest host response. The presence of each additional VF gene associated with host immune evasion increased the odds of having a more severe immune response by 0.074 (P = 0.003) (having one more host immune evasion gene made the isolate 1.07 times more likely to cause a more severe inflammation [i.e., increased SCC and/or CM]). Other types of VFs, however, were not associated with an increased risk of having a more severe immune response (Table 1). To investigate whether any unique pattern of the presence and absence of VF could predict disease outcome, a decision tree was generated. Although this revealed many unique patterns of VF distributions, none of these patterns were clearly associated with the level of host immune responses (low SCC, medium SCC, high SCC, and CM). In dendrograms, created by Ward clustering and complete-linkage clustering methods, placement of isolates (from low, medium, or high SCC, and CM) was polyphyletic. Isolates from all categories were randomly distributed over dendrograms (Fig. S3 and S4). No single clusters dominated by isolates from a single category were identified. Rather, most of the clustering in these dendrograms was according to species. Isolates from the same species, regardless of their isolation stage (low, medium, or high SCC or CM) were grouped together. Similarly, based on distribution patterns (generated using t-SNE algorithms), most NAS isolates clustered according to their respective species (Fig. 11A). For example, isolates from SAG, SHY, SAR, SEP, SEQ, SCI, and SSI were grouped in distinct clusters according to respective species (Fig. 11A). However, isolates from SCH were grouped in two clusters, whereas intermixing of isolates was also present for some isolates of various NAS species (Fig. 11A). Similarly, a t-SNE plot was also generated after grouping isolates into low, medium, or high SCC or CM. However, there were no unique clusters or unique patterns distinct for isolates from each of these categories (Fig. 11B). Interestingly, a comparison of clusters between (Fig. 11A and B) revealed that clustering in Fig. 11B was based on species.

TABLE 1

Regression model results showing the association between numbers of VFs and log SCC and inflammation severity

VF gene type	No. of VF genes	Linear regressiona			Logistic regressionb
		Coefficient	Standard error	P value	Coefficient	Standard error	P value
Adherence	28	0.067	0.036	0.065	−0.020	0.039	0.607
Exoenzyme	21	0.098	0.055	0.078	0.076	0.060	0.199
Host immune evasion	20	0.045	0.023	0.054	0.074	0.025	0.003
Iron uptake	29	0.011	0.029	0.714	−0.006	0.031	0.858
Toxin	93	0.067	0.034	0.045	0.019	0.040	0.638
Total	191	0.024	0.010	0.017	0.015	0.011	0.168

a

Mixed-effects linear regression using log SCC as the outcome. The coefficient represents the increase in log SCC (log number of cells/milliliter) for each additional VF detected.

b

Mixed-effects ordinal logistic regression using inflammation severity as the outcome (measured as low, medium, and high SCC and clinical mastitis in order of increasing severity). The coefficient represents the increase in odds of the inflammation being more severe for each additional VF detected.

FIG 11

Dimensionality reduced clustering for NAS isolates determined using t-SNE. (A) Clusters labeled by 25 NAS species. (B) Clusters labeled by the inflammation severity determined in the original milk sample: low SCC, somatic cell count of ≤150,000 cells/ml; medium SCC, somatic cell count between 150,000 and 250,000 cells/ml; high SCC, somatic cell count of >250,000 cells/ml; CM, clinical mastitis.

DISCUSSION

In this study, we determined the distribution of 191 Staphylococcus VF genes, across 441 isolates from all 25 known NAS species isolated from multiple Canadian dairy herds. Most earlier studies on VFs in NAS species were limited by the following: (i) using only a few NAS species (36), (ii) testing for a few and different VF genes (37–39), or (iii) collecting samples from very few herds (28). The first study identifying a large number of VFs in bovine-associated staphylococci was published recently (31). However, this study had only three NAS species. The 191 VFs tested in this study were grouped into five functional categories, and their distributions were determined in 441 isolates from 25 NAS species.

In this study, 28 VF genes involved in adherence and biofilm formation were tested in 25 bovine NAS species. Of these 28 adherence- and biofilm-related genes, the ica operon, encoding the polysaccharide intercellular adhesins (PIA), is the earliest recognized and most widely distributed genetic determinant of biofilms, mostly in human-associated NAS species (16, 40, 41). Biofilm formation starts with adhesion, initiated mostly by atl, followed by PIA production from the ica operon, which consists of icaA, icaB, icaC, icaD, and a promoter icaR. The icaA gene encodes a transferase enzyme, and icaB promotes deacetylation of PIA. Furthermore, icaC synthesizes N-acetylglucosamine, and icaD complements icaA action (41–43). In previous studies, there were contradictory results regarding the distribution of ica genes in bovine NAS species. For example, icaA was detected in ∼50% of NAS isolates by Tremblay et al. (37), whereas a much lower prevalence (∼5%) of icaA was reported by Piessens et al. (44). A difficulty in determining the role of the ica gene in biofilm formation by bovine NAS is that various studies have tested for different ica genes (37, 45, 46). We tested all genes and determined that icaC (17/25) was the most frequent gene, followed by icaA, icaD (8/25), and icaB (7/25). In contrast, most previous studies have either used icaA or icaD (or both) in biofilm studies (28, 46, 47), which may have underestimated the true prevalence of ica-based biofilm formation ability of NAS, especially in Canada.

None of the 25 NAS species tested in this study contained aap (encodes the accumulation-associated protein Aap) (48) which contributes to the primary attachment and establishment of intercellular connections in biofilms. Biofilms formed by Aap are proteinaceous in nature, in contrast to ica-dependent polysaccharide-based biofilms (48, 49). The absence of aap in our analysis, consistent with previous studies, indicated the ica-dependent biofilm-forming potential of our isolates (48, 49). However, distribution of bap (encoding biofilm-associated protein), the second most important biofilm-related gene, was sporadic. Previously, bap was described as a cattle-specific pathogenic factor of biofilm formation (50–52). However, in our study, this gene was detected only in species of clade D2 (SPA and SWA) and some species of clade E (SSU, SGA, SCO, and SEQ) with modest distribution frequencies (44% to 53%). Interestingly, bap was not detected in SCH, SSI, SHA, SXY, and SEP, the five most prevalent species causing IMI in Canada (1). However, our results are in contradiction to a recent study of three bovine NAS species, which reported the presence of bap in SAG and in a few isolates from SCH (31).

After adhering to the host surface, bacterial pathogens often produce a variety of exoenzymes to neutralize the immune system, damage host tissue, and degrade complex macromolecules (53) to be used as nutrients. Consistent with the results of a recent study (31), nuc (encodes thermonuclease) was the most frequent exoenzyme gene in our study, followed by aur and sspA. Adenosine synthase A (adsA) was detected mostly in species of clades B, C, and D; it converts AMP to adenosine and is considered a critical VF of SAU (54). Among proteases, cysteine proteases (ssp locus encoded) were more widespread than serine proteases (spl encoded). In SAU pathogenesis, both serine and cysteine proteases have important roles in cleaving the activity of neutrophil serine proteases, which have key roles in immune defense (55–57). Lipases were frequently distributed among NAS, except in species of clade A (Fig. 2A). The exact role of lipase in NAS pathogenicity is not understood; however, lipases are involved in release of free fatty acids (octadecanoic acid) in blood plasma (58–60). These free fatty acids affect several immune system functions and thus indirectly enhance pathogenic potential (58, 61, 62). Interestingly, SAG, SHY, and SCH contain vWbp, which along with coa and clfA were considered sufficient to convert commensal SSI into an invasive pathogen (63). Importantly, coa and clfA were not detected in these species in our study. However, the presence of coa in SAG was reported recently (31). The variable results in the coagulase test of SAG, SHY, and SCH strains (64, 65) may be due to the presence and expression of vWbp in these species.

Apart from exoenzyme production, encapsulation is another strategy of pathogenic bacteria to avoid host immune responses. Staphylococci, especially SAU strains, are equipped with genes enabling production of capsular polysaccharide or capsule, to protect them from phagocytosis, enhancing virulence and persistence (66–68). Among 11 known capsular polysaccharide serotypes, cap5A-P and cap8A-P are the most widespread in SAU isolated from human and bovine infections (35, 69). In our study, capM was the most frequent cap gene and was detected in all 25 NAS species, followed by capP and capC. In other studies, depending on the animal model of infection, cap genes either enhanced or decreased virulence of SAU. For example, capsular production enhanced virulence of SAU in murine models of bacteremia (70) and surgical wound infection (71). In contrast, there was a lower virulence of capsular gene-containing SAU in IMI (72, 73) and in catheter-induced endocarditis (74) compared to corresponding capsule gene mutants. Interestingly, based on recovery of more capsule-negative isolates from human patients with osteomyelitis, mastitis, or cystic fibrosis, the lack of a capsule was suggested to be advantageous for SAU during chronic infections (72, 75, 76). In our study, 12 capsular genes (capA to capL), in agreement with the results of a recent study (31), were present in very low frequencies (7 to 11%) among SCH isolates, the most common species of NAS in IMI worldwide. The absence of capsular genes may explain the persistence of SCH in Canadian IMI (5, 9, 77).

Apart from capsular genes, other genes common in the immune evasion cluster (IEC), such as scn, chp, spa, and sbi, have important roles in host immune evasion (78, 79). Staphylococcal complement inhibitor (scn product) and the chemotaxis inhibitory protein (chp product) are thought to be highly specific for staphylococcal isolates of human origin (80, 81). Consistent with these studies, chp was not detected in any of the 441 bovine NAS isolates in this study. However, scn was detected in species of clade B (SAG, SHY, and SCH), which is in contrast with the results of a recent study (31) who did not detect scn in any NAS isolates.

All 25 NAS species contained at least one gene of the iron-responsive surface determinant (isd) operon, with isdI the most frequently distributed among all NAS species, which is consistent with the results of a recent study (31). Similar to most bacteria, staphylococci must acquire iron to replicate and to sustain infections (82–84). Iron becomes scarce during infections, as body fluids are actively depleted of free iron by the host to prevent bacterial growth, a process called “nutritional immunity” (82, 83, 85). Various iron acquisition mechanisms have key roles in the pathogenesis of SAU (82). Two principal mechanisms of iron uptake and metabolism are well studied in SAU (82–84). The first mechanism involves direct uptake of iron from molecules such as heme, using isd genes, whereas a second mechanism involves production of siderophores called staphyloferrin A and staphyloferrin B, along with surface transporters (82, 86, 87). Both isdI and isdG are necessary genes for SAU pathogenesis (88, 89). Most ABC transporter genes (siderophore receptors) were detected in NAS species.

Staphyloferrin A synthesis genes were more widespread than staphyloferrin B and isd genes in NAS species; therefore, we inferred that staphyloferrin A production was the predominant mechanism of iron acquisition among NAS isolates. Staphyloferrins have also been linked to increased virulence in human NAS infections (90). However, the frequent distribution of iron uptake and metabolism-related genes among bovine NAS species (Fig. 4A) may or may not be linked to their pathogenesis, and these genes may just be required for commensalism and survival in the host.

Virulence of staphylococci, especially SAU, can be related to their ability to produce a variety of toxins. These toxins can broadly be classified as cytotoxins (hemolysins, leukotoxins, and leukocidins) and superantigenic toxins (enterotoxins, exfoliative toxins, and toxic shock syndrome toxin [TSST]). In our study of cytotoxins, beta hemolysin (hlb), in accordance with the results of a recent study (31), was the most frequent and mainly detected in species of clade B (SAG, SHY, and SCH) and clade D3 (SEP, SCR, and SCI), whereas delta hemolysin (hld) was only detected in 16% of SWA isolates. Alpha and gamma hemolysins were not detected in any of our NAS isolates. A similar study in NAS of human and bovine origin, conducted in Iran, detected the presence of hla, hlb, and hld and the absence of hlg in bovine species (91). Similarly, both hla and hlb were not detected in 76 bovine NAS isolates, although that study was limited to one herd in China (28) and genes for delta and gamma hemolysin were not included. Recently, in India, gamma hemolysin gene was present in 10 bovine NAS species; however, other hemolysin genes (hla, hlb, and hld) were not tested (39). None of the leukocidins (lukM, lukF-like, lukS-PV, and lukF-PV) or leukotoxins (lukDE) were detected in our study. However, the presence of lukD in one isolate of bovine SSI was reported recently (31). Recently, the pvl gene was reported in a few NAS isolates (4/62) in India (39). It is noteworthy that Panton-Valentine leukocidin was strongly active against human neutrophils but only weakly active against bovine neutrophils (92–94); therefore, its absence in the present study, in agreement with the results of a recent study (31), was not surprising. Inconsistencies with other studies might have been due to methodical differences and/or inclusion of isolates from different geographical niches/areas. Of four exfoliative toxin genes, eta was detected in all isolates of SAG, SHY, and SEQ, whereas etb was only in SAG (8%) and etc and etd were not detected in any NAS species in our study. Previously, etb was detected in SAG (31); however, eta, etc, and etd were not tested in this study. The etd was detected in SXY (95), and in SAG and SHA (39) isolated from milk from cows with mastitis. We also detected the presence of six (esxA, esaA, esaB, essA, essB, and essC) type VII secretion system (T7SS) genes in multiple species. T7SS is a protein secretion pathway in Gram-positive bacteria (96–98). In SAU, T7SS is dispensable for laboratory growth, but it is essential for virulence (97, 99, 100). The six genes detected in our study code for core components of the secretion apparatus (97, 100). On the basis of the presence of all core genes of T7SS system, we inferred that these genes may have important roles, although further characterization is needed.

Alpha-type phenol-soluble modulins (PSMs) were not detected in any NAS species, consistent with previous findings, as α-type PSMs are considered more aggressive, and have mostly been associated with more virulent SAU strains (101). According to their length, PSMs can be classified as α-type (∼20 to 25 amino acids) or β-type (43 to 45 amino acids). In our study, most NAS species contained β-type PSMs, considered less aggressive forms of PSMs. Similarly, in previous studies, SEP produced β-type but not α-type peptides (102–104). PSMs have recently emerged as a novel toxin family of staphylococci and are considered major determinants of SAU virulence (94, 102, 105). PSMs have multiple roles in staphylococcal pathogenesis, including lysis of red and white blood cells, development of biofilm, and stimulation of inflammatory responses. Since PSM genes are encoded within the core genome, they are present in virtually all Staphylococcus species (102, 106, 107). In contrast to PSMs being present in all staphylococci, PSMβ were not detected in three species of clade A (SVI, SFL, and SSC) and two species of clade E (SAR and SKL). However, the absence of these genes may have been due to limitations of similarity search methods. Although not true for β-type PSMs, as they have been detected by similarity searches (107), PSMα, because of their short size, often do not yield meaningful results in similarity searches (102, 107).

Much of the toxicity of staphylococci can be attributed to superantigens. Among the superantigens, toxic shock syndrome toxin (tsst) gene was not present in any NAS species. This was consistent with several other studies (28, 108). Although tsst is mainly detected in SAU, the presence of this gene in NAS has been reported (39, 109). Enterotoxins and staphylococcal exotoxins, identified originally from SAU, have been studied extensively in staphylococcal isolates originating from humans (110, 111) and animals (108, 112, 113) and in animal-derived food products (114–116). Enterotoxins produced by some staphylococcal species, apart from interrupting host immune responses, also cause food poisoning. Therefore, the presence of these toxins in milk is of great concern for the dairy industry (115, 117, 118). In most studies, prevalence of enterotoxins was determined for NAS as a group, ignoring species-specific prevalences (108, 112, 113). In our study, 22 NAS species did not contain any enteroexotoxin genes. Interestingly, all enterotoxins containing species belong to clade B of NAS (32). Unlike most other bovine NAS studies (39, 108, 113), sea, one of the classical enterotoxins, was not detected in any of our isolates. The absence of sea in NAS and in bovine NAS specifically has been reported previously (28, 119).

All isolates of each NAS species contained on average 30 or more VF genes, with the highest virulence potential (defined by total number of VFs), assigned to SAG, SHY, and SCH (clade B), largely due to exotoxins, host evasion and capsular genes, followed by SPA and SEP (clade D) and SAR and SXY (clade E). Virulence potential was lowest for SFL (clade A) and species of clade E (SAC, SSU, and SNE), which all contained <21 VF genes. However, when associations between the total number of genes and disease severity (SCC category or presence of CM) were assessed, there was no significant association of the total number of genes with disease severity (Table 1). Notwithstanding, the mere presence of genes does not guarantee their expression (120). Additionally, many other factors (e.g., host environment, nutritional status, presence of other competing microbes, and host genetics) have crucial roles in successful colonization, persistence, and pathogenicity of mammary pathogens. For example, conversion of nonpathogenic bacteria into invasive pathogens in immunocompromised hosts has been reported (13, 121). Pathogenesis is complex and often involves an organized and systematic participation of various VFs to establish disease. Often VFs complement each other to promote pathogen colonization and persistence of disease (13).

NAS isolates clustered according to their distinct species when analyzed by t-SNE and in dendrograms created by Ward clustering and complete-linkage clustering methods, indicating that isolates from each NAS species are distinct and may represent distinct pathogens. Thus, effects of NAS on udder health should not be assessed as a single group; rather, they should be characterized according to individual species. Intermixing of various NAS species observed in t-SNE plot corresponded with phylogenetic relationships of these species. For example, an intermixed cluster of SAG and SHY and of SVI, SFL, and SSC reflects their close evolutionary relationship, as evident in phylogenetic trees (32).

We also computed the difference in gene associations among NAS species and for isolates from low, medium, and high SCC and CM. Differences in associations for individual NAS species and isolates from various inflammatory responses suggest complex interplay among virulence genes in causing disease. Unraveling these interactions will be important to elucidate distinctive pathogenic mechanisms of individual NAS species and assessing species-specific effects on udder health. However, caution should be exercised in predicting synergistic functions of genes, solely based on genetic studies, as expression of one gene could have antagonistic effects on other genes (13), as demonstrated for luk-PV and lukED (122), where expression of luk-PV inhibited expression of lukE and lukD genes.

Although we completed the largest VF screening of all known NAS species isolated from dairy cow’s milk, our study also had some limitations. Identification of VFs based on genetic similarity could be problematic for two reasons. First, we assumed a high level of sequence conservation signified conserved function. Identification and prediction of functions of NAS VFs were extrapolated from well-characterized analogues in SAU or from NAS of human origin. However, SAU VF genes may have niche-adapted functions, impacting their roles in virulence. Second, similarity-based identification between two genes may change over time. For example, if another gene with higher similarity was established as a separate VF gene, it would alter frequency distributions in our results. Furthermore, the presence/absence of virulence-associated genes may not directly correlate with pathogenesis and severity of disease, since in addition to expression of these genes, host environment and host genetics also have major roles in disease development and progression (13). Additionally, apart from the presence or absence of VF genes, mutations in regulators, such as two-component systems, are often involved in increased virulence (13), which were not investigated in this study.

MATERIALS AND METHODS

Selection of isolates.

A total of 441 NAS isolates were selected for whole-genome sequencing (WGS) from a collection of 5,507 isolates as described previously (32). The milk samples were not collected specifically for this study, rather the isolates were later selected from an already existing collection. Originally, milk samples were collected in 2007 to 2008 as described previously (123), and they were stored at the Mastitis Pathogen Collection of the Canadian Bovine Mastitis and Milk Quality Research Network (CBMQRN). The SCC was measured using the Fossomatic method (Fossomatic 4000 series; Foss Electric, Hillerød, Denmark) within 5 days after sample collection (123). Milk samples were collected from 87 herds distributed over four geographic regions encompassing Atlantic Canada (Nova Scotia, Prince Edward Island, New Brunswick), Québec, Ontario, and Western Canada (Alberta) (123). Isolates included 68 NAS isolates from CM cases (defined as abnormal milk or severe clinical signs), 26 isolates with a multidrug-resistant (MDR) phenotype (124), 1 isolate per cow of uncommon species (defined as species isolated from fewer than 20 cows), and 1 randomly selected isolate per cow for all other species (Table 2).

TABLE 2

NAS isolates grouped according to species and inflammation groups

NAS species	No. of isolatesa	% of isolates	No. of isolates in the following group:
			Low SCC	Medium SCC	High SCC	Clinical mastitis
S. agnetis	13	2.9	2	1	8	2
S. arlettae	15 (1)	3.4	9	1	4	1
S. auricularis	2	0.5	2
S. capitis	22	5.0	11		10	1
S. caprae	1	0.2		1
S. chromogenes	83 (13)	18.8	33	7	20	23
S. cohnii	24 (2)	5.4	16		7	1
S. devriesei	8	1.8	4	2	1	1
S. epidermidis	27 (1)	6.1	14	2	7	4
S. equorum	17	3.9	15	1	1
S. fleuretti	2	0.5	2
S. gallinarum	21	4.8	12	1	7	1
S. haemolyticus	28	6.3	12	2	10	4
S. hominis	11	2.5	7	1	3
S. hyicus	3	0.7	1		2
S. kloosii	1	0.2	1
S. nepalensis	2	0.5	2
S. pasteuri	6	1.4	2	1	3
S. saprophyticus	16	3.6	9	1	6
S. sciuri	29	6.6	13	2	6	8
S. simulans	42 (3)	9.5	16	2	8	16
S. succinus	15	3.4	11	3	1
S. vitulinus	6	1.4	4		2
S. warneri	19 (1)	4.3	10	5	4
S. xylosus	28 (5)	6.3	12	7	7	2
Total	441	100	220	40	117	64

a

Values in parentheses indicate the numbers of MDR isolates sequenced.

DNA extraction and whole-genome sequencing, assembly, and annotation.

Genomic DNA extraction, WGS, and genome assembly and annotation for 441 NAS isolates were performed as described previously (32, 125). Briefly, DNA was extracted with DNeasy Blood & Tissue kit (Qiagen, Toronto, ON, Canada). Quality and DNA concentrations were determined with a Qubit 2.0 fluorometer (Invitrogen, Burlington, ON, Canada). Sequencing libraries were prepared with an Illumina Nextera XT DNA library preparation kit (Illumina, San Diego, CA, USA). Paired-end sequencing (2 × 250 bp) was performed on the Illumina MiSeq platform (Illumina). An in-house pipeline, implemented in Snakemake workflow engine (126), was used to assemble and annotate genomes. Adapter sequences from the raw reads were removed using cutadapt (127), implemented in Trim Galore! 0.4.0. De novo assembly of these genomes was performed with Spades version 3.6.0 (128, 129). The quality of each assembled genome was assessed using Quality Assessment Tool for Genome Assemblies (QUAST) (130). Gene annotations for all genomes and gene predictions for contigs of >200 bp were done with Prokka 1.11 (131).

Identification of virulence factors in NAS genomes.

To ensure comprehensive analysis of all VFs reported for staphylococci, a comprehensive VF data set of staphylococci (CVFS) was created (see Data Set S1 in the supplemental material). For this data set, VF sequences were obtained from publicly available databases, including the VFDB database (132), the Victors database (http://www.phidias.us/victors/), the PATRIC database (133), and phenol-soluble modulin sequences from the UniProtKB database (134). To identify VFs in NAS genomes, local “blastdbs” of individual NAS species as well as a combined blast database of all 441 NAS genomes were created with the “makeblastdb” application from BLAST+ 2.5.0 (135). All VF sequences from CVFS were blasted against NAS genomes by conducting BLASTp searches, and a single best hit for each VF query from each NAS genome was selected as a final hit. Homology between query protein sequences and blast hits was determined by calculating H scores (124, 136). The H scores between protein sequences, labeled Ha (where a represents amino acid) were calculated using the following formula: Ha = VFid × Lm/Lq (124, 136), with VFid representing the level of BLASTp identities between the VF query sequence and identified protein sequence, articulated as a frequency ranging from 0 to 1, Lm representing the length of the matching sequence from the hit, and Lq denoting the length of the query sequence.

A cutoff 30% sequence similarity and 50% query length coverage were used for initial searches (124, 137–139). All genomic hits that met the minimum cutoff for each individual query were selected at this stage. A final blast hit table containing all possible hits for all query sequences from all NAS genomes was imported into R v.3.4.2 (140). Hits from each query sequence were then arranged according to Ha score, using “dplyr” version 0.7.2 in R (141). From this list, only hits with the highest Ha score (highest sequence similarity and query length coverage) were selected as potential VFs in NAS genomes and the rest were discarded. In order to prevent one VF query returning hits to two different genes within a particular genome, only the highest scoring hit selected based on Ha score was retained (124). After identification of putative NAS VFs, to confirm orthology between identified putative NAS VF and CVFS sequences, reciprocal blast searches between putative NAS VFs and the CVFS database were performed (142). Putative NAS VF sequences that failed to match with corresponding VFs from the CVFS database as the best hit in reciprocal blast searches were not considered true orthologues and were excluded from further analyses.

Classification and distribution of virulence factors.

A total of 191 VFs were classified into five functional categories: adherence (n = 28; Fig. 1A), exoenzymes (n = 21; Fig. 2A), host immune evasion (n = 20; Fig. 3A), iron uptake and metabolism (n = 29; Fig. 4A), and toxins (hemolysin, leukocidins, leukotoxins, toxic shock syndrome toxin, exfoliative toxins, type VII secretion system genes, phenol-soluble modulins, enterotoxins, and exotoxins) (n = 93; Fig. 5 and 6). Distributions of each of these VFs in all NAS isolates were determined by calculating the proportion of isolates in which they were identified.

For ease of reporting, 25 non-aureus staphylococcal species were abbreviated as follows: Staphylococcus agnetis = SAG, Staphylococcus arlettae = SAR, Staphylococcus auricularis = SAC, Staphylococcus capitis = SCI, Staphylococcus caprae = SCR, Staphylococcus chromogenes = SCH, Staphylococcus cohnii = SCO, Staphylococcus devriesei = SDE, Staphylococcus epidermidis = SEP, Staphylococcus equorum = SEQ, Staphylococcus fleurettii = SFL, Staphylococcus gallinarum = SGA, Staphylococcus haemolyticus = SHA, Staphylococcus hominis = SHO, Staphylococcus hyicus = SHY, Staphylococcus kloosii = SKL, Staphylococcus nepalensis = SNE, Staphylococcus pasteuri = SPA, Staphylococcus saprophyticus = SSA, Staphylococcus sciuri = SSC, Staphylococcus simulans = SSI, Staphylococcus succinus = SSU, Staphylococcus vitulinus = SVI, Staphylococcus warneri = SWA, and Staphylococcus xylosus = SXY.

In the figures (Fig. 1A, 2A, 3A, 4A, 5, and 6), these 25 NAS species are grouped into five distinct clades, according to phylogenetic placement (32).

Associations between virulence factors.

Plots were generated to determine associations among genes of five VF functional categories and between VF functional categories. Matrices of associations among VFs were computed using the phi coefficient (vcd package) in R v.3.4.2 (143). Before generating plots, VFs not identified in any NAS isolates were excluded from the data set. For visualization, resulting matrices were plotted using the “corrplot” package in R (144). Association plots were also generated, after classifying isolates into four SCC classes defined according to the inflammation severity score from 1 to 4, where 1 corresponds to low SCC (< 150,000 cells/ml), 2 is average SCC (150,000 – 250,000 cells/ml), 3 is high SCC (>250,000 cells/ml), and 4 corresponds to isolates from clinical cases of mastitis. In all association plots, positive associations are shown in blue, and negative associations are shown in red. Positive and negative associations were defined as the chance that if one gene was identified, then the other gene would be more or less likely identified in the same isolates, respectively.

Association between the presence of virulence factors and mastitis.

Relationships between four SCC classes and VFs were examined after dichotomizing genes by presence or absence in all isolates. Mixed-effects ordinal logistic regression models were used to compare the odds of having an increased inflammation severity score with an increase in the total number of VFs. Ordinal logistic regression models the log odds of having more severe inflammation compared to a reference category. The models had the form:

ln (oddss)=κs−βX−μiwhere odds_s is the odds of interest and represents the probability of score ≤ s divided by the probability of score > s, κ_s is the odds of having score = s divided by the odds of having score ≠ s, β is the regression coefficient associated with either the total number of virulence factors or the number of virulence factors within a given functional category, and μ_i is the random-effects term to account for within-herd correlations.

A mixed-effects linear regression model was used to determine the association between the total number of virulence factors or the number of virulence factors within a given functional category and the natural logarithm of the SCC of the sample. This model had the form:

ln (SCC)=βX+ε+μiwhere β is the regression coefficient associated with either the total number of virulence factors or the number of virulence factors within a given functional category, μ_i is the random-effects term to account for within-herd correlations, and ε is the residual error.

Significance was assessed using a cutoff P of <0.05, and the strength of the association and model fit were assessed using the log likelihood for the ordinal logistic regression models, and the adjusted multiple-R² for linear regression models. All statistical analyses were conducted using the “ordinal” package (145) in R v.3.4.2 (140). To determine whether VF distributions were associated with four SCC classes, agglomerative dendrograms were generated with the “AgglomerativeClustering” module, specifying four clusters (low SCC, medium SCC, high SCC, and CM) using either Ward clustering (based on analysis of within-cluster variances) or complete-linkage clustering (based on maximum within-cluster distances) methods. After dendrograms were generated, labels of SCC level were applied, and their distribution among four clusters was visually examined. Clustering was assessed using the “Scikit-Learn” module in Python 3.6 (146). A decision tree based on the presence or absence of VFs was also constructed using the “DecisionTreeClassifier” module in “Scikit-Learn” (146). A tree of depth 9 was used to determine whether specific combinations of VFs could be used to predict CM. Model generalizability was assessed by calculating the classification accuracy of the model for the training set and comparing it to the accuracy of the model for a previously unseen validation set. To visualize similarities in VF distributions between isolates, t-distributed stochastic neighbor embedding (t-SNE) (147) was conducted using the “manifold.t-SNE” module within “Scikit-Learn” (146). Distributions were reduced to two dimensions and plotted such that distances between points were proportional to relative differences between isolates. Plots were labeled with CM or NAS species and visually examined to find patterns in clusters in relation to either label.

Conclusions.

On the basis of whole-genome sequencing (WGS) data of 441 distinct NAS isolates, we established a comprehensive VF gene (n = 191) profile of 25 NAS species and determined associations among VF genes. The VF profiles of various NAS species and their VF gene associations differed, indicating their existence as distinct pathogens. Isolates from the same species usually had the same virulence potential, with variation within species being smaller than variation between species. The variation of VF genes among isolates of the same species may represent evolution and adaptation to distinct niches or environments within the host. We also investigated the link between VF genes and inflammatory response of the mammary gland but failed to detect a clear link between VF genes and mastitis. Regardless, comprehensive studies such as ours provide foundational genetic knowledge to design and conduct further experiments, focusing on understanding the synergy between VF and roles of individual NAS species in IMI and characterizing species-specific effects on udder health. To the best of our knowledge, this was the first study to determine the distributions of 191 virulence genes in 25 NAS species isolated from bovine IMI.

Data availability.

All whole-genome sequencing data used in this study are available (with no restrictions) from NCBI under BioProject ID PRJNA342349.