INTRODUCTION
Antimicrobial resistance (AMR) is a leading public health concern, causing significant morbidity and mortality, and increased healthcare-related costs [Reference Smith and Coast1]. Antimicrobial use (AMU) in food animals has been posited as a driver of AMR by national and international governing agencies [2]; however, studies in this area have produced ambiguous results [Reference Rao3–Reference Gow7]. The significance of AMU in animals on AMR outcomes in humans is difficult to quantify [Reference Hurd8–Reference Hurd12]. Despite best efforts to address this knowledge gap, substantial logistical challenges remain with respect to designing and conducting commercial field studies to quantify associations between AMU and AMR in food animals. Such challenges include uniquely identifying and tracking animals over time, maintaining and accessing detailed health and antimicrobial drug (AMD) exposure records for individual animals, and obtaining biological samples longitudinally from a sufficient number of these uniquely identified animals.
Even large field studies conducted in commercial feedlots that have overcome these challenges have produced inconsistent results [Reference Checkley5–Reference Gow7], including differences in detecting associations between tetracycline use and resistance in beef cattle feedlots [Reference Rao3, Reference Morley4], and differences in the association between tetracycline resistance and use of chlortetracycline and oxytetracycline in swine [Reference Alali13]. One reason for such differential results may stem from differences in laboratory and analytical methods. Potential bench-side differences include the use of multiple resistance testing methods that may produce different results [e.g. disk diffusion (DD) vs. broth microdilution (BM)] [Reference Benedict14], a lack of host- and agent-specific breakpoints for interpretation of resistance results, and debate as to which indicator bacterial species best represent resistance dynamics in food animals. Analytically, the datasets resulting from large representative studies of AMU-AMR in commercial settings present several challenges including accounting for multiple levels of non-independence (i.e. clustering), accounting for time of sampling, validating quantification of antimicrobial drug exposures, accounting for AMU exposures that occurred prior to feedlot entry, and classifying resistance outcomes when there are many isolates within a sample, and one can use multiple testing methods for resistance characterization.
This report is part of a larger study that aimed to overcome many of these obstacles with the goal of developing a longitudinal AMR surveillance system for use in beef feedlot cattle [Reference Benedict15]. The purpose of this investigation was to estimate resistance prevalence in isolates of non-type-specific E. coli (NTSEC) recovered from faeces of beef feedlot cattle, to model associations between AMU and AMR in these isolates using various modelling strategies, and to compare these strategies in order to understand how they may differentially impact measures of association between AMU and AMR.
MATERIALS AND METHODS
Study overview
Pens of cattle (N = 300) in four feedlots in Alberta, Canada were enrolled in the study and composite pen-floor faecal samples were collected at the beginning of the feeding period, and at least once again during the remainder of the feeding period. Samples were cultured for NTSEC and isolates were susceptibility tested using BM and DD. Antimicrobial exposures for all cattle in enrolled pens were tracked throughout the study period. The prevalence of recovery of AMR NTSEC was estimated, various methods of multivariable logistic regression were used to investigate associations between AMU and AMR in NTSEC isolates, and model results were compared.
Sample collection
Bacterial isolates were collected as part of a project intended to develop and evaluate methods for surveillance of AMR in beef feedlot cattle, the details of which have been reported elsewhere [Reference Benedict15]. Sample collection ran from September 2007 to January 2010, during which time feedlot collaborators enrolled 30% of all arriving pens using a randomization table. Pens were sampled when filled to capacity, and at least once again later during the feeding period. Twenty fresh faecal samples of ~14 g each were collected in a standard spatial pattern from each pen floor and then mixed together for 1 min, and 10 g was removed for further processing. All cattle handling and sampling procedures were approved by the Animal Care Committee of Feedlot Health Management Services Ltd (FHMS) and the Institutional Animal Care and Use Committee of Colorado State University.
Susceptibility testing
Faecal samples were cultured for NTSEC, and up to five isolates from each sample were selected for susceptibility testing, which was performed using both automated DD (BioMIC, (Giles Scientific USA)) and BM (Sensititre, Thermo Fisher Scientific, USA, panel type: CMV1AGNF). The antimicrobial panels differed between the two test methods, as they were selected independently for surveillance purposes. An automated visual imaging system recorded zone diameter from DD tests, while laboratory personnel recorded the minimum inhibitory concentration (MIC) from BM tests. All susceptibility testing was conducted in accordance with standards established by the Clinical Laboratory Standards Institute [16–20]. Details about quality control assessments and interpretive criteria have been previously described in detail [Reference Benedict15].
Antimicrobial use data
Feedlot personnel used a chute-side, customized health information system (iFHMS software, Feedlot Health Management Services, Canada) to record all parenteral and in-feed AMD treatments for all cattle in enrolled pens from the day of arrival until the day that the last sample was collected (Table 1). Records were subsequently exported for analyses, and included individual animal and pen identification, AMD type, dosage, and route and date(s) of administration. For analysis, antimicrobial doses were converted to units of animal defined daily dose (ADD), a metric that defines the number of days that a single dose remains in the recipient's target tissue(s) based on labelling claims for therapeutic use in treatment of respiratory disease, and therefore approximates ‘exposure days’ for one standard treatment [Reference Benedict21]. This standardization enables comparison of AMD exposures across varying sizes of cattle with varying antimicrobial dosages. To aggregate AMD exposures at the group (pen) level, each treatment was converted to ADDs, multiplied by both the number of animals exposed and the duration of treatment (in days), and then summed by AMD class (β-lactams, phenicols, quinolones, sulfonamides, tetracyclines and macrolides). The AMD categories were further subdivided into the route of administration as being in-feed or parenteral.
Antimicrobial drugs used in this population
AMD, Antimicrobial drug; ADD, animal defined daily dose; BW, body weight; BRD, bovine respiratory disease.
* Assuming 9 kg of dry matter intake per individual animal per day.
Multivariable modelling to estimate adjusted prevalence of resistance
Adjusted prevalence of resistance in individual isolates was estimated from marginal (adjusted) means obtained from Poisson regression modelling. Estimates were stratified by the number of days that cattle had been in the feedlot when faecal samples were collected: 0–3 days-on-feed (DOF), 4–70, 71–120, 121–180, and >180 DOF. This enabled us to adjust resistance prevalence relative to DOF, and to determine if DOF was a statistically significant predictor of AMR levels. Cut-offs for DOF categorization were chosen based on the likelihood of disease occurrence and therefore AMU at different phases of the feeding period, as well as the distribution of collected samples. Generalized estimating equations (GEE) with an exchangeable correlation structure were used to control for clustering of isolates within samples. For isolates that were susceptibility tested using both BM and DD (repeated measures), we included only BM results in order to meet the assumption of independence between observations. Pen size and feedlot could potentially confound resistance prevalence, and thus were included as fixed effects.
Multivariable modelling to test for associations between use and resistance
A variety of multivariable modelling techniques were used to analyse potential associations between AMU and AMR in NTSEC isolates. The primary outcome of interest in all models was the antimicrobial susceptibility status of NTSEC isolates, defined dichotomously as either resistant or non-resistant, the latter of which included intermediate and susceptible classifications. In order to maintain temporal logic for associations between AMU and AMR, only AMDs given prior to sample collection were included in these analyses, and therefore only isolates recovered from non-arrival samples were interrogated. Resistance status for each of the 19 AMDs included on the test panels was evaluated in parallel in separate models.
The primary independent variables of interest were exposures to AMDs, which were categorized and summed as described above. Despite the inherently continuous nature of pen-level AMD data, class-specific distributions for AMDs were strongly right-skewed and zero-inflated. We attempted to assess linearity using quadratic terms for AMD exposures, as well as by modelling deciles and quintiles of AMD exposures in order to compare parameter estimates with the logit of the resistance outcome. However, these models either would not converge or the Hessian matrix was not positive definite, and linearity could not be formally assessed. Therefore, AMD exposures were modelled as both continuous and categorical variables and model results were compared (Table 2). For categorization, parenteral AMD exposures were dichotomized (‘no exposure’ vs. ‘any exposure’), while in-feed exposures were categorized into four levels: no exposure, low exposure (<25th percentile), medium exposure (25th–89th percentiles) and high exposure (⩾90th percentile). To investigate the impact of temporality on the association between AMD exposures and AMR, we grouped AMD exposures based on the time of sample collection as either recent or non-recent, i.e. AMDs administered ⩽6 days prior to sample collection and those administered ⩾7 days prior to sample collection.
Model specifications for all 10 models used to analyse associations between AMU and AMR in NTSEC isolates
AMU, Antimicrobial use; AMR, antimicrobial resistance; NTSEC, non-type-specific Escherichia coli; AMD, antimicrobial drug; GEE, generalized estimating equations; GEE/ALR, generalized estimating equations with alternating logistic regression; GLMM, generalized linear mixed modelling; BM, broth microdilution; DD, disk diffusion.
Two variables were included in all models as potential confounders: number of animals housed in the pen from which the composite faecal sample was obtained (‘pen size’) and DOF. Pen size was modelled as a five-level ordinal variable due to nonlinearity with the logit of the outcome (<101, 101–200, 201–300, 301–400, and >400 animals). The DOF variable exhibited a linear relationship with the logit of the outcome and was modeled as both a continuous and an ordinal variable (0–3, 4–70, 71–120, 121–180, and >180 DOF) for comparison purposes.
Population and study design factors created numerous issues related to data hierarchy, clustering and repeated measures. Resistance outcomes could be clustered within feedlot (n = 4 feedlots), within pens (n = 275 pens), and within samples, as multiple isolates were collected from each sample (n = 564 samples). In addition, repeated measures were present at two levels: first, pens were sampled multiple times throughout the feeding period; and second, a subset of NTSEC isolates were tested by two different resistance testing methods.
Because there are a number of equally valid modelling approaches to analyze such data, we used and compared a variety of methods. In all model types, feedlot was included as a categorical fixed effect due to the small number of feedlots and the fact that none of the model predictors were considered to be feedlot-level effects. GEE with alternating logistic regression (ALR) was used to account for clustering of isolates within samples as well as repeat susceptibility testing on some isolates [Reference Carey, Zeger and Diggle22]. For these models, sample identification number was specified as the repeated subject with an exchangeable correlation structure and susceptibility testing method was specified as the subcluster with a 1-nested log odds ratio structure. When data sparseness did not support GEE with ALR modelling (i.e. the correlation matrix was not positive definite or parameter estimates were unrealistically large), GEE without ALR was used, and only BM results were analysed to avoid issues of non-independence. Finally, generalized linear mixed modelling (GLMM) with random effects and Laplace estimation [Reference Wolfinger23] was used to control for multiple levels of clustering. When model convergence was achieved, isolates, samples and pens were specified as random effects. When convergence could not be attained on all three levels of clustering, only susceptibility results from BM were used in the analysis and separate models specifying pens and samples as random effects were compared. To facilitate comparison of model results, subject-specific (SS) parameter estimates from GLMM models were converted to population-averaged (PA) parameter estimates using the equation [Reference Zeger, Liang and Albert24]:
$$\beta ^{{\rm PA}} \, \approx \,\displaystyle{{\beta ^{{\rm SS}}} \over {\sqrt {1 + 0.346\sigma _h^2}}}. $$
The modelling decisions made at the outset of this analysis, combined with the limitations imposed by model convergence and stability, yielded 10 permutations of modelling methods (models A–J) (Table 2). These included various combinations of model type (GEE with ALR vs. GEE vs. GLMM), clustering specification, and methods for quantification of DOF and AMD exposures (continuous vs. categorical).
The same approach to model development was used for all 10 modelling methods (A–J). First, univariable screening models were used to analyse associations between each class of AMD exposure (n = 8) and every resistance outcome (n = 19 AMDs for models using both BM and DD, and n = 15 AMDs for models using only BM results). The AMD exposures were split into recent and non-recent exposures as described above. Each screening model included fixed-effect variables for feedlot, DOF, and pen size. The AMD exposure variables with a P value ⩽0·20 for either recent or non-recent exposures were included in the initial candidate multivariable model. Backward elimination was used to refine multivariable models using a critical alpha for retention of 0·05 with recent and non-recent AMD exposures considered independently. Variables exhibiting confounding upon stepwise removal (defined as >20% change in any parameter estimate) were added back to the model. Once all variables met the critical alpha value, confounding was reassessed and variables were removed if confounding was no longer present.
Models for each resistance outcome were assessed independently using each of the modelling methods (A–J), resulting in development of 139 different multivariable models (Tables 2 and 3). Model fit was evaluated using Akaike's Information Criterion (AIC) for GLMM models, and Quasi-Information Criterion (QIC) for GEE models.
Results of modelling process for all antimicrobial resistance outcomes
✓ Model converged and results are presented.
✓* Model converged, but no results are presented because no antimicrobial drug exposures were statistically significant in final model.
– Model would not converge.
‡ Model not needed because models with random effects at all levels did converge.
† Model not needed because resistance was only tested with one susceptibility test (i.e. no repeated measures on isolates).
RESULTS
A total of 300 pens were sampled at least once. Twenty-five pens were excluded from analysis due to missing exposure or resistance information, resulting in 275 pens for use in multivariable analyses that ranged in size from 27 to 555 head (median 168). Of these 275 pens, 23% (64/275) were sampled once during the feeding period, 48% (132/275) twice, and 29% (79/275) three times, resulting in collection of 564 composite faecal samples, from which a total of 2911 NTSEC isolates were cultured. Almost all isolates were tested using BM (2903/2911), while 41% (1192/2911) were tested by DD, for a total of 4095 susceptibility test results for analysis. Of the 564 composite faecal samples, 40% (226/564) were collected between 0 and 3 DOF, another 39% (214/564) between 4 and 120 DOF, and 22% (124/564) later in the feeding period (median 152 DOF, range 121–244 DOF).
Prevalence of resistance to most AMDs was too low for valid estimation of least-square means for adjusted prevalence estimates (Table 4). Only six of the 19 resistances tested attained crude prevalence >2%, namely ampicillin, chloramphenicol, streptomycin, sulfisoxazole, tetracycline and florfenciol (Fig. 1). DOF was significantly associated with resistance levels for all modelled AMDs except ampicillin (Fig. 1).
Adjusted prevalence of resistance in non-type-specific Escherichia coli isolates for six antimicrobial drugs (AMDs) (adjusted for days-on-feed at sample collection, as well as pen size and feedlot). Clustering of isolates within samples was controlled for using multivariable generalized estimating equations with a Poisson distribution. * Indicates AMDs for which there were statistically significant (P < 0·05) differences in resistance prevalence by days-on-feed (DOF).
Crude prevalence of antimicrobial resistance in NTSEC isolates (for antimicrobials with prevalence that would not support multivariable modelling of adjusted prevalence of estimates)*
NTSEC, Non-type-specific Escherichia coli; CI, confidence interval.
* Adjusted prevalence of resistance for other drugs can be found in Figure 1.
† Tested by broth microdilution only (n = 2903).
‡ Tested by broth microdilution and disk diffusion (n = 4095).
§ Tested by disk diffusion only (n = 1192).
|| Adjusted CI for binomial proportions (adding two successes and two failures) were estimated as previously described [Reference Agresti and Coull25].
Antimicrobial drug use
The most commonly used AMD was in-feed tetracycline, which was administered to all pens of cattle in this study population; however, because some of the pens were only sampled at arrival, and not again later in the feeding period, only 83% of pens were exposed prior to collection of the ‘last’ sample (Table 5). Pen-level prevalence of exposure to many parenterally administered AMD classes was <50% (Table 5), and even in exposed pens, distribution of AMD exposures had a strong right skew because most pens contained only one or two exposed animals, while a handful of pens received pen-wide metaphylactic treatment. Over 96% of all AMD exposures occurred ⩾7 days prior to sample collection and were classified as ‘non-recent’. This was consistent across all AMD classes with the exception of in-feed macrolides, of which 31·3% of ADDs were classified as ‘recent’, i.e. administered ⩽6 days with respect to sample collection. This disparity stems from the AMD protocols utilized in the four study feedlots.
Total pen-level exposure to AMDs, by class
AMD, Antimicrobial drug; ADD, animal defined daily dose.
* Among pens receiving at least 1 ADD of AMD prior to collection of last composite faecal sample.
Results of multivariable modelling
For AMDs with resistance prevalence sufficient to support multivariable modelling (Table 3) we found several statistically significant associations between AMU and AMR. However, most of these associations were not consistent between different model specifications. Details of comparisons are given below.
Effect of quantifying AMD exposures as continuous vs. categorical variables
The quantification of AMD exposures as continuous or categorical variables exerted the strongest influence over model results, best illustrated in the results for chloramphenicol resistance (Table 6). When AMD exposures were quantified as continuous variables (models A, B and I), only non-recent parenteral macrolide exposures were significantly associated with an increase in the odds of isolating chloramphenicol-resistant NTSEC. However, this association was not seen when modelling the AMD exposures as categorical variables (models F and G), and instead there was a positive association with non-recent in-feed tetracycline and parenteral quinolone exposures, as well as a negative association with recent parenteral quinolone exposures. Quantification of AMD exposures produced a similar effect when modelling streptomycin resistance, as models in which exposures were treated continuously contained a different significant association than models in which exposures were treated categorically (Table 7), models B, D and I vs. E, F and G, respectively.
Comparison of multivariable model results [population-averaged odds ratios (OR) and 95% confidence intervals (CI)] for the outcome of chloramphenicol resistance
AMD, Antimicrobial drug; Ref., reference category.
* OR and 95% CI represent the change in odds of resistance for a median increase in animal defined daily dose for the pertinent exposure.
† OR and 95% CI have been converted from a subject-specific to a population-averaged estimate.
‡ Recent exposures are those that occurred ⩽6 days prior to sample collection.
§ Non-recent exposures are those than occurred ⩾7 days prior to sample collection.
|| Variable was not significant in the final multivariable model, and therefore results are not shown.
Comparison of multivariable model results [population-averaged odds ratios (OR) and 95% confidence intervals (CI)] for the outcome of streptomycin resistance
AMD, Antimicrobial drug; Ref., reference category.
* OR and 95% CI represent the change in odds of resistance for a median increase in animal defined daily dose for the pertinent exposure.
† OR and 95% CI have been converted from a subject-specific to a population-averaged estimate.
‡ Recent exposures are those that occurred ⩽6 days prior to sample collection.
§ Non-recent exposures are those than occurred ⩾7 days prior to sample collection.
|| Predictor was not significant in the final multivariable model, and therefore results are not shown.
Effects of different modelling strategies
While GEE, GEE/ALR and GLMM are all valid methods for regression modelling of clustered data, our results show that choice of modelling strategy can result in substantial differences with respect to which subset of AMD exposures exhibit significant association with resistance outcomes, as seen in the case of ampicillin resistance (Table 8). Models A and B indicate that recent parenteral exposures to macrolides, recent in-feed exposures to tetracyclines and recent in-feed exposures to macrolides were associated with decreased odds of recovering ampicillin-resistant NTSEC, while non-recent exposures to in-feed tetracycline were associated with increased odds. However, model I suggests that only the association between recent parenteral macrolide exposures and ampicillin resistance was statistically significant. The only difference between these models was that models A and B used GEE for clustered data, while model I used GLMM.
Comparison of multivariable model results [population-averaged odds ratios (OR) and 95% confidence intervals (CI)] for the outcome of ampicillin resistance
AMD, Antimicrobial drug; Ref., reference category.
* OR and 95% CI represent the change in odds of resistance for a median increase in animal defined daily dose for the pertinent exposure.
† OR and 95% CI have been converted from a subject-specific to a population-averaged estimate.
‡ Recent exposures are those that occurred ⩽6 days prior to sample collection.
§ Non-recent exposures are those than occurred ⩾7 days prior to sample collection.
|| Predictor was not significant in the final multivariable model, and therefore results are not shown.
Effects of clustering
In a few instances, the hierarchical level of non-independence accounted for in model specifications exerted a strong influence on model results, for example in the results of models F and G for sulfisoxazole resistance (Table 9). These models were identical except that model F included a random effect at the pen level and model G at the sample level. This change resulted in vastly different final multivariable models, with model G showing only one significant association between AMD exposure and sulfisoxazole resistance, and model F showing four significant associations.
Comparison of multivariable model results [population-averaged odds ratios (OR) and 95% confidence intervals (CI)] for the outcome of sulfisoxazole resistance
AMD, Antimicrobial drug; Ref., reference category.
* OR and 95% CI represent the change in odds of resistance for a median increase in animal defined daily dose for the pertinent exposure.
† OR and 95% CI have been converted from a subject-specific to a population-averaged estimate.
‡ Recent exposures are those that occurred ⩽6 days prior to sample collection.
§ Non-recent exposures are those than occurred ⩾7 days prior to sample collection.
|| Predictor was not significant in the final multivariable model, and therefore results are not shown.
Effect of modelling DOF as a categorical vs. continuous variable
It is often considered best practice to model inherently continuous data as continuous [Reference Royston, Altman and Sauerbrei26]. However, this may not be possible due to violation of the assumption of linearity, or categorization may be preferred in order to interpret model results in context of real-world practices [Reference Dohoo, Martin and Stryhn27]. We modelled DOF as both a continuous and categorical variable and found that this difference, while less impactful than other modelling decisions, can substantially change model results (Table 10). Models A and B for tetracycline resistance are specified identically except DOF is categorical in model A and continuous in model B. Model A identified only one statistically significant association between AMD exposures and resistance, while model B identified four (Table 10).
Comparison of multivariable model results [population-averaged odds ratios (OR) and 95% confidence intervals (CI)] for the outcome of tetracycline resistance
AMD, Antimicrobial drug; Ref., reference category.
* OR and 95% CI represent the change in odds of resistance for a median increase in animal defined daily dose for the pertinent exposure.
† OR and 95% CI have been converted from a subject-specific to a population-averaged estimate.
‡ Recent exposures are those that occurred ⩽6 days prior to sample collection.
§ Non-recent exposures are those than occurred ⩾7 days prior to sample collection.
|| Predictor was not significant in the final multivariable model, and therefore results are not shown.
Model comparisons
Despite large variability in models results, three AMU-AMR relationships were detected in all models. Recent exposures to parenteral tetracycline and sulfonamide were positively associated with odds of isolating NTSEC resistant to tetracycline and sulfisoxazole, respectively (Tables 9 and 10), while recent parenteral macrolide exposures were associated with decreased odds of ampicillin-resistant NTSEC (Table 8).
The point estimates for odds ratios for tetracycline resistance with respect to the median increase of recent parenteral tetracycline exposures ranged from 1·12 to 1·19 for models in which exposures were modeled as continuous variables, and from to 3·12 to 3·23 for models in which exposures were categorized (Table 10). Different model types (GEE/ALR, GEE, or GLMM) produced very similar odds ratios if the method for quantifying AMD exposures was the same (models A–D and J, respectively), despite the fact that only BM test results were used in GEE models and BM and DD susceptibility test results were used in GEE/ALR and GLMM models. Similar patterns emerged from model results for the association between recent parenteral sulfonamide exposure and resistance to sulfisoxazole (Table 9).
Results differ somewhat for the association between recent parenteral macrolide exposure and ampicillin resistance (Table 8). Across all models, recent parenteral macrolide administration demonstrated a negative association with ampicillin resistance, with the magnitude of this effect strongest in GEE models in which macrolide exposures were quantified as a continuous variable (models A, B and I).
Model fit statistics were not useful in identifying an optimal modelling strategy across all resistance outcomes, as no single model type consistently exhibited better goodness-of-fit (Table 11). In GEE models, model E – in which AMD exposures were categorized – had the lowest QIC for three of five resistances analysed, suggesting that categorization improved model fit. However, in GLMM models, treatment of AMD exposures as a continuous variable (i.e. model I) resulted in the lowest AIC for three of five resistances analysed. Interestingly, model G consistently exhibited lower AIC than model F, suggesting that all else being equal, sample-level clustering resulted in better model fit than pen-level clustering. Surprisingly, models H and J had a higher AIC than models F, G and I; the former models specified random effects at the pen, sample and isolate levels, while the latter models specified random effects at only one level. This finding could stem from decreased model parsimony due to inclusion of additional random effects, as marginal AIC statistics do not account for random effects and therefore their inclusion can decrease goodness-of-fit [Reference Greven and Kneib28].
Comparison of final multivariable model goodness-of-fit statistics
QIC, Quasi-Akaike's Information Criterion; AIC Akaike's Information Criterion.
* Final multivariable model results not available due to lack of model convergence.
DISCUSSION
These results highlight the inherent complexity in modelling AMU-AMR associations from real-world data. Despite this complexity, we have uncovered several important relationships between AMU and AMR outcomes in faecal NTSEC. The direction and magnitude of these relationships remained consistent despite changes in modelling technique, suggesting that they are likely true drivers of NTSEC resistance in this study population. Specifically, use of parenteral tetracyclines and sulfonamides increased the odds of recovering NTSEC resistant to tetracyclines and sulfisoxazoles, respectively, while use of parenteral macrolides decreased the odds of recovering ampicillin-resistant NTSEC. Importantly, these three relationships apply only to AMD given within 1 week prior to sample collection, suggesting that shifts in the most predominant resistance phenotypes are short-lived (i.e. reverts back to less resistant populations), or that other changes to the microbiota eventually ‘crowd out’ these resistant populations. In addition, the shortest withdrawal times for injectable tetracyclines, sulfonamides and macrolides are currently 18, 5 and 18 days, respectively. In addition, most AMDs are given early in the feeding period. Given these factors and the short-lived relationship observed in this study between exposures to these AMDs and AMR, it is questionable whether or not such use plays a significant role in dissemination of resistance through the food chain.
The three associations described above are consistent and easily interpretable; the rest of our results, however, highlight the crucial role that modelling decisions play in model results, the most obvious example being quantification of AMD exposures. In some instances, categorization of exposures produced models with wholly different subsets of significant predictors when compared to models in which exposures were not categorized. Given the highly skewed distribution of AMD exposures across pens, we believe that categorization is the more legitimate method for modelling AMU data. We attempted to evaluate linearity both by adding the quadratic term for AMD exposures and by modelling AMD exposures as quintiles and graphing resulting parameter estimates against the logit of the outcome. Unfortunately, these models either would not converge or lacked of a positive definite Hessian matrix, and therefore we were unable to formally evaluate linearity. Ideally, goodness-of-fit statistics would provide generalizable guidelines for model selection. However, our results show that the ‘best fit’ model depends greatly on the resistance outcome being modeled (Table 11).
The issue of quantification applied to confounders as well as variables of primary interest. Our analysis accounted for DOF as a known confounder of the AMU-AMR relationship [Reference Morley4]. The question is how to best account for DOF, i.e. as a categorical or continuous variable, and does this make a difference? Our results uncovered instances in which even this seemingly minor decision produced differential model results with respect to the set of significant predictors in the final model (Table 10, models A and B).
Clustering of subjects and samples occurs frequently in agriculture production settings [Reference Dohoo and Stryhn29] and modelers can choose from several techniques to account for this. Ideally, these techniques would result in similar model results; however, we have shown this is not always the case. GLMM may produce a final model with a wholly different subset of significant predictors than GEE, as seen in the final model results for tetracycline resistance in which AMDs were categorized (Table 10). While this is not surprising, it is an important finding given that model selection is often based on necessity rather than choice, for instance when data are sparse and mixed models do not converge. Furthermore, choice of modelling technique can be based on subjective criteria such as researcher preference.
The role of timing of AMD exposures in driving or mitigating resistance is made clear in these results, and in some cases, there is a strong dichotomy in this regard. Recent parenteral quinolone exposures decrease the odds of recovering chloramphenicol-resistant NTSEC, while non-recent parenteral quinolone exposures strongly increase these odds (Table 6). Multiple factors could be driving this finding, including co-selection of resistance genes on plasmids. The interlinked nature of AMR outcomes has been previously observed [Reference Wagner30], and future analyses may wish to account for this with multivariate modelling approaches. Alternatively, this finding could be explained by microbial population dynamics in which the short-term alteration of the microbiome caused by exposure to quinolone creates a long-term niche for chloramphenicol-resistant NTSEC bacterial subpopulations, perhaps due to the presence of unrelated traits. Other studies have found similar interactions between timing of antimicrobial exposure and development of resistance [Reference Lowrance31], and together these results highlight the importance of accounting for timing in models of AMU/AMR. We divided AMD exposures at 7 days relative to sampling; however, other schemes may be just as legitimate and should be based on treatment protocols employed in the study and/or target population.
The veterinary community has long recognized the difficulties of researching AMU/AMR in production settings. The analysis presented here is founded on a unique database in which many of these difficulties were overcome, including linkage of pen-level AMD exposures with susceptibility test results from pooled faecal samples. The pooling of faeces from the pen floor may have introduced sampling bias, particularly for AMDs used in low frequency, when faeces from exposed individuals were less likely to be randomly included in the 25 individual faecal samples that were pooled. However, the enrolment of 300 pens and the collection of multiple pooled samples per pen mitigated this potential bias, as well as provided sufficient power to uncover associations that may have remained hidden in smaller studies. In addition, the longitudinal sampling design allowed exploration of resistance changes over time. Despite these advancements, the complex nature of AMU and AMR renders interpretation of results challenging.
The findings presented here add another layer of complexity by demonstrating that modelling decisions greatly impact model results, especially with respect to the set of significant predictors in the final model. Which results are we to believe? For example, some AMU-AMR associations were seen across all model types, while others were present in <25% of models assessed. Does this mean that the latter associations have less real-world applicability? Hill's criteria of causation state that consistency of results across different studies, times, circumstances and populations strengthens the evidence for a causal relationship [Reference Hill32], suggesting that those AMU-AMR associations identified in all models are likely causally linked. In addition, Hill's criterion of plausibility is readily apparent in the associations between sulfonamide exposures and sulfisoxazole resistance and tetracycline use and tetracycline resistance, lending more weight to the causal nature of these relationships. The negative association between macrolide use and ampicillin resistance is less directly understood, and further research is needed to understand microbiome-wide population genetic dynamics that may be influencing this relationship.
Importantly, however, most model results were not consistent, and the question then becomes how to determine which inconsistently identified AMU-AMR associations to consider valid. From a technical perspective, the most valid model is that which most appropriately handles the data structure and provides best data fit. However, there is often a gap between theory and application, and choice of modelling technique is not always dictated solely by mathematical appropriateness. Often, factors such as study design, the ability of the model to converge given sparse data, the need for a population-averaged estimate, and even researcher preference will weigh heavily in the decision [Reference Dohoo, Martin and Stryhn33]. Therefore, while perhaps unsurprising, the finding that modelling technique can substantially change the set of significant predictors in the final model is disconcerting, particularly when the results have potential ramifications for public health and policy decisions, such as the case with AMU and AMR in livestock production. It is therefore our opinion that stakeholders in the AMU-AMR debate deserve to see results of all legitimate models, and that decisions about which model results to act on should be based on a careful assessment of the relative risk, the consistency and the plausibility of each AMU-AMR relationship being analysed.
ACKNOWLEDGEMENTS
We gratefully acknowledge Dr Trevor W. Alexander, Shaun R. Cook, and Lorna J. Selinger for their technical assistance regarding bacterial culture, isolation, and susceptibility testing at the Lethbridge Research Center, and for Dr Danielle Daignault's technical assistance regarding susceptibility testing at the Public Health Agency of Canada. We also gratefully acknowledge Chelsea Flaig at Feedlot Health Management Services Ltd for coordinating sampling. Thanks are also extended to the commercial feedlot operators for providing us access to their facilities to obtain the field isolates of Escherichia coli that were an integral part of the study.
This study was supported by grants from: the Advancing Canadian Agriculture and Agri-Food (ACAAF) Program, 2006–2011; the Alberta Beef Producers, 2006–2008 (project no. 0007-038RDB); the Canadian Cattlemen's Association, Beef Cattle Research Council, 2006–2008 (project no. BCRC 6·41); and the Surveillance Division, Center for Food-borne, Environmental Zoonotic Infectious Diseases, Public Health Agency of Canada, Canadian Integrated Program for Antimicrobial Resistance Surveillance (CIPARS).
DECLARATION OF INTEREST
Dr Booker is an owner at Feedlot Health Management Services, Okotoks, Alberta. This is a private company that provides expert consultation regarding management of feedlot cattle, including medical treatments such as use of antimicrobial drugs. They also conduct research on a fee-for-service basis regarding a variety of topics including the efficacy of different antimicrobial drugs. Dr Booker is an owner of two other veterinary companies, Lethbridge Animal Clinic (LAC), Lethbridge, Alberta and Southern Alberta Veterinary Services (SAVS), Okotoks, Alberta. LAC is a private company that provides veterinary services to beef cattle operations, which would include recommendations for medical treatments such as use of antimicrobial drugs. SAVS in a private company that sells pharmaceuticals and consumables, including prescription antimicrobial drugs, to end-users such as cow/calf producers or commercial feedlots.
