Metagenomic Analysis of Antibiotic Pressure Antibiotic selection pressure determination through ...

43 Pages·2015·3.15 MB·English

AAC Accepted Manuscript Posted Online 14 September 2015
Antimicrob. Agents Chemother. doi:10.1128/AAC.01504-15
Copyright © 2015, American Society for Microbiology. All Rights Reserved.
Metagenomic Analysis of Antibiotic Pressure
1 Antibiotic selection pressure determination through sequence-based metagenomics
2
3 Matthias Willmanna,b#, Mohamed El-Hadidic, Daniel H. Husonc, Monika Schütza,b,
4 Christopher Weidenmaiera,b, Ingo B. Autenrietha,b and Silke Petera,b
5
6 Institute of Medical Microbiology and Hygiene, University of Tübingen, Tübingen,
D
o
7 Germanya; German Center for Infection Research (DZIF), partner site Tübingen, Tübingen, w
n
lo
8 Germanyb; Center for Bioinformatics, University of Tübingen, 72076 Tübingen, Germanyc. ad
e
d

9 fro
m

h
10 #Address correspondence to Matthias Willmann, [email protected]. tt
p
:
/
/
a
11 a
c
.
a
s
12 Running Title: Metagenomic Analysis of Antibiotic Pressure m
.
o
r
g
/

o
n

J
a
n
u
a
r
y

3
1
,

2
0
1
8

b
y

g
u
e
s
t
1
Metagenomic Analysis of Antibiotic Pressure
13 Abstract
14 The human gut forms a dynamic reservoir of antibiotic resistance genes (ARGs). Treatment
15 with antimicrobial agents has a significant impact on the intestinal resistome, and leads to
16 enhanced horizontal transfer and selection of resistance. We have followed the development
17 of intestinal ARGs over a six-days course of ciprofloxacin (Cp) treatment in two healthy
D
18 individuals by using sequenced-based metagenomics and different ARG quantification o
w
n
19 methods. Fixed- and random-effects models were applied to determine the change in ARG lo
a
d
20 abundance per defined daily dose of Cp as an expression of the respective selection pressure. e
d

f
r
21 Among various shifts in the composition of the intestinal resistome we found in one o
m

h
22 individual a strong positive selection for class D beta-lactamases which were partly located on t
t
p
:
/
23 a mobile genetic element. Furthermore, a trend to a negative selection has been observed with /a
a
c
24 class A beta-lactamases (-2.66 hits per million sample reads / defined daily dose; p = 0.06). .a
s
m
25 By 4 weeks after the end of treatment, the composition of ARGs returned toward their initial .o
r
g
26 state but to a different degree in both subjects. We present here a novel analysis algorithm for o/
n

J
27 the determination of antibiotic selection pressure which can be applied in clinical settings to a
n
u
28 compare therapeutic regimes regarding their effect on the intestinal resistome. This a
r
y

3
29 information is of critical importance for clinicians to choose antimicrobial agents with a low 1
,

2
30 selective force on their patients’ intestinal ARGs, likely resulting in a diminished spread of 0
1
8

31 resistance and a reduced burden of hospital-acquired infections with multidrug resistant b
y

g
u
32 pathogens. e
s
t
2
Metagenomic Analysis of Antibiotic Pressure
33 Background
34 Antibiotic resistance is at present one of the most serious public health issue across the globe,
35 threatening the achievements of modern medicine. The occurrence of hospital-acquired
36 infections has increased steadily, and this leaves only a limited number of therapeutic options
37 particularly in immunocompromised patients who are reliant on effective antimicrobial
D
38 substances for treatment and prevention of infection (1, 2). The human intestinal microflora is o
w
n
39 known to harbor a vast reservoir of antibiotic resistance genes (ARGs) with geographically lo
a
d
40 variable composition and abundance as shown in samples from Spain that exhibited a higher e
d

f
r
41 relative resistance potential than samples from Denmark or the U.S. (3, 4). The intestinal o
m

h
42 resistome is a dynamic and open system which constantly alters due to a transfer of resistance t
t
p
:
/
43 determinants within and between bacterial species via mobile genetic elements, such as /a
a
c
44 transposons, prophages and plasmids (5). Horizontal transfer of ARGs from the intestinal .a
s
m
45 flora to pathogenic microorganisms poses a notable danger. This process is enhanced by .o
r
g
46 treatment with broad-spectrum antibiotics, which also have an impact on non-targeted o/
n

J
47 organisms like commensals, and it presumable enriches the pool of resistance genes a
n
u
48 accessible to pathogens (6-9). Furthermore, antibiotic treatment might enable already resistant a
r
y

3
49 pathogens to gain dominance over other bacterial species, leading to an increased risk of 1
,

2
50 infection with the resistant strain (10). This indicates the relevance of a customized 0
1
8

51 administration of antimicrobial agents, particularly in a hospital environment due to a possibly b
y

g
u
52 massive spread of resistance. e
s
t
53 Antibiotic stewardship programs (ABS) are one way to improve antibiotic usage in hospitals
54 (11). However, ABS strategies should be based on a detailed understanding of the impact of
55 various antimicrobial agents on the human intestinal resistome in different patient populations
56 and hospital settings. This would allow a more efficient allocation of resources and enable
57 correct decision making about which antibiotics should be primarily used or restricted. Most
3
Metagenomic Analysis of Antibiotic Pressure
58 of our knowledge about the association of antibiotic selection pressure and resistance spread
59 is derived from studies of cultured microorganisms (12, 13), which reflect only a minor
60 fraction of the overall microbial population in the gut. Thus, such studies do not represent the
61 immense diversity and variation of the intestinal resistome. Alternatively, in a sequence-based
62 metagenomics approach a comprehensive sequence collection, which includes unculturable
63 and rare taxa in addition to the usually culturable pathogens, is generated from microbial D
o
w
64 DNA of the human faeces and mapped to a resistance gene database, allowing a deeper n
lo
a
65 exploration of ARGs (14, 15) and determination of a more representative selection pressure. d
e
d

f
r
66 In this study, we investigated the impact of a six-days course of ciprofloxacin administration o
m

h
67 on the intestinal resistome and microbiota of two healthy volunteers by applying a sequence- t
t
p
:
/
68 based metagenomics time series analysis and different ARG quantification methods. We /a
a
c
69 propose a novel algorithm for the determination of treatment-induced selection pressure in .a
s
m
70 various groups of ARGs and discuss the merit of such an approach for clinical studies that .o
r
g
71 aim at identifying antibiotic regimes which impose a low selective force on clinically relevant o/
n

J
72 ARGs and ARG groups. a
n
u
a
73 ry

3
1
,

74 Materials and Methods 2
0
1
8

b
75 Participants and sampling y

g
u
e
76 Two healthy male adult volunteers were recruited in Tuebingen, Germany. Participants were s
t
77 excluded when treated with antibiotics within the previous 12 months, age under 18 years,
78 pregnancy and past reaction to fluoroquinolones. Study approval was given by the local ethics
79 committee (662/2013BO1); and written consent was obtained. A six-day course of
80 ciprofloxacin (Cp) was administered with 500 mg orally twice daily corresponding to one
81 defined daily dose (DDD) (16). The two participants were asked to document any symptoms
4
Metagenomic Analysis of Antibiotic Pressure
82 during Cp treatment but none were reported. Stool samples from both participants were
83 collected before treatment (day 0), during treatment (day 1, 3 and 6), and 2 and 28 days after
84 treatment.
85 Sample collection and screening for resistant bacteria
86 Stool samples were collected in sterile plastic devices (nerbe plus GmbH, Germany) and
D
o
87 stored within 30 min at 4°C. Samples were plated within three hours of collection on selective w
n
lo
88 agars for the detection of extended-spectrum beta-lactamase producing Enterobacteriaceae a
d
e
89 (chromID ESBL, bioMérieux, France), vancomycin-resistant enterococci (chromID VRE, d
f
r
o
90 bioMérieux, France), methicillin-resistant Staphylococcus aureus (MRSA) (chromID MRSA, m

h
t
91 bioMérieux, France) and Pseudomonas aeruginosa (Cetrimide Agar, BectonDickinson, tp
:
/
/
a
92 Germany). Plates were incubated aerobically at 37°C for 48 h and were inspected for growth a
c
.
a
93 at 24 h and 48 h. s
m
.
o
r
g
94 DNA extraction and sequencing /

o
n

J
95 DNA extraction was performed within 3 h of sample collection according to the Human a
n
u
a
96 Microbiome Project protocol with minor modifications r
y

3
97 (http://hmpdacc.org/doc/HMP_MOP_Version12_0_072910.pdf, accessed November 2013) 1
,

2
0
98 and using the Power Soil DNA Isolation Kit (MO BIO Laboratories, CA, USA). In brief, 1
8

b
99 subsamples of ̴ 3 g stool were suspended in 7 ml lysis buffer, shaken for 30 sec and y

g
u
100 centrifuged at 3000 rpm for 5 min. 1000 µl supernatant was transferred into a bead tube e
s
t
101 (garnet bead tubes, MO BIO Laboratories, CA, USA) and incubated at 65°C for 10 min
102 followed by an incubation at 95°C for 10 min, shaking at 300 rpm. Three tubes per sample
103 were stored immediately at -80°C prior further steps. After unfreezing, 500 µl of tube content
104 was transferred into a new garnet tube after removal of 500 µl lysis buffer from the new tube.
105 DNA extraction was subsequently performed according to the manufacturers’ instructions.
5
Metagenomic Analysis of Antibiotic Pressure
106 DNA was eluted in 100 µl water (HyClone water, molecular biology grade, GE Healthcare,
107 Germany). DNA quantification was carried out using Qubit Analyzer (Invitrogene, Life
108 Technologies, Singapore) following standard protocols. Metagenomic shotgun sequencing
109 was performed at GATC Biotech AG (Konstanz, Germany) on an Illumina HiSeq 2000
110 platform using a paired-end sequencing approach with a targeted read length of 100 bp and an
111 insert size of 180 bp. D
o
w
n
112 Real-time PCR for detection of tetQ and cblA lo
a
d
e
d
113 The tetQ and cblA genes were detected using the primers tetQ-fw
f
r
o
m
114 (ATCTGCTGTTTGCCAGTGGA), tetQ-rv (TGTATGCCTTCCTTTGCGGA), cblA-fw

h
t
t
115 (CACTTCCCCTTGCTCAGTGT), and cblA-rv (CGTGAAATCCTGGTCGGGAA) p
:
/
/
a
116 purchased from tib-molbiol, Germany. The rtPCR was performed in a total volume of 20 µl a
c
.
a
117 using the SYBR green Jump Start mix (Sigma-Aldrich, Germany), 10 µM of each primer, 25 sm
.
o
118 mM MgCl2 and 1 ng DNA template. HyClone water (HyClone water, molecular biology rg
/

o
119 grade, GE Healthcare, Germany) was used as a negative control. The following cycling n

J
a
120 conditions were applied: 95°C 30 s, 40 cycles of 95°C 15 s, 55°C 30 s and 72°C 30 s. The n
u
a
r
121 specificity of the PCR product was confirmed by agarose gel electrophoresis following y

3
1
122 Sanger sequencing (GATC Biotech AG, Germany). ,
2
0
1
123 DNA abundance calculation was based on the quantification cycle c defined as the number of 8

b
y
124 cycles at which the fluorescence crosses the threshold. It was assumed that the quantity of the g
u
e
125 DNA target would double every cycle. DNA target quantity T was calculated as follows: s
t
1
(cid:1846)=
2(cid:3030)
126 T was normalized to the value of 1 for the sample with the highest value with all other values
127 as fractions of 1. The same was applied to the metagenomic quantification values relative
128 abundance (RA) and omega (ω). Normalized values were visually compared in terms of their
6
Metagenomic Analysis of Antibiotic Pressure
129 quantity and progression over time. A high concordance between rtPCR and sequenced-based
130 metagenomics was assumed when an a priori maximum deviation of 0.25 in each time point
131 between all measures of quantification was not exceeded.
132 Sequence data analysis and database construction
133 For 12 samples, a total number of 820,599,346 raw reads with a length of 101 bp was
D
o
134 generated. Quality control was performed using FastQC (available from w
n
lo
135 http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) with all sequenced libraries a
d
e
d
136 passing commonly applied quality measures. Raw reads were trimmed to contain a minimum
f
r
o
m
137 quality phred score of 30 using the FASTX-Toolkit (available from

h
t
t
138 http://hannonlab.cshl.edu/fastx_toolkit/index.html). Read pairs were placed into one file p
:
/
/
a
139 before further processing. Reads were mapped against an assembly of the human genome a
c
.
a
140 (release CRCh37, downloaded from NCBI) using an in-house Smith-Waterman based DNA sm
.
o
141 aligner called SASS aligner and any read that matched human sequence was removed. Raw rg
/

o
142 reads without human sequences were uploaded to the European Nucleotide Archive (study n

J
a
143 accession number: PRJEB10391). n
u
a
r
y

144 Metagenomic assembly of each sample was performed using Ray Meta version 2.3.1 with a 3
1
,

2
145 K-mer length of 21 (17). MetaGeneMark version 2.8 was applied to predict proteins from the 0
1
8
146 assembled contigs using the default setting (18). In order to generate a comprehensive study b
y

g
147 protein catalogue (CSPC), all predicted proteins from 12 samples were pooled together into u
e
s
148 one file. Sequences with a length shorter than 140 amino acids (aa) were removed. The t
149 threshold of 140 aa is equivalent to the 10th percentile length of NCBI RefSeq bacterial
150 sequences. Redundant sequences were removed using the cluster algorithm of CD-HIT
151 version 4.6 (19). Proteins with an alignment coverage greater or equal to 90% of the full
152 length of the shorter sequence and with more than 95% sequence similarity were removed as
153 redundancies. NCBI RefSeq bacterial protein sequences were downloaded from NCBI FTP
7
Metagenomic Analysis of Antibiotic Pressure
154 site (access date: 31.12.2014) and merged with those sequences of the subject specific
155 database. An in-house DNA to amino acid alignment program called MALT (available from
156 http://ab.inf.uni-tuebingen.de/software/malt/) was used to translate reads in all six readings
157 frames and then to perform a semi-global alignment against the given database. Sequences
158 were removed from the subject specific database if they were at least 99% identical to RefSeq
159 proteins, thus keeping the RefSeq annotation in the merged database. This database was D
o
w
160 subsequently merged with the comprehensive antibiotic resistance database (CARD; access n
lo
a
161 date: 05.02.2015; http://arpcard.mcmaster.ca). Redundant proteins were removed as described d
e
d
162 above, keeping the CARD annotation. Generation of the CSPC as well as the analytical fr
o
m
163 flowchart are illustrated in figure 1. h
t
t
p
:
/
164 Taxonomic profiling /a
a
c
.
a
165 Quality controlled reads were mapped against NCBI nr database (access date: 31.12.2014) s
m
.
o
166 using DIAMOND (20) and imported in MEGAN version 5 for taxonomic assignment (21). r
g
/

o
167 Lowest common ancestor parameters were used as default with the exception of a maximum n

J
a
168 expected value of 0.001. Phylogenetic parameters of all samples are displayed in table S1. n
u
a
169 Principle coordinate analysis of species distribution was conducted with the cluster analysis ry

3
1
170 function in MEGAN. Diversity index values (Shannon and Simpson index) were determined ,
2
0
1
171 with MEGAN on species level. Before comparison all samples were normalized to an equal 8

b
y
172 number of reads to account for unequal sampling effort. Stream plots were generated by the
g
u
e
173 streamgraph htmlwidget R package, horizon plots by the lattice extra R package. The s
t
174 horizonscale was not specified as by default, and the Δ value was calculated for each phylum
175 by subtracting the lowest from the highest sample read count followed by a division by three
176 (http://lmdvr.r-forge.r-project.org/figures/figures.html).
177 Quantification of antibiotic resistance genes and groups
8
Metagenomic Analysis of Antibiotic Pressure
178 Antibiotic resistance genes were defined as such genes with the potency to directly mediate
179 resistance to an antimicrobial agent. Genes that mediate resistance only in the presence of
180 certain mutations or those with regulatory functions were excluded. This was done because
181 short reads do not always cover the variable gene regions which are involved in mediating
182 resistance. Such reads cannot be used for metagenomic quantification due to the lack of
183 information whether they are derived from a wild type or a resistant variant of the gene. D
o
w
184 Regulatory elements on the other hand do not directly mediate antimicrobial resistance and n
lo
a
185 were excluded for this reason. d
e
d

f
r
186 Three different methods of metagenomic quantification were applied and compared. A best o
m

h
187 hit algorithm was used for all methods: Quality controlled reads were aligned against the t
t
p
:
/
188 CSPC database or CARD using DIAMOND (20) with an either 90% or 100% similarity /a
a
c
189 threshold. Mapped reads were counted once respecting the match with the highest bit score. .a
s
m
190 The maximum e-value to report alignments was 0.001. .o
r
g
/

o
191 The calculation of the relative abundance (RA) has been reported previously (22). For any n

J
a
192 sample S, the copy number of the antibiotic resistance gene i is given by n
u
a
r
y
(cid:1854) = (cid:1876)(cid:3036) 3
(cid:3036) (cid:1838) 1
(cid:3036) ,
2
0
1
193 and the relative abundance of the antibiotic resistance gene i is given by 8

b
y
(cid:1854)(cid:3036) g
(cid:1844)(cid:1827)(cid:3036) = ∑(cid:1854) ue
s
t
194 where
195 RAi is the relative abundance of antibiotic resistance gene i in sample S,
196 Li is the length of antibiotic resistance gene i,
197 xi is the number of reads that were mapped against the antibiotic resistance gene i in sample S,
198 bi is the copy number of the antibiotic resistance gene i in sample S, and
199 ∑b is the sum of the copy number of all detected genes in sample S.
9
Metagenomic Analysis of Antibiotic Pressure
200
201 We considered this calculation as reference standard in our study. However, one drawback is
202 that the results depend on the total number of mapped reads which might considerably
203 fluctuate between samples. Therefore, we used two alternative metagenomic quantification
204 methods that are independent of mapping fluctuations.
D
205 To define the first alternative quantification method, for each sample S and each o
w
n
206 corresponding count xi of reads assigned to an antibiotic gene or group i, we set lo
a
d
e
(cid:2033) = (cid:1876)(cid:3036) 10(cid:2874) d
(cid:3036) (cid:1865) fr
o
m

h
207 where m denotes the number of reads in the sample S. This is the normalized hit count ωi for tt
p
:
/
208 the antibiotic resistance gene i in sample S. The unit of ωi is hits per million sample reads /aa
c
.
209 (hmr). Gene groups can be quantified in the same way. a
s
m
.
o
210 The third metagenomic quantification method ωCi for the antibiotic gene i in sample S is rg
/

o
211 calculated in the same way than ωi. The only difference between both approaches is the n
J
a
212 database used for the initial mapping, which is the CSPC for ωi and CARD for ωCi. nu
a
r
y

213 Calculation of selection pressure 3
1
,

2
0
214 We collected cross-sectional time-series data to investigate the amount of ARGs across time 1
8

b
215 and Cp administration. Two techniques are usually used to analyse such data, both accounting y

g
u
216 for individual heterogeneity: fixed- and random-effects models. These models can calculate e
s
t
217 the change in ARGs per defined daily dose of an antibiotic, and coefficients constitute
218 appropriate estimates of the antibiotic selection pressure over the period of treatment. The
219 preferred model was chosen on the base of the Hausman test (23, 24). Furthermore, models
220 were tested for time-fixed effects and heteroskedasticity, and robust standard errors were
221 calculated when the latter was present. Since random-effects models are based on a high
10