In China, bovine mastitis is a major health problem in dairy herds, causes substantial economic losses and remains one of the main indications for antimicrobial treatment (He et al., Reference He, Ma, Lei, He, Li, Tao, Wang, Song, Wang, Wang, Shen, Cai and Wu2020; Fan et al., Reference Fan, Qadeer, Asiri, Alzahrani, Alzahrani, Alsharif, Khan and Jiang2025). Surveillance and farm-based investigations across multiple provinces consistently identify Staphylococcus aureus as a leading mastitis pathogen in China (Wang et al., Reference Wang, Cha, Liu, Deng, Yang, Xu, Wang, Zhang, Gu, Huang and Qu2022b; Saeed et al., Reference Saeed, Kamaruzzaman, Gahamanyi, Nguyen, Hossain and Kahwa2024). A meta-analysis of studies conducted between 2000 and 2020 estimated that approximately one-third of bovine mastitis cases involved S. aureus, although estimates varied widely between studies (Wang et al., Reference Wang, Cha, Liu, Deng, Yang, Xu, Wang, Zhang, Gu, Huang and Qu2022b). The emergence of multidrug-resistant lineages, including occasional livestock-associated methicillin-resistant S. aureus (LA-MRSA) in dairy settings, further heightens concerns about udder health and the spread of resistance determinants along the dairy food chain (Bal et al., Reference Bal, Coombs, Holden, Lindsay, Nimmo, Tattevin and Skov2016; Cui et al., Reference Cui, Li, Ali, Kalim, Wang, Song, Li, Ren, Ma, Zou, Shen and Xu2020; Zhang et al., Reference Zhang, Wang, Wang, Zhang, Shang, Li, Song, Li, Cheng, Zhang, Zhao, Shen and Cui2023).
Genomic epidemiology provides a high-resolution framework for describing circulating lineages and for linking phenotypic resistance with underlying resistance determinants and mobile genetic elements (Quainoo et al., Reference Quainoo, Coolen, Van Hijum, Huynen, Melchers, Van Schaik and Wertheim2017; Partridge et al., Reference Partridge, Kwong, Firth and Jensen2018). Whole-genome sequencing (WGS), coupled with MLST and spa typing and core-genome SNP phylogenetics, improves comparability between studies and strengthens inference about transmission, persistence and the spread of resistance-associated elements (Schürch et al., Reference Schürch, Arredondo-Alonso, Willems and Goering2018; Naushad et al., Reference Naushad, Nobrega, Naqvi, Barkema and De Buck2020).
The overall objective of this study was to characterise the clonal structure and antimicrobial resistance profiles of S. aureus recovered from bulk-tank milk in Nanping, a subtropical dairy region in south-eastern China. Dairy production in Nanping has expanded since the 1950s, and the area is now Fujian's principal milk-sourcing region, with large-scale Chinese Holstein herds supplying regional processors under structured herd health and mastitis control programmes. Despite this development, genomic data on dairy-associated S. aureus from this subtropical region remain limited. Specifically, this study aimed to: (i) describe the distribution of major MLST/spa lineages and their phylogenetic relationships; (ii) summarise antimicrobial non-susceptibility patterns and associated resistance determinants; and (iii) explore an unsupervised multivariate grouping framework to support risk-tier stewardship for mastitis control.
Materials and methods
Sample collection
We collected 162 bulk-tank raw milk samples (≥40 mL) from all 14 large-scale Chinese Holstein dairy herds in in Nanping, Fujian Province, between March 2023 and November 2024. Each herd was sampled at approximately 2-month intervals. Together, these herds housed approximately 20,000 lactating cows and were milked in rotary parlours. This was an exploratory pilot study based on a census of all eligible large-scale herds in the study area; no additional exclusion criteria were applied. Bulk-tank milk was aseptically collected from the tank outlet immediately after thorough agitation by trained laboratory staff using sterile 50 mL tubes. Samples were kept at 4°C on farm and transported to the laboratory in insulated coolers within 4 hour of collection. Bulk-tank somatic cell counts (SCC) were routinely monitored by the dairy processor for quality assurance under the Chinese raw-milk standard; however, detailed SCC records were not available for inclusion in the present analysis.
Isolation and phenotypic characterisation
For S. aureus isolation, 25 mL of each sample was enriched in 225 mL of 7.5% NaCl broth at 36°C for 24 h and streaked onto Baird–Parker agar and 5% sheep blood agar. Colonies with typical S. aureus morphology were subcultured and confirmed as Gram-positive cocci in clusters, with a positive tube coagulase (rabbit plasma) test, and identified as S. aureus using the VITEK 2 Compact system (China NMPA & NHFPC Reference China2017). Antimicrobial susceptibility testing was performed by disc diffusion according to CLSI M100 (35th edition), using 14 antimicrobial agents (listed in Fig. 2). Inhibition zone diameters were interpreted as susceptible (S), intermediate (I) or resistant (R) using CLSI breakpoints (CLSI, 2025). Multidrug resistance (MDR) was defined as non-susceptibility to at least one agent in three or more antimicrobial classes (Magiorakos et al., Reference Magiorakos, Srinivasan, Carey, Carmeli, Falagas, Giske, Harbarth, Hindler, Kahlmeter, Olsson-Liljequist, Paterson, Rice, Stelling, Struelens, Vatopoulos, Weber and Monnet2012). CLSI-recommended quality-control strains were included in each batch.
DNA extraction, whole-genome sequencing and assembly
Genomic DNA was extracted from overnight cultures, quantified and used to prepare paired-end (2 × 150 bp) libraries for sequencing on an Illumina NovaSeq 6000 platform. Raw reads were demultiplexed with ‘bcl2fastq’, quality assessed with ‘FastQC’, trimmed/filtered to ensure that >90% of bases had quality scores ≥ Q30 and mean sequencing depth exceeded 100 × . Reads were de novo assembled with ‘SPAdes’ v3.15.5. All genome assemblies have been deposited in NCBI under BioProject PRJNA1258060.
In silico genomic analyses, phylogenomics, and functional-risk stratification
Typing and annotation: Multilocus sequence typing (MLST) was performed using the PubMLST scheme (Enright et al., Reference Enright, Day, Davies, Peacock and Spratt2000; Jolley et al., Reference Jolley, Bray and Maiden2018), and spa typing was assigned using the Ridom SpaServer (Friedrich et al., Reference Friedrich, Witte, Harmsen, De Lencastre, Hryniewicz, Scheres and Westh2006). Antimicrobial resistance genes were predicted with RGI v6.0.0 against CARD v4.0.0 using Perfect and Strict criteria (minimum identity ≥95%, high-quality/coverage) (McArthur et al., Reference McArthur, Waglechner, Nizam, Yan, Azad, Baylay, Bhullar, Canova, De Pascale, Ejim, Kalan, King, Koteva, Morar, Mulvey, O'Brien, Pawlowski, Piddock, Spanogiannopoulos, Sutherland, Tang, Taylor, Thaker, Wang, Yan, Yu and Wright2013; Alcock et al., Reference Alcock, Huynh, Chalil, Smith, Raphenya, Wlodarski, Edalatmand, Petkau, Syed, Tsang, Baker, Dave, McCarthy, Mukiri, Nasir, Golbon, Imtiaz, Jiang, Kaur, Kwong, Liang, Niu, Shan, Yang, Gray, Hoad, Jia, Bhando, Carfrae, Farha, French, Gordzevich, Rachwalski, Tu, Bordeleau, Dooley, Griffiths, Zubyk, Brown, Maguire, Beiko, Hsiao, Brinkman, Van Domselaar and McArthur2023). Virulence factors were annotated with ‘ABRicate’ v1.0.1 against VFDB (core set; accessed April 2025) using identity ≥75% and coverage ≥80% (Liu et al., Reference Liu, Zheng, Jin, Chen and Yang2019, Reference Liu, Zheng, Zhou, Chen and Yang2022; Seemann, Reference Seemann2020). Plasmid replicons were identified with ‘ABRicate’ against PlasmidFinder (accessed April 2025) using an identity ≥75% (default coverage) (Camacho et al., Reference Camacho, Coulouris, Avagyan, Ma, Papadopoulos, Bealer and Madden2009; Carattoli et al., Reference Carattoli, Zankari, García-Fernández, Voldby Larsen, Lund, Villa, Møller Aarestrup and Hasman2014). SCCmec types were inferred using ‘ABRicate’ against SCCmecFinder (accessed April 2025; default settings) (Kaya et al., Reference Kaya, Hasman, Larsen, Stegger, Johannesen, Allesøe, Lemvigh, Aarestrup, Lund and Larsen2018). Immune-evasion cluster (IEC) types were inferred from the presence/absence patterns of scn, chp/sak and sea/sep (van Wamel et al., Reference van Wamel, Rooijakkers, Ruyken, van Kessel and van Strijp2006).
Phylogenomic reconstruction: Trimmed reads were mapped to S. aureus reference genome NCTC 8325 (RefSeq assembly GCF_000013425.1; NCBI, n.d.) using Snippy v4.6.0 with BWA-MEM (Seemann, Reference Seemann2015 [15 May 2014]; Shi et al., Reference Shi, Nayfach and Pollard2023; Yoshimura et al., Reference Yoshimura, Kajitani, Gotoh, Katahira, Okuno, Ogura, Hayashi and Itoh2019). Core-genome SNPs were extracted and used to construct a maximum-likelihood phylogenetic tree in IQ-TREE v2.1.2 (Nguyen et al., Reference Nguyen, Schmidt, Von Haeseler and Minh2015; Minh et al., Reference Minh, Schmidt, Chernomor, Schrempf, Woodhams, Von Haeseler and Lanfear2020) with 1,000 bootstrap replicates. Trees were visualized and annotated in iTOL v7 (Letunic and Bork, Reference Letunic and Bork2024).
Unsupervised functional-risk stratification: Phenotypic resistance was coded as S = 0, I = 0.5 and R = 1 and integrated with presence/absence data for resistance genes, virulence factors and plasmid replicons to generate a 59 × P feature matrix (after excluding invariant features). The remaining features were z-standardised, reduced by principal component analysis (PCA) and clustered by K-means (scikit-learn 1.4.2) (Pedregosa et al., Reference Pedregosa, Varoquaux, Gramfort, Michel, Thirion, Grisel, Blondel, Prettenhofer, Weiss, Dubourg, Vanderplas, Passos, Cournapeau, Brucher, Perrot and Duchesnay2011). The elbow method was used to select the optimal number of clusters (K), and cluster centroids were ranked along PC1 to label low-, medium- and high-risk groups.
Statistical analyses
Statistical analyses were performed in Python 3.11 using SciPy and StatsModels. Fisher's exact tests (expected counts <5) or χ2 tests (expected counts ≥5), with Benjamini–Hochberg correction, were used to evaluate associations; q <0.05 was considered statistically significant. Full 2 × 2 contingency tables, P-values, q-values and odds ratios (ORs) with 95% confidence intervals are provided in the Supplementary Tables. When a 2 × 2 table contained a zero cell, only the exact two-tailed P value was reported and the odds ratio was denoted as not estimable (NE).
A complete laboratory and bioinformatic workflow is summarised in Fig. 1.
Workflow for isolation, phenotypic characterisation and genomic analysis of Staphylococcus aureus from bulk-tank milk, and subsequent functional-risk stratification.
Figure 1 Long description
The infographic illustrates a workflow divided into three main sections: Sample Collection, Isolation & Phenotypic Characterisation; Whole-Genome Sequencing & In-Silico Genomic Analyses; and Phylogenomics & Functional-Risk Stratification. 1. Sample Collection, Isolation & Phenotypic Characterisation: - 162 bulk-tank milk samples from 14 dairy farms in Nyanza, Northern Highlands, Kenya. - Enrichment in 7.5 percent NaCl broth. - Streaking on Baird-Parker & blood agar plates. - Identification using VITEK 2 Compact. - Antimicrobial Susceptibility Testing. - n = 59 Staphylococcus aureus. 2. Whole-Genome Sequencing & In-Silico Genomic Analyses: - DNA Extraction. - Illumina NovaSeq 6000 Sequencing. - Assembly of Sequences. - Detection of Resistance Genes. - Detection of Virulence Factors. - Detection of Plasmids. - SNP Phylogenetic Tree. - MLST & spa Typing. - PanMLST & Pan Typing. 3. Phylogenomics & Functional-Risk Stratification: - Functional-Risk Stratification using 2-D Principal-Component Analysis (PCA). - A scatter plot showing Principal Component 1 and Principal Component 2, with points labeled as 'Pathogenic' and 'Non-pathogenic'.
Results
Isolation rate, molecular typing, and core-genome phylogeny
Staphylococcus aureus was recovered from 59 of 162 bulk-tank milk samples (36.4%). Eight spa types and seven MLST sequence types were identified, and two lineages predominated: ST1-t114/t127 (40.7%) and ST9784-t189 (20.3%). ST9784 is a newly assigned sequence type following submission to PubMLST and is defined by a novel gmk allele (gmk-750) (full allelic profile in Supplementary Table S1). Most isolates belonged to CC1 (91.5%), with minor contributions from CC398 and CC59 (Fig. 3; full metadata in Supplementary Table S1).
Core-genome SNP analysis (Fig. 3) produced a well-resolved phylogeny, with 48 of 57 internal nodes supported by bootstrap values ≥70%. Clustering corresponded closely to spa types, whereas ST1 isolates separated into multiple subclades, consistent with within-ST1 genomic diversification. Lineage counts are summarised in Supplementary Table S2.
Antimicrobial susceptibility profiles
Figure 2 summarises non-susceptibility (NS; I + R) rates for the 14 tested antimicrobial agents, consistent with international antimicrobial-resistance surveillance reporting practice (Centers for Disease Control and Prevention (U.S.), 2019; European Centre for Disease Prevention and Control and World Health Organization, 2022). Overall, 27/59 isolates (45.8%) were non-susceptible to at least one agent; non-susceptibility was most frequent for penicillin and ampicillin (28.8% each). Only one isolate (ST59-t437) was cefoxitin non-susceptible, consistent with MRSA. All isolates remained susceptible to Penicillin/Neomycin, ceftaroline, cefazolin, meropenem, and rifampin.
Antimicrobial susceptibility and non-susceptibility rates of 59 Staphylococcus aureus isolated from bulk-tank raw milk.
Figure 2 Long description
The stacked bar graph lists antimicrobial agents with horizontal stacked bars. Columns shown: Antimicrobial agent, Concentration, Abbreviation, Susceptible rate left parenthesis percent right parenthesis, Non-susceptible rate left parenthesis percent right parenthesis. Horizontal axis tick labels: 0 percent, 20 percent, 40 percent, 60 percent, 80 percent, 100 percent. Legend under Non-susceptible rate left parenthesis percent right parenthesis: S, I, R. Rows and values: Penicillin, 10 I U, P: Susceptible rate 71.2; Non-susceptible rate 28.8. Ampicillin, 10 microgram, AMP: Susceptible rate 71.2; Non-susceptible rate 28.8. Penicillin slash Neomycin, 10 I U slash 30 microgram, P N V: Susceptible rate 100. Cefoxitin, 30 microgram, FOX: Susceptible rate 98.3; Non-susceptible rate 1.7. Ceftaroline, 30 microgram, CPT: Susceptible rate 100. Cefazolin, 30 microgram, CZ: Susceptible rate 100. Meropenem, 10 microgram, MEM: Susceptible rate 100. Gentamicin, 10 microgram, CN: Susceptible rate 96.6; Non-susceptible rate 3.4. Erythromycin, 15 microgram, ERY: Susceptible rate 74.6; Non-susceptible rate 25.4. Tetracycline, 30 microgram, TET: Susceptible rate 96.6; Non-susceptible rate 3.4. Ciprofloxacin, 5 microgram, CIP: Susceptible rate 98.3; Non-susceptible rate 1.7. Trimethoprim slash Sulfamethoxazole, 23.5 slash 1.5 microgram, SXT: Susceptible rate 84.7; Non-susceptible rate 15.3. Rifampin, 5 microgram, RIF: Susceptible rate 100. Clindamycin, 2 microgram, DA: Susceptible rate 83.1; Non-susceptible rate 16.9.
Lineage-level differences in non-susceptibility were evident (Fig. 3; Supplementary Table S4-1). The ST1 and ST8396 group (n = 32) was broadly susceptible; ST1 isolates showed significantly lower non-susceptibility to penicillin, ampicillin, and erythromycin (all q < 0.05). In contrast, the newly identified ST9784-t189 lineage (n = 12) exhibited the most pronounced multidrug-resistant (MDR) profile, with universal non-susceptibility to penicillin and ampicillin, and markedly elevated resistance to erythromycin and clindamycin (OR ≥ 20; q < 0.001). ST81-t5100 (n = 7) was uniquely associated with sulfamethoxazole-trimethoprim (SXT) resistance (q = 0.002). Notably, resistance within ST398 differed by spa type: t34 isolates were fully susceptible, whereas t11 isolates showed broad MDR, indicating within-lineage phenotypic heterogeneity.
Core-genome SNP phylogeny of 59 Staphylococcus aureus isolates with MLST/spa types and binary antimicrobial non-susceptibility profiles.
Figure 3 Long description
The image displays a circular phylogenetic tree representing the core-genome SNP phylogeny of Staphylococcus aureus isolates. The tree is divided into colored segments, each representing different MLST types: ST1 in light blue, ST59 in blue, ST81 in purple, ST9784 in pink, ST398 in orange, ST2990 in yellow and ST8396 in green. Surrounding the tree are boxes indicating spa types: t111 in light green, t114 in light blue, t127 in red, t189 in yellow, t437 in pink, t1119 in dark green and t5100 in dark blue. The branches of the tree radiate outward from the center, showing evolutionary relationships among the isolates. Each branch is labeled with specific identifiers and the tree is organized to show the distribution of MLST and spa types among the isolates. The outer ring of boxes provides a visual representation of the spa types associated with each isolate, corresponding to the colors in the legend.
Resistome, virulome, and plasmidome
Figure 4 presents a binary heatmap of antimicrobial resistance genes (ARGs), virulence factors and plasmid replicons across all isolates. Lineage-level counts are provided in Supplementary Table S3, and statistical comparisons are summarised in Supplementary Table S4-2.
Core-genome SNP phylogeny of 59 Staphylococcus aureus isolates with MLST/spa types and binary heatmaps of antimicrobial resistance genes, virulence factors and plasmid replicons.
Figure 4 Long description
The binary heat map displays the presence or absence of antimicrobial resistance genes, virulence factors and plasmid replicons across isolates. Rows represent individual isolate identifiers, including sequence types such as ST1, ST8, ST30, ST398 and others. Columns are grouped into three main blocks: resistance genes, virulence factors and plasmid replicons. Within the virulence factors block, sub-categories are labeled as adhesins and surface proteins, capsule and biofilm, toxins and immune evasion. Filled circles indicate gene or factor presence. Circles within the resistance genes block appear in one distinct shade, circles within the adhesins and surface proteins columns appear in a second distinct shade, circles within the capsule and biofilm columns appear in a third distinct shade, circles within the toxins columns appear in a fourth distinct shade and circles within the immune evasion columns appear in a fifth distinct shade. Each distinct shade corresponds to its respective column category block. Stars mark the presence of plasmid replicons in the plasmid block. Risk group labels are shown to the right of the rows. Four groups are labeled: low-risk group, high-risk group and medium-risk group, with low-risk group appearing more than once. Isolates assigned to the high-risk group show a broader spread of filled circles across resistance gene and virulence factor columns. Isolates in the low-risk groups show fewer filled circles overall. The medium-risk group shows an intermediate pattern of filled circles and stars.
Resistome: All 59 genomes shared a small core resistome: seven efflux-pump genes plus a low-identity vanT-like hit. Accessory ARGs were highly lineage-specific. ST9784-t189 carried most of the accessory resistance genes; blaZ, msrA and lnuA were detected in ≥58% of ST9784-t189 isolates and were enriched compared with other sequence types (e.g., blaZ q < 0.001; OR = NE). In contrast, ST1 isolates rarely carried these genes (e.g., blaZ q = 0.048, OR = NE), consistent with their largely susceptible phenotypes. The two MDR outliers – MRSA ST59-t437 and ST398-t11 – carried dense β-lactam, macrolide, tetracycline and aminoglycoside resistance gene complements, in line with their broad phenotypic non-susceptibility.
Virulome: Adhesin genes (fnbA, clfB, sdrC) were ubiquitous, whereas clfA was comparatively infrequent in ST9784-t189 (q < 0.02). IEC type B (sak-chp-scn) was exclusive to ST9784-t189, while IEC type D (sea-sak-scn) and four classical enterotoxin genes (sea, seb, selk, selq) were detected only in ST59-t437.
Plasmidome: Overall, 79% of genomes carried at least one plasmid replicon. ST1 was strongly associated with MSSA476 (rep7C) (OR = 36, q = 0.002). ST9784-t189 was enriched for pLNU4 and pSJH101/901 (all q < 0.001), consistent with its elevated clindamycin non-susceptibility.
Full cross-lineage comparisons of ARGs, virulence factors and plasmid prevalence are provided in Supplementary Table S4.
Functional-risk stratification
Unsupervised K-means clustering (K = 3, elbow criterion) of 49 integrated genotypic and phenotypic features partitioned the 59 isolates into three groups: high risk (n = 2), medium risk (n = 12) and low risk (n = 45) groups (Fig. 5a–b). A descriptive summary is provided in Table 1. For statistical comparisons, high- and medium-risk isolates were pooled (n = 14) and compared with the low-risk group (n = 45), given the small size of the high-risk cluster (n = 2).
Functional-risk stratification of 59 dairy-derived Staphylococcus aureus isolates.
Figure 5 Long description
The image A showing (a). The x-axis label is K, with values 1, 2, 3, 4, 5, 6, 7, 8, 9. The y-axis label is Inertia, with values 0, 200, 400, 600, 800, 1000, 1200, 1400, 1600. A line with circular markers passes through points at K equals 1 near 1600, K equals 2 near 900, K equals 3 near 200, K equals 4 near 100, K equals 5 near 50, K equals 6 near 20, K equals 7 near 10, K equals 8 near 10, K equals 9 near 10. One marker near K equals 3 is circled. The image B showing (b). The x-axis label is PC1 (log2), with values 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5. The y-axis label is PC2, with values negative 7.5, negative 5.0, negative 2.5, 0.0, 2.5, 5.0, 7.5. A legend titled Risk Level lists Low, Medium, High. Scatter points include a cluster of triangle markers around PC1 from about 0.0 to 1.2 and PC2 from about negative 2.5 to 2.5 and a cluster of square markers around PC1 from about 1.2 to 2.0 and PC2 from about 0.0 to 5.0. Text labels near points include ST9784-t189 near about PC1 1.7 and PC2 4.0, ST398-t11 near about PC1 3.0 and PC2 negative 5.0, ST59-t437 near about PC1 3.4 and PC2 negative 5.5 and MRSA near about PC1 3.2 and PC2 negative 5.2.
Functional-risk stratification of raw-milk Staphylococcus aureus isolates from northern Fujian
Table 1 Long description
The table compares three K-means risk groups of raw-milk Staphylococcus aureus isolates by dominant lineages, multidrug resistance, common non-susceptibility patterns, and enriched resistance genes, virulence factors, and plasmid replicons. The high-risk group had two isolates, both multidrug resistant, with lineages ST59-t437 (including one MRSA) and ST398-t11 and broad non-susceptibility patterns spanning beta-lactams and multiple other drug classes. High-risk isolates were enriched for blaZ, mecA, ermB, tet(M), tet(45), and an aminoglycoside gene cluster, and their virulence profile emphasized toxins and immune-evasion, including sea, seb, selk, selq, plus sak and scn; IEC type D was noted for ST59-t437. The medium-risk group had twelve isolates, all ST9784-t189, with multidrug resistance in six and non-susceptibility commonly involving penicillin and ampicillin with erythromycin and clindamycin; blaZ, msrA, and lnuA were enriched. Medium-risk virulence factors were mainly immune-evasion, with sak, chp, scn, and lukF-PV, and IEC type B. The low-risk group had forty-five isolates across several lineages dominated by ST1 variants and had no multidrug resistance; non-susceptibility was mainly penicillin and ampicillin, with some enrichment of erythromycin with clindamycin and trimethoprim-sulfamethoxazole. Low-risk resistance genes were limited to the core efflux gene lmrS, while virulence factors emphasized adhesins and capsule or biofilm, including clfA and clfB and sbi, with lukF-PV described as sporadic and no IEC type reported. Multiple-plasmid carriage was enriched in high- and medium-risk groups but not detected in low-risk, with different key replicons by group. Arrows indicate features more common in one group than the others, and the high-risk group is very small, so its patterns should be interpreted cautiously.
Notes: Summary of predominant MLST-spa lineages, multidrug resistance (MDR) prevalence, common phenotypic non-susceptibility (NS) patterns, and cluster-associated antimicrobial resistance genes (ARGs), virulence factors (VFs) and plasmid replicons across the three K-means-defined risk groups. Arrows (↑) denote features enriched in a given risk group relative to the other groups; ‘0’ indicates not detected/0% in that group; and ‘–’ indicates absent/not applicable. Antimicrobial abbreviations are as defined in Fig. 2. Full 2 × 2 contingency tables, odds ratios (ORs) with 95% confidence intervals, Fisher's exact-test P-values, and Benjamini–Hochberg–adjusted q values are provided in Supplementary Table S4.
High-risk group (n = 2): The two isolates (MRSA ST59-t437 and ST398-t11) were multidrug-resistant and showed broad phenotypic non-susceptibility across multiple drug classes. Both carried multiple resistance determinants (including blaZ, mecA, macrolide resistance genes, tet family genes and aminoglycoside-modifying enzymes) and multiple plasmid replicons. The ST59-t437 isolate carried IEC type D (sea–sak–scn) and additionally harboured enterotoxin genes (seb, selk, selq).
Medium-risk group: This cluster consisted exclusively of the newly identified lineage ST9784-t189. Half of the isolates were MDR. The group was significantly enriched for blaZ, msrA, lnuA, IEC type B (sak–chp–scn) and the clindamycin-associated plasmid pLNU4 (all q < 0.001; Supplementary Table S4-2).
Low-risk group: The low-risk cluster was dominated by ST1 and ST8396 and had a sparse accessory resistome. MDR prevalence was 0% (0/45), and plasmid carriage was limited, with rep7C (MSSA476) predominating.
Pooling the high/medium groups versus the low-risk group confirmed a clear gradient. MDR (8/14 vs 0/45, q < 0.001; OR = NE), blaZ (5/14 vs 9/45, q = 1.09 × 10−⁴; OR = 2.22) and msrA (7/14 vs 7/45, q = 1.20 × 10−3; OR = 5.43) were enriched in the high/medium category. IEC type B remained a distinctive marker (12/14 vs 0/45, q = 9.85 × 10−⁸, OR = NE). Taken together, the enrichment of IEC type B and plasmid-associated resistance determinants in the ST9784-t189 (medium-risk) cluster suggests a mobile-element–rich lineage with elevated resistance potential, warranting targeted surveillance.
Discussion
In this study, S. aureus was detected in 36.4% of bulk-tank milk samples from northern Fujian. This prevalence lies towards the upper end of estimates reported in China (Supplementary Table S8) and exceeds several international reports, although direct comparison is limited by heterogeneity in sampling strategies and laboratory methods (Zhang et al., Reference Zhang, Wang, Jin, Li, Zhang, Shi and Zhao2022). Using WGS together with unsupervised K-means clustering, we stratified 59 isolates into three exploratory functional-risk tiers (high, medium and low). This grouping provides a basis for prioritising surveillance and on-farm control alongside routine culture and susceptibility testing and should not be interpreted as a stand-alone clinical decision tool.
The low-risk tier comprised several sequence types – dominated by ST1-t114/t127 and ST8396-t114 – and accounted for most isolates. These strains showed little accessory resistance, typically carrying only a small MSSA476-type replicon alongside core efflux-associated genes, and MDR was not observed in this group. Within ST1, the core-genome phylogeny resolved multiple subclades, indicating substantial within-lineage genetic heterogeneity. Similar patterns have been interpreted as compatible with recurrent human-to-bovine spillover and ongoing host adaptation (Greening et al., Reference Greening, Zhang, Midwinter, Wilkinson, McDougall, Gates and French2021; Huang et al., Reference Huang, Zhang, Sun, Jiang, Cao, Shang, Zhang, Gu, Lv, Guo, Li, Li, Guo, Zhu, Huang and Li2023). In practice, routine mastitis control measures (including post-milking teat disinfection and standard dry-cow management) are likely to be adequate for this tier, alongside culture and susceptibility testing and local veterinary guidance.
In contrast, the medium-risk tier comprised exclusively the newly identified lineage ST9784-t189 (approximately 20% of isolates). This lineage exhibited a characteristic accessory gene profile, including blaZ, msrA and lnuA together with immune-evasion genes, and half of isolates were multidrug-resistant. ST9784 is a newly assigned sequence type, distinguished from livestock-associated ST188 by a novel gmk allele (gmk-750), and clustered closely with ST188 in the core-genome phylogeny, suggesting recent emergence and possible local establishment.
Within ST9784-t189, clindamycin non-susceptibility was confined to isolates carrying both pLNU4 and lnuA (6/6 vs 0/6; Fisher's exact test, P = 0.002). Among non-susceptible isolates, lnuA co-localised with the pLNU4 replicon on the same assembled contig in 5/6 (83.3%), consistent with a plasmid-associated lincosamide resistance determinant (Lüthje et al., Reference Lüthje, von Köckritz-blickwede and Schwarz2007); however, short-read assemblies do not resolve plasmid architecture or transmissibility, which warrants confirmation using long-read sequencing and transfer assays (Huisman et al., Reference Huisman, Vaughan, Egli, Tschudin-Sutter, Stadler and Bonhoeffer2022). Taken together, these findings suggest that macrolide–lincosamide regimens (including pirlimycin) should be guided by susceptibility testing for ST9784-t189, and merits further evaluation (More et al., Reference More, McAloon, Silva Boloña, O'Grady, O'Sullivan, McGrath, Buckley, Downing, Kelly, Ryan and McCoy2022). Where supported by susceptibility data, β-lactamase-stable options and strengthened environmental measures (e.g., bedding management and fly control) may therefore be prioritised for this tier.
The high-risk tier contained two isolates: one LA-MRSA ST59-t437 (detected in raw milk in Fujian for the first time in this study) and one MDR ST398-t11 isolate. Both carried broader resistance gene complements, including macrolide resistance and tetracycline resistance determinants.
The ST59-t437 isolate remained susceptible to ceftaroline and meropenem, supporting the fitness-cost hypothesis typically associated with community-associated MRSA (David and Daum, Reference David and Daum2010; Wang et al., Reference Wang, Xu, Zhao, Wang, Rao, Guo, Yi, Hu, Chen, Han, Zhou, Xiang, Hu, Chen and Yu2022a). Furthermore, IEC type D (sea–sak–scn) and multiple enterotoxin genes (sea, seb, selk and selq) were confined to ST59-t437. This contrasts with the IEC loss often described in livestock-adapted lineages and is compatible with relatively recent human-associated ancestry (Cuny et al., Reference Cuny, Abdelbary, Layer, Werner and Witte2015; Wardyn et al., Reference Wardyn, Forshey, Farina, Kates, Nair, Quick, Wu, Hanson, O'Malley, Shows, Heywood, Beane–Freeman, Lynch, Carrel and Smith2015). Its sak⁺/PVL− profile aligns with Asia-Pacific ST59 lineages, and its phylogenetic proximity to hospital-derived strains is consistent with possible human-to animal transmission (Bai et al., Reference Bai, Chen, Lin, Ye, Fan, Wen, Zeng, Huang, Mo, Lei and Liao2021; Wang et al., Reference Wang, Xu, Zhao, Wang, Rao, Guo, Yi, Hu, Chen, Han, Zhou, Xiang, Hu, Chen and Yu2022a). The presence of a compact SCCmec IV element may facilitate persistence and spread across host contexts (Katayama et al., Reference Katayama, Ito and Hiramatsu2000; Wang et al., Reference Wang, Xu, Zhao, Wang, Rao, Guo, Yi, Hu, Chen, Han, Zhou, Xiang, Hu, Chen and Yu2022a).
ST398 displayed striking spa-type heterogeneity. Although based on few isolates, the ST398-t11 (high-risk tier) showed MDR across seven antimicrobial classes; it was the only isolate resistant to ciprofloxacin and carrying fexA, and it harboured multiple plasmid replicons. By contrast, ST398-t34 (low-risk tier) were fully susceptible, with no plasmid replicons or major acquired resistance genes detected. This contrast may reflect diversification within CC398 under differing host niches and antimicrobial selection pressures in the dairy environment (Cuny et al., Reference Cuny, Abdelbary, Layer, Werner and Witte2015).
The detection of the ancient Enterococcus-derived plasmid pAMalpha1 carrying tet(L) in ST398-t11, together with pKH1 carrying tet(K) in ST59-t437, underscores the circulation of distinct tetracycline-resistance elements and is consistent with ongoing interspecies exchange of mobile genetic material on local dairy farms (Yagi and Clewell, Reference Yagi and Clewell1977; Francia and Clewell, Reference Francia and Clewell2002). Given reports of MRSA carriage among dairy workers of up to 14% in some settings (Wardyn et al., Reference Wardyn, Forshey, Farina, Kates, Nair, Quick, Wu, Hanson, O'Malley, Shows, Heywood, Beane–Freeman, Lynch, Carrel and Smith2015), strengthened biosecurity and risk-proportionate periodic screening may be considered in accordance with local regulations and farm-level risk assessment.
Despite long-standing β-lactam use in dairy production (Gustafson and Bowen, Reference Gustafson and Bowen1997), overall β-lactam non-susceptibility in this study (28.8%) was lower than that reported in several other regions of China (Kou et al., Reference Kou, Cai, Huang, Ni, Luo, Qian, Ji and Wang2021; Huang et al., Reference Huang, Zhang, Sun, Jiang, Cao, Shang, Zhang, Gu, Lv, Guo, Li, Li, Guo, Zhu, Huang and Li2023; Ning et al., Reference Ning, Zhou and Li2023; Cui et al., Reference Cui, Wang, Li, Han, Li, Peng, Zhang, Zhao and Wang2024). As blaZ carriage only partly accounted for β-lactam phenotypes, other mechanisms may contribute (point mutations in intrinsic pbp1–3 genes or overexpressed PBP4 that diminish β-lactam binding, and synergistic pbp4–gdpP mutations that elevate intracellular c-di-AMP, thicken the cell wall and confer broad-spectrum resistance) (Ba et al., Reference Ba, Harrison, Edwards, Holden, Larsen, Petersen, Skov, Peacock, Parkhill, Paterson and Holmes2014; Da Costa et al., Reference Da Costa, De Oliveira, Chambers and Chatterjee2018; Lai et al., Reference Lai, Satishkumar, Cardozo, Hemmadi, Marques, Huang, Filipe, Pinho, Chambers and Chatterjee2024). Targeted molecular work will be needed to test these hypotheses.
Overall, the risk-tier approach presented here is intended to support discussion of relative hazards and antimicrobial stewardship priorities, rather than to prescribe treatment. To aid interpretation, we provide an illustrative intervention matrix (Table 2) that should be applied in conjunction with culture and susceptibility testing and within prevailing stewardship guidance (including restrictions on critically important antimicrobials in food-producing animals). Within this framework, our findings may support risk-tier–aware mastitis control and milk-quality assurance in subtropical dairy regions, and help to prioritise surveillance for S. aureus in China.
Illustrative risk-tier options (to be tailored with veterinary oversight)
Table 2 Long description
The table groups mastitis-associated Staphylococcus aureus clones into high, medium, and low tiers and links each tier to lactation treatment, dry-off strategy, and biosecurity or zoonotic precautions. High tier examples include ST59-t437 and ST398-t11, with culture- and susceptibility-guided therapy, possible segregation or milking last, and removal from the milking herd for chronic or recurrent cases. High tier dry-off emphasizes selective dry-cow therapy with internal teat sealant when indicated and avoiding critically important antimicrobials unless clearly justified, alongside reinforced milking hygiene, equipment sanitation, and discouraging raw-milk consumption. The medium tier example ST9784-t189 recommends penicillinase-stable intramammary beta-lactams when susceptible, avoiding macrolides or lincosamides when non-susceptible, and strengthening bedding and fly control. Medium tier dry-off remains selective and guided by somatic cell count history and culture, with monitoring of bulk-tank somatic cell count and periodic culture-based surveillance. Low tier examples ST1-t114 or t127 and ST8396-t114 focus on treating only when supported by clinical assessment and culture, following standard herd dry-off protocols, and maintaining routine hygiene procedures with periodic surveillance as resources allow. The tiers are illustrative and should be adapted with veterinary oversight and local policy, especially for staff screening and antimicrobial choices.
Several limitations warrant consideration. Sampling was confined to a single province and the high-risk tier was represented by few isolates; accordingly, inferences beyond this setting should be drawn with caution. Bulk-tank SCC values and some herd-level metadata (e.g., mastitis history and antimicrobial use) were not available for analysis, which limited our ability to relate isolate profiles to udder-health status and to adjust for potential confounding when comparing herds. Finally, the proposed functional-risk tiers are exploratory and depend on feature encoding and unsupervised clustering choices; larger, multi-region datasets integrating SCC and treatment metadata will be needed to assess robustness and practical utility.
Conclusion
This pilot, single-province study provides a genomic snapshot of S. aureus in bulk-tank raw milk from northern Fujian and introduces an exploratory, unsupervised functional-risk stratification. High-risk isolates were uncommon (3.4%), but the detection of LA-MRSA ST59-t437 in raw milk raises a potential zoonotic concern and is compatible with recent human-to-animal transmission, though this could not be confirmed in the present study. The medium-risk lineage ST9784-t189, characterised by moderate MDR and a distinct plasmid/IEC profile, warrants continued surveillance for possible local establishment and spread. Across tiers, resistance features co-segregated with mobile genetic elements (e.g., pLNU4, pAMalpha1 and SCCmec), emphasising the role of gene mobility in shaping functional risk. Larger, multi-region studies with richer metadata will be required to assess generalisability and to determine how best to translate this framework into practical, risk-based stewardship under veterinary oversight.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0022029926102106.
Acknowledgements
We are grateful to Dr Junxiang Chen (Fujian Institute of Research on the Structure of Matter, Chinese Academy of Sciences) for generous discussions and technical advice on implementing and interpreting the functional-risk stratification workflow.
Competing interests
The authors declare no competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Funding
This work was supported by the Science & Technology Program of Fujian Administration for Market Regulation (Grant No. FJMS2022008).
Data availability statement
The whole-genome assemblies generated in this study have been deposited in the NCBI GenBank/Assembly under BioProject PRJNA1258060. Individual BioSample accessions are listed in Supplementary material Table S1 (SAMN48309129 – SAMN48309187).
Author contributions
S.Z. conceived and designed the study, obtained funding, carried out the statistical analyses and functional-risk stratification, and drafted the manuscript. D.F. led herd enrolment, bulk-tank sampling and laboratory testing, and contributed to data curation and quality control. P.M. assisted with sample processing, antimicrobial susceptibility testing and data curation. Y.H. performed the whole-genome sequencing–based bioinformatic analyses, including typing, phylogenomics and multivariate clustering, and assisted with figure preparation. T.H. provided overall supervision, contributed to study design and interpretation of the results, and critically revised the manuscript. All authors contributed to data interpretation, commented on successive drafts and approved the final version of the manuscript.
AI-assisted tools
The authors used an AI-based language tool only to assist with grammar and phrasing. All scientific content, including study design, data generation, analysis and conclusions, was developed and checked by the authors.
Ethical statement
This study based on bulk-tank milk did not involve experimental procedures on animals or humans and complied with local regulations for food safety surveillance.
