Introduction
Antimicrobial-resistant (AMR) bacterial infections constitute a major public health crisis, causing an estimated 700,000 deaths annually worldwide, with this number projected to rise to 10 million by 2050. [1] Although methicillin-resistant Staphylococcus aureus (S. aureus) (MRSA) continues to be prevalent in the United States (US), its rates have declined in several patient populations, with a dramatic reduction observed among patients in Veterans Health Administration (VHA) acute care and nursing home facilities. [Reference Evans2–Reference Hidron8] Between 2005 and 2017, MRSA rates in VHA inpatients decreased by 55%, driven by the nationwide implementation of the MRSA bundle in 2007. [Reference Jain3, Reference Jones5] A recent study found similar trends in the US veteran outpatient population, with peak MRSA prevalence of 53.6% in 2010, dropping to 38.8% in 2019. [Reference Carrel4]
Despite its overall decline, increased resistance to additional antibiotics in MRSA, such as tetracyclines and trimethoprim–sulfamethoxazole (TMP–SMX), has heightened concerns of multi-drug-resistant strains both domestically and globally. [Reference Lynch and Zhanel9–Reference Diekema11] In the US veteran outpatient population, tetracycline- and TMP–SMX-resistant community-associated MRSA (CA-MRSA) infections have increased from 2010 to 2019 by 9.2% and 6.6%, respectively. [Reference Carrel4] Increases in resistance to non-beta-lactam antibiotics among MRSA have also been seen in the US inpatient and paediatric populations, most notably for clindamycin, mupirocin, and TMP–SMX. [Reference Carrel4, Reference Sutter12–Reference Sader16] This is particularly concerning, as many clinicians prescribe non-beta-lactam antibiotics for empirical coverage of MRSA in infections that do not warrant hospitalization.
Current research on S. aureus resistance rates to non-beta-lactam antibiotics is limited, largely due to data availability. Nevertheless, geographic variability in temporal trends has been observed across the US for both the prevalence of MRSA and S. aureus resistant to non-beta-lactam antibiotics. Regional differences have been reported, with the greatest reductions in MRSA percentage and increases in tetracycline and TMP–SMX resistance found in the south and northeast. [Reference Carrel4, Reference Ham13] Finer-scale analysis revealed that while regional tetracycline and TMP–SMX resistance rates may be rising, these trends often vary among counties within the region. [Reference Carrel4] This intra-regional variation suggests that factors beyond patient-level clinical characteristics, such as ecological and socioeconomic differences, shape S. aureus resistance patterns, emphasizing the need for analysis at finer spatial scales. Leveraging more spatially representative data, such as VHA electronic health records, can help address this need, providing insights into how geographically variable factors may influence resistance rates.
The aim of this study was to examine the spatiotemporal patterns of the relative risks of tetracycline- and TMP–SMX-resistant MRSA among VHA outpatients in eastern US commuting zones (CZs) from 2018 to 2022. By (1) identifying hotspots of resistant MRSA relative risk and (2) areas where local trends are increasing at a faster rate relative to the overall trend, this study can generate hypotheses about place-based sociodemographic and ecological factors that may be influencing the timing and rate of increase in these areas.
Methods
Study population and data source
We obtained data from the VHA’s electronic health record repository, VA Corporate Data Warehouse, for VHA outpatients at 127 medical centres and over 1,400 outreach clinics with positive S. aureus cultures from 1 January 2010 to 30 September 2023. The data included patient-level residential information and were linked to microbiology results containing antimicrobial susceptibility data for several antibiotic classes, including tetracyclines and TMP–SMX. Approximately 90% of facilities had on-site microbiology laboratories, required to perform routine quality control and use Food and Drug Administration (FDA)–approved methods in accordance with VHA-designated accreditation organizations. [17, 18]
Inclusion criteria were applied for patients aged 18 years or older, residing in contiguous states of the US or the District of Columbia (DC), with address-level geographic data, and at least one recorded antimicrobial susceptibility result. The data were further filtered to exclude observations with missing MRSA classification and retain only the first observation per 30-day period per patient. The data were stratified into MRSA and methicillin-susceptible S. aureus (MSSA) based on resistance, including intermediate resistance, to beta-lactam agents with antistaphylococcal activity, specifically cephalosporins, antistaphylococcal penicillin, or carbapenem. To examine overall trends in resistance, the monthly prevalence of tetracycline and TMP–SMX resistance across the entire US from 2010 to 2023 was calculated and plotted using the ggplot2 version 3.5.1 package in R version 4.4.1.
For the Bayesian analyses, MRSA isolates from 2018 to 2022 were aggregated to annual total counts and counts of cultures resistant to tetracyclines and TMP–SMX at the CZ level to prioritize more recent trends. Focusing on this time period allowed the model to capture recent net changes in relative risk while minimizing the influence of any larger prevalence shifts from earlier in the study period. CZs use commuting data to capture local economics, linking areas to their nearest economic centres. [Reference Tolbert and Killian19] Due to insufficient sample sizes in many western CZs and prior evidence of greater increases in tetracycline and TMP–SMX resistance in the south and northeast, Bayesian analyses were restricted to the eastern US. [Reference Carrel4, Reference Ham13] Of the 593 CZs defined in the 2020 delineation, 331 of the eastern CZs were included in the Bayesian model. [Reference Fowler20]
Statistical analysis
To account for variations across space and time, a hierarchical spatiotemporal Bayesian Poisson model with a log-link function was applied to cases of tetracycline- and TMP–SMX-resistant MRSA within each CZ i and year j:
$$ \log \left({\mu}_{ij}\right)=\log \left({\mathrm{offset}}_{ij}\right)+\alpha +{s}_i+{\beta}_j+{\gamma}_{i,j} $$
The offset term was the log of the expected count, calculated by multiplying the proportion of tetracycline- or TMP–SMX-resistant cases by the adjusted number of MRSA cases within each CZ and year. The intercept term (
$ \alpha $
) represents the mean log risk for the response from 2018 and 2022. The spatially correlated random effect (
$ {s}_i $
) follows a Besag–York–Mollié (BYM2) model, including an intrinsic conditional autoregressive (ICAR) component and an unstructured random effect, with the overall variance parameter assigned a weakly informative N(0,1) prior. [Reference Besag, York and Mollié21, Reference Simpson22] A temporal effect term (
$ {\beta}_j $
) estimates the average time trend across all CZs, with year coded as an indexed effect, while the space–time interaction (
$ {\gamma}_{i,j} $
) estimates any deviation of a CZ’s temporal trend from this average. A weakly informative N(0, 10) prior was assigned to the intercept, temporal effect, and space–time interaction term.
The model was run in Stan version 2.32.2, utilizing the RStan version 2.32.6 R package. [23] Three chains were run for 50,000 iterations, removing the first 25,000 iterations as burn-in and keeping every fifth sample for a total of 15,000 posterior samples. To assess model convergence, the R-hat and number of effective samples (
$ {\mathrm{N}}_{\mathrm{eff}} $
) were reviewed. The linear predictor (
$ {\eta}_{ij} $
) was generated, representing the log-relative risk of each CZ and year. Posterior samples for the linear predictor, temporal effect, and space–time interaction terms were extracted, and the mean yearly, CZ, and CZ-year absolute relative risks were calculated. The mean yearly differences in the space–time interaction term were calculated to infer CZ-specific temporal trends.
All CZs were categorized based on their relative risk and temporal trend using a two-step classification procedure. [Reference Li24] A CZ was classified as high risk if the posterior probability of its relative risk exceeding one was >0.8, moderate risk if between 0.2 and 0.8, and low risk if <0.2. Temporal trends were classified as increasing if the posterior probability of a positive yearly change in the interaction effect was >0.8, stationary if between 0.2 and 0.8, and decreasing if <0.2. An alternative procedure using more lenient cut-offs for temporal trends (>0.6 for increasing; <0.4 for decreasing) was applied to assess whether less stringent criteria improved the sensitivity of CZ classification. The CZ-specific mean relative risks, temporal trends, and classifications were mapped in ArcGIS Pro with the 20 most populated cities based on the 2020 US Census (Esri, Redlands, CA). The yearly absolute relative risks were plotted for the CZs with the three highest posterior probabilities for increasing trend and compared to the average. Model and classification details, additional Stan settings, and R code are provided in the Supplemental Methods.
To evaluate model robustness, we conducted several sensitivity analyses. Potential confounding between the temporal effect term and the space–time interaction term was assessed by comparing the posterior distributions of the temporal effect from the original model and a model fitted without the space–time interaction. The original model specified weakly informative N(0, 10) priors for both the temporal effect and space–time interaction terms. To assess sensitivity to these priors, the model was re-estimated using a non-centred first-order random walk prior (RW1) for the temporal effect and a structured non-centred RW1 prior over time within each CZ for the space–time interaction term. Robustness was evaluated by examining correlations between posterior mean relative risks and the space–time interaction estimates, comparing posterior differences in temporal effects, and reassessing the two-step CZ risk and trend classification relative to the original model. Ethical approval was given by the institutional review board at the University of Iowa and the Research and Development Committee at the Iowa City Veterans Affairs Health Care System with a waiver of informed consent because the study was a retrospective analysis of health records with no direct contact with patients. This study followed the Strengthening the Reporting of Observational Studies in Epidemiology (STROBE) reporting guideline.
Results
From the 558,737 S. aureus culture results compiled from the VHA’s electronic health record repository, 469,624 were included in the analysis. The first result per patient per month was retained (excluding n = 75,842), and patients were excluded if aged less than 18 years (n = 3), resided outside the 48 contiguous states or DC (n = 5,277), were missing geographic data and all antibiotic resistance results (n = 1,517), or lacked MRSA classification (n = 6,474) (Supplementary Table S1).
Prevalence of tetracycline (Figure 1a) and TMP–SMX resistance (Figure 1b) from 2010 to 2023 was plotted using 205,887 MRSA isolates from 149,422 outpatients and 263,737 MSSA isolates from 199,871 outpatients across the entire United States. Tetracycline resistance in MRSA increased from a 2010 low of 4.57% to a 2023 high of 20.98%. TMP–SMX resistance in MRSA increased from a 2010 low of 2.48% to a 2023 high of 11.77%. This trend was not seen in MSSA isolates. Tetracycline resistance in MSSA only varied by approximately 5.66%, with a 2010 low of 3.78% and a 2023 high of 8.35%. TMP–SMX resistance in MSSA had a smaller range of 1.93%, with a 2010 low of 0.83% and a 2023 high of 2.76%.
Figure 1. Monthly prevalence from January 2010 to September 2023 for VHA outpatients for (a) tetracycline-resistant MRSA and MSSA and (b) TMP–SMX-resistant MRSA and MSSA.
A hierarchical Bayesian model was fit using aggregated data of 58,901 MRSA isolates collected from 45,325 outpatients between 2018 and 2022. The estimated mean absolute relative risk of tetracycline-resistant MRSA rose from 0.77 (95% confidence interval (95% CI) = 0.58, 1.26) in 2018 to 1.10 (95% CI = 0.85, 1.66) in 2022 (Supplementary Figure S1A). The estimated mean absolute relative risk of TMP–SMX-resistant MRSA also increased from 0.68 (95% CI = 0.51, 1.20) in 2018 to 0.96 (95% CI = 0.73, 1.59) in 2022 (Supplementary Figure S1B).
The CZ-specific average absolute relative risk for tetracycline (Figure 2a) and TMP–SMX (Figure 2b) resistance in MRSA was geographically heterogeneous. For both resistance types, higher-risk CZs were distributed throughout the study area, in both densely populated metropolitan and rural areas. High-risk hot spots for tetracycline-resistant MRSA were evident along the eastern states and near central metropolitan cities. High-risk CZs for TMP–SMX-resistant MRSA were located near metropolitan cities, with a notable cluster in southern Florida. To quantify the certainty in the mean relative risk estimates, the posterior probabilities of the relative risks being greater than 1 were calculated. High-risk CZs with posterior probabilities greater than 0.8 are visualized with a hatched pattern in Figure 2. For tetracycline-resistant MRSA, CZs categorized as high risk were concentrated along the east coast states, with two centrally located CZs containing or adjacent to metropolitan cities. CZs categorized as high risk for TMP–SMX-resistant MRSA were largely located in South Florida, with additional high-risk CZs containing major metropolitan cities. All CZ-specific posterior probabilities for relative risk being greater than 1 can be visualized in Supplementary Figure S2.
Figure 2. Average absolute relative risk in eastern US commuting zones between 2018 and 2022 of (a) tetracycline resistance and (b) TMP–SMX resistance in MRSA. Hatched commuting zones are classified as high risk based on the first step of the classification procedure described in the methods section, where
$ p\left(\exp \left({\eta}_{s_i}\right)>1| data\right)>0.8 $
.
The temporal trend of each CZ over the study period is represented by the mean yearly change in the spatiotemporal interaction effect. Positive values suggest increasing trends, while negative values indicate decreasing trends. These spatiotemporal trends show high heterogeneity for both tetracycline-resistant (Figure 3a) and TMP–SMX-resistant MRSA (Figure 3b).
Figure 3. Spatiotemporal trends of eastern US commuting zones between 2018 and 2022. The mean yearly change in the interaction effect for (a) tetracycline resistance and (b) TMP–SMX resistance in MRSA indicates the presence of an overall increasing (positive) or decreasing (negative) trend. The hatched commuting zone in Figure 5a is classified as having an increasing trend based on the second step of the classification procedure described in the methods section, where
$ p({\beta}_{s_i,{t}_j-}{\beta}_{s_i,{t}_{j-1}}>0|data)>0.8 $
.
To assess the certainty of the spatiotemporal trends, posterior probabilities of the yearly change in the interaction effect being greater than zero were calculated. Visualization of these probabilities was stratified by risk group, determined by the first step of the classification procedure outlined in the methods (Figure 4). Low- and moderate-risk CZs for tetracycline-resistant (Figure 4a,b) and TMP–SMX-resistant (Figure 4d,e) MRSA show posterior probabilities primarily between 0.2 and 0.8. One moderate-risk CZ for tetracycline-resistant MRSA, located along the North Carolina and Virginia border, had a posterior probability of greater than 0.8, indicating a likely increasing trend. CZs classified as high risk for tetracycline-resistant (Figure 4c) and TMP–SMX-resistant (Figure 4f) MRSA had posterior probabilities between 0.4 and 0.6, suggesting stable trends.
Figure 4. Posterior probabilities of increasing spatiotemporal trends in eastern US commuting zones between 2018 and 2022 stratified by risk group. The darker the colour, the greater the likelihood of that CZ having a positive yearly change in the spatiotemporal interaction term, indicating an overall increasing trend. (a) Low-risk, (b) moderate-risk, and (c) high-risk commuting zones for tetracycline resistance in MRSA. (d) Low-risk, (e) moderate-risk, and (f) high-risk commuting zones for TMP–SMX resistance in MRSA.
CZ risk and trend were categorized based on the classification procedure (Table 1). For tetracycline-resistant MRSA, 9 (3%) were categorized as high risk, 189 (57%) as moderate risk, and 132 (40%) as low risk. Nearly all CZs (99.7%) were categorized as having a stationary trend, with one moderate-risk CZ showing an increasing trend. For TMP–SMX-resistant MRSA, 13 (4%) CZs were categorized as high risk, 140 (42%) as moderate risk, and 178 (54%) as low risk. All CZs were categorized as having a stationary temporal trend. Alternative classifications using more lenient cut-offs for temporal trends (<0.4 and > 0.6) altered the categorization of moderate- and low-risk CZs. However, no high-risk CZs with increasing or decreasing trends were identified, further supporting the finding that high-risk CZs remained stable throughout the study period (Supplementary Table S2).
Table 1. Two-step classification of commuting zones based on their absolute relative risk and spatiotemporal trend
Although categorization facilitates straightforward interpretation, aggregating to discrete categories obscures finer differences, particularly in temporal trends. To capture this variation, the continuous distribution of the posterior probabilities was plotted for tetracycline-resistant (Figure 5a) and TMP–SMX-resistant (Figure 6a) MRSA. Low-risk CZs for both resistance types were concentrated around a posterior probability of a positive change in the spatiotemporal interaction term of 0.5, consistent with a stationary trend. In contrast, moderate-risk CZs, as well as low-risk CZs near the 0.2 threshold, demonstrated greater variability in their posterior probabilities of increasing trends. High-risk CZs also remained near a posterior probability of 0.5, suggesting a likely stationary pattern.
The yearly mean absolute relative risks for the three CZs with the highest posterior probability of increasing trends were plotted. For tetracycline-resistant MRSA, these included CZ 401 (Virginia–North Carolina), CZ 380 (New York), and CZ 248 (Maine) (Figure 5b). All showed increasing trends from 2018 to 2022, with a sharp increase between 2020 and 2021. CZs 401 and 248, categorized as moderate risk, showed trends rising above the average, whereas CZ 380 remained mostly below the average, consistent with its categorization of low risk despite an increasing trend.
Figure 5. (a) Continuous joint distribution of commuting zone (CZ) mean posterior probabilities of relative risk for tetracycline resistance in MRSA being greater than one (x-axis) and the change in the spatiotemporal interaction term being greater than zero (y-axis). (b) Yearly mean absolute relative risks of tetracycline-resistant MRSA for CZ 401 (Virginia–North Carolina), CZ 380 (New York), and CZ 248 (Maine), plotted alongside the average temporal trend (dotted line) of the study area.
For TMP–SMX-resistant MRSA, the three CZs with the highest posterior probabilities of increasing trends included CZ 163 (Indiana), CZ 398 (North Carolina), and CZ 464 (Pennsylvania) (Figure 6b). All exhibited increasing trends from 2018 to 2022, though the rate and timing of change varied. CZ 163 exhibited a steady increase, remaining below the average trend until an almost fivefold increase between 2021 and 2022. CZ 398 increased at a gradual rate, with a steep increase between 2019 and 2020, though it remained below the average trend, consistent with its low-risk categorization. In contrast, CZ 464 showed a fluctuating trend, with sharp increases elevating it above the average trend in 2020 and 2022, and a decrease between 2020 and 2021, temporarily bringing it below the average.
Figure 6. (a) Continuous joint distribution of commuting zone (CZ) mean posterior probabilities of relative risk for TMP–SMX resistance in MRSA being greater than one (x-axis) and the change in the spatiotemporal interaction term being greater than zero (y-axis). (b) Yearly mean absolute relative risks of TMP–SMX-resistant MRSA for CZ 163 (Indiana), CZ 398 (North Carolina), and CZ 464 (Pennsylvania), plotted alongside the average temporal trend (dotted line) of the study area.
Sensitivity analyses indicated that the primary findings were robust to alternative model specifications. Posterior distributions of the main temporal effect were nearly identical between models fitted with and without the space–time interaction, with posterior mean differences centred near zero and 95% CIs overlapping zero across all years, indicating minimal confounding between terms (Supplementary Table S3). Re-estimation of the models using non-centred RW1 priors resulted in CZ-level posterior mean relative risks that were highly correlated with those from the original model (Spearman’s ρ = 0.87–0.92; Pearson’s r = 0.85–0.88), with slightly lower but positive correlations for the space–time interaction term (Spearman’s ρ = 0.74–0.81; Pearson’s r = 0.73–0.78). Posterior differences in the temporal effect under RW1 prior were small and again included zero in the 95% CIs across all years (Supplementary Table S4). Although reclassification resulted in more CZs categorized as high and low risk, none were reclassified as having an increasing or decreasing temporal trend, and the overall pattern of temporal stability remained unchanged (Supplementary Table S5).
Discussion
Electronic health record data from a nationwide sample of VHA outpatients demonstrated that the prevalence of tetracycline- and TMP–SMX-resistant MRSA increased by 16.41% and 9.29%, respectively, from 2010 to 2023, consistent with previous studies. [Reference Carrel4, Reference Ham13, Reference Khamash15] This rise was not mirrored in MSSA, suggesting fundamental differences in their genetic compositions and epidemiology. These patterns, along with the emergence of new CA-MRSA strains and the spread of healthcare-associated (HA-MRSA) strains into the general population, suggest the presence of a large reservoir of MRSA outside healthcare facilities, with antimicrobial resistance patterns being shaped by a multitude of selective pressures. [Reference David and Daum25] At the same time, the long-dominant CA-MRSA clone USA300 has declined in certain regions since its peak in the mid-to-late 2000s, potentially facilitating clonal replacement by more resistant strains during the mid-2010s. [Reference Planet26] Changes in prescribing practices may have further amplified these trends, as tetracycline use has increased in recent years among US outpatients. [Reference Kim27] These findings highlight the need for antibiotic stewardship strategies that address evolving resistance patterns, particularly for antibiotics commonly used in outpatient settings, as well as infection control strategies implemented in both community and healthcare settings.
The geographic distributions of relative risk and temporal trend were highly heterogeneous, with variation observed within states. This suggests that the risk of resistant MRSA infections is more likely driven by local, rather than broader statewide factors, emphasizing the importance of examining underlying place-based factors. While the use of antibiotics remains the primary driver of resistance, sociodemographic factors such as household crowding and poor health literacy have been shown to contribute to racial and economic health disparities in AMR bacterial infections across the US. [Reference Endale, Mathewos and Abdeta28–Reference Immergluck32] These disparities have been linked to neighbourhood characteristics, such as education and rurality, with significant clustering of resistant infections observed in areas of high deprivation. [Reference See33, Reference Cooper34] Climate factors can further influence patient behaviour and their interactions with their environment, likely contributing to the higher rates of resistance observed in the south. [Reference Sun, Klein and Laxminarayan35, Reference Klein36] The identification of CZs with elevated risk raises questions about which local-level sociodemographic, environmental, or public health factors may be driving these differences, underscoring the importance of local-level surveillance to inform future studies.
Although high- and low-risk CZs exhibited relatively stable temporal trends, CZ-specific relative risk plots revealed distinct temporal patterns in the rate and timing of change, highlighting the importance of localized surveillance, as larger spatial scales may obscure these differences. Several notable changes in relative risk were observed, particularly during the post-COVID-19 period from 2020 to 2022. In all three CZs examined for trends in tetracycline-resistant MRSA, increases were observed between 2020 and 2021. The timing and magnitude of these increases may reflect local factors or COVID-19-related disruptions, including shifts in prescribing practices and changes in healthcare delivery. While several studies reported decreased antibiotic prescriptions during the pandemic’s first wave, the sharp decline in in-person visits, except in patients with more severe infections, may have contributed to the underdiagnosis of bacterial infections, inflated resistance rates, and suboptimal antibiotic use. [Reference King37–Reference Baker41]
This study has several limitations. Although the VHA provides a nationwide sample with standardized microbiology testing, it is not representative of the broader US population, as VHA patients are more likely to be male gender, older in age, and White race. Aggregating to finer geographic units, such as CZs, results in low or zero counts, particularly in the Mountain West Region, requiring the analysis to be restricted to the eastern US. Additionally, use of a classification procedure for categorizing temporal trends, while conceptually advantageous, is better suited to detecting linear trends. High posterior probabilities are most likely in CZs with consistent increases over time, reducing the method’s sensitivity to fluctuating or non-linear trends. As demonstrated in the results, CZ-specific trends often exhibit non-linearity, which is better captured by examining posterior probabilities on a continuous scale rather than through discrete categorization.
Despite the decline of MRSA in the US, it remains a prominent cause of infection, with growing concerns over rising resistance to additional, commonly utilized antibiotics. This study characterizes the spatiotemporal patterns of relative risk for tetracycline- and TMP–SMX-resistant MRSA in the eastern US, identifying CZs with elevated risk and increasing local trends. In doing so, the findings illustrate how finer spatial scales can capture local variation that may be obscured at more broadly aggregated geographic units. The identification of stable yet spatially heterogenous patterns in tetracycline- and TMP–SMX-resistant MRSA points towards the influence of local, place-based factors, rather than broader statewide determinants. This highlights the need for a deeper understanding of how sociodemographic, ecological, and healthcare characteristics, particularly in the wake of the COVID-19 pandemic, shape the geographic distribution and temporal trends of AMR.
Supplementary material
The supplementary material for this article can be found at http://doi.org/10.1017/S0950268826101216.
Data availability statement
The data that support the findings of this study are from the US VHA electronic health record repository. Restrictions apply to the availability of the data, which contain protected health information and are not publicly available. Access may be granted with appropriate VHA approval. The R code used in this study is available in the Supplementary Methods.
Author contribution
Conceptualization: M.B., M.C., M.G.; Data curation: Q.S.; Supervision: M.C., M.G., S.H., J.O.; Formal analysis: M.B., J.O.
Competing interests
The authors declare none.
Open access
