Dairy management practices associated with multi-drug resistant fecal commensals and Salmonella in cull cows: a machine learning approach

Farms surveys and sampling

The study was approved by the University of California, Davis’s Institutional Animal Care and Use Committee (protocol number 18019). Six dairy farms located in the San Joaquin Valley of California were reenrolled as a convenience sample and followed up for a second year after being part of an earlier study (Abu Aboud et al., 2016). Briefly, cull cows were identified for fecal sampling once during each season between 2015 and 2016, specifically during summer (July 1–September 30, 2015), fall (October 1–December 31, 2015), winter (January 1–March 31, 2016) and spring (April 1–June 30, 2016). The choice of sampling week to collect fecal samples from the cull cows during any of the four seasons was also by convenience. From the list of cows selected by the dairy farms for culling, 10 cows were randomly selected for fecal sampling on the day of their removal from the herd using a random number generator (Excel; Microsoft Corp., Redmond, WA, USA). Several lists of random numbers were prepared in multiples of 10 ranging from 11–20 to 91–100 cows and the respective list was selected depending on the number of cows presented for culling on the day of sampling. If a producer had less than 11 cows presented for culling, all the cows were sampled. An individual disposable sleeve was used to manually collect fecal samples from the rectum of each selected cow. The fecal samples were transported to the Aly Lab (Dairy Epi Lab) on wet ice for processing on the same day.

A survey was also completed with the help of the herd manager on the same sampling day. The survey questions were described in an earlier report (Pereira et al., 2019). Briefly, the questions targeted management of the herd in the previous 4 months and collected information on herd size, breed, rolling herd average, cull rate, frequency of culling per month, the proportion of cows sold for beef (compared to as dairy), proportion and reasons for condemnation of culled cows. The survey also collected information on the proportion of culled cows that received injectable medical treatments in the 3 weeks prior to culling, the role of dairy staff allowed to treat cows on the dairy, practices to avoid drug residue violations (use of specific drug types, following withdrawal times, milk and/or urine testing prior to the cow being culling, or other practices), tracking of drug withdrawal intervals, having a drug inventory system in place, and extra-label drug use (ELDU, frequency and familiarity). Finally, a backup of the herd record software was obtained to collect information on the culled cows’ milk production and health events. A relational database was used to house and merge data from the surveys, dairy records, and test results (Access; Microsoft Corp., Redmond, WA, USA).

Bacteriological culture

The California Animal Health and Food Safety (CAHFS) laboratory conducted all the study sample testing for Salmonella as described by Adaska et al. (2020). Briefly, 1 g of feces was inoculated into tetrathionate selective enrichment broth and incubated at 37 ± 2 °C. The next morning a cotton swab was used to inoculate the overnight broth onto XLD and XLT-4 plates and these were incubated overnight at 37 ± 2 °C. H₂S positive, Salmonella suspect colonies from each set of plates were subcultured onto individual bi-plates (5% sheep blood agar-MacConkey agar) and incubated overnight at 37 ± 2 °C. One colony from each bi-plate was used for biochemical testing which included triple sugar iron, urea, motility indole ornithine, citrate, O-nitrophenyl-beta-D-galactopyranoside, and lysine iron agar slants. Serogrouping and serotyping were performed, using the White–Kauffmann–Le Minor scheme (Grimont & Weill, 2007) on colonies with biochemical test results consistent with Salmonella (Quinn et al., 2002).

Fecal samples were also cultured for E. coli and Enterococcus spp. isolation as described previously (Li et al., 2014). Briefly, a 40 mL solution of buffered peptone water containing 5 g of feces in a 50 mL polypropylene tube was homogenized using a mechanical shaker for 15 min before filtering using gauze. A total of 1000, 100 and 10 μl of the filtered solution were then streaked on CHROMAgar ECC (Chromagar, Paris) and Enterococcus Indoxyl-β-D-Glucoside agar (mEI) plates (Becton, Dickinson and Company, Franklin Lakes, NJ, USA) both incubated at 37 °C for 24 h. Reference strains ATCC 25922 (E. coli) and ATCC 29212 (Enterococcus faecalis) were plated on agar plates as positive controls. Two pure colonies were isolated of each species after presumptive colonies were confirmed using biochemical tests (E. coli were confirmed using urea, indole, triple sugar iron, Methyl Red–Voges–Proskauer and Simmons citrate; Enterococcus were confirmed using bile esculin, brain heart infusion agar, and growth in broth with and without 6.5% NaCl).

Antimicrobial susceptibility testing

Salmonella and E. coli bacterial resistance was evaluated with a broth microdilution method using a Gram-negative Sensititre plate (CMV2AGNF) and Enterococcus spp. evaluated on gram-positive Sensititre plate (CMV3AGPF) (Sensititre, Thermo Fisher Scientific, MA, USA) according to the manufacturer’s instructions and as described in previous studies (Li et al., 2018; Pereira et al., 2019). The minimum inhibitory concentration (MIC) values were the lowest concentrations of AMD that inhibited visible growth of bacteria. Interpretations of antimicrobial resistance were based on breakpoints recommended by the National Antimicrobial Resistance Monitoring System (https://www.cdc.gov/narms/antibiotics-tested.html; and https://www.ars.usda.gov/ARSUserFiles/60400520/NARMS/ABXEntero.pdf) and the Clinical and Laboratory Standards Institute (CLSI, 2014; CLSI, 2018). Due to the inherent resistance of Salmonella to cephalosporins, aminoglycosides, lincosamides, oxazolidinones and glycopeptides, drugs from these classes were excluded from the antimicrobial resistance analysis. In addition, the following drug classes to which E. coli is inherently resistant were excluded from the analysis: lincosamides, oxazolidinones, penicillins, streptogramins, glycopeptides. Similarly, the drug classes exclude due to the inherent resistance of Enterococcus spp. included cephalosporins, lincosamides, fluoroquinolones, aminoglycosides, aminocyclitols, sulfonamides, folate pathway antagonists. Isolates from any of the three species (Salmonella, E. coli, Enterococcus spp.) were identified as multi-drug resistant if resistance to at least one AMD in each of three or more drug classes was observed (Magiorakos et al., 2012).

Development of classification algorithms

Classification tree models were developed to test if herd management practices and features related to dairy cows can predict multi-drug resistance phenotype in Salmonella and fecal commensal E. coli, and Enterococcus spp. isolated from the same cows. A cow was considered to shed MDR bacteria if any of its Salmonella, Enterococcus spp. and/or E. coli isolates showed resistance for three or more antimicrobial classes (regardless of whether the species with resistance was Salmonella, Enterococcus spp. and/or E. coli), in which case the cow was labeled as ‘shedding bacteria with multi-drug resistance (MDR) phenotype’ (numerical label = 2). In contrast, if a cow shed bacteria that were resistant to only one or two antimicrobial classes (regardless of species) the bacteria isolated were labeled as ‘shedding bacteria with antimicrobial resistance (AMR) phenotype’ (numerical label = 1). If the bacteria isolated from a cow did not show any resistance across all three bacterial species to any AMD, the bacteria were labeled as non-resistant (numerical label = 0). The mutually exclusive definitions were necessary to develop a single model that predicts one of three resistance states MDR, AMR, or no resistance. Similarly, cows were also classified based on resistance phenotypes separately observed in each bacterial species isolated (three separate classification labels based on resistance phenotypes of Salmonella, Enterococcus spp. and E. coli isolates). Finally, classification labels were also generated based on resistant phenotypes observed in either commensal bacteria (Enterococcus spp. and E. coli) shed in feces collectively. Using features from herd surveys, classification tree models were trained to predict the MDR phenotype of bacteria shed in the feces of the study cows.

Three machine learning algorithms, specifically decision tree classifier (DTC), random forest (RF) and gradient boosting (GB) were developed to explore the risk factors for resistance phenotypes in the study isolates using herd survey data, specifically to predict the multilabel outcome based on the resistance phenotype of bacteria shed by cows (non-resistant, AMR, and MDR). For each algorithm, data from the entire study cohort described by Pereira et al. (238 cows) were used, except models specific to predicting AMR phenotypes in Salmonella were restricted to the cohort of 58 Salmonella positive cows. Table 1 describes all classification algorithms trained and developed for various bacteria-specific outcomes.

For all three classification algorithms, 25 features related to herd management and 12 features related to individual cows collected from the survey were used as predictive features (Table 2). The DTC generates an optimum tree based on attributes by recursive selections to split data into classes and it was only used to visualize the optimum tree and as a contrast to the remaining classification tree models (RF and GB) that prevent overfitting, unlike DTC. The RF and GB algorithms both generate a series of recursive trees of binary splits for randomly sampled predictor variables. While all tree classification algorithms handle interaction effects between predictors, within GB, boosting builds and combines collective models improving the predictive performance of many weak models substantially, and fits complex nonlinear relationships (Elith, Leathwick & Hastie, 2008). For the validation of each algorithm, the data was split into training and validation datasets. To identify the best hyperparameters of the classification algorithms (hyper-tuning), a grid search was implemented on the training dataset. Grid search is a tuning technique that computes optimum values of model parameters automatically by an exhaustive search performed on a set of parameter values. Training datasets, composed of 80% of the data, were randomly selected for gird search, maintaining the outcome proportional to the original dataset, with three-fold cross-validation (internal validation). Model parameters tested for hyper-tuning of DTC, RF, and GB are given in Table 1. The best performing model parameters were chosen based on the accuracy of the model. For each algorithm, the external validation of the best performing model was done on the validation dataset (20% remaining random sample of the original dataset) to quantify the performance of the model on the completely independent validation dataset. The DTC and RF models were implemented using the Scikit-learn machine learning package (Pedregosa et al., 2011) and GB was implanted using XGBoost in python (Chen & Guestrin, 2016).

Validated models were eventually fit on the complete dataset to produce the final predictions. The relative influence (importance) of features for the random forest model was estimated using average gini importance, permutation, and feature drop methods. The gini importance for a feature is defined as the sum over the number of splits (across all tress) that include the feature, proportionally to the number of samples it splits (Breiman, 2001). The gradient boosting model with the XGboost platform was evaluated using Shapely Additive Explanations (SHAP) that assigned each predictive feature an additive feature unifying six existing methods (Lundberg & Lee, 2017). Partial dependence of gradient boosted model prediction on model features (expected output response trend as a function of feature) was explored to understand the associations of the herd and cow-related features with predictions (Friedman, 2001). The Python code used for pre-processing the data, training and validating models and generating figures can be found here in the Zenodo repository: DOI 10.5281/zenodo.4387017.

Distribution of antimicrobial resistance

Seasonality of multi-drug resistance

The overall prevalence of shedding MDR bacteria (Salmonella, E. coli, and Enterococcus spp.) in cows was 30.54%, (±2.97, n = 238). The prevalence of cows shedding bacteria that are resistant to <=2 antimicrobial drug classes (AMR) was 43.93% (±3.21, n = 238). The highest prevalence of MDR was seen in the summer (50.00% ± 11.18) followed by fall (34.00%, ±4.73). The highest prevalence of shedding MDR bacteria was seen in herd 4 (52.00%, ±7.89) and the lowest was seen in herd 3 (10.00% ± 4.74). Herds 1 and 2 each showed the lowest prevalence of shedding AMR bacteria (37.50% ± 7.65), while herd 4 showed the highest prevalence of shedding AMR (50.00%, ±7.90). The seasonal prevalence of MDR and AMR resistance across three bacterial species is presented in Fig. S1.

The annual prevalence of cows shedding MDR Salmonella was 8.62% (±3.65, n = 58). The highest prevalence for MDR was detected in the winter season (10.53%, ±7.04, n = 19) and the highest AMR prevalence was detected in fall (42.31% ± 9.68, n = 26). Salmonella isolated from cows in the summer season did not show any AMR or MDR phenotypes (Fig. 1A). The annual prevalence for MDR phenotypes within E. coli isolates was 2.92% (±1.09, n = 238). The highest prevalence for MDR phenotypes in E. coli was detected in winter (4.08% ± 2.82, n = 49) and the highest AMR prevalence in E. coli was detected in spring (34.28% ± 5.67, n = 70, Fig. 1B). The annual prevalence for MDR phenotypes within Enterococcus spp. isolates was 26.35% (±2.84, n = 238). The highest prevalence for MDR Enterococcus spp. was detected in the summer season (50.00% ± 11.18, n = 20) and the highest AMR prevalence in Enterococcus spp. was detected in the winter season (48.97% ± 7.14, n = 49, Fig. 1C).

Model performance and tuning results

Grid search of model parameters with 3-fold cross-validation yielded satisfactory results in classifying MDR cows based on herd management practices and cow-related features. The sensitivity of models ranged from 0.47 to 0.74 with precision ranging from 0.46 to 0.66. Details related to hyperparameters, best performing decision tree classifiers, random forest, and XGboost models and the performance of the selected models are given in Table 1.

Association between herd management practices and cow-related features with shedding multi-drug resistant bacteria

Decision tree classification models

For all DTC models, the impurity (gini), a measure of the heterogeneity of the outcome in a subset of samples resulting from a split in a decision tree, was reduced in samples the most by the rolling herd average milk production, denoting its highest position in decision trees generated by all five models (Fig. 1, Figs. S2–S5). Other management features that formed nodes showing high measures of split quality (gini) including the proportion of Jersey cows in the herd, administration of tetracyclines, monthly culling percentage, culling frequency, and the number of culled individuals. Of the five DCT models, all but the model for overall resistance across Salmonella and commensal bacterial species were able to generate nodes that classify cows based on their fecal shedding of bacteria to a single class of AMD, or not resistant (Figs. S2–S5) with 100% purity (homogenous subsets). While only the DCT model for commensal bacteria, resulted in pure nodes for MDR cows, none of the other models were able to classify cows into pure MDR nodes (Fig. S2). Figure 2 shows the optimum decision tree for classifying cows into MDR, AMR, or not resistant based on phenotypes of all bacteria (Salmonella, Enterococcus spp., E. coli).

Random forest models

Results of random forest models indicated common herd management practices that influence the shedding of MDR and AMR phenotypes of Salmonella, Enterococcus spp., and E. coli collectively, as well as individually (Fig. 3). The number of culled animals, monthly culling frequency and percentage, herd size, and proportion of Holstein cows in the herd were found to be influential herd characteristics (top ten features by relative influence) predicting MDR phenotypes in all algorithms. Random forest algorithms for predicting AMR phenotype in Enterococcus spp., commensal bacterial species combined, and in all bacterial species showed the same top ten most influential herd management features. Individual-level features such as culling due to milk production, mastitis, reproductive and other reasons appeared important for predicting resistance phenotypes in Salmonella. The use of the chalk method for withdrawal determination was in the top ten most influential features and for the E. coli and Salmonella models.

Gradient boosting classification models

Herd characteristics that showed higher variable importance in SHAP values for predicting cows shedding MDR resistance bacteria based on all bacterial species were herd size, the proportion of Jersey cows, sampling season, the frequency of extra-label drug use, rolling herd average, and culling related features. These features also showed high marginal contributions in predicting AMR phenotypes. For predicting AMR phenotypes, features such as the number of monthly veterinary treatments, antibiotic dose tracking, and the proportion of Holstein cows also showed higher SHAP values (Fig. 4).

Partial dependence of gradient boosting model prediction on important features showed the possible relationships of these herd and cow-related features in predicting the AMR of their fecal bacteria as MDR. Higher culling frequency and monthly culling percentages were associated with cows with MDR phenotypes from all bacteria, whereas an overall increase in the number of culled animals from the herd showed a negative trend with classifying a cow as shedding MDR bacteria. Winter season was negatively associated with MDR phenotype bacteria shed by cows compared to cows sampled in Spring. Similarly, herds with more than 10,432 kgs of rolling average milk production showed a positive trend with MDR positive cows. Cows from herds with a higher proportion of Jersey cows were negatively associated with being classified as shedding MDR bacteria by the gradient boosting algorithm (Fig. 5). Within other herd characteristics that were identified as important, herd size showed a varying trend with classifying cows shedding MDR bacteria, with some herd sizes showing a positive association of classifying cows as shedding MDR bacteria (Fig. 5).

For predicting MDR phenotypes in Salmonella shed by cows, rolling herd average milk production, sampling season, chalk methods for tracking withdrawal, monthly culling frequency, and Salmonella vaccine were found to be important predictive features with high SHAP values (Fig. 6). The exploration into partial dependence of these features gave insights into relationships of feature values with classifying a cow as shedding MDR phenotype Salmonella (Fig. 7). Increasing rolling herd average milk production and monthly culling percentage was positively associated with MDR phenotypic Salmonella in cow feces. Similarly, cows sampled in spring were more likely to be classified as shedding MDR Salmonella.

Models predicting phenotypes in commensals E. coli, and Enterococcus spp.; the number of culled animals in the previous year was always in the top ten most important features to classify cows as shedding MDR bacteria (Figs. S5–S8). Rolling herd average milk production features describing culling practices, and herd size, consistently featured as important in classifying cows as shedding MDR phenotypic bacteria for all the other three models (Commensals, E. coli and Enterococcus spp.). The frequency of extra-label drug use was an important feature for models separately predicting MDR in Enterococcus spp. (Fig. S7) and E. coli shed by cows.