Introduction

Effective antimicrobial stewardship is necessary to limit the emergence and spread of antimicrobial resistance (AMR),¹ and this includes the preferential use of agents least likely to select for drug resistant pathogens among the commensal flora. This typically consists of avoiding “broad” spectrum antimicrobials, as well as those to which resistance is currently rare, and these principles underlie the World Health Organisation AWaRe classification (Access, Watch and Reserve).² AWaRe is reflected in the current UK National Health Service standard contract, which requires hospitals to increase the proportion of “Access” antibiotics used, replacing previous requirements to reduce overall antibiotic use.³ However, these classifications are largely based on activity against pathogens rather than direct measures of AMR gene selection or microbiome disruption, and stewardship could be improved if such measures were available. For example, two antibiotics being considered for use may have similar spectra against a target pathogen, but very different impacts on AMR selection, either because the amounts reaching the commensal flora differ, or because they have differing spectra against non-pathogenic commensals that protect against colonisation with drug-resistant pathogens. Metagenomic sequencing provides a direct measure of AMR genes and microbiome composition in a sample, and its increasing availability in recent years is starting to allow the fuller impacts of antimicrobials to be measured.

The large intestine contains the vast majority of human commensal bacteria and is the primary reservoir for several clinically important commensal pathogens, particularly the Enterobacteriaceae (including Escherichia coli and Klebsiella pneumonia) and Enterococcus faecium. Increasing multi-drug resistance in these organisms represents a major global public health challenge.^4,5 Existing evidence linking antibiotic exposure to individual-level selection for AMR in the gut flora comes predominantly from small, healthy volunteer studies, which have shown that antibiotics can cause rapid microbiome disruption, but provide limited comparative data between antibiotics. They also may also have limited applicability to real patients, who often have recent exposure to several different agents and are at high risk of colonisation with drug resistant organisms. Randomised trials are a robust method to assess the impacts of different treatment approaches, but few have reported microbiome outcomes.^6–9 Another approach is to exploit variation in routine use of antibiotics in groups of patients at high risk of AMR to understand the nature and extent of differences between agents which cannot easily be achieved in other designs. Here we report results from a prospective observational study assessing the impact of antimicrobial use on the gut microbiome and selection of AMR genes. This study included two different sampling frames to produce independent and complimentary estimates, one cross-sectional, analysing a single stool sample from each participant, and one longitudinal, analysing changes in serial samples taken from the same participant admitted for haematopoietic stem cell transplant to enrich for broad spectrum antimicrobial exposure.

Methods

Study design and participants

The Antibiotic Resistance in the Microbiome – Oxford (ARMORD) study was an observational study that recruited healthy individuals living in Oxfordshire, and patients at the Oxford University Hospitals NHS Foundation Trust (OUH). The study was coordinated by the Nuffield Department of Medicine, University of Oxford, and was approved by the East Midlands-Leicester Central Research Ethics Committee (15/EM/0270).

The study involved two sampling strategies:

1. Cross-sectional sampling. Participants provided a single stool sample, and measures of the microbiome and AMR gene abundance were related to exposures recorded at the time of sampling.

2. Longitudinal sampling. Participants provided serial stool samples, and changes in the microbiome and AMR gene abundance between serial samples from the same individual were related to exposures between samples. Longitudinal sampling was only performed in patients admitted to the OUH haematology ward for haematopoietic cell transplant (HCT). The initial sample collected from these patients was also used in the cross-sectional analysis.

Participants were eligible if they were ≥18 years old, had no stoma or active inflammatory bowel disease, and were able to give informed consent and provide a history of recent antibiotic use. Inpatients and outpatients at OUH were approached by a member of the study team during routine care, and healthy individuals responded to articles in local media. After providing written informed consent, participants were interviewed and their medical notes were reviewed to collect information about antimicrobial exposures in the past year, recent travel, diet, alcohol and tobacco use, animal exposures, healthcare exposure and use of specific drugs (case report form in Appendix p8). Electronic patient records available at OUH included i) inpatient, emergency department and outpatient attendances, ii) inpatient and discharge antimicrobial prescriptions iii) microbiology, haematology and biochemistry results, iv) inpatient clinical observations (used to calculate the National Early Warning Score 2 [NEWS2] summary score of physiological abnormality),¹⁰ and v) discharge coding (including Charlson co-morbidity index).

Sample collection

In the cross-sectional stratum, participants were asked to provide the first stool sample passed after recruitment. This was stored at ambient temperature for a maximum of 24 hours before being frozen at −80°C. In the longitudinal stratum, participants were asked to provide a stool sample every other day until discharge. These were stored at 4-8°C for up to 72 hours before being frozen at −80°C.

DNA extraction, sequencing and bioinformatic analysis

DNA extraction was performed by bead beating in Lysing Matrix E tubes (MP Biomedicals) followed by QIAGEN Fast DNA Stool MiniKit (QIAGEN) (details in Appendix p1). Samples were sequenced in batches of 56-114 at the Oxford Centre for Genomics, and had automated normalisation and library preparation using either NEBNext Ultra or NEBNext Ultra II FS kits. All samples from the same batch were pooled and sequenced across 2-8 Illumina HiSeq4000 lanes using 150bp paired-end sequencing. Following human read removal all samples were subsampled to a depth of 3.5 million paired reads, and samples with fewer reads were excluded. Taxonomic classification was performed with MetaPhlAn2 (for diversity indices) and Kraken2 (for abundance of specific taxa). AMR gene detection in metagenomic sequence data was performed with the ARIBA software package, using the CARD database and ontology (v3.0.2) (details in Appendix p2).

Outcomes

Sequence data was used to derive three types of outcome:

1. Shannon diversity index - a single metric of diversity for each sample (calculated as the sum of p*ln(p) for all species, where p is proportional abundance). This index incorporates abundance and the number of species detected (richness).

2. Log relative abundance of specific bacterial taxa. The taxa of interest were two major groups of opportunistic pathogens (family Enterobacteriaceae and genus Enterococcus), and three major groups of anaerobes (phyla Bacteroidetes and Actinobacteria, and class Clostridia). Together these non-overlapping taxa accounted for 89% of the total abundance of organisms across all samples. If a particular taxon was not detected, its relative abundance was imputed as 10⁻⁶ (the approximate lower limit of detection).

3. Log relative abundance of specific classes of AMR genes. Five classes of clinically important resistance mechanisms were of interest: Clinically significant beta-lactamases (CTX-M, OXA, TEM, SHV), tetracycline ribosomal protection proteins, aminoglycoside transferases (AAC, ANT, APH), macrolide/clindamycin resistance genes (erm and mef), and the vanA vancomycin resistance gene. If a particular gene was not detected its relative abundance was imputed as 10⁻⁵ (the approximate lower limit of detection). Further details are in Appendix p3.

Statistical analysis

Cross-sectional sampling frame

Multivariable linear regression was used to estimate the effects of specific antimicrobial exposures on the outcomes above. Covariates were: age, sex, participant category (healthy, general medical, autologous stem cell transplant, allogeneic stem cell transplant), days of chemotherapy received (0 for non-HCT participants, truncated at 14), maximum Charlson comorbidity score in the year before sampling (identified from electronic health records), and the following physiological abnormalities in the fortnight prior to sample collection; maximum NEWS2 score, C-reactive protein (CRP) >50 mg/dL, white cell count (WCC) >11×10⁹/L, and WCC <0.5×10⁹/L. Healthy volunteers lived in the OUH catchment area but most had no previous activity at OUH, and in this case normal values were imputed (i.e. Charlson Index 0, NEWS2 score 0, CRP <50 mg/dL, WCC <11 and >0.5×10⁹/L). Other covariates such as diet and travel were not included in the final model, as they were not significantly associated with Shannon diversity in multivariable analysis and their inclusion did not materially affect other estimates. We did not make formal adjustment for multiple testing, as many of the antibiotic exposures were correlated, but interpreted level of evidence supporting findings in the context of the number of comparisons performed.

All antimicrobial exposures observed in >5 participants were included in the model. Individual agents were categorised separately if given by different administration routes (e.g. oral vancomycin and intravenous vancomycin), but route of administration was ignored in analyses of antimicrobial class (e.g. glycopeptides). In categorising antimicrobial class, beta-lactams were divided into ‘narrow spectrum’ (defined as penicillin, amoxicillin, flucloxacillin and first generation cephalosporins) and ‘broad spectrum’ (all others). Exposure to each antimicrobial was included as a separate variable on a scale of 0 (no exposure) to 1 (full exposure). In order to reflect both recency and total duration of antimicrobial use, this exposure was modelled as the area under an exponential decay curve of the form y = 2^−x/λ, in which λ is the microbiome disruption half-life, and x is time before sample collection. A single value of λ was used for all analyses, chosen as the common value across all antimicrobial exposures with the lowest Akaike Information Criterion across 1 to 14 days in the cross-sectional model for microbial diversity (6 days; Supplementary Table 1). The disruption half-life of 6 days means that after 6 days of an antibiotic course a patient would have an exposure of 0.5 to that agent, and after 12 days they would have an exposure of 0.75. Details of the exposure calculation, including graphical depictions are in the Appendix (p3 and Supplementary Figures 1-2).

Longitudinal sampling frame

Serial samples collected from participants undergoing HCT were used for the longitudinal analysis. The unit of analysis was a pair of consecutive samples collected from the same individual. For patients with >2 samples, each consecutive pair was used (i.e. sample 2 in pair 1 was sample 1 in pair 2, and so on, so the total number of pairs per participant is one less than the number of samples). Only pairs collected within 2-30 days were used. A multivariable linear regression model was used that was analogous to the cross-sectional model described above, except that the outcome was the change between the first and second samples in a pair, rather than absolute values (i.e. change in Shannon diversity, or change in log relative abundance of bacterial taxa or AMR genes). Because pairs of samples from the same participant may not be independent of one another, robust standard errors were used to allow for possible clustering.

Covariates were: age, sex, type of transplant (allogeneic or autologous), days of chemotherapy received at collection of sample 1 (0 in non-HCT patients, and truncated at 14), number of days between sample 1 and sample 2, change in NEWS2 score between sample 1 and sample 2, and the presence of the following new physiological abnormalities recorded after collection of sample 1 and before collection of sample 2; CRP >50 mg/dL, WCC >11×10⁹/L, WCC<0.5×10⁹/L. For each outcome, the value for the first sample in the pair was also included as a covariate (i.e. baseline diversity or log relative abundance of taxa/AMR genes).

Antimicrobial exposures were calculated for samples as in the cross-sectional analysis above, and the exposure for each pair was the difference between the first and second samples. This has the following implications: i) if a patient starts an antimicrobial after the first sample is collected then the exposure for that pair is the same as for sample 2, ii) if a patient is on long-term antimicrobial treatment then the exposure for a pair is zero (as one would not expect this to lead to a difference between samples), iii) if a patient stops antimicrobial treatment shortly after sample 1 then the exposure to that agent will be negative (as one expects gut microbiome diversity to increase after stopping an antimicrobial). Truncating the small number of negative values at zero had little impact on results (data not shown).

Analysis was performed in R v4.2.3 with the ontologyIndex package (v2.4).

Results

Between July 2015 and November 2018, 225 participants were recruited and had at least one sequenced stool sample, all of whom were included in the cross-sectional analysis (Table 1 and Supplementary Figure 3 [CONSORT diagram]). Thirty-three (15%) were healthy volunteers, 91 (40%) were general medical patients, and 101 (45%) were HCT patients. The healthy volunteers were on average younger, more likely to be female, and rarely had recent antibiotic exposure. In contrast, 184/192 (96%) of medical and haematology patients had received antibiotics in the past year, and 97 (51%) of these were receiving antibiotics at the time of sampling. Of those with antibiotic exposure in the past month, 99/148 (67%) had received >1 type of antimicrobial. Of 101 HCT participants in the cross-sectional analysis, 79 had >1 sequenced sample and so also contributed to the longitudinal analysis, and 173 sample pairs were included in this analysis.

Table 1

Impact of antimicrobial exposures on gut microbiome diversity

A half-life of 6 days was identified as the best fit for the antimicrobial exposure in the cross-sectional model of microbiome diversity, and this value was used for all subsequent analyses (Supplementary table 1). The independent effects of antimicrobial exposure on gut microbiome diversity in the cross-sectional and longitudinal models are shown in Figure 1. The cross-sectional model (Figure 1A) provided more precise estimates of microbiome disruption than the longitudinal model, and four antimicrobial exposures were associated with a significant reduction in gut Shannon diversity: iv co-amoxiclav (−3.0 [95% Confidence Interval −4.7 to −1.4; p = 0.0005), iv piperacillin-tazobactam (−2.6, [−3.4 to −1.8]; p <10⁻⁹), oral clindamycin (−2.2 [−3.5 to −0.8]; p = 0.002), and iv meropenem (−1.6 [−2.6 to −0.5]; p = 0.003). There was no evidence that the impact of iv co-amoxiclav was greater than iv meropenem (p=0.13). In contrast, co-trimoxazole and doxycycline were associated with minimal reductions in gut microbiome diversity, as were azole antifungals and acyclovir (lower 95% CI above −1.0). In the longitudinal analysis (Figure 1B), only five exposures had data from more than 20 sample pairs, and these were consistent with the cross-sectional results, including large reductions in diversity associated with exposure to piperacillin-tazobactam (−3.9 [−5.0 to −2.9]; p <10⁻¹⁰) and meropenem (−3.0 [−4.1 to −1.9]; p <10⁻⁶) (data on iv co-amoxiclav and oral clindamycin insufficient for comparison). In the longitudinal model, gentamicin was associated with marginally significant increased gut microbiome diversity (+6.6 [+1.3 to +11.8]; p = 0.02).

Figure 1 -

Analyses by antimicrobial class were also largely consistent between cross-sectional and longitudinal analyses (Supplementary Figure 4). In both analyses ‘narrow’ spectrum beta-lactams were associated with substantially lower microbiome disruption than ‘broad’ spectrum beta-lactams. In the cross-sectional analysis, several other classes of antibiotics also had little to no impact on the microbiome diversity; antifolates, macrolides, aminoglycosides and tetracyclines (lower 95% CI above −1.0).

Effects of antimicrobials on the abundance of specific taxa and AMR genes

Piperacillin-tazobactam, meropenem, and ciprofloxacin were associated with significant decreases in the relative abundance of Enterobacteriaceae in both models, as were ceftriaxone and ceftazidime in the cross-sectional model (>1,000 fold decreased relative abundance, p<0.01, Figure 2 & Supplementary Figure 5). There was no evidence of effects of iv co-amoxiclav (p=0.90) or oral clindamycin (p=0.31) on relative abundance of Enterobacteriaceae. Piperacillin-tazobactam and meropenem were also associated with large increases in the relative abundance of Enterococcus in both models (>100 fold increased relative abundance, p<0.01), as were iv co-amoxiclav (p=0.04), iv ceftriaxone (p=0.02) and iv vancomycin (p=0.03) in the cross-sectional model. Many antibiotics were associated with large reductions in the relative abundance of major anaerobe groups, particularly Actinobacteria (reductions in which were associated with exposure to ceftazidime, oral ciprofloxacin, oral clindamycin, iv co-amoxiclav, meropenem, piperacillin-tazobactam, and iv vancomycin). When analysed by antimicrobial class, broad spectrum beta-lactams were associated with reductions all anaerobe groups (Supplementary Figure 5).

Figure 2 -

Several antimicrobial exposures were also associated with increases in the relative abuindance of AMR genes (Figure 3, Supplementary Figure 6). In the cross-sectional analysis, iv piperacillin-tazobactam and meropenem were both associated with increases in several AMR mechanisms, including aminoglycoside transferases and macrolide resistance genes (all p<0.02) but there was no evidence that these antibiotics had any effect on overall relative abundance of beta-lactamases (p≥0.14), potentially reflecting lower abundance of relative species (see above) but greater AMR gene carriage in those surviving. In the cross-sectional model, exposure to several classes of antimicrobials were associated with an increase in the abundance of corresponding resistance genes, including glycopeptides, tetracyclines, macrolides and clindamycin (all p≤0.02). Broad-spectrum beta-lactam use was not associated with an increase in the abundance of beta-lactamases, but it was associated with an increase in the abundance of glyopeptide and aminoglycoside AMR genes (both p≤0.0001). Increased aminoglycoside AMR gene abundance was also associated with clindamycin and glycopeptide use (both p≤0.02), but not aminoglycoside use.

Figure 3 -

Discussion

Our findings demonstrate that simultaneous modelling of multiple antimicrobial exposures in a heterogeneous and heavily antibiotic exposed population can produce direct, quantitative comparisons of the impact of different agents on multiple aspects of the gut microbiome. Several broad-spectrum antibiotics were associated with large and rapid reductions in gut microbiome diversity, including decreased abundance of several major anaerobe groups, and increased abundance of Enterococcus species, along with glycopeptide and aminoglycoside resistance mechanisms often found in Enterococcus faecium.⁵ The plausibility of these results is reinforced by the consistency between the two independent analysis frameworks that were used; cross-sectional and longitudinal. Microbiome disruption was clearest withclindamycin and broad spectrum beta-lactams including iv co-amoxiclav, piperacillin-tazobactam, and meropenem, consistent with some previous microbiome studies of these drugs.^11–14 This is also in keeping with estimates of risk of Clostridioides difficile infection (CDI) which is highest with clindamycin, fluoroquinolones, carbapenems and third-generation cephalosporins.^15,16 The limited impact of doxycycline on diversity is also consistent with the lower risk of CDI observed with tetracyclines. The decreased relative abundance of Enterobacteriaceae seen in this study with piperacillin-tazobactam, ceftriaxone, meropenem, and ciprofloxacin correspond with reductions seen in culture based studies.¹⁷ The absence of any clear effect of beta-lactams on the abundance of beta-lactamase gene abundance may be because there were too few patients in this study to create robust categories of beta-lactamse by resistance spectrum, so opposite impacts on different beta-lactamase genes would have been combined.

Observational studies of the gut microbiome allow large numbers of patients to be recruited much more easily than interventional studies such as clinical trials. However, using data from these to inform antibiotic usage is complicated by potential biases from confounding, and because of difficulties in accurately modelling multiple exposures in patients receiving many different antibiotics. In ARMORD, the apparent impact of vancomycin and clarithromycin on diversity was substantially reduced when adjusting for other antibiotic exposures (Figure 1A), in keeping with the fact that these drugs are typically given alongside broad-spectrum beta-lactams in the UK. Nevertheless, the complexity of antibiotic exposure captured in observational data also more closely reflects real life, whereas clinical trials generally manipulate one single antibiotic, and restrict background antibiotic exposure to produce a cleaner, but potentially less generalisable, comparison.

The large degree of inter-individual variation in the gut microbiome introduces an additional source of variation to cross-sectional studies compared to longitudinal sampling. Despite this, cross-sectional sampling in ARMORD generally gave more precise estimates than longitudinal sampling, despite using a similar number of samples. Recruiting any hospitalised patients for longitudinal sampling before they have received antibiotics is challenging, and this was true in ARMORD, even among patients admitted electively for HCT. In ARMORD the quasi-experimental approach of longitudinal sampling starting before antimicrobial exposure had no advantage over cross-sectional sampling, and complicated recruitment.

The ARMORD study has important limitations. Along with age and sex, several markers of acute illness and comorbidity were included in the model to adjust for potential confounding, but there may be residual and unmeasured confounding related to factors not adequately represented in the model. The microbiome disruption half-life used was 6 days, as this value best fitted the overall data, which implies that the majority of the disruption and recovery of the bowel flora diversity occurs rapidly after starting and stopping antimicrobials (Supplementary Figure 2). This is in keeping with interventional studies in volunteers,^13,18 but is a simplification that does not account for some longer term impacts of antibiotic use that can be detected months after treatment.^19,20 A larger study may have allowed derivation of exposure models that take this into account, or could test multiple different ways to represent the impact of antibiotic exposure on the gut microbiome, or that allow for different disruption dynamics for different antimicrobials. Also, because the study was observational, no data were available on drugs that are not commonly used at OUH, including imipenem, cefepime and aztreonam. Finally, although we estimate the independent effects of multiple exposures, the uncertainty is still large making it challenging to use these results to adequately inform clinical use, as this would require clear distinction between antibiotics that might be given for the same indication. For example, the cross-sectional data are consistent with an identical average reduction in diversity with piperacillin-tazobactam, meropenem, ceftazidime, ceftriaxone, intravenous co-amoxiclav, ciprofloxacin and clindamycin, but they are also consistent with important differences between these drugs. A larger study is required to identify or rule out such differences.

Overall, however, simultaneous estimation of the impact of over 20 antimicrobials on the gut microbiome and AMR gene abundance highlighted important differences between individual drugs, and that some drugs in the WHO Access group (co-amoxiclav, clindamycin) had similar magnitude impact on microbiome diversity to those in the Watch group (meropenem, piperacillin-tazobactam) with potential implications for acquisition of other resistant organisms. The consistency of the ARMORD results between sampling frames and with previous studies supports the wider use of observational metagenomic studies to compare the impact of antimicrobials on the gut microbiome. Although some challenges remain, such as identifying an optimal measure of antimicrobial exposure, this is a practical approach to inform future research and stewardship.

Data Availability

All data produced in the present study are available upon reasonable request to the authors

Acknowledgements

We would like to thank all the participants who took part in ARMORD, and also the nurses, healthcare support workers, and doctors on the OUH Haematology ward for their help with this study. This work was supported by the National Institute for Health Research Health Protection Research Unit (NIHR HPRU) in Healthcare Associated Infections and Antimicrobial Resistance at Oxford University in partnership with the UK Health Security Agency (NIHR200915), and the NIHR Biomedical Research Centre, Oxford. ASW is an NIHR Senior Investigator. The views expressed are those of the authors and not necessarily those of the NHS, the NIHR, the Department of Health or the UK Health Security Agency. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Declaration of interests

No authors have no conflict of interest to declare.