Analysis of multiple bacterial species and antibiotic classes reveals large variation in the association between seasonal antibiotic use and resistance.

Understanding how antibiotic use drives resistance is crucial for guiding effective strategies to limit the spread of resistance, but the use-resistance relationship across pathogens and antibiotics remains unclear. We applied sinusoidal models to evaluate the seasonal use-resistance relationship across 3 species (Staphylococcus aureus, Escherichia coli, and Klebsiella pneumoniae) and 5 antibiotic classes (penicillins, macrolides, quinolones, tetracyclines, and nitrofurans) in Boston, Massachusetts. Outpatient use of all 5 classes and resistance in inpatient and outpatient isolates in 9 of 15 species-antibiotic combinations showed statistically significant amplitudes of seasonality (false discovery rate (FDR) < 0.05). While seasonal peaks in use varied by class, resistance in all 9 species-antibiotic combinations peaked in the winter and spring. The correlations between seasonal use and resistance thus varied widely, with resistance to all antibiotic classes being most positively correlated with use of the winter peaking classes (penicillins and macrolides). These findings challenge the simple model of antibiotic use independently selecting for resistance and suggest that stewardship strategies will not be equally effective across all species and antibiotics. Rather, seasonal selection for resistance across multiple antibiotic classes may be dominated by use of the most highly prescribed antibiotic classes, penicillins and macrolides.

Introduction

Antibiotic resistance is a growing threat to society, with important public health [1] and economic consequences [2]. Antibiotic use is considered a primary driver of resistance not only in the pathogen targeted by the antibiotic but also in host-associated bacteria subject to “bystander selection” [3]. As such, stewardship programs to reduce overall antibiotic prescribing have become a popular strategy for broadly reducing the burden of resistance. However, the efficacy of stewardship efforts has varied widely and, in some cases, shown a limited impact on reducing rates of resistance [4,5]. These findings reflect the complexity of the antibiotic use–resistance relationship, underscoring the need to characterize this relationship across a wide range of bacterial species and antibiotics and identify factors that influence the strength of this association.

Temporal studies have shown that an association between population-level antibiotic use and resistance can be detected on rapid timescales, where seasonal fluctuations in use have been accompanied by seasonal fluctuations in resistance with up to a few months lag [6-9]. To interpret the lag between seasonal use and resistance, Blanquart and colleagues proposed a model for the relationship between short-term sinusoidal fluctuations in antibiotic use and resistance [10]. This model predicts that antibiotic use determines the rate of change of resistance, such that the derivative of the prevalence of resistance should depend on the level of use. Thus, if antibiotic use follows a sine function over a 12-month period, then peak resistance should lag peak use by a quarter period or 3 months. The lag can be shortened by including a “stabilizing force” in the model to account for forces that counteract the effect of use and drive fluctuations in resistance toward equilibrium.

Findings from previous seasonality studies have been largely consistent with this model. Studies that focused on antibiotics with wintertime peaks in use (e.g., penicillins, macrolides, and quinolones) have identified positive associations with winter/spring peaks in resistance lagged by 0 to 3 months in Streptococcus pneumoniae [6], Escherichia coli [7,8], Staphylococcus aureus [7], and Neisseria gonorrhoeae [9]. In contrast, a study from the Netherlands that analyzed antibiotics with summer/autumn peaks in use (e.g., nitrofurantoin, fosfomycin, and trimethoprim) found that resistance in E. coli and Klebsiella pneumoniae still peaked in the winter/spring and lagged use by 3 to 6 months [11], inconsistent with the Blanquart and colleagues’ model. The authors of this study attribute the longer lag time to the weaker seasonal fluctuations and lower overall rates of antibiotic use in their study population. However, it is unclear whether resistance to other antibiotics with winter peaks in use would exhibit similarly long lag times in this population or whether these findings reflect a broader phenomenon where despite different seasonal patterns of use, resistance always peaks in the winter/spring due to other ecological factors.

We aimed to characterize the seasonal relationship between antibiotic use and resistance across antibiotic classes with winter, summer, and biannual peaks in use [6,7,9,12] in Boston, Massachusetts. We studied 3 clinically relevant species—S. aureus, E. coli, and K. pneumoniae—which represent skin/nasal and gut colonizing bacteria [13,14] that cause a diversity of infections types and are subject to strong bystander selection [3]. We obtained antibiotic use data from a centralized statewide insurance claims database and resistance data from 2 major Boston area hospitals. Given the near-universal health insurance coverage in Massachusetts [15], this analysis provided a unique opportunity to characterize the antibiotic use–resistance relationship in a dataset that captures nearly all antibiotic use in a population.

Results

Seasonality in antibiotic use varies across classes

The 5 antibiotic classes included in this study each displayed statistically significant seasonal patterns of use (). Penicillins and macrolides were most frequently prescribed, with year-round averages of 4.8 and 4.1 daily claims per 10,000 people, respectively. Quinolones, tetracyclines, and nitrofurans were prescribed with year-round averages of 1.8, 1.0, and 0.54 daily claims per 10,000 people, respectively.

Seasonal patterns of antibiotic use by class.

(A) Average daily antibiotic claims per 10,000 people by calendar month in Boston, Massachusetts from 2011 to 2015. Lines indicate LOESS smoothing curves and shaded regions indicate 95% CIs. (B) Sinusoidal model fits for monthly prescribing rate. Points indicate monthly mean seasonal deviates in average daily antibiotic claims per 10,000 people by calendar month and error bars indicate the standard error of the mean. Lines indicate the point estimate for the amplitude and phase of the sinusoidal model. Shaded regions indicate the 95% CIs for the amplitude. Asterisks indicate the amplitude of seasonality is statistically significant (FDR < 0.05). FDR, false discovery rate. Underlying data are available at https://github.com/gradlab/use-resistance-seasonality/tree/master/figure_data/Fig1 [16].

Penicillins had the greatest magnitude change in prescribing rate across seasons, with the seasonal component having an amplitude of 1.1 additional daily claims per 10,000 people (peak to mean) (95% CI, 0.96 to 1.3). This was followed by macrolides (amplitude, 0.74; 95% CI, 0.59 to 0.89), quinolones (amplitude, 0.081; 95% CI, 0.044 to 0.12), nitrofurans (amplitude, 0.044; 95% CI, 0.024 to 0.064), and tetracyclines use (amplitude, 0.031; 95% CI, 0.011 to 0.052) ().

The timing of peak prescribing varied by antibiotic class (). Macrolide and penicillin use peaked in the winter, around late January (phase, 1.7 months; 95% CI, 1.3 to 2.1; note that phase is indexed to 1.0 representing January 1) and early February (phase, 2.2 months; 95% CI, 2.0 to 2.5), respectively. Tetracycline and nitrofuran use peaked in the summer, around mid-June (phase, 6.5 months; 95% CI, 4.6 to 8.5) and late August (phase, 8.8 months; 95% CI, 8.2 to 9.4), respectively. Finally, quinolone use peaked twice a year in early January and early July (phases, 1.0 (95% CI, 0.6 to 1.5) and 7.0 (95% CI, 6.6 to 7.5) months).

Seasonality in antibiotic resistance is prevalent across species and antibiotic classes

Resistance was seasonal for 9 out of 15 species–antibiotic combinations (Figs ), with statistically significant amplitudes of seasonality (false discovery rate (FDR) < 0.05) ranging from a peak log2 (minimum inhibitory concentration, MIC) increase of 0.028 to 0.063 above the yearly average (). Ciprofloxacin resistance and nitrofurantoin resistance were seasonal in all 3 species with a 12-month period (Figs ). Resistance to erythromycin in S. aureus was also seasonal with a 12-month period (amplitude, 0.048; 95% CI, 0.012 to 0.083). Conversely, tetracycline resistance was not seasonal in any of the 3 species (). The seasonal patterns of resistance to penicillin class antibiotics were variable across species. Oxacillin resistance in S. aureus was seasonal with a 12-month period (amplitude, 0.031; 95% CI, 8.8e-3 to 0.054), while both penicillin resistance in S. aureus (amplitude, 0.010; 95% CI, −6.6e-3 to 0.027) and amoxicillin/clavulanate resistance in E. coli (amplitude, 0.010; 95% CI −5.1e-4 to 0.021) and K. pneumoniae (amplitude, 0.034; 95% CI, 1.2e-3 to 0.067) did not meet our criterion for seasonality. Ampicillin resistance in E. coli was the only species–antibiotic combination with a statistically significant amplitude (0.034; 95% CI, 0.019 to 0.049) that showed a 6-month period in seasonality. However, despite having a slightly worse fit, the 12-month period model of ampicillin resistance in E. coli also indicated seasonality (amplitude, 0.041; 95% CI, 0.019 to 0.062) ().

Seasonality of antibiotic use and resistance by class in Staphylococcus aureus.

Solid lines indicate point estimates of the amplitude and phase from the best-fitting sinusoidal model of resistance (comparing 6- and 12-month periods) to each antibiotic, colored by class. Dashed gray lines indicate point estimates of the amplitude and phase of sinusoidal models for use of the corresponding antibiotic class. Shaded regions indicate the 95% CIs for the amplitude. Points indicate the monthly mean seasonal deviates in resistance, and error bars indicate the standard error of the mean. Asterisks indicate the amplitude of seasonality in resistance is statistically significant (FDR < 0.05). FDR, false discovery rate; MIC, minimum inhibitory concentration. Underlying data are available at https://github.com/gradlab/use-resistance-seasonality/tree/master/figure_data/Fig2 [16].

Amplitudes and phases of seasonality by species and antibiotic class.

(A) Comparison of amplitudes estimated from best-fitting sinusoidal models of resistance (comparing 6- and 12-month periods) across antibiotics in Staphylococcus aureus, Escherichia coli, and Klebsiella pneumoniae. Error bars indicate 95% CIs of the amplitude. Point color indicates the antibiotic class. (B) Comparison of phases of seasonality of use and resistance across species and antibiotic classes. Points indicate peak month(s) of seasonal resistance estimated by the best-fitting sinusoidal model (comparing 6- and 12-month periods) for each species–antibiotic combination, and error bars indicate the 95% CIs. Included are species–antibiotic combinations for which the amplitude of seasonality of resistance was statistically significant (FDR < 0.05). Vertical lines indicate the peak month(s) of seasonal use estimated by the best-fitting sinusoidal model (comparing 6- and 12-month periods) for each antibiotic class, and shaded regions indicate the 95% CIs. AMC, amoxicillin-clavulanate; AMP, ampicillin; CIP, ciprofloxacin; ERY, erythromycin; FDR, false discovery rate; NIT, nitrofurantoin; OXA, oxacillin; PEN, penicillin; TET, tetracycline. Underlying data are available at https://github.com/gradlab/use-resistance-seasonality/tree/master/figure_data/Fig3 [16].

Resistance peaked in the winter to spring months in all 9 seasonal species–antibiotic combinations, with peaks ranging from early December to mid-April (). Comparing across species, resistance in E. coli (median phase, 3.2; range, 2.5 to 4.5) tended to peak slightly later in the year than resistance in S. aureus (median phase, 1.6; range, 0.98 to 2.1) and K. pneumoniae (median phase, 1.1; range, 0.0 to 2.2). Peak resistance to macrolides and penicillins in S. aureus and the first peak in resistance to ampicillin in E. coli occurred around the same time of year as peak use of macrolides and penicillins, with lags of −1.2 to 2.3 months. However, resistance to ampicillin in E. coli also peaked a second time during the year in October, although use did not. Resistance to nitrofurans in all 3 species peaked between 3.2 and 5.7 months after peak nitrofuran use. Finally, resistance to quinolones in all 3 species peaked once a year about 0.79 to 2.2 months after the first peak in quinolone use.

Since the antibiotic use dataset was restricted to outpatient prescribing for people under 65 years of age, we repeated the resistance analysis on data subset to isolates from outpatients under 65 years old, representing 53%, 31%, and 47% of the E. coli, K. pneumoniae, and S. aureus isolates, respectively. Resistance to ampicillin and nitrofurantoin in E. coli and nitrofurantoin in S. aureus remained seasonal with statistically significant seasonal amplitudes (FDR < 0.05) and showed the same periods and phases of seasonality as in the analyses with the full dataset including all isolates (). Resistance to ciprofloxacin in E. coli also remained seasonal with a statistically significant amplitude in both the 12- and 6-month period model; however, a 6-month period model performed marginally better than the 12-month model (Akaike information criterion (AIC) difference, 0.3). In contrast, resistance to all antibiotics in K. pneumoniae and ciprofloxacin, erythromycin, and oxacillin resistance in S. aureus no longer met our criterion for seasonality (amplitude FDR > 0.05) after restricting our analysis to outpatients under 65 years old. For ciprofloxacin resistance with a 12-month period model and erythromycin resistance with a 6-month period model in S. aureus, the amplitude p-values were <0.05, but did not remain significant after multiple testing correction.

To further explore whether the observed seasonality of resistance could be attributable to seasonally varied sampling of patient demographics and sites of infection, we repeated the resistance analysis on the full dataset with all isolates after including covariates in our model to adjust for patient age and sex (Eq 3 in Materials and methods) and patient age, sex, and site of infection (Eq 4 in Materials and methods). Resistance remained seasonal for all 9 species–antibiotic combinations after adjusting for patient age and sex, with statistically significant amplitudes of seasonality (FDR < 0.05), although the magnitude of the estimated amplitudes decreased by 0% to 32% compared to the unadjusted model (). After adjusting for site of infection in addition to age and sex, resistance to ciprofloxacin, erythromycin, and oxacillin in S. aureus no longer met our criterion for seasonality (amplitude FDR > 0.05), and the estimated amplitudes decreased by 29% to 62% compared to the unadjusted model, while resistance in E. coli and K. pneumoniae remained seasonal (amplitude FDR < 0.05), with a 0% to 14% decrease in amplitude compared to the unadjusted model (). The amplitude p-values were <0.05 for ciprofloxacin and oxacillin resistance in S. aureus after adjusting for age, sex, and site of infection, but did not remain significant after multiple testing correction. In this model (Eq 4 in Materials and methods), the coefficient for sex was significantly positive (FDR < 0.05; where a positive β indicates that being male is associated with higher MICs) across all antibiotics in E. coli (median, 0.30; range, 0.053 to 0.54) and K. pneumoniae (median, 0.25; range, 0.18 to 0.40) and largely nonsignificant in S. aureus (median, −2.3e-3; range, −0.079 to 0.020) (). The coefficient for age was significantly positive (FDR < 0.05; where a positive β indicates that older ages are associated with higher MICs) in E. coli (median, 1.9e-3; range, 7.6e-4 to 0.013) and significantly negative in all antibiotics except ciprofloxacin in K. pneumoniae (median, −1.0e-3; range, −2.7e-3 to 8.0e-4) (). In S. aureus, the coefficient for age was significantly positive for ciprofloxacin, erythromycin, and oxacillin resistance, significantly negative for penicillin resistance, and nonsignificant for nitrofurantoin and tetracycline resistance (median for all antibiotics in S. aureus, 2.3e-3; range, −9.4e-4 to 0.016). Finally, at least one of the coefficients for site of infection was significant (FDR < 0.05) in all 15 species–antibiotic combinations, indicating that the site of infection was also an important determinant of MIC ().

Seasonal resistance is positively correlated with use of winter peaking antibiotic classes

Spearman correlation coefficients between use–resistance antibiotic pairs varied widely across antibiotics, species, and lag times, ranging from −0.91 to 0.92 (Figs ). The number of statistically significant correlations between use–resistance pairs was maximized when the lag was 0 months in S. aureus and K. pneumoniae and 1 month in E. coli. Resistance across multiple antibiotics was most positively correlated with use of winter peaking classes, penicillins, and macrolides (median Spearman’s ρ across all lags, 0.45; interquartile range (IQR), 0.058 to 0.76). Resistance to most antibiotics also showed a negative correlation with use of summer peaking classes, tetracyclines and nitrofurans (median Spearman’s ρ across all lags, −0.35; IQR, −0.63 to −0.068). Finally, resistance was not significantly correlated with use of quinolones, which peaked twice a year, for almost all antibiotics and species (median Spearman’s ρ across all lags, −0.21; IQR, −0.32 to −0.051).

Spearman correlations between seasonal use and resistance with 0 to 3 months lag in Staphylococcus aureus.

Spearman rank correlation coefficients were calculated between the monthly mean seasonal deviate in resistance (in log2(MIC)) and the monthly mean seasonal deviate in use (in average daily claims per 10,000 people) with 0, 1, 2, or 3 months lag between use and resistance, for each pairwise combination of antibiotics and classes. Error bars indicate the 95% CIs. Colors indicate the use antibiotic class. CIP, ciprofloxacin; ERY, erythromycin; Mac, macrolide; MIC, minimum inhibitory concentration; Nit, nitrofuran; NIT, nitrofurantoin; OXA, oxacillin; Pen, penicillin; Qui, quinolone; Tet, tetracycline. Underlying data are available at https://github.com/gradlab/use-resistance-seasonality/tree/master/tables/correlations.csv [16].

Discussion

Under a model in which antibiotic use drives resistance, seasonal variation in antibiotic consumption is expected to be associated with variation in population-level resistance that is in phase with or lagged up to a quarter period behind use [10]. However, we found that resistance to all antibiotics, including those with summer and biannual peaks in use, best correlated temporally with use of winter peaking antibiotics—penicillins and macrolides—at a 0- to 1-month lag.

The observed patterns of use and resistance for penicillins and macrolides were mostly consistent with previous findings [6,7,9] and with model predictions [10]. Use of penicillins and macrolides peaked in the winter, likely due to increased wintertime prescribing for respiratory infections [12]. In S. aureus, resistance to oxacillin and erythromycin peaked in the winter and was most correlated with penicillins and macrolides use with no lag. This was consistent with a study that compared seasonal macrolide use and resistance in methicillin-resistant Staphylococcus aureus (MRSA) [7] and findings in other species–antibiotic combinations that have shown winter peaks in use and resistance with little to no lag [6-8]. In E. coli, the first peak in ampicillin resistance occurred in the spring, lagging penicillins use by about 2.3 months, but also showed a second peak in resistance 6 months later, in the absence of a second peak in penicillins use. However, we note that a 12-month period model for ampicillin resistance in E. coli also met our criterion for seasonality and showed a single winter peak, which is more consistent with previous findings in E. coli [7,8].

Antibiotic classes with different seasonal patterns of use, such as nitrofurans and quinolones, showed seasonal patterns of resistance inconsistent with model predictions and not readily explained based on their patterns of use. Nitrofurans, which are almost exclusively used to treat urinary tract infections (UTIs) [17], showed summer peaks in use during the same season as peak UTI incidence [18-21]. However, resistance to nitrofurantoin in all 3 species peaked in the winter and lagged use by 3.2 to 5.7 months. This was consistent with a report that nitrofurantoin resistance lagged use by 3 to 6 months in E. coli and K. pneumoniae urinary tract isolates [11]. Quinolones, which are used to treat both respiratory infections and UTIs [22], showed both winter and summer peaks in use, but only a single winter/spring peak in ciprofloxacin resistance in all 3 species.

Several species–antibiotic combinations did not show seasonality in resistance. In each case, this may be due to a lack of association between use and resistance, a lack of strong seasonal variation in use of some antibiotic classes, or other factors, such as a signal too small to be identified or dampened by the combination of inpatient and outpatient samples in our dataset. The lack of seasonality in S. aureus resistance to penicillin may be explained by the high prevalence of penicillin resistance in S. aureus [23] (83% in this dataset), which may limit the observable effect of increased wintertime use of penicillins. For resistance to amoxicillin-clavulanate in E. coli and K. pneumoniae, the estimated amplitude of seasonality missed the 5% FDR cutoff for statistical significance, but additional data may narrow the CIs around the estimated amplitude. The absence of observed seasonality in tetracycline resistance in all 3 species is consistent with a previous study that found no significant correlation between tetracycline use and resistance in E. coli [7] and may be explained by the lack of strong seasonal variation in tetracycline use ().

Our finding that resistance to all antibiotics most correlated with use of winter peaking antibiotic classes suggests that the simple model in which use of a given antibiotic independently selects for resistance is insufficient to explain the full seasonal use–resistance landscape. Below, we discuss 4 factors that may contribute to this result, while acknowledging that there may be additional unknown seasonally varying determinants of resistance.

First, the lack of association between seasonal use and resistance in some antibiotics may be an artifact of comparing between use and resistance in overlapping but not identical populations. Data availability limited us to comparing between antibiotic use in outpatients under age 65 and resistance measured in inpatients and outpatients at 2 hospitals with patient populations that skew toward older ages. Comparisons between population-level community use and hospital resistance have frequently been utilized in previous ecological [24-28] and seasonal [6-8,11] studies, due to common difficulties in both obtaining within-hospital use data and tracking community infections that do not result in a healthcare visit. The volume of antibiotic use in the community is much greater than in the hospital [29] and thus can have a strong impact on resistance in both settings [30,31]. This impact appears to vary across species–antibiotic combinations and by demographics. For quinolones, we observed biannual peaks in use and winter peaks in resistance in all 3 pathogens; it may be that the populations receiving quinolones in the summer and winter are not sampled equally in our resistance dataset. Including only outpatients under 65 years old in our resistance analysis resulted in the loss of statistically significant seasonality in S. aureus and K. pneumoniae for some antibiotics. While this may be attributable to reduction in signal (this subset represents only 30% to 50% of the total isolates), this result could also suggest that the observed seasonal trends in resistance for some species–antibiotic combinations are disproportionately driven by older populations with infections associated with hospitalization. In contrast, the seasonal patterns of resistance among outpatients under 65 largely remained the same across antibiotics in E. coli and for nitrofurantoin in S. aureus. Thus, the disparity in community and hospital use is unlikely to fully explain the lack of association in seasonal patterns of use and resistance.

Second, the observed seasonal peaks in resistance could be driven by seasonal variation in the incidence of infection with a given species, although no specific mechanism for this association has been proposed to our knowledge. In our dataset, isolate counts peaked in the summer for all 3 species considered (), while most resistance peaked in the winter. Therefore, for seasonal variations in incidence to drive seasonality in resistance, there would need to be an inverse association between incidence and resistance; however, since there is no proposed mechanism for this association, we did not explore it statistically.

In contrast, a third mechanism by which resistance could vary seasonally is that certain patient demographics, or certain sites of infection, are associated with resistance and themselves vary seasonally. Rates of resistance have been shown to vary by age, sex, and site of infection [32-34], and the incidence of infections from these groups have been shown to vary by season in our dataset () and others [21,35]. If these factors fully explained seasonal variation in resistance, incorporating them as covariates in our seasonal regression (assuming the association was correctly specified) should have accounted for the seasonal signal and led to near-zero estimates for the sinusoidal component of the regression. However, resistance remained seasonal with amplitudes decreasing by 0% to 32%, but remaining statistically significant, in all 9 species–antibiotic combinations after accounting for age and sex (). After also accounting for site of infection, the seasonal amplitudes of resistance in 6 of those 9 combinations remained significant, decreasing by 0% to 14% compared to the unadjusted model. However, for the 3 antibiotics in S. aureus where resistance was no longer seasonal, the amplitudes decreased substantially by 29% to 62% compared to the unadjusted model. Therefore, seasonally varied sampling of isolates from different demographic groups and sites of infection likely contributes to but does not fully explain the observed seasonality in resistance.

Fourth, the winter peaks in resistance to antibiotics with different seasonal peaks in use could be explained by coselection, where use of one antibiotic can indirectly select for resistance to a second antibiotic in bacteria that are coresistant to both antibiotics [36]. Coresistance between penicillins/macrolides and other antibiotics is common across many bacterial species, including those in our study [37]. Therefore, the winter peaks in resistance to other antibiotics may be driven by coselection by winter peaking use of penicillins and macrolides. We might expect that selection by use of pencillins and macrolides dominates over selection by use of other antibiotics because they are prescribed at substantially higher rates and show greater seasonal variations in use [12]. Antibiotics with higher rates of use showed stronger correlations between seasonal use and resistance [7]. In addition, use of macrolides and penicillins have been shown to be more strongly correlated with resistance than less frequently prescribed antibiotics [24].

There were several limitations to this study. First, we measured antibiotic use in the population by the number of claims per capita, rather than daily doses, making the assumption that the average dose and duration do not vary greatly within the short timescales in which we are measuring seasonal variations in use and that there were not major selective differences between the effects of one prescription for different members of an antibiotic class. Second, we were unable to link antibiotic prescriptions to specific pathogens, and, thus, we could not assess the extent to which antibiotic resistance in a given species is attributable to antibiotic use for the treatment of infections caused by that species. However, given that bystander selection has been predicted to account for over 80% of the total antibiotic selection experienced by S. aureus, E. coli, and K. pneumoniae [3], the analysis we performed comparing total antibiotic use to resistance in each of these species may in fact yield more relevant interpretations of the use–resistance relationship [38]. Third, we determined antibiotic exposure as the total population-level use of each antibiotic class across the individuals included in our dataset, thus making the simplifying assumption that population-level use applies a uniform selective pressure. Finally, the antibiotic use and resistance datasets that were available to us for each species and antibiotic often spanned overlapping but different year ranges. Therefore, we aggregated monthly use and resistance data across years to perform our correlation analyses. We accounted for variability in use and resistance between years by adjusting for annual trends in use and resistance in our model.

In conclusion, this work contributes to describing the complexity of the antibiotic use–resistance relationship. Our finding that resistance to all antibiotics peaked in the winter/spring, regardless of patterns of use, is consistent with studies from a range of geographic scales and regions, including South Carolina [8], the United States of America overall [7,9], Israel [6], and the Netherlands [11] (note that these regions are all in high-income countries in the northern hemisphere). This suggests a general phenomenon in which selection or coselection by those antibiotics with high-volume, winter peaking use and/or other ecological factors results in wintertime peaks of resistance for all antibiotics and indicates that the simplest model of antibiotic use independently driving resistance to the same antibiotic is inadequate. We show that additional factors, such as seasonal variations in the contribution of isolates from different demographic groups, may also contribute to the observed winter peaks in resistance. This study lays the groundwork for future work to further identify and describe the factors that shape the use–resistance landscape across diverse pathogens and antibiotics, with important implications for informing and monitoring the outcome of efforts to reduce antibiotic resistance.

Materials and methods

Antibiotic use data

Outpatient antibiotic use data were obtained from the Massachusetts All Payer Claims Database [39], which covers >94% of outpatient prescriptions claims for Massachusetts residents under the age of 65 [40]. Rates of use for each antibiotic class were measured as the average daily number of antibiotic claims per 10,000 people during each calendar month from January 2011 to May 2015. These data were subset to include only individuals residing in “Boston City” census tracks, as defined by the US Census Bureau [41], to capture the antibiotic use patterns in the communities served by the hospitals in our resistance dataset. We aggregated antibiotic use data by class according to the World Health Organization’s Anatomical Therapeutic Chemical Classification System [42] (). We included 5 antibiotic classes in our analysis, which together make up 74% of the total outpatient antibiotic claims in Boston: penicillins, macrolides, quinolones, tetracyclines, and nitrofurans. Given that bystander selection likely accounts for most of the antibiotic selection experienced by S. aureus, E. coli, and K. pneumoniae [3], we included use data for all antibiotics within each class, regardless of the target pathogen for which they were prescribed.

Antibiotic resistance data

Clinical microbiology data were obtained for S. aureus, E. coli, and K. pneumoniae isolates collected at 2 tertiary care hospitals in Boston, Massachusetts: Brigham and Women’s Hospital (BWH) and Massachusetts General Hospital (MGH), from 2007 to 2019 and 2007 to 2016, respectively. Included in this analysis were all nonsurveillance isolates from inpatients and outpatients of all demographics, collected from the 5 most common sites of infection across the 3 species: blood, skin and soft tissue, abscess/fluid, respiratory tract, and urinary tract (). Isolates of the same species that were collected from the same patient within 2 weeks were assumed to represent a single infection and thus treated as a single isolate. Our final dataset comprised of 47,374 S. aureus, 130,407 E. coli, and 27,178 K. pneumoniae isolates.

Antibiotic susceptibility testing was performed on each isolate either by automated broth microdilution (Vitek 2, bioMérieux, Marcy-l’Étoile, France) or by gradient diffusion. Resulting MIC values were log2 transformed. When MICs were reported with an inequality sign, we used only the numerical value in our quantitative analyses. Due to variations in hospital testing guidelines across the years, we excluded tests on isolates that did not report an MIC value, either because a different test method was used (e.g., disk diameter) or due to missing data. We excluded years/months from our analysis for each species–antibiotic combination in each hospital where MIC values were reported for fewer than 80% of isolates or only a subset of isolate types (e.g., only testing nitrofurantoin resistance in urinary tract isolates). lists date ranges and percent resistance in each hospital, calculated as the percentage of nonsusceptible isolates out of the total number of isolates from each hospital with a reported MIC value, for each species–antibiotic combination included in our analysis. In , we determined antibiotic susceptibility by applying the Clinical & Laboratory Standards Institute 2017 breakpoints [43] to the reported MIC values of each isolate; the determined susceptibility was then adjusted based on β-lactamase screen and cefoxitin screen results if available for penicillin and oxacillin, respectively. For all other analyses other than in , we use the log2-transformed MIC as the unit of resistance. This study was approved by the Mass General Brigham Institutional Review Board (protocol number: 2016P001671).

Statistical methods

We quantified the extent of seasonality in antibiotic use and resistance by fitting the use and MIC data to a pair of mathematical models, based on a previously described method [9,10]. Both models consist of (a) a sinusoidal component to describe seasonal deviations from average year-round use and MICs; and (b) a linear component to adjust for secular trends, such as declines in use and resistance across years [23,40]. This model makes no assumptions about the underlying mechanism of resistance and is generalizable to any species and antibiotic [9,10]. As in the Olesen and colleagues’ model [9], we chose to substitute MIC for the original outcome (proportion resistance) of the Blanquart and colleagues’ model [10] to allow for detection of seasonal variations in the quantitative level of resistance, as measured by MIC, even if that variation occurs without crossing a defined breakpoint for antibiotic susceptibility.

To describe the seasonality of use, monthly claims data for each antibiotic class were fit to where u is the mean daily reported claims per 10,000 people during calendar month t, A is the amplitude of use seasonality, ω is the frequency of seasonality, where , P is the phase of use seasonality, and B and C are the within-year slope and intercept terms. To describe the seasonality of resistance, MICs for each isolate were fit to where y is the log2-transformed MIC and t is the calendar month of collection of the i isolate, A is the amplitude of resistance seasonality, ω is the frequency of seasonality, where , P is the phase of resistance seasonality, and B and C are the within hospital/year slope and intercept terms.

We further fit the MIC data to 2 additional models to account for the effect of seasonally varied sampling of patient demographics and sites of infection on the observed seasonality of resistance. To describe the seasonality of resistance while adjusting for patient demographics, age and sex, we fit the MIC data to where β and β are the coefficients for age and sex, respectively, a is the age (in years) of the patient from which the i isolate is collected, and is an indicator variable for whether the patient is male. To describe the seasonality of resistance while adjusting for patient demographics and site of infection, we fit the MIC data to where β, β, β, and β are the coefficients for blood, respiratory tract, skin and soft tissue, and abscess/fluid, respectively, and , , , and are the indicator variables for whether the i isolate was collected from the blood, respiratory tract, skin and soft tissue, or abscess/fluid, respectively.

The amplitude, phase, slope, intercept, and demographic coefficient terms in each model were estimated by nonlinear regression, using the nls function in R (version 3.6.2) [44]. We examined periods of both 12 and 6 months to account for annual or biannual cycles in use and resistance. We justified using these fixed periods by performing a wavelet analysis, using the WaveletComp package [45] in R, on the raw antibiotic use data to show that the dominant periods of variations in use across the included years are at 12 and 6 months (). To determine whether to use a 12- or 6-month period for each species–antibiotic combination, we performed model comparisons using the AIC and used the period that resulted in the lower AIC (). We determined that there was seasonality in use or resistance if the amplitude was statistically significant after accounting for multiple comparisons by applying the Benjamini–Hochberg correction with a 5% FDR.

We quantified the association between the observed seasonal patterns of use and resistance using Spearman rank correlations. To eliminate the impact of annual trends, we calculated correlations between the average monthly seasonal deviates in use and resistance, aggregated across years, rather than the raw use and MIC data. We define a “seasonal deviate” as the deviation in use or MIC at a given time of year from the year-round average, which we estimated by the linear component of the models. Seasonal deviates in use for each year and month were calculated as where is the seasonal deviate of the mean reported daily claims per 10,000 people during calendar month t, u is the mean reported daily claims per 10,000 people during calendar month t, and and are the within-year slope and intercept terms estimated from the model fit (Eq 1) for the corresponding year. For resistance, we calculate the seasonal deviate for each isolate as where is the seasonal deviate of the log2-transformed MIC of the i isolate, y is the log2-transformed MIC of the i isolate, t is the calendar month of collection of the i isolate, and and are the hospital/year slope and intercept estimated from the model fit (Eq 2) for the hospital and year that the i isolate was collected in.

Since the working model for the use–resistance relationship predicts that seasonal fluctuations in resistance can lag use by up to 3 months [10], we calculated Spearman correlations between use and resistance seasonal deviates with no lag and lags of 1, 2, and 3 months. In addition, because we observed some seasonal patterns of resistance that better aligned with use of noncognate antibiotic classes, we calculated use–resistance correlations between each pairwise combination of target antibiotics and use classes. We only included use–resistance pairs in this analysis for which both use and resistance met our criterion for seasonality.

All analyses were performed in R version 3.6.2 [44]. Data and code are available at https://github.com/gradlab/use-resistance-seasonality [16].

Seasonality of antibiotic use and resistance by class in Escherichia coli.

Solid lines indicate point estimates of the amplitude and phase from the best-fitting sinusoidal model of resistance (comparing 6- and 12-month periods) to each antibiotic, colored by class. Dashed gray lines indicate point estimates of the amplitude and phase from sinusoidal models of use of the corresponding antibiotic class. Shaded regions indicate the 95% CIs for the amplitude. Points indicate the monthly mean seasonal deviates in resistance, and error bars indicate the standard error of the mean. Asterisks indicate the amplitude of seasonality in resistance is statistically significant (FDR < 0.05). Amox/Clav, amoxicillin-clavulanate; FDR, false discovery rateS. Underlying data are available at https://github.com/gradlab/use-resistance-seasonality/tree/master/figure_data/S1_Fig [16].

(TIFF)

Click here for additional data file.

Seasonality of antibiotic use and resistance by class in Klebsiella pneumoniae.

Solid lines indicate point estimates of the amplitude and phase from the best-fitting sinusoidal model of resistance (comparing 6- and 12-month periods) to each antibiotic, colored by class. Dashed gray lines indicate point estimates of the amplitude and phase from sinusoidal models of use of the corresponding antibiotic class. Shaded regions indicate the 95% CIs for the amplitude. Points indicate the monthly mean seasonal deviates in resistance, and error bars indicate the standard error of the mean. Asterisks indicate the amplitude of seasonality in resistance is statistically significant (FDR < 0.05). Amox/Clav, amoxicillin-clavulanate; FDR, false discovery rate. Underlying data are available at https://github.com/gradlab/use-resistance-seasonality/tree/master/figure_data/S2_Fig [16].

(TIFF)

Click here for additional data file.

Seasonality of penicillins use and ampicillin resistance in Escherichia coli with a 12-month period model.

Solid line indicates point estimates of the amplitude and phase from a 12-month period sinusoidal model of resistance to ampicillin in E. coli. Dashed gray line indicates point estimates of the amplitude and phase from a 12-month period sinusoidal model of use of penicillin class antibiotics. Shaded regions indicate the 95% CIs for the amplitude. Points indicate the monthly mean seasonal deviates in resistance, and error bars indicate the standard error of the mean. Asterisk indicates the amplitude of seasonality in resistance is statistically significant (FDR < 0.05). FDR, false discovery rate. Underlying data are available at https://github.com/gradlab/use-resistance-seasonality/tree/master/figure_data/S3_Fig [16].

(TIFF)

Click here for additional data file.

Spearman correlations between seasonal use and resistance with 0 to 3 months lag in (A) E. coli and (B) K. pneumoniae. Spearman rank correlation coefficients were calculated between the monthly mean seasonal deviate in resistance (in log2 (MIC)) and the monthly mean seasonal deviate in use (in average daily claims per 10,000 people) with 0, 1, 2, or 3 months lag between use and resistance, for each pairwise combination of antibiotics and classes. Error bars indicate the 95% CIs. Colors indicate the use antibiotic class. AMP, ampicillin; CIP, ciprofloxacin; Mac, macrolide; MIC, minimum inhibitory concentration; Nit, nitrofuran; NIT, nitrofurantoin; Pen, penicillin; Qui, quinolone; Tet, tetracycline. Underlying data are available at https://github.com/gradlab/use-resistance-seasonality/tree/master/tables/correlations.csv [16].

(TIF)

Click here for additional data file.

Seasonal incidence of infection by demographic group or site of infection.

Bars show the total number of isolates by month included in the resistance dataset for each species, colored by (A) age group, (B) sex, and (C) site of infection. NOS, not otherwise specified. Underlying data are available at https://github.com/gradlab/use-resistance-seasonality/tree/master/figure_data/S5_Fig [16].

(TIF)

Click here for additional data file.

Wavelet analysis of antibiotic use by class.

Dotted lines show 12-month (upper line) and 6-month (lower line) periods. Solid lines indicate regions where the amplitude p-value is less than 0.05. Shaded areas indicate the “cone of influence” where edge effects are important. Underlying data are available at https://github.com/gradlab/use-resistance-seasonality/tree/master/figure_data/S6_Fig [16].

(TIFF)

Click here for additional data file.

Amplitudes and phases of seasonality of resistance in outpatients under 65 years old.

Amplitudes and phases were estimated from the best-fitting sinusoidal model of resistance (comparing 6- and 12-month periods) for each species–antibiotic combination, using data that was subset to include only outpatients under age 65. Asterisks indicate that the amplitude is significant after Benjamini–Hochberg multiple testing correction (FDR < 0.05). AMC, amoxicillin-clavulanate; AMP, ampicillin; CIP, ciprofloxacin; ERY, erythromycin; FDR, false discovery rate; NIT, nitrofurantoin; OXA, oxacillin; PEN, penicillin; TET, tetracycline.

(DOCX)

Click here for additional data file.

Age, sex, and site of infection coefficients for adjusted sinusoidal model of seasonal resistance.

Demographic and site of infection coefficients were estimated from the sinusoidal model of resistance that was adjusted for patient age, sex, and site of infection (Eq 4 in Materials and methods). In parentheses are the 95% CIs on the coefficient estimates. Asterisks indicate that the coefficient is significant after Benjamini–Hochberg multiple testing correction (FDR < 0.05). β, coefficient for patient age; β, coefficient for patient sex; β, coefficient for if the isolate was a blood isolate; β, coefficient for if the isolate was a respiratory tract isolate; β, coefficient for if the isolate was a skin/soft tissue isolate; β, coefficient for if the isolate was an abscess/fluid isolate; AMC, amoxicillin-clavulanate; AMP, ampicillin; CIP, ciprofloxacin; ERY, erythromycin; FDR, false discovery rate; NIT, nitrofurantoin; OXA, oxacillin; PEN, penicillin; TET, tetracycline.

(DOCX)

Click here for additional data file.

Percent of claims by individual antibiotics within each class.

Antibiotic claims data for Boston, Massachusetts residents was obtained from the Massachusetts All Payers Claims Database and subset by antibiotic class according to the following WHO ATC codes [42]: penicillins, J01C; macrolides, J01F; quinolones, J01M; tetracyclines, J01A; and nitrofurans, J01XE. All antibiotics under each ATC header that were present in this dataset were included in the analysis. The percentage of the total claims within each antibiotic class made up by each individual antibiotic is listed in the last column of the table.

(DOCX)

Click here for additional data file.

Total number of isolates by demographics.

Table of the total number of isolates of each species within each clinical or demographic category that was used in this analysis. In parentheses is the percent of the total number of isolates for that species. BWH, Brigham and Women’s Hospital; MGH, Massachusetts General Hospital; NOS, not otherwise specified.

(DOCX)

Click here for additional data file.

Antibiotics included in analysis and percent resistance by hospital.

Percent resistance was calculated as the percentage of nonsusceptible isolates out of the total number of isolates from each hospital with a reported MIC for that antibiotic. BWH, Brigham and Women’s Hospital; MGH, Massachusetts General Hospital; MIC, minimum inhibitory concentration.

(DOCX)

Click here for additional data file.

Comparison of the AIC values between 6- and 12-month period models for antibiotic use.

In parentheses is the difference in AIC from the model with the lower AIC. AIC, Akaike information criterion.

(DOCX)

Click here for additional data file.

Comparison of the AIC values between 6- and 12-month period models for antibiotic resistance.

In parentheses is the difference in AIC from the model with the lower AIC. AIC, Akaike information criterion.

(DOCX)

Click here for additional data file.

Table 1

Comparison of estimated amplitudes of seasonality across 3 sinusoidal models for resistance.

Species	Abx	Period (months)	Amplitude (95% CI)
			Model A (unadjusted)	Model (adjusted for age and sex)	Model C (adjusted for age, sex, and site of infection)
E. coli	AMC	6	0.01 (−5.1e-04, 0.021)	0.011 (4.3e-04, 0.022)	0.011 (6.3e-04, 0.022)
E. coli	AMP	6	0.034 (0.019, 0.049)*	0.034 (0.019, 0.048)*	0.034 (0.019, 0.048)*
E. coli	CIP	12	0.051 (0.031, 0.072)*	0.045 (0.026, 0.065)*	0.044 (0.025, 0.064)*
E. coli	NIT	12	0.028 (0.02, 0.037)*	0.028 (0.02, 0.036)*	0.028 (0.02, 0.036)*
E. coli	TET	6	0.013 (−3.6e-03, 0.029)	0.013 (−3.4e-03, 0.029)	0.013 (−3.3e-03, 0.029)
K. pneumoniae	AMC	12	0.034 (1.2e-03, 0.067)	0.034 (1.3e-03, 0.066)	0.032 (−2.1e-04, 0.065)
K. pneumoniae	CIP	12	0.053 (0.023, 0.083)*	0.05 (0.02, 0.081)*	0.048 (0.018, 0.078)*
K. pneumoniae	NIT	12	0.035 (9.6e-03, 0.061)*	0.034 (7.9e-03, 0.06)*	0.034 (7.6e-03, 0.06)*
K. pneumoniae	TET	6	0.021 (−9.8e-03, 0.051)	0.019 (−0.011, 0.049)	0.019 (−0.012, 0.049)
S. aureus	CIP	12	0.063 (0.034, 0.093)*	0.043 (0.016, 0.07)*	0.024 (5.0e-04, 0.048)
S. aureus	ERY	12	0.048 (0.012, 0.083)*	0.042 (7.7e-03, 0.077)*	0.032 (−1.0e-03, 0.065)
S. aureus	NIT	12	0.033 (0.022, 0.043)*	0.034 (0.023, 0.044)*	0.038 (0.027, 0.048)*
S. aureus	OXA	12	0.031 (8.8e-03, 0.054)*	0.027 (5.5e-03, 0.049)*	0.023 (8.2e-04, 0.044)
S. aureus	PEN	6	0.01 (−6.6e-03, 0.027)	9.5e-03 (−7.3e-03, 0.026)	9.7e-03 (−7.1e-03, 0.027)
S. aureus	TET	12	0.013 (−7.2e-03, 0.033)	0.012 (−7.3e-03, 0.032)	0.013 (−6.9e-03, 0.033)

Model A does not adjust for patient demographics or site of infection, Model B adjusts for patient age and sex, and Model C adjusts for patient age, sex, and site of infection. In parentheses are the 95% CIs on the amplitude estimates. Asterisks indicate that the amplitude is significant after Benjamini–Hochberg multiple testing correction (FDR < 0.05).

AMC, amoxicillin-clavulanate; AMP, ampicillin; CIP, ciprofloxacin; ERY, erythromycin; FDR, false discovery rate; NIT, nitrofurantoin; OXA, oxacillin; PEN, penicillin; TET, tetracycline.