Skip to main content

Bovine milk somatic cell transcriptomic response to Staphylococcus aureus is dependent on strain genotype



Mastitis is an economically important disease of dairy cows with Staphylococcus aureus a major cause worldwide. Challenge of Holstein-Friesian cows demonstrated that S. aureus strain MOK124, which belongs to Clonal Complex (CC)151, caused clinical mastitis, while strain MOK023, belonging to CC97, caused mild or subclinical mastitis. The aim of this study was to elucidate the molecular mechanisms of the host immune response utilising a transcriptomic approach. Milk somatic cells were collected from cows infected with either S. aureus MOK023 or MOK124 at 0, 24, 48, 72 and 168 h post-infection (hpi) and analysed for differentially expressed (DE) genes in response to each strain.


In response to MOK023, 1278, 2278, 1986 and 1750 DE genes were found at 24, 48, 72 and 168 hpi, respectively, while 2293, 1979, 1428 and 1544 DE genes were found in response to MOK124 at those time points. Genes involved in milk production (CSN1, CSN10, CSN1S2, CSN2, a-LACTA and PRLR) were downregulated in response to both strains, with a more pronounced decrease in the MOK124 group. Immune response pathways such as NF-κB and TNF signalling were overrepresented in response to both strains at 24 hpi. These immune pathways continued to be overrepresented in the MOK023 group at 48 and 72 hpi, while the Hippo signalling, extracellular matrix interaction (ECM) and tight junction pathways were overrepresented in the MOK124 group between 48 and 168 hpi. Cellular composition analysis demonstrated that a neutrophil response was predominant in response to MOK124, while M1 macrophages were the main milk cell type post-infection in the MOK023 group.


A switch from immune response pathways to pathways involved in maintaining the integrity of the epithelial cell layer was observed in the MOK124 group from 48 hpi, which coincided with the occurrence of clinical signs in the infected animals. The higher proportion of M1 macrophages in the MOK023 group and lack of substantial neutrophil recruitment in response to MOK023 may indicate immune evasion by this strain. The results of this study highlight that the somatic cell transcriptomic response to S. aureus is dependent on the genotype of the infecting strain.

Peer Review reports


Bovine mastitis is a disease of the udder which affects milk yield and quality and causes significant economic losses to the dairy industry [1,2,3]. Mastitis can present as clinical or subclinical. Clinical mastitis presents with overt signs such as clots in milk or swelling of the quarter. Clots in milk are a result of the influx of substantial numbers of immune cells into the udder and are formed by aggregates of leukocytes and blood clotting factors [4]. Subclinical mastitis presents with no obvious signs but typically causes a moderate influx of immune cells into the udder, leading to an increase in somatic cell count (SCC) in milk [5]. Individual cow SCC of < 100,000 cells/ml is typically considered to be reflective of a healthy mammary gland, while an SCC > 200,000 cells/ml may be indicative of bacterial infection [6].

Staphylococcus aureus is a major cause of bovine mastitis in Ireland [7, 8] and worldwide [9,10,11,12] and is also one of the aetiological agents of human mastitis [13, 14]. S. aureus is a clonal pathogen which evolves slowly compared to many other bacteria [15]. Nevertheless, substantial genetic diversity exists among S. aureus strains, with many lineages identified to date, often adapted to different hosts [16]. Four main bovine-adapted lineages have been identified among Irish clinical mastitis isolates: Clonal Complex (CC)71, CC97, ST136 and CC151 [17]. CC97 and CC151 are also the most common mastitis-associated S. aureus lineages in Europe [18]. CC97 is a predominantly bovine-associated lineage although host jumps into humans with subsequent global spread have been demonstrated [16], while CC151 appears to be a uniquely bovine lineage which has not been found in other ruminants, avian species or in humans [16, 19,20,21].

During bovine IMI, pathogens enter the mammary gland through the teat canal. The immune response in the mammary gland is presumably initiated by bovine mammary epithelial cells (bMEC) in the alveoli of the udder and by resident macrophages [22]. Pathogen-Associated Molecular Patterns such as peptidoglycan and lipoteichoic acid on the surface of Gram-positive bacteria are recognised by Toll-like receptor 2 (TLR2), which stimulates activation of the NF-κB signalling pathway in bMEC as well as macrophages [23, 24]. NF-κB and other immune pathways stimulate the production of cytokines and chemokines, which cause inflammation in the udder and stimulate neutrophil recruitment [25]. Neutrophils are the first cell type recruited into the udder to engulf the bacteria and attempt to clear the infection [26]. Previous studies have reported that S. aureus IMI induces a weak inflammatory response with inhibition of the TLR2 and NF-κB pathways and low levels of pro-inflammatory cytokine production [27,28,29]. However, these studies used a variety of S. aureus strains and the magnitude of the immune response to S. aureus has been reported to be strain-dependent [30]. Therefore, the general supposition that S. aureus causes a low immune response in bovine IMI [31] may only be true for specific strains of S. aureus and there is a need to further study the immune response to different strains of S. aureus during IMI.

The transcriptomic response to intra-mammary challenge with S. aureus has been previously studied in bovine gland cistern tissue using microarrays [32, 33] and RNA sequencing [28]; however, studies were often short in duration [28, 32] or involved only one time point [33]. Characterising the transcriptomic response to S. aureus over several time points would be desirable, since it was previously demonstrated in vitro and in vivo that the response to S. aureus IMI changes over time [30, 34]. Furthermore, previous transcriptomic studies utilised only single strains of S. aureus for intra-mammary challenge, and none used a strain belonging to CC151, a globally distributed bovine-adapted lineage. Immune gene expression in bMEC in response to S. aureus differs according to bacterial strain or lineage [25, 34, 35]. Intramammary challenge of Holstein-Friesian cows with strains MOK023 (CC97) or MOK124 (CC151) demonstrated that S. aureus strain MOK023 caused mild or subclinical mastitis, while strain MOK124 caused severe clinical mastitis [30]. Differences were observed in secretion of cytokines IL-8 and IL-1β into milk between the MOK023 and the MOK124 infected cows, which further suggested that the immune response is strain-specific. Due to the differences observed between the immune response to MOK023 and MOK124 in bMEC in vitro and also in vivo, a further study to elucidate the molecular mechanisms of the host immune response in vivo by utilising a transcriptomic approach was performed. The objective of this study was to identify patterns of gene expression associated with the response to infection with two strains of S. aureus. Expression of cell markers was also examined in order to explore the proportions of immune cells in milk somatic cell samples in response to each strain.


Summary of in vivo disease presentation

Milk somatic cells collected during a previously described S. aureus IMI challenge trial were used. The influence of bacterial strain on clinical signs, bacterial load, SCC and the concentration of select cytokines and anti-S. aureus antibodies in milk were previously reported [30]. Briefly, cows infected with MOK124 generally developed clinical mastitis, while cows in the MOK023 group developed mild or subclinical mastitis. There was a greater drop in milk yield in the MOK124 group and SCC was significantly higher between 24 and 216 hpi in the MOK124 group compared to the MOK023 group. The MOK124 bacterial load was significantly higher than that of MOK023 at 24 and 31 hpi and subsequently decreased, while MOK023 bacterial load slowly increased and became significantly higher than MOK124 later in the infection (216 to 504 hpi). The levels of interleukin 8 (IL-8) and IL-1β were also significantly higher in milk of cows challenged with MOK124 compared to cows challenged with MOK023 during the first week of infection. Post-infection milk somatic cells were primarily neutrophils, which comprised 80–89% of somatic cells in the MOK124 group and 41–69% of somatic cells in the MOK023 group. Infection with MOK124 resulted in a significantly higher proportion of neutrophils in milk than infection with MOK023 [30].


The severity of bovine IMI has been previously associated with secretion of the bi-component leukocidin LukMF′ by S. aureus [36, 37]. As the genes encoding this toxin are found in MOK124, but not MOK023, LukM′ in milk from infected animals was quantified. No LukM′ was detected in either MOK023 or MOK124 group milk samples at any time point.

Sequence and mapping quality

Somatic cells collected pre-infection and at 24, 48, 72 and 168 hpi were subjected to 100 bp stranded paired-end sequencing. On average, 95% of the bases across all samples had a phred score of ≥30 and 98% had a score of ≥20. The mean number of reads per sample was 83.5 M; however, less than 60 M reads were generated from 3 samples from time 0 due to the low number of somatic cells in the samples. The distribution of read counts, phred scores and GC content per sample is available in Supplementary File 1. Quality control, as assessed by FastQC, yielded satisfactory results for all samples. On average, 93% of reads mapped uniquely to the bovine genome UMD 3.1.93.

Data exploration

Raw read counts were used for DE gene analysis using DESeq2. After filtering for low counts (sum of reads from all samples < 1) and non-protein coding genes, 17,334 genes were included in further analysis. The expression of the 500 most variable genes was used for principal component analysis (PCA) within the DESeq2 software package. Time 0 (pre-infection) samples clustered separately from infected samples (Fig. 1A). Expression patterns from cows infected with either strain did not separate on the PCA plot at 24 hpi. However, at 48 and 72 hpi, separation between the response to the 2 strains was evident. At 48 hpi MOK124 infected animals’ response had progressed to a state that was maintained until 168 hpi, while the response in MOK023 infected animals progressed more slowly to the same state and samples clustered in distinct groups at 48, 72 and 168 hpi. The difference in progression of gene expression in response to each strain could also be seen in the expression profile of a number of core immune genes across the experiment time course. This was illustrated by the NF-κB subunit 1 (NFKB1) gene, which is a key gene involved in switching on production of immune mediators through the NF-κB signalling pathway. A rapid increase in expression of this gene at 24 hpi was followed by a rapid decrease in MOK124 infected animals. On the other hand, infection with MOK023 caused a more sustained NFKB1 induction (Fig. 1B).

Fig. 1
figure 1

Results of exploration of raw counts data. A PCA plot generated from milk somatic cell gene expression data of the 500 most variable genes following infection with MOK023 (circles) or MOK124 (triangles). Samples from each time point are represented by a different colour: 0 hpi - orange, 24 hpi - purple, 48 hpi - green, 72 hpi - blue, 168 hpi - red. HPI - hours post infection. B Normalised expression (mean ± SEM) of the NF-κB subunit 1 (NFKB1) gene in milk somatic cells in response to infection with MOK023 (blue) or MOK124 (red)

Differentially expressed genes

To identify highly up and downregulated genes, a log2 fold change threshold of 1 (indicating a 2 fold increase or decrease in expression) was applied to the DE gene results in DESeq2 and genes with an adjusted P value < 0.05 were subsequently selected as significant. In response to MOK023, 1278, 2278, 1986 and 1750 significant DE genes were found at 24, 48, 72 and 168 hpi, respectively, while 2293, 1979, 1428 and 1544 significant DE genes were found in response to MOK124 at those time points (Fig. 2A). At 24 hpi the majority of genes with altered expression were downregulated in response to both strains, which could indicate a suppression of normal cell function; while at 72 and 168 hpi, upregulated DE genes predominated in response to both strains. Details of all significant DE genes are available in Supplementary File 2.

Fig. 2
figure 2

Differentially expressed genes in response to S. aureus MOK023 and MOK124. A The number of significant DE genes (Benjamini-Hochberg adjusted P-value < 0.05, log2 fold change threshold > 1) which were up- (red) and downregulated (blue) in milk somatic cells in response to strains MOK023 and MOK124 relative to pre-infection. B Top 5 most significant up- and downregulated DE genes in response to infection with MOK023. C Top 5 most significant up- and downregulated DE genes in response to infection with MOK124. Genes with lowest adjusted P value were selected at each time point. In the case of genes upregulated at 24 hpi in response to infection with MOK023 two genes had the same adjusted P value and so both are included

Among the most significant upregulated DE genes in response to MOK023, the gene for oncostatin M (OSM) was in the top 5 at all time points post-infection (Fig. 2B). OSM encodes a cytokine involved in inflammation and regulation of production of other cytokines. Oncostatin is also secreted by neutrophils during their rolling in blood vessels in order to enhance neutrophil adhesion to endothelial cells [38]. In response to MOK124, the top 5 upregulated genes included those involved in the immune response such as interleukin genes (IL17F, IL18) and interleukin-2 receptor (CD25) gene (Fig. 2C). The top 5 downregulated genes in response to MOK124 included genes encoding milk proteins lactalbumin alpha (a-LACTA), casein alpha subunit 2 (CSN1S2), casein alpha subunit 1 (CSN1) and casein kappa (CSN10). In order to examine the effect of strain on transcription of major genes involved in milk production, log2 fold changes of relevant significantly DE genes were examined (Table 1). The main genes involved in milk production were reviewed by Strucken et al. (2015) [39] and include leptin, lactalbumin, caseins, prolactin and several genes involved in fat transport, energy production and growth. No significant change in expression of the selected milk production genes was observed in MOK023 infected animals at 24 hpi; however, all examined genes with the exception of the gene coding for 1-acylglycerol-3-phosphate O-acyltransferase 6 (AGPAT6) were significantly downregulated at all other time points in this group. In the MOK124 group, downregulation of the expression of milk production genes was already evident by 24 hpi and continued until 168 hpi. Log2 fold changes in expression of genes encoding milk-related proteins (a-LACTA, CSN1, CSN10, CSN1S2, CSN2) were 2–3 log lower in the MOK124 group than in the MOK023 group, indicating a more severe reduction of protein synthesis in the first group.

Table 1 Expression of genes involved in milk production

A total of 402 genes in the MOK023 group were DE at all time points post-infection, while 362 DE genes were common in the response to MOK124 at all time points (Fig. 3 A-B). Of these, 130 were shared between MOK023 and MOK124 at all time points and are proposed as strain-independent biomarkers of S. aureus IMI during the first week of infection (Fig. 3C). Among the common genes, 42 were involved in immune response processes.

Fig. 3
figure 3

Common DE genes in milk somatic cells across time points. A Common DE genes in response to MOK023. B Common DE genes in response to MOK124. C Overlap of the DE genes common at all time points between the two strains. HPI - hours post infection

Pathway analysis

KEGG pathways

Pathway analysis of KEGG metabolic processes revealed 67 unique significantly overrepresented pathways in response to MOK023 and 74 significantly overrepresented pathways in response to MOK124. The pathways were enriched at one or more time points. Overall, 84 unique pathways were enriched in response to either strain. The top 10 most significant KEGG pathways at each time point are presented in Fig. 4, while a complete list of pathways at all time points is in Supplementary File 3.

Fig. 4
figure 4

KEGG pathways in response to infection with S. aureus MOK023 or MOK124 over time. Significant (P < 0.05) and unique (filtered for shared genes) pathways are shown. Colour scale represents P value, grey squares represent non-significance. Up to 10 of the most significant pathways at each time point are presented. Hpi - hours post infection

KEGG pathway analysis illustrated that there was a sustained somatic cell immune response to MOK023. At 24 hpi, the top 2 most significant pathways in response to both MOK023 and MOK124 were NF-κB and TNF signalling pathways, both of which are involved in the immune response (Fig. 4). Genes involved in these two pathways displayed similar patterns of expression change in response to both strains at 24 hpi (Fig. 5). Immune response pathways continued to be significantly upregulated in the MOK023 group at 48 and 72 hpi, but not in the MOK124 group, even though mastitis was more severe in the MOK124 group (Fig. 5). Instead, a specific response to MOK124 was evident at later time points, characterised by overrepresentation of the Hippo signalling, extracellular matrix receptor interaction (ECM) and tight junction pathways. The Hippo signalling pathway is involved in regulating cell apoptosis or proliferation [40], while ECM receptor and tight junction pathways both represent groups of genes that influence communication between epithelial cells and maintaining the integrity of the epithelial cell layer [41, 42]. When examining genes involved in the Hippo signalling pathway, the majority of genes were downregulated relative to time 0 in response to MOK124 (Fig. 6A).

Fig. 5
figure 5

Log2fold change of genes involved in selected KEGG pathways in response to MOK023 and MOK124. A NF-κB signalling pathway. B TNF signalling pathway. Grey tiles show genes that were not significantly DE

Fig. 6
figure 6

Log2fold change of genes involved in selected KEGG pathways in response to MOK023 and MOK124. A Hippo signalling pathway. B ECM receptor interaction pathway. C Tight junction pathway. Grey tiles show genes that were not significantly DE

The ECM receptor interaction pathway was enriched in the MOK124 group at 48 hpi. Among ECM receptor interaction genes, collagen genes COL4A1, COL4A2, COL6A1, COL6A2, COL8A2, COL12A1, COL14A1 and COL18A1 were downregulated in response to MOK124 between 48 and 168 hpi (Fig. 6B). The collagen genes encode proteins that are components of basement membranes. Several integrin genes were upregulated in response to MOK124 at 48 hpi (ITGA5, ITGA9, ITGA2B, ITGB7), while some were also downregulated between 48 and 168 hpi (ITGA2, ITGA3, ITGB4, ITGB6).

For the tight junction pathway, in response to MOK124, most paracellular space genes were downregulated, including genes encoding claudin (CCLN), occludin (OCLN), Marvel domain-containing protein 3 (MARVELD3) and immunoglobulin superfamily member (IGSF5) (Fig. 6C). Several genes from this pathway were upregulated at 24 hpi; however, the pathway was not significantly enriched in the MOK124 group at that time point. Only two genes were upregulated at later time points: Rho associated coiled-coil containing protein kinase 1 gene ROCK1 was upregulated at 48 hpi and myosin heavy chain 1 gene MYH1 was upregulated at 72 hpi, both of which are involved in regulation of actin cytoskeleton. All other genes, involved in cell polarity (PARD6B, TJP3, PATJ, AMOT, PRKCZ, LLG2, TJP1, DLG3), decreased paracellular permeability (TJP1, CGN, CGNL1, CLDN3), reduced cell proliferation (YBX3, CCND1), cell differentiation (YBX3, ERBB2) and tight junction assembly (RAB13, TJP1, CGN) were downregulated between 48 and 168 hpi.

While pathways specific to MOK124 emerged from 48 hpi, the tuberculosis pathway was solely overrepresented in response to MOK023, at 24, 48 and 72 hpi. Notably, among the upregulated genes in this pathway were CTSD, ATP6A1 and LAMP1, which encode cathepsin D, ATPase H+ transporting accessory protein 1 and lysosomal associated membrane protein 1, respectively. All of these proteins are involved in internalisation into macrophages and arrest of phagosome maturation. Other pathways specific to MOK023 were inflammation related pathways, such as inflammatory bowel disease (IBD) pathway and rheumatoid arthritis pathway.

Gene ontology

Enrichment analysis of GO biological processes revealed 36 different significantly overrepresented terms in response to MOK023 and 36 significantly overrepresented terms in response to MOK124. Overall, 66 different GO terms were enriched in response to either strain. The top 10 enriched terms in response to each strain at each time point are presented in Fig. 7, while a complete list at each time point is attached in Supplementary File 3.

Fig. 7
figure 7

GO biological processes in response to infection with S. aureus MOK023 or MOK124 over time. Top 10 most significant (P < 0.05) and unique (filtered for shared genes) are shown. Colour scale represents P value, grey squares represent non-significance. Hpi - hours post infection

Immune system processes were overrepresented in the MOK023 group at 24 and 48 hpi, while they were only significant in MOK124 group at 24 hpi. While the GO terms in the MOK023 group varied across time points, a small number were consistently overrepresented in the MOK124 group between 48 and 168 hpi. These included extracellular matrix organisation, wound healing and cell junction organisation. Additionally, cell death was an enriched process in both groups at 168 hpi.

Immune cell populations

A significantly higher proportion of neutrophils was observed in milk from cows in the MOK124 group than in the MOK023 group [30]. In order to assess whether immune cell proportions in the milk somatic cell samples could be inferred using transcriptomic data, a cell population analysis based on gene expression of cell markers was performed using quanTIseq [43]. The cell proportions at 0 hpi were similar between the two groups, with most cells classified as “other”, likely representing mammary epithelial cells (Fig. 8A). Genes for epithelial cell markers such as cytokeratins 14, 18 and 19 (KRT14, KRT18, KRT19), E-cadherin (CDH1) and epithelial cell adhesion molecule (EPCAM) [44, 45] showed the highest expression at 0 hpi and decreased post-infection (Fig. 8B). Immune cell populations predominated at later time points and immune cell proportions differed between the two groups (Fig. 8A). At 24 hpi, neutrophils represented the largest proportion of milk somatic cells in both MOK023 and MOK124 infected animals; however, in the MOK124 group neutrophils represented a greater proportion of somatic cells than in the MOK023 group. A major difference in cell proportions was observed between 48 and 168 hpi. In the MOK023 group, the largest proportion of cells between 48 and 168 hpi were M1 macrophages, whereas in the MOK124 group neutrophils were the predominant cell type. A slightly larger proportion of M2 macrophages could also be observed at each time point post-infection in the MOK124 group when compared to the MOK023 group. M2 macrophages seemed to increase in proportion in both groups over time.

Fig. 8
figure 8

Cell composition analysis. A Changes in cellular composition in milk somatic cells over time illustrated by cell proportions. B Expression of epithelial cell marker genes in milk somatic cells over time in response to infection with MOK023 (blue) or MOK124 (red). KRT14 - cytokeratin 14, KRT18 - cytokeratin 18, KRT19 - cytokeratin 19, CDH1 - E-cadherin, EPCAM - epithelial cell adhesion molecule, HPI - hours post infection

Milk secretion of chemokines and cytokines

Cytokine concentration was measured in milk at 0, 7, 24, 48, 72, 168, 336 and 504 hpi in order to reflect both the acute and chronic response to infection. There were no differences in the concentration of any of the examined cytokines or chemokines between the MOK023 and MOK124 groups prior to intramammary challenge or after 168 hpi. Protein concentrations increased post-infection in both groups with a peak expression at 168 hpi, except for TNFα and IL-8 which demonstrated a sharp increase at 24 hpi. A significant group x time interaction (P < 0.05) was observed for IL-1α, IL-1β, IL-1RA, and IL-6, which are typical pro-inflammatory cytokines, but also IL-2, IL-4, IL-10 and CCL3 (Fig. 9). In the case of all the above cytokines, MOK124 caused significantly higher protein secretion at 168 hpi, which is also a time point at which the difference in milk yield was low.

Fig. 9
figure 9

Concentration of selected immune markers in milk. The immune markers are divided into: pro-inflammatory cytokines: IL-1α, IL-1β, IL-1RA, IL-6, TNFα; cytokines involved in adaptive immunity: IFNγ, IL-2, IL-4, IL-10, IL-17A; and chemokines: CCL2, CCL3, CCL4, CXCL10, IL-8. Colour scale represents mean fluorescence intensity (MFI). Significant time points where cytokine concentration in MOK124 group was higher than in the MOK023 group are marked with *


In this study, the transcriptional response of milk somatic cells to two genotypically distinct strains, representing two important globally distributed bovine-adapted lineages of S. aureus, was assessed. Differences in the response to the 2 strains were found in gene expression patterns and in overrepresented pathways. An initial immune response to both strains was evident but the transcriptional immune response was sustained in the MOK023 group between 24 and 72 hpi, while immune genes were only significantly upregulated in the MOK124 group at 24 hpi. This may indicate a difference in disease progression or a response to different virulence strategies employed by the two strains. The expression of genetic markers of various immune cells also differed in response to each strain, indicating differing immune cell proportions among somatic cells recruited in response to each strain, a result that was confirmed by a differential milk somatic cell count using microscopy [30]. More neutrophils were recruited in response to MOK124 from 48 hpi and there were proportionally more M1 macrophages in the MOK023 group.

Exploration of normalised gene expression data yielded an interesting insight into the dynamics of the transcriptional response to the two infecting strains, with a strain-specific transcriptional response evident by 48 hpi. At this time point, gene expression in the MOK023 group was similar to that at 24 hpi while in the MOK124 group gene expression was more similar to that at 72 and 168 hpi. Examination of normalised counts of NFKB1 and other key genes suggested a difference in the pro-inflammatory response, with a more rapid and higher response to MOK124, which later subsided; in contrast, more sustained upregulation of immune gene expression was evident in response to MOK023. The sustained transcriptional upregulation of immune genes in the MOK023 group was supported by pathway analysis – the overrepresentation of both the NF-κB and TNF signalling pathways was sustained for longer after infection with MOK023. This difference in the transcriptional response may be related to a rapid switch from activation of the immune response to a disruption in cellular proliferation and the integrity of the epithelial cell layer in the MOK124 group. The 168 h time point could represent the end of the acute innate immune response in both groups, since immune pathways were not significantly enriched at this time point and few genes involved in the NF-κB and TNF signalling pathways were upregulated in either group.

The transcriptional response observed in this study supported several findings of the corresponding in vivo infection, such as the larger reduction in milk yield, higher inflammatory response at 24 hpi and higher SCC in the MOK124 group, as well as the increase in the proportion of neutrophils in milk observed in the MOK124 group. The downregulation of expression of milk production genes in the MOK124 group coincided with a significant reduction in milk yield, which was observed in response to MOK124 between 24 and 144 hpi when compared to MOK023 [30]. This downregulation also indicates that infection with MOK124 caused a greater disruption to the epithelial layer of the mammary gland, possibly due to epithelial cell dysfunction and sloughing, when compared to MOK023 infection. Overall, the cellular composition analysis indicated that the strain-specific response was mediated by different cell populations. M1 macrophages were mainly involved in the response to MOK023, while the response to MOK124 led to a higher neutrophil recruitment. There was no significant difference between the MOK023 and MOK124 groups in the proportion of large mononuclear cells by microscopy; however, large mononuclear cells would have been comprised of M1 and M2 macrophages as well as monocytes, and shifts in those respective populations could be observed when looking at expression of cellular markers. The transcriptomic data therefore has the potential to provide more detail on the proportions of different cell populations in milk, especially in species where antibodies to cell markers are not available for the fine characterization of these cell subpopulations.

A larger proportion of M2 macrophages in somatic cells from the MOK124 group may indicate a trend towards tissue repair. M2 macrophages synthesise polyamines and proline which stimulate cell growth, collagen formation and tissue repair [46]. M2 macrophages are also capable of producing ECM components [47] and may have contributed to the ECM related pathways which were overrepresented in response to MOK124. Therefore, it is possible that MOK124 creates high epithelial dysfunction, which exposes the stroma to bacteria and bacterial products and induces a strong inflammation that the host attempts to repair. This leads to a higher recruitment of neutrophils and more plasma leakage due to failure of the epithelium to function as a barrier.

Pathway analysis revealed that MOK124 elicited a strain-specific response indicative of tissue injury, as illustrated by the downregulation of Hippo signalling, ECM receptor interaction and tight junction pathways between 48 and 168 hpi. Moreover, MOK124-specific GO terms identified at those time points, such as wound healing and ECM organisation, corresponded to the KEGG pathways identified. These pathways are well studied in humans. The key genes in the Hippo signalling pathway are YAP1 and WWTR which encode transcription factors YAP and TAZ. A successful activation of this pathway depends on translocation of the YAP/TAZ proteins to the nucleus, followed by binding of the transcription factors to a TEAD protein. The YAP/TAZ-TEAD complex causes activation of transcriptional responses which result in cell proliferation by controlling the expression of genes directly involved in the cell cycle and by upregulation of anti-apoptotic genes [40]. The genes involved in the Hippo signalling pathway were downregulated, which could indicate a decrease in cell proliferation. However, protein interactions and localisation in the cell would need to be studied in order to confirm this hypothesis.

The ECM receptor interaction pathway is involved in maintenance of cell shape and the interaction of epithelial cells with the basement membrane. Within the ECM receptor interaction pathway, several integrin genes were upregulated in response to MOK124 at 48 hpi (ITGA5, ITGA9, ITGB7, ITGA2B). The integrin proteins mediate cell to matrix and cell to cell adhesion [48]. ITGA5 encodes the light and heavy chains that comprise the α5 subunit. This subunit associates with the β1 subunit to form a fibronectin and fibrinogen receptor. Both fibronectin and fibrinogen can mediate adhesion of S. aureus to host cells [49], and increase of fibronectin levels in whey of Holstein-Friesian cows is associated with mastitis progression [50]. The integrin α9 subunit encoded by the ITGA9 gene, when bound to the β1 chain, forms an integrin that is a receptor for vascular cell adhesion molecule 1 (VCAM1), cytotactin and osteopontin, all of which are involved in immune response. ITGB7 encodes integrin subunit β7 which forms dimers with an α4 chain or an αE chain and plays a role in leukocyte adhesion [48]. These two upregulated integrin genes could indicate leukocyte adhesion and migration was taking place in MOK124 infected animals at 48 hpi. Integrins can also act with cytoskeleton proteins from the tight junction pathway and clustering of integrins increases cellular tension, which can lead to the activation of the Hippo signalling pathway [48].

Tight junctions form barriers between cells to the flux of ions and molecules. Pore and leak pathways can be differentiated within the action of the tight junction, and while claudins are involved in the pore pathway, zonula occludens 1 (ZO-1) and occludins are involved in the leak pathway [51]. In the MOK124 group, claudin (CCLN), occludin (OCLN), Marvel (MARVELD3) and immunoglobulin superfamily member (IGSF5) genes were downregulated as part of the tight junction pathway, which suggests a decrease in cell to cell interaction and cell detachment. TNFα has been shown to induce barrier loss in epithelial cells during Crohn’s disease in vitro and in mice. The loss of integrity of intestinal epithelium in those studies was mediated by TNF-induced occludin internalisation [51]. Therefore, upregulation of the TNF pathway in the MOK124 group at 24 hpi could have contributed to loss of cell integrity at later time points. On the other hand, dysregulation of cytoskeleton genes in the tight junction pathway could be indicative of S. aureus internalisation into host cells. Modification of the actin cytoskeleton in cells infected by S. aureus has previously been described [28] and was attributed to host cell internalisation. A CC151 strain, RF122, has been hypothesised to be adapted to an intracellular niche [19]. In vitro, strain MOK023 was a significantly higher internaliser in bMEC than strain MOK124 [34]. The cytoskeleton reorganisation in response to MOK124 is therefore more likely to be related to cell damage, loss of cellular structure and an attempt to repair the tissue rather than due to host cell internalisation. The presence of clinical signs in MOK124 challenged animals [30] also confirms that cell damage likely took place and a healing response would be necessary to rebuild the damaged tissue. The ability of MOK124 to cause tissue damage may be associated with the toxin encoding genes present in this lineage. MOK124 encodes enterotoxin genes seg, sei, selm, seln, seln, selo, egc, selu and ORFCM14 as well as the leukocidins lukF-PV and lukM, none of which are present in MOK023 [17]. Staphylococcal enterotoxin genes seg, sei, selm and sei as well as agrII (which is the agr type of MOK124) have recently been found to be associated with a higher likelihood of causing clinical mastitis in a collection of European S. aureus strains isolated from natural infections [18]. While LukM′ was not detected in milk of MOK124 infected cows in this study, the other toxins encoded by this strain could play a role in the virulence of this strain. Other CC151 strains have been found to produce low amounts of LukM′ in previous studies [52], therefore the mechanism of neutrophil killing by CC151 strains remains to be determined. Since immune genes tended to be upregulated for longer in response to MOK023 than in response to MOK124, the greater influx of neutrophils in response to MOK124 is unlikely entirely due to immune signalling by milk somatic cells. Instead, a larger influx of neutrophils could be related to cellular and tissue damage.

Although MOK023 caused immune gene induction both in bMEC in vitro and in milk somatic cells in vivo [30, 34], immune signalling from bMEC was followed by a temporary recruitment of neutrophils at 24 hpi, and M1 macrophages became the predominant cell type from 48 hpi. Since macrophages are normally present in a healthy mammary gland, while neutrophils tend to be recruited to a site of intramammary infection [53], the difference in cellular composition in milk post-infection in response to the two strains may indicate that MOK023 subverted the immune response which failed to recruit neutrophils to a significant extent. This is also supported by the fact that even though immune gene expression was sustained in response to MOK023, bacterial load increased in the MOK023 infected animals, to reach its maximum at 168 hpi. Lack of inflammation and relatively low neutrophil numbers in response to MOK023 could be due to host cell internalisation. It was demonstrated that MOK023 can internalise into bMEC in vitro [34]. It has also been reported that human-associated S. aureus can internalise and survive within human [54] and bovine macrophages [55]. S. aureus SH1000 (CC8) internalised within murine or human macrophages can exhibit a persister phenotype and evade antimicrobial treatment [56], and strain ATCC 29213 (CC5 [57];) was found to promote its survival in bovine macrophages by blocking autophagic flux [55]. Internalisation of MOK023 into MEC and/or macrophages could aid in the evasion of immune response. Interestingly, among significantly DE genes, LAMP1 was found to be overexpressed solely in response to MOK023 (Supplementary File 2). This gene has previously been associated with S. aureus survival inside a mature phagosome in macrophage cell lines [58, 59]. Upregulation of the tuberculosis pathway in the MOK023 group further supports the hypothesis that MOK023 may internalise into host cell macrophages and be immune evasive. The tuberculosis pathway represents a process where tubercle bacilli enter macrophages in the lung, and following host cell internalisation they interfere with phagosomal maturation, antigen presentation, apoptosis and host bactericidal responses to establish persistent or latent infection [60]. Persistence of infection of MOK023 was illustrated by a consistent bacterial load in infected animals over 30 days; however, more research is needed to determine the precise mechanisms by which this strain persists within the host. Notably, a recent study associated high internalisation into bMEC, biofilm formation and low cytotoxicity with high within-herd prevalence of a S. aureus CC9 strain [61].

The transcriptomic response to S. aureus mastitis in dairy cows has been studied previously with samples taken at single time points: 3 hpi [28, 32], 24 hpi [33] and 48 hpi [62]. Udder biopsies were used in those studies to analyse gene expression in response to a S. aureus IMI, and single S. aureus strains were used. This study characterises later events occurring in IMI, and represents a response of different cell populations, as it would be expected that an udder biopsy would be mostly comprised of bMEC. Transcriptomic analysis of S. aureus mastitis in naturally infected dairy cows has also been performed on post-mortem mammary tissue and peripheral blood leukocytes [63, 64], with the extent of time post-infection being unknown, and the infecting strains not identified. This is the first transcriptomics study of S. aureus mediated IMI where the response to two different strains of S. aureus are compared, with full strain genotype information available. Genotyping of bacterial strains used, use of multiple strains of S. aureus as well as examination of the cellular populations of studied samples should be considered in future studies. All of the above can ensure better understanding of strain-specific responses and an improved ability to compare different studies of IMI.

The concentration of selected cytokines and chemokines in milk throughout the course of infection confirmed and extended the finding that MOK124 induced higher IL-1β and IL-8 concentrations in milk. A discrepancy was found between the gene expression results and the cytokine profiles since cytokine and chemokine protein levels remained elevated in response to MOK124 at later time points, while most of the cytokine genes were not DE post 24 h in this group. While correlation between transcriptomic and proteomic data is often poor [65,66,67] due to protein turnover, RNA degradation or post-translational modifications, there could be other reasons for the discrepancy in this study. Firstly, the protein levels were measured in milk, while the transcript levels were measured in milk somatic cells. The mammary epithelial cells lining the bovine mammary gland could be contributing to cytokine secretion into milk, and would not be captured by the RNA sequencing. Therefore, mammary epithelial cells and other cells of the epithelial lining such as resident macrophages and lymphocytes could be major contributors to cytokine and chemokine levels in the milk. Secondly, the number and the type of milk cells likely influenced absolute protein concentration with higher neutrophil counts in MOK124 infected cows. Finally, a systemic response could have contributed to increased milk cytokine levels as exudation of plasma was likely in the MOK124 group with higher milk Ig concentrations, but not evident in the MOK023 group [30]. The milk cytokine data highlights the need to implement proteomic studies alongside transcriptomic studies to better understand the molecular mechanisms of host immune response. A study where each milk somatic cell type was assessed separately during an in vivo trial, along with udder biopsies or post-mortem examinations at one pre-defined time point, would be beneficial to further understanding of bovine immune response to S. aureus. Furthermore, the pattern of cytokine expression measured in vitro on MEC culture infected with either strain [34] does not correlate with the secretion profile detected in milk upon mammary gland infection of dairy cows, showing that bacteria-MEC interactions do not summarize the complex mechanisms occurring during mastitis and may not be predictive of the pathogenicity of S. aureus strains.


This study confirms the ability of bovine S. aureus strains to cause a strain-specific immune response during IMI, with each strain resulting in a unique signature of infection in terms of the number and type of immune cells recruited, the biological function of the recruited cells and the resulting cytokine and chemokine concentration in milk. MOK124 caused more inflammation and udder damage, as illustrated by the disruption of tight junction and Hippo signalling pathways between 48 and 168 hpi, coupled with clinical signs observed in the infected animals. On the other hand, MOK023 caused less inflammation and may be able to limit neutrophil recruitment and hence persist in the mammary gland, as suggested by the sustained immune gene expression but low increase in SCC, sustained bacterial load and mild clinical signs in MOK023 infected cows. Since differences in early transcriptional response between the two strains were observed, followed by either a severe or a mild, persistent infection, it is possible that the magnitude of this early response may be crucial in determination of the course of staphylococcal IMI. The response to each strain also involved recruitment of specific immune cells, with M1 macrophages primarily recruited in response to MOK023, while neutrophils were the main immune cell type recruited in response to MOK124. This indicates that different cell types as well as the transcriptional activity of each cell type may influence gene expression differences between groups when analysing mixed cell population tissues such as blood or milk. Further studies are necessary to further understand the molecular mechanisms utilised by these strains to cause inflammation or evade the immune response.


Study design

The aim of this study was to elucidate the molecular mechanisms of the host immune response in vivo by utilising a transcriptomic approach. Fourteen first lactation Holstein-Friesian cows were purchased from two commercial farms which were leptospirosis, salmonellosis and Johne’s disease free. Animals were selected based on composite milk SCC < 100,000 cells/ml. Milk bacteriology was performed prior to challenge to confirm lack of contaminating organisms or mastitis pathogens. Cows were infected with approximately 500 cfu of S. aureus strain MOK023 (ST3170, CC97) or MOK124 (ST151, CC151) in 400 μl of PBS as described previously [30]. Animals in which a S. aureus infection was not established (two cows challenged with MOK023 and one cow challenged with MOK124) were excluded from the study. Therefore, the MOK023 infected group included 5 animals, while the MOK124 infected group included 6 animals. Due to severe clinical mastitis developed by 2 cows in the MOK124 group, these animals were treated with antibiotic and milk samples for transcriptomic analysis were not taken after antibiotic treatment: cow 1776 was removed from the trial at 31 hpi, while cow 1920 was removed at 120 hpi. A further sample from the MOK124 group at 48 hpi was removed from the analysis for quality control purposes as a sample from an uninfected quarter was mistakenly sequenced instead of the sample from the infected quarter. The number of animals in each experimental group at each time point is presented in Table 2.

Table 2 The numbers of animals (N) included in the study at each time point

RNA extraction from milk somatic cells

Milk from the infected quarter was collected into a sterile universal tube (Sarstedt, Nümbrecht, Germany). In order to obtain a sufficient quantity of RNA, 100 ml of milk was used for somatic cell isolation. Milk was centrifuged at 800 x g, at 15 °C for 10 min. The fat layer and supernatant were removed and tubes were placed upside down on tissue paper to drain for 10 min. Cell pellets were subsequently resuspended in 900 μl of Qiazol (Qiagen, Manchester, UK) and vortexed for 20 s to lyse the cells. Lysed samples were stored in Qiazol at − 20 °C until RNA extraction. Samples were rapidly thawed at 37 °C and maintained at room temperature (RT) for the duration of the extraction. RNA was extracted using an RNeasy Plus Universal kit (Qiagen, Manchester, UK) as per the manufacturer’s instructions. RNA was eluted in 30 μl of RNase-free water (Sigma Aldrich Ireland Ltd., Wicklow, Ireland).

RNA quality was assessed using 2100 Bioanalyser RNA Nano chips (Agilent Technologies Ireland Limited, Cork, Ireland) as well as with a NanoDrop™ 1000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). RNA quantity was determined with a Qubit RNA broad range assay (Invitrogen: Bio-Sciences, Dun Laoghaire, Ireland). Sample RIN values ranged between 5.9 and 9.5 with a mean RIN value of 7.82 (Supplementary Table S1), and concentrations ranged from 4.1–930 ng/μl (mean concentration of 154 ng/μl).

Library preparation and sequencing

RNA was sent on dry ice to Macrogen (South Korea) for library preparation and sequencing. RNA quality was confirmed on arrival using a TapeStation HSRNA Screen Tape and analysed on a 2200 TapeStation (Agilent Technologies Korea Ltd., Seoul, Korea). Library preparation was performed with 100 ng of RNA using a TruSeq Stranded mRNA LT Sample Prep Kit (Illumina, San Diego, USA) according to the manufacturer’s instructions (Part #15031047 Rev. E). Library size was verified with a 2100 Bioanalyzer using a DNA 1000 chip (Agilent, Korea) and library quantity was determined with qPCR according to the Illumina qPCR Quantification Protocol Guide. The mean library fragment size was 348 bp while the mean library concentration was 92.63 nM. Libraries were diluted and pooled with 300 pM of each library and sequenced (100 bp, paired-end) on a NovaSeq (Illumina, USA). Samples were divided across two sequencing runs.

Quality control, assembly and DE gene analysis

All the bioinformatics pipeline bash, Perl and R scripts used for computational analyses were deposited in a GitHub repository at RStudio v1.1.463 and R v3.5.2 were used for the analyses. Sequence quality was assessed using FastQC v11.5 [68]. Trimming was performed using fastp v12.1 [69] with default settings as well as enabling base correction for paired end data and enabling overrepresented sequence analysis; on average 0.5% of the reads were trimmed per sample. Trimmed fastq files were re-checked with FastQC to confirm sequence read quality. Trimmed sequences were mapped to the bovine genome (UMD3.1) using the STAR aligner v2.5.2 [70], with gene counts generated simultaneously with the STAR software using the UMD3.1 v93 of the B. taurus reference genome and annotation. Genes with zero read counts across all samples as well as non-protein coding genes were removed in a filtering step [71]. DeSeq2 v1.18.1 [72] was then used to visualise gene expression data and perform differentially expressed (DE) gene analysis. Data was normalised with a vst function and visualised with a principal component analysis (PCA) plot of the samples based on the 500 genes with the highest inter-sample variation. DE gene lists were generated using a negative binomial generalized linear model. Wald tests were performed to obtain DE genes separately for each group at each time point using time 0 as a reference, with P values adjusted for multiple comparisons using a Benjamini and Hochberg (B-H) method. DE genes with an adjusted P value < 0.05 and a log2 fold change threshold of 1 were used for further DE gene data exploration and pathway analysis.

Pathway and gene correlation analysis

Significant (P adj < 0.05) and highly overexpressed (|log2 fold change| > 1) DE genes were converted to their 1 to 1 human orthologs using BioMart, ordered by adjusted P value and analysed for enrichment of Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways using gProfiler R package v0.6.6 [73]. When a group of DE genes matched to multiple pathways, the pathway that was represented by the highest number of DE genes was selected in a filtering step.

Cell composition analysis

Transcripts per million (TPMs) were generated from raw gene counts. Bovine Gene IDs converted to 1 to 1 human orthologs were used for the analysis. To estimate the proportion of the different immune cell populations in each somatic cell sample, a quanTIseq algorithm which is based on expression of cell specific surface marker genes was used [43].

Quantification of a panel of cytokines and chemokines in milk

Milk was defatted by centrifugation at 800 x g, 15 °C for 10 min. The top fat layer was discarded and the skim milk supernatant was stored in 1 ml aliquots at − 20 °C until use. Prior to protein quantification milk was diluted 1:2 in assay buffer (Merck-Millipore, Burlington, USA). Concentrations for 15 cytokines (IL-1α, IL-1β, IL-1RA, IL-2, IL-4, IL-6, IL-17a, IFNγ, CCL-2, CCL-3, CCL-4, CXCL8, CXCL10, and TNFα) were determined using a custom bovine cytokines MilliPlex xMAP assay (Merck Millipore, USA). Data were recorded on a MAGPIX flow cytometer using Xponent software (Luminex, Austin, USA).


ELISA to detect the LukM subunit of the LukMF′ Staphylococcal toxin in milk samples was performed as described previously [37]. Briefly, anti-LukM polyclonal bovine IgG (32 μg/ml) was used to coat the wells of an ELISA plate at 4 °C overnight. The plate was emptied and washed twice with 300 μl of PBS-T (DBPS: Gibco, UK; 0.05% Tween-20: Sigma-Aldrich, St. Louis, MO), blocked with 200 μl of 4% skimmed milk powder (Marvel, DSDelta Ltd) in PBS-T for 30 min with shaking and washed once with 300 μl PBS-T. A 100 μl aliquot of each sample (skimmed milk heated to 95 °C for 10 min to inactivate any existing antibodies against LukM) or standard (recombinant LukM; kindly provided by Utrecht University) was added to the wells and incubated for 60 min at room temperature (RT) with shaking, followed by 2 washes with 300 μl PBS-T. Mouse anti-LukM monoclonal antibody (KoMa43, Podiceps, The Netherlands) was diluted to 3 μg/ml in 1% skimmed milk-PBS-T and 50 μl was added to each well. The plate was incubated for 60 min at RT, with shaking, and washed 3 times. Polyclonal HRP conjugated goat anti-mouse IgG (BioLegend, London, UK) was diluted 1 in 1000 and 100 μl was added to each well. The plate was incubated for 60 min at RT, with shaking, and washed 3 times. TMB (100 μl; Thermo Fisher Scientific, Waltham, USA) was subsequently added and incubated in the dark at RT for 15 min, followed by addition of 100 μl of the stop solution (0.05 M H2SO4; Thermo Fisher Scientific, Waltham, USA). Absorbance was read at 450 nm using a microplate reader (xMark, Bio-Rad, Watford, UK).

Statistical analysis

Data generated from cytokine and LukM′ protein quantification was analysed with a repeated measures ANOVA using a MIXED® procedure in SAS (SAS Version 9.4, SAS Institute, Cary NC, USA) to compare strain-dependent differences. Strain, time and the strain x time interaction were included in the model, with time as a repeated measure. Time 0 measurements for each cytokine were included as a covariate, and all variable interactions were explored. Log10 transformations were used for all cytokine variables to satisfy the distributional requirements of ANOVA. Correlations among repeated measurements across time between strains were modelled using the compound symmetry covariance structure for each parameter analysed. Tukey’s post-hoc multiple comparison test was applied as appropriate.

Availability of data and materials

Raw fastq files used during the current study are available in the European Nucleotide Archive (ENA) repository with a project number PRJEB43443 .

Code used to analyse the data is available at:

Any other data from the current study are available from the corresponding author on reasonable request.



bovine mammary epithelial cells


clonal complex


differentially expressed


extracellular matrix


hours post-infection


Gene Ontology


intramammary infection


Kyoto Encyclopedia of Genes and Genomes


principal component analysis


room temperature


sequence type. Gene symbol abbreviations – Supplementary Table S2.


  1. Geary U, Lopez-Villalobos N, O'Brien B, Garrick DJ, Shalloo L. Estimating the impact of somatic cell count on the value of milk utilising parameters obtained from the published literature. J Dairy Res. 2014;81(2):223–32.

    Article  CAS  PubMed  Google Scholar 

  2. Halasa T, Huijps K, Osteras O, Hogeveen H. Economic effects of bovine mastitis and mastitis management: a review. Vet Q. 2007;29(1):18–31.

    Article  CAS  PubMed  Google Scholar 

  3. Rollin E, Dhuyvetter KC, Overton MW. The cost of clinical mastitis in the first 30 days of lactation: an economic modeling tool. Preventive Veterinary Medicine. 2015;122(3):257–64.

    Article  CAS  PubMed  Google Scholar 

  4. Harmon RJJods. Physiology of mastitis and factors affecting somatic cell counts. 1994;77(7):2103–12.

  5. Bradley AJ. Bovine Mastitis: An Evolving Disease. Vet J. 2002;164(2):116–28.

    Article  CAS  PubMed  Google Scholar 

  6. Bradley A, MJIp G. Use and interpretation of somatic cell count data in dairy cows. 2005;27(6):310–5.

  7. Barrett DJ, Healy AM, Leonard FC, Doherty ML. Prevalence of pathogens causing subclinical mastitis in 15 dairy herds in the Republic of Ireland. Ir Vet J. 2005;58(6):333–7.

    Article  PubMed  PubMed Central  Google Scholar 

  8. Keane OM, Budd KE, Flynn J, McCoy F. Pathogen profile of clinical mastitis in Irish milk-recording herds reveals a complex aetiology. Veterinary Record. 2013;173(17).

  9. Barkema H, Schukken Y, Lam T, Beiboer M, Wilmink H, Benedictus G, et al. Incidence of clinical mastitis in dairy herds grouped in three categories by bulk milk somatic cell counts. 1998;81(2):411–9.

  10. Østeras O, Sølverød L, Reksen O. Milk culture results in a large Norwegian survey—effects of season, parity, days in Milk, resistance, and clustering. J Dairy Sci. 2006;89(3):1010–23.

    Article  PubMed  Google Scholar 

  11. Persson Y, Nyman AK, Gronlund-Andersson U. Etiology and antimicrobial susceptibility of udder pathogens from cases of subclinical mastitis in dairy cows in Sweden. Acta Vet Scand. 2011;53(1):36.

    Article  PubMed  PubMed Central  Google Scholar 

  12. Gianneechini R, Concha C, Rivero R, Delucci I, JMJAVS L. Occurrence of Clinical and Sub-Clinical Mastitis in Dairy Herds in the West Littoral Region in Uruguay. 2002;43(4):221.

  13. Delgado S, García P, Fernández L, Jiménez E, Rodríguez-Baños M, del Campo R, et al. Characterization of Staphylococcus aureus strains involved in human and bovine mastitis. FEMS Immunology & Medical Microbiology. 2011;62(2):225–35.

    Article  CAS  Google Scholar 

  14. Holmes MA, Zadoks RN. Methicillin resistant S. aureus in human and bovine mastitis. J Mammary Gland Biol Neoplasia. 2011;16(4):373–82.

    Article  PubMed  Google Scholar 

  15. Lindsay JA, Holden MT. Understanding the rise of the superbug: investigation of the evolution and genomic variation of Staphylococcus aureus. Funct Integr Genomics. 2006;6(3):186–201.

    Article  CAS  PubMed  Google Scholar 

  16. Fitzgerald JR. Livestock-associated Staphylococcus aureus: origin, evolution and public health threat. Trends Microbiol. 2012;20(4):192–8.

    Article  CAS  PubMed  Google Scholar 

  17. Budd KE, McCoy F, Monecke S, Cormican P, Mitchell J, Keane OM. Extensive genomic diversity among bovine-adapted Staphylococcus aureus: evidence for a genomic rearrangement within CC97. PLoS One. 2015;10(8):e0134592.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Hoekstra J, Zomer AL, Rutten VPMG, Benedictus L, Stegeman A, Spaninks MP, et al. Genomic analysis of European bovine Staphylococcus aureus from clinical versus subclinical mastitis. Sci Rep. 2020;10(1):18172.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Herron-Olson L, Fitzgerald JR, Musser JM, Kapur V. Molecular correlates of host specialization in Staphylococcus aureus. PLoS One. 2007;2(10):e1120.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Merz A, Stephan R, Johler S. Staphylococcus aureus Isolates from Goat and Sheep Milk Seem to Be Closely Related and Differ from Isolates Detected from Bovine Milk. Front Microbiol. 2016;7(319).

  21. Schlotter K, Ehricht R, Hotzel H, Monecke S, Pfeffer M, Donat K. Leukocidin genes lukF-P83 and lukM are associated with Staphylococcus aureus clonal complexes 151, 479 and 133 isolated from bovine udder infections in Thuringia, Germany. Vet Res. 2012;43(1):42.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Le Gall A, Plommet M, editors. Observations sur la croissance des staphylocoques et la réaction leucocytaire au cours des premières heures de la mammite expérimentale de la brebis. Annales de Biologie Animale Biochimie Biophysique; 1965: EDP Sciences.

  23. Yang W, Zerbe H, Petzl W, Brunner RM, Günther J, Draing C, et al. Bovine TLR2 and TLR4 properly transduce signals from Staphylococcus aureus and E. coli, but S. aureus fails to both activate NF-κB in mammary epithelial cells and to quickly induce TNFα and interleukin-8 (CXCL8) expression in the udder. Mol Immunol. 2008;45(5):1385–97.

    Article  CAS  PubMed  Google Scholar 

  24. Gilbert FB, Cunha P, Jensen K, Glass EJ, Foucras G, Robert-Granié C, et al. Differential response of bovine mammary epithelial cells to Staphylococcus aureus or Escherichia coli agonists of the innate immune system. Vet Res. 2013;44(1):1–23.

    Article  CAS  Google Scholar 

  25. Lahouassa H, Moussay E, Rainard P, Riollet C. Differential cytokine and chemokine responses of bovine mammary epithelial cells to Staphylococcus aureus and Escherichia coli. Cytokine. 2007;38(1):12–21.

    Article  CAS  PubMed  Google Scholar 

  26. Oviedo-Boyso J, Valdez-Alarcon JJ, Cajero-Juarez M, Ochoa-Zarzosa A, Lopez-Meza JE, Bravo-Patino A, et al. Innate immune response of bovine mammary gland to pathogenic bacteria responsible for mastitis. J Inf Secur. 2007;54(4):399–409.

    Article  Google Scholar 

  27. Bannerman DD, Paape MJ, Lee JW, Zhao X, Hope JC, Rainard P. Escherichia coli and Staphylococcus aureus elicit differential innate immune responses following intramammary infection. Clin Diagn Lab Immunol. 2004;11(3):463–72.

    Article  PubMed  PubMed Central  Google Scholar 

  28. Günther J, Petzl W, Bauer I, Ponsuksili S, Zerbe H, Schuberth HJ, et al. Differentiating Staphylococcus aureus from Escherichia coli mastitis: S. aureus triggers unbalanced immune-dampening and host cell invasion immediately after udder infection. Sci Rep. 2017;7(1):4811.

    Article  PubMed  PubMed Central  Google Scholar 

  29. Riollet C, Rainard P, Poutrel B. Differential induction of complement fragment C5a and inflammatory cytokines during Intramammary infections with Escherichia coli and Staphylococcus aureus. Clin Diagn Lab Immunol. 2000;7(2):161–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Niedziela DA, Murphy MP, Grant J, Keane OM, Leonard FC. Clinical presentation and immune characteristics in first-lactation Holstein-Friesian cows following intramammary infection with genotypically distinct Staphylococcus aureus strains. J Dairy Sci. 2020;103(9):8453–66.

    Article  CAS  PubMed  Google Scholar 

  31. Barkema HW, Schukken YH, Zadoks RN. Invited review: the role of cow, pathogen, and treatment regimen in the therapeutic success of bovine Staphylococcus aureus mastitis. J Dairy Sci. 2006;89(6):1877–95.

    Article  CAS  PubMed  Google Scholar 

  32. Petzl W, Günther J, Mühlbauer K, Seyfert H-M, Schuberth H-J, Hussen J, et al. Early transcriptional events in the udder and teat after intra-mammary Escherichia coli and Staphylococcus aureus challenge. Innate Immun. 2016;22(4):294–304.

    Article  CAS  PubMed  Google Scholar 

  33. Fang L, Hou Y, An J, Li B, Song M, Wang X, et al. Genome-Wide Transcriptional and Post-transcriptional Regulation of Innate Immune and Defense Responses of Bovine Mammary Gland to Staphylococcus aureus. Frontiers in Cellular and Infection Microbiology. 2016;6(193).

  34. Murphy MP, Niedziela DA, Leonard FC, Keane OM. The in vitro host cell immune response to bovine-adapted Staphylococcus aureus varies according to bacterial lineage. Sci Rep. 2019;9(1):6134.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Zbinden C, Stephan R, Johler S, Borel N, Bu J, Bruckmaier RM, et al. The Inflammatory Response of Primary Bovine Mammary Epithelial Cells to Staphylococcus aureus Strains Is Linked to the Bacterial Phenotype. PLoS ONE. 2014;9(1):e87374.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Hoekstra J, Rutten V, Sommeling L, van Werven T, Spaninks M, Duim B, et al. High production of LukMF' in Staphylococcus aureus field strains is associated with clinical bovine mastitis. Toxins (Basel). 2018;10(5):200.

    Article  CAS  Google Scholar 

  37. Vrieling M, Boerhout EM, van Wigcheren GF, Koymans KJ, Mols-Vorstermans TG, de Haas CJC, et al. LukMF′ is the major secreted leukocidin of bovine Staphylococcus aureus and is produced in vivo during bovine mastitis. Sci Rep. 2016;6(1):37759.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Setiadi H, Yago T, Liu Z, McEver RP. Endothelial signaling by neutrophil-released oncostatin M enhances P-selectin-dependent inflammation and thrombosis. Blood Adv. 2019;3(2):168–83.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Strucken E, Laurenson Y, Brockmann G. Go With the Flow - Biology and Genetics of the Lactation Cycle2015.

  40. Totaro A, Panciera T, Piccolo S. YAP/TAZ upstream signals and downstream responses. Nat Cell Biol. 2018;20(8):888–99.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Anderson JM, Van Itallie CM. Physiology and function of the tight junction. Cold Spring Harb Perspect Biol. 2009;1(2):a002584.

    Article  PubMed  PubMed Central  Google Scholar 

  42. Dedhar S. Cell–substrate interactions and signaling through ILK. Curr Opin Cell Biol. 2000;12(2):250–6.

    Article  CAS  PubMed  Google Scholar 

  43. Finotello F, Mayer C, Plattner C, Laschober G, Rieder D, Hackl H, et al. Molecular and pharmacological modulators of the tumor immune contexture revealed by deconvolution of RNA-seq data. Genome Medicine. 2019;11(1):34.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Polisetti N, Agarwal P, Khan I, Kondaiah P, Sangwan VS, Vemuganti GK. Gene expression profile of epithelial cells and mesenchymal cells derived from limbal explant culture. Mol Vis. 2010;16:1227–40.

    CAS  PubMed  PubMed Central  Google Scholar 

  45. Arévalo Turrubiarte M, Perruchot M-H, Finot L, Mayeur F, Dessauge F. Phenotypic and functional characterization of two bovine mammary epithelial cell lines in 2D and 3D models. Am J Phys Cell Phys. 2015;310(5):C348–C56.

    Article  Google Scholar 

  46. Martinez FO, Sica A, Mantovani A, Locati M. Macrophage activation and polarization. Front Biosci. 2008;13(1):453–61.

    Article  CAS  PubMed  Google Scholar 

  47. Novak ML, Koh TJ. Macrophage phenotypes during tissue repair. J Leukoc Biol. 2013;93(6):875–81.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Ivaska J, Heino J. Adhesion receptors and cell invasion: mechanisms of integrin-guided degradation of extracellular matrix. Cellular and Molecular Life Sciences CMLS. 2000;57(1):16–24.

    Article  CAS  PubMed  Google Scholar 

  49. Hauck CR, Ohlsen K. Sticky connections: extracellular matrix protein recognition and integrin-mediated cellular invasion by Staphylococcus aureus. Curr Opin Microbiol. 2006;9(1):5–11.

    Article  CAS  PubMed  Google Scholar 

  50. Maity S, Das D, Ambatipudi K. Quantitative alterations in bovine milk proteome from healthy, subclinical and clinical mastitis during S. aureus infection. J Proteome. 2020;223:103815.

    Article  CAS  Google Scholar 

  51. Shen L, Weber CR, Raleigh DR, Yu D, Turner JR. Tight junction pore and leak pathways: a dynamic duo. Annu Rev Physiol. 2011;73(1):283–309.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. Vrieling M, Koymans KJ, Heesterbeek DAC, Aerts PC, Rutten VPMG, de Haas CJC, et al. Bovine Staphylococcus aureus Secretes the Leukocidin LukMF′ To Kill Migrating Neutrophils through CCR1. mBio. 2015;6(3):e00335–15.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  53. Leitner G, Shoshani E, Krifucks O, Chaffer M, Saran A. Milk leucocyte population patterns in bovine udder infection of different Aetiology. J Veterinary Med Ser B. 2000;47(8):581–9.

    Article  CAS  Google Scholar 

  54. Kubica M, May R, Potempa J, Koziel J, Guzik K, Zarebski M, et al. A Potential New Pathway for Staphylococcus aureus Dissemination: The Silent Survival of S. aureus Phagocytosed by Human Monocyte-Derived Macrophages. PLoS ONE. 2008;3(1):e1409.

    Article  PubMed  PubMed Central  Google Scholar 

  55. Cai J, Li J, Zhou Y, Wang J, Li J, Cui L, et al. Staphylococcus aureus facilitates its survival in bovine macrophages by blocking autophagic flux. J Cell Mol Med. 2020;24(6):3460–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  56. Peyrusson F, Varet H, Nguyen TK, Legendre R, Sismeiro O, Coppée J-Y, et al. Intracellular Staphylococcus aureus persisters upon antibiotic exposure. Nat Commun. 2020;11(1):2200.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  57. Suzuki M, Matsumoto M, Takahashi M, Hayakawa Y, Minagawa H. Identification of the clonal complexes of Staphylococcus aureus strains by determination of the conservation patterns of small genomic islets. J Appl Microbiol. 2009;107(4):1367–74.

    Article  CAS  PubMed  Google Scholar 

  58. Lacoma A, Cano V, Moranta D, Regueiro V, Domínguez-Villanueva D, Laabei M, et al. Investigating intracellular persistence of Staphylococcus aureus within a murine alveolar macrophage cell line. Virulence. 2017;8(8):1761–75.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  59. Tranchemontagne ZR, Camire RB, O'Donnell VJ, Baugh J, Burkholder KM. Staphylococcus aureus strain USA300 perturbs Acquisition of Lysosomal Enzymes and Requires Phagosomal Acidification for survival inside macrophages. Infect Immun. 2015;84(1):241–53.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  60. Koul A, Herget T, Klebl B, Ullrich A. Interplay between mycobacteria and host signalling pathways. Nat Rev Microbiol. 2004;2(3):189–202.

    Article  CAS  PubMed  Google Scholar 

  61. Grunert T, Stessl B, Wolf F, Sordelli DO, Buzzola FR, Ehling-Schulz M. Distinct phenotypic traits of Staphylococcus aureus are associated with persistent, contagious bovine intramammary infections. Sci Rep. 2018;8(1):1–10.

    Article  CAS  Google Scholar 

  62. Wang X, Fan Y, He Y, Han Z, Gong Z, Peng Y, et al. Integrative analysis of miRNA and mRNA expression profiles in mammary glands of Holstein cows artificially infected with Staphylococcus aureus. Pathogens. 2021;10(5):506.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  63. Wang XG, Ju ZH, Hou MH, Jiang Q, Yang CH, Zhang Y, et al. Deciphering transcriptome and complex alternative splicing transcripts in mammary gland tissues from cows naturally infected with Staphylococcus aureus mastitis. PLoS One. 2016;11(7):e0159719.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  64. Wang D, Liu L, Augustino SMA, Duan T, Hall TJ, MacHugh DE, et al. Identification of novel molecular markers of mastitis caused by Staphylococcus aureus using gene expression profiling in two consecutive generations of Chinese Holstein dairy cattle. J Anim Sci Biotechnol. 2020;11(1):98.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  65. Hack CJ. Integrated transcriptome and proteome data: the challenges ahead. Briefings in Functional Genomics. 2004;3(3):212–9.

    Article  CAS  Google Scholar 

  66. Haider S, Pal R. Integrated analysis of transcriptomic and proteomic data. Curr Genomics. 2013;14(2):91–110.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  67. Nie L, Wu G, Culley DE, Scholten JCM, Zhang W. Integrative analysis of transcriptomic and proteomic data: challenges, solutions and applications. Crit Rev Biotechnol. 2007;27(2):63–75.

    Article  CAS  PubMed  Google Scholar 

  68. Andrews S. FastQC: A quality control tool for high throughput sequence data. 2018 [Available from:

  69. Chen S, Zhou Y, Chen Y, Gu J. Fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics. 2018;34(17):i884–i90.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  70. Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29(1):15–21.

    Article  CAS  PubMed  Google Scholar 

  71. Smedley D, Haider S, Ballester B, Holland R, London D, Thorisson G, et al. BioMart–biological queries made easy. BMC Genomics. 2009;10(1):22.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  72. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  73. Reimand J, Arak T, Adler P, Kolberg L, Reisberg S, Peterson H, et al. G: profiler—a web server for functional interpretation of gene lists (2016 update). Nucleic Acids Res. 2016;44(W1):W83–W9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references


We would like to acknowledge the farm staff in UCD Lyons Research Farm, particularly Eddie Jordan, Michael Clarke and Joseph Callanan for invaluable assistance with sample collection for this study. Thank you to Dr. Matthew McCabe (Teagasc) and Dr. Elaine Kenny (ELDA Biotech) for advice on RNA extraction. The LukM′ ELISA was performed by Margaret Murray in Teagasc Grange. Antibodies and protocol for this ELISA were generously provided by Manouk Vrieling, Kok van Kessel and Lindert Benedictus (Utrecht University).


This study was funded by a grant from the Department of Agriculture, Food and the Marine of Ireland (14/S/802).

Author information

Authors and Affiliations



DAN optimised and conducted the laboratory-based work, performed bioinformatics analysis of the RNAseq data and drafted the manuscript. OMK conceived the initial idea of the study and participated in its design and coordination. FCL conceived the initial idea of the study and participated in its design and coordination. PC provided mentoring and supervision during the analysis of RNAseq data and performed the cell composition analysis. GF conducted/supervised protein laboratory-based work and statistical analysis of protein data. All authors helped to draft the manuscript and read and approved the final manuscript.

Corresponding author

Correspondence to Orla M. Keane.

Ethics declarations

Ethics approval and consent to participate

All procedures involving animals were performed in accordance with the relevant guidelines and regulations of the European Directive 2010/63/EU and received approval from the University College Dublin Animal Research Ethics Committee (protocol number Leonard AREC 16 44), and a project authorisation from the Health Product Regulatory Authority (HPRA) (project authorisation number AE18982/P108). The study was carried out in compliance with the ARRIVE guidelines and has been described in detail previously (Niedziela DA, et al. 2020. J Dairy Sci. 103 [9]:8453–66).

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1. Supplementary Table S1.

Sample metadata and parameters of sequencing and alignment quality.

Additional file 2. Supplementary Table S2.

Log2 fold changes and adjusted P values of all differentially expressed genes (Benjamini‐Hochberg adjusted P‐value < 0.05, log2 fold change threshold > 1) in response to MOK023 and MOK124 at each time point.

Additional file 3.

 All significantly overrepresented gene ontology (Table S3) and KEGG (Table S4) pathways derived from DE gene lists at every time point post-infection in response to MOK023 and MOK124.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Niedziela, D.A., Cormican, P., Foucras, G. et al. Bovine milk somatic cell transcriptomic response to Staphylococcus aureus is dependent on strain genotype. BMC Genomics 22, 796 (2021).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: