Characterization of antibiotic resistomes by reprogrammed bacteriophage-enabled functional metagenomics in clinical strains

This research complies with all relevant ethical regulations approved by the Human Investigation Review Board of Albert Szent-Györgyi Clinical Centre of the University of Szeged and the National Biodiversity Authority (NBA) of India. Permission for the faecal sample collection was obtained from the Human Investigation Review Board of Albert Szent-Györgyi Clinical Centre, University of Szeged (registered under 72/2019-SZTE). Volunteer participants were selected on the basis of strict criteria that (1) they did not take any antibiotics for at least one yr before sample donation and (2) they are in a good health. These requirements are standard in the field and secure a bias-free comparison of the antibiotic resistomes in the healthy human gut microbiome. Informed consent was obtained from all participants. Soil and river sediment sample collection from around the city of Hyderabad and Lucknow was approved by the National Biodiversity Authority (NBA), India (application number: NBA/Tech Appl/9/1822/17/18-19/3535). No statistical methods were used to pre-determine sample sizes, but our sample sizes are similar to those reported in previous publications^18,51,52. Samples were not allocated to groups. Samples for each individual experiment were handled by one person in charge. Data collection and analysis were not performed blind to the conditions of the experiments. No data were excluded from the analysis. Unless otherwise stated, when using a kit, we followed the manufacturer’s instructions.

Plasmid construction for DEEPMINE

A custom plasmid was created from pZE21 expression vector (Supplementary Table 11) for compatibility with the T7 transduction and the sequencing pipelines. Specifically, the replication origin was switched from ColE1 to p15A, and the packaging signal of the T7 bacteriophage was introduced (enzymes and primers used are listed in Supplementary Table 11). Subsequently, the pZE21_p15A vector was amplified by PCR using a mixture of primers containing 10-nt-long random barcodes (Supplementary Table 11), followed by digestion and self-ligation.

Sample collection and construction of metagenomic libraries

For the gut microbiome library, we collected faecal samples from 10 unrelated, healthy individuals with no history of taking antibiotics in the year before sample donation. For the anthropogenic soil microbiome, samples were collected from highly antibiotic-contaminated industrial areas in India⁵³. Metagenomic DNA from the gut and soil samples was extracted using DNeasy PowerSoil kit (Qiagen, 47016). Genomic DNA of clinical bacterial isolates (Supplementary Table 1) was isolated using the Sigma GenElute bacterial genomic DNA kit (Sigma, NA2110-1KT).

From each sample, 40 µg of extracted DNA was digested with MluCI enzyme (NEB, R0538L) (10 min, 37 °C), followed by inactivation (20 min, 85 °C). The quantity of the MluCI enzyme was varied to obtain DNA in the target size range of 1–5 kbp. DNA was isolated with pulsed field gel electrophoresis (Sage Science, PB02901) with a 0.75% agarose gel cassette and low-voltage 1–6 kbp marker S1 cassette definition. The metagenomic DNA fragments were ligated into the pZE21_p15A plasmid at the EcoRI site using a 3:1 mass ratio of insert:vector. Pure ligation mixture was electroporated into 40 µl of either E. coli MegaX (Invitrogen, C640003) or E. coli 10G ELITE (Lucigen, 60080-2) cells. Following one h of incubation at 37 °C, transformants were plated onto 50 µg ml⁻¹ kanamycin containing Luria Bertani (LB) agar plates in 10¹×, 10²× and 10³× dilutions for colony forming unit determination. The rest of the recovered cells were grown overnight on LB agar plates supplemented with kanamycin. The next day, plasmids were isolated. Insert size distribution was estimated by PCR amplification of relevant plasmid regions from 10–20 randomly selected clones. The average insert size was determined to be 2–3 kbp.

Transducing hybrid bacteriophage particle preparation

Transducing hybrid bacteriophage preparation was adapted from ref. ²¹. In brief, E. coli BW25113 cells containing phage-tail-encoding plasmids (Supplementary Table 11) were grown to optical density (OD)_600nm ~0.7 (250 r.p.m. at 37 °C), then placed on ice for 15 min. Next, cultures were centrifuged (2,200 × g, 4 °C, 10 min), supernatant was discarded and the cells resuspended in the same amount of medium (LB or Terrific Broth (TB)). Afterwards, T7 bacteriophages lacking T7 fibre-encoding regions (T7∆(gp11-12-17)) were used to infect cells at multiplicity of infection (MOI) 2–3. Following 2 h of incubation (100 r.p.m., 37 °C), cells were treated with 2% chloroform and vortexed. The mixture was then centrifuged with the same parameters as above. Finally, the supernatant containing phage particles was collected.

Measuring transduction efficiency

Transduction efficiencies were measured as previously described²¹. In brief, target bacterial cells were grown to OD₆₀₀ ~0.5 (250 r.p.m. at 37 °C), followed by 15-min-long incubation on ice, during which dilutions of the transducing phage particles were prepared with tenfold dilution steps. Then, 50 µl of target cells were mixed with 50 µl of phage particles from each dilution. Plates were incubated at 37 °C at 180 r.p.m. for 1 h. Samples then were spotted on antibiotic-supplied agar plates. Transductant forming units per ml (t.f.u. ml⁻¹) were calculated on the basis of colony counts.

Assembly of transducing particles containing the metagenomic libraries

E. coli K12 BW25113 strain containing phage-tail-encoding plasmids were electroporated with 30 ng of each plasmid library in five parallels to achieve suitable colony numbers, then plated on antibiotic-containing LB agar plates and grown overnight. Following growth, cells were stored in 20% glycerol at −80 °C. Next, frozen cells containing the library were grown in 40 ml LB supplemented with kanamycin 50 and streptomycin 100 by shaking at 230 r.p.m. at 37 °C until OD₆₀₀ 0.7. Cells were cooled down on ice, centrifuged at 2,000 × g (4 °C, 10 min) and resuspended in LB medium. Then, the T7∆(gp11-12-17) bacteriophage was used to infect cells at MOI 2–3. Following 2 h of incubation (100 r.p.m. at 37 °C), cells were treated with 2% chloroform and vortexed. The mixture was then centrifuged and supernatant was collected.

Delivery of the metagenomic libraries by transducing phage particles and by electroporation

Overnight cultures of the corresponding bacterial strains were diluted to OD₆₀₀ 0.1 in 50 ml LB medium to grow at 230 r.p.m. at 37 °C until OD₆₀₀ 0.5. Next, we added 20 ml of library containing transducing particles to the cells, followed by one h incubation at the same parameters. Next, cells were centrifuged at 2,200 × g for 10 min at 4 °C, resuspended in 1–5 ml LB medium, plated on LB + kanamycin 50 and grown overnight. The next day, cells were collected and stored with glycerol at −80 °C. Of each library, 50 ng was electroporated into E. coli K12 BW25113 in five parallels. Cells were recovered in SOC medium for one h at 37 °C and plated on LB + kanamycin50 plates and grown overnight. The next day, cells were collected and stored in 20% glycerol at −80 °C.

Phage tail mutagenesis

To locate the HRDRs of the tail fibre genes, we used pairwise sequence alignment, where the recently identified HRDRs of gp17 of T3 coliphage²⁹ were aligned to the tail fibre sequences of Escherichia phage T7 gp17, Salmonella phage ΦSG-JL2 gp17 and Salmonella phage Vi06 gp43. The determined sites and the proximal regions were then subjected to targeted mutagenesis by DIvERGE²², a technique based on the targeted incorporation of mutational load carrying 90-mer oligos. In brief, E. coli BW25113 cells carrying the phage-tail-encoding plasmid to be mutated and the plasmid mediating the mutagenesis²² were grown to ~OD₆₀₀ 0.3–0.4 in TB (250 r.p.m. at 37 °C) supplied with appropriate antibiotics. Next, m-toluic acid was added (1 mM final concentration) to induce gene expression and after one h incubation, cells were transferred to ice for 15 min. Cell culture was made electrocompetent by repeated washing and centrifuging (2,200 × g, 4 °C, 10 min, three times), then electroporated with 2.5 µM oligos (Supplementary Table 11). Following recovery in TB (250 r.p.m., 37 °C, one h), cells were transferred to 19 ml TB supplied with appropriate antibiotics and left to grow overnight. Mutagenesis cycle was repeated if it was deemed necessary.

Selection of mutant phage tails with improved transduction efficiency

To select for tail mutants with improved delivery capacity, we applied a transduction optimization protocol. In brief, we chose three pathogenic bacterial strains (Enterobacter cloacae ATCC 23355, Shigella sonnei HNCMB 25021 and E. coli NCTC 13351) based on initial weak T7 bacteriophage infectivity. These target bacterial cells were grown to ~OD₆₀₀ 0.5 (250 r.p.m. at 37 °C) in LB, cells were placed on ice for 15 min, mixed with 2 ml of phage particles in a 1:1 volume ratio, and incubated at 37 °C and 100 r.p.m. for one h. The mixture was then plated and placed at 37 °C to grow overnight. The same protocol was carried out with non-mutagenized wild-type phage-tail-carrying particles. Colonies were pooled the next day and plasmid DNA was isolated using GeneJET plasmid miniprep kit (Thermo Fisher), then further purified using DNA Clean and Concentrator-5 (Zymo Research kit, D4004). Of the plasmids, 100 ng were electroporated into E. coli BW25113 cells. After recovery, cells were supplied with appropriate antibiotics, spread onto agar plates after one h of incubation and left to grow overnight. The following day, the cells were pooled in 4 ml LB, 250 µl were transferred into 40 ml TB supplied with appropriate antibiotics and grown to ~OD₆₀₀ 0.7 (250 r.p.m. at 37 °C). After growth, cells were placed on ice for 15 min, centrifuged (2,200 × g, 4 °C, 10 min) and resuspended. Next, cell cultures were infected with T7∆(gp11-12-17) bacteriophages. After two h (100 r.p.m. at 37 °C), cells were treated with 2% chloroform and vortexed. After centrifugation at the above parameters, phages present in the supernatant were collected. The transduction of the investigated bacterial strain was repeated until saturation in the number of transduced cells (~two or three rounds) was observable. Finally, plasmids from single colonies were sequenced to reveal tail mutations.

Quantifying replicative phage contamination

E. coli cells containing MGP4240 or MGP4240_gp17^V544G and pZE21_p15A plasmids were infected with T7Δ(gp11-12-17) phage to package the pZE21_p15A plasmid. The resulting phage particles were used to generate phage lysates in E. coli BW25113 and S. sonnei HNCMB 25021 harbouring either MGP4240 or MGP4240_gp17^V544G. The presence of the phage-tail-encoding plasmids in the target cells was necessary for replicative phage contamination to form plaques. For the plaque assays, 4 ml top agar was prepared and supplemented with 100 µg ml⁻¹ streptomycin (Sigma, S6501-25G) and 400 µl of the overnight cultures. Finally, from each phage stock, 10 µl was dropped onto the top agar in 1–10¹⁰ times dilutions.

Site-directed mutagenesis of phage-tail-encoding plasmids

For functional metagenomic library delivery, the mutation identified in the T7 gp17^V544G phage tail variant was introduced into plasmid MGP4240²¹ by using whole plasmid amplification with primers carrying the corresponding mutation, followed by DpnI (Thermo Fisher, ER1701) treatment to eliminate the original, methylated template plasmid DNA and subsequent gel electrophoresis, gel extraction and self-ligation. The plasmids were then electroporated into E. coli BW25113 cells. Transformants carrying the desired constructs were identified by PCR and validated via sequencing.

Functional selection of antibiotic resistance

Functional selections for resistance were performed on Mueller Hinton Broth (Sigma, 90922) agar plates containing a concentration gradient of a given antibiotic (adapted from ref. ⁵⁴). Antibiotics were purchased from Sigma or MedChem Express. The number of plated cells covered at least 10× the size of the corresponding metagenomic library. Plates were incubated at 37 °C for 24 h. For each functional selection, a control plate was prepared with the same number of cells containing the empty plasmid (that is, the plasmid without a cloned DNA fragment in the multiple cloning site) that showed the inhibitory zone of the antimicrobial compound for the cells without any resistance plasmid. The resistant clones from the libraries were isolated by washing together the sporadic colonies from the plate region (distal to the inhibition zone and containing higher antibiotic concentration), defined by visual inspection in comparison to the inhibition zone from the control plate. Half of the culture suspended in LB was used for plasmid isolation (GeneJET plasmid miniprep kit; Thermo Fisher, PLN70-1KT), and the rest was frozen with glycerol and stored at −80 °C.

Sample preparation for sequencing

The obtained resistance-conferring plasmids were sequenced with a hybrid sequencing pipeline (Extended Data Fig. 5) based on ref. ³⁴. Long-read sequencing identifies the metagenomic DNA fragments (inserts) and the two 10-nt-long random barcodes pre-cloned up- and down-stream (uptag and downtag, respectively) of each metagenomic DNA fragment. Aliquots of plasmid DNA preparations obtained from each screen were pooled in an equimolar ratio. Genomic DNA contamination was removed from the mixture by Lambda-exonuclease and Exonuclease-I double digestion. The resulting sample was cleaned (DNA Clean and Concentrator-5, Zymo Research kit) and quantified. Next, the plasmid mixture was linearized by adding 5 U of SrfI restriction endonuclease (NEB, R0629S) for every 1 µg of plasmid DNA (one h at 37 °C, followed by inactivation at 65 °C for 20 min), and DNA was quantified using Qubit dsDNA broad-range assay kit (Thermo Fisher,Q33266) before applying to Oxford Nanopore long-read sequencing. Parallel, multiplexed short-read deep sequencing was applied on each functional metagenomic plasmid DNA preparation (previous pooling) to associate nanopore contigs with screening samples (Extended Data Fig. 5). To this end, we amplified the up- and downtag barcodes on the plasmid preparations of each selection experiment separately, using Illumina specific forward and reverse primer pairs. Each primer pair contained P5 and P7 adapter sequences, respectively, and 8-nt-long barcodes for multiplexing and plasmid annealing sites (Supplementary Table 11). We performed PCR using Phusion high-fidelity DNA polymerase (Thermo Fisher, F530S) using the following reaction mixture: 15 ng of template plasmid DNA, 4 µl 5× GC buffer, 0.2 µl Phusion high-fidelity DNA polymerase, 0.6 µl DMSO (dimethyl sulfoxide), 0.2 mM dNTPs, 0.5–0.5 µM forward and reverse primers and water in a final volume of 20 µl. The following thermocycler conditions were used: 95 °C for five min, 30 cycles of 95 °C for 30 s + 59 °C for 30 s + 72 °C for 5 s, 72 °C for seven min. Following concentration measurement of each PCR reaction, we mixed the samples in a 1:1 mass ratio. Next, we isolated the 137-bp-long fragment mixture from 0.75% agarose gel.

Nanopore sequencing

Libraries were prepared by using a ligation sequencing kit (Oxford Nanopore Technologies, SQK-LSK109) with 1 µg plasmid DNA. The DNA was end-prepped with the NEBNext FFPE Repair (M6630S) and Ultra II End Prep kit (E7546S), purified using Agencourt AMPure XP (Beckman Coulter, A63882) and then the adapter ligated using NEBNext Quick T4 DNA ligase (E6056S). Finally, the adapted library was purified by Agencourt AMPure, quantified using Qubit 3.0, mixed with ONT running buffer and loading beads, primed with FLO-MIN106 9.4.1 SpotON flow cell attached to a MinION device and run for 72 h. Guppy algorithm (v8.25) with high-accuracy config settings was used for basecalling. Raw reads were filtered on the basis of quality value (QC ≥ 7) and length (4,000–8,000 bp) using NanoFilt v2.7.1⁵⁵. Reads were mapped to the reference sequence with minimap2 (v2.17)⁵⁶; SAM files were converted to sorted BAMs; the insert sequences were exctracted, and barcodes were identified and added to the read/insert names applying samtools tview (1.11-9-ga53817f) subcommand⁵⁷; individual FASTQ files were created using SEQTK (v0.13.2)⁵⁸; consensus sequences were generated using SPOA (v4.0.2)⁵⁹ with the following parameters: -l 0 -r 0 -g -2. Finally, the raw consensus inserts were polished using the relevant set of insert sequences by minimap2 and racon (v1.4.19)⁵⁶ to create the final consensus inserts with at least 100× coverage. Delivered metagenomic DNA fragment lengths and diversities were determined by using long-read deep sequencing right after electroporation into E. coli BW25113 and transduction into Salmonella enterica subsp. enterica serovar Typhimurium str. LT2, K. pneumoniae NCTC 9131 and S. sonnei HNCMB 25021. Shannon alpha diversity indices (H) were calculated on the basis of the frequency of each of the contigs of all hosts using the vegan R package (2.5-7)⁶⁰.

Illumina sequencing

Pooled sequencing libraries were denatured with 0.1 M NaOH, diluted to 12 pM with HT1 hybridization buffer (Illumina) and mixed with 40% PhiX Control v3 (Illumina) sequencing control library. Denatured sequencing pools were loaded onto MiSeq Reagent kit V2-300 (Illumina) and 2 × 70 bp sequence reads were generated with an Illumina MiSeq instrument with custom read 1, read 2 and index 1 sequencing primers spiked in the appropriate cartridge positions (12, 14 and 13, respectively) at a final concentration of 0.5 µM.

Host ranges of the ARGs encoded by the functional metagenomic DNA contigs

Resistant plasmid pools collected from the metagenomic screen were mixed and re-transformed or re-electroporated into the four hosts. Selection experiments were performed on gradient agar plates as described previously (see ‘Functional selection of antibiotic resistance’ above). Resistant colonies were collected and following plasmid preparation, barcodes on the plasmids were sequenced by Illumina sequencing (Supplementary methods). For calculating the overlaps between functional ARG sets across species, we first estimated the accuracy of the screen by comparing the results to that of the MIC measurements of the 13 selected resistance-conferring DNA fragments. On the basis of these comparisons, we estimated the true positive, false positive, true negative and false negative rates of the screen. Next, we calculated an adjusted Jaccard index for each species pair, which takes into account the screen’s accuracy as follows. For each species, we replaced the original vector of presence/absence of detected resistance instances with a new vector where the original presence (absence) values were randomly kept with a probability equal to the positive (negative) predictive value (that is, the proportion of true positives among all positive cases and the proportion of true negatives among all negative cases). The procedure was repeated 50,000 times, and the medians and 95% confidence intervals of the Jaccard indices between pairs of species were calculated.

Resistance levels in the bacterial hosts

We measured how DNA fragments that provide antibiotic resistance to E. coli influence susceptibility in Shigella sonnei HNCMB 25021, K. pneumoniae NCTC 9131 and Salmonella enterica subsp. enterica serovar Typhimurium str. LT2. For this purpose, we used a representative set of 13 plasmids that were isolated in our antibiotic selection screens. For each strain, the provided resistance levels (that is, the MIC) were measured with a standard 12-step microdilution method in 96-well plates, and the MIC fold change was determined by comparing them to the MIC of the corresponding empty vector harbouring control strains. MICs were determined on the basis of cell growth (OD₆₀₀) after 24 h incubation (37 °C, 180 r.p.m.).

Sequencing data analysis and functional annotation of ARGs

Each consensus insert sequence from nanopore sequencing was associated with screening samples (host, resistome, antibiotic) by combining the Nanopore and Illumina datasets through the unique uptag and downtag barcodes with a custom R script. To identify ARGs in the metagenomic contigs, two parallel approaches were used: (1) Open Reading Frame (ORF) prediction with prodigal ⁶¹, followed by annotation with BLASTP search against CARD³⁵ and ResFinder³⁶ databases, with coverage >50 bp at e-value < 10⁻⁵ and (2) BLASTX search with the same parameters but without ORF prediction to decrease the risk of truncated ORFs due to frame-shifting sequencing errors. To remove low-fidelity sequencing data from the dataset, metagenomic DNA fragments supported by <10 consensus insert sequences in the nanopore dataset and <9 reads in the Illumina uptag and downtag barcode dataset were filtered out.

If a metagenomic DNA fragment contained more than one predicted ARG, ARGs known to act on an antibiotic class (based on CARD and ResFinder reference databases) other than the one we used in the selection experiment were filtered out. ARG sequences having at least 95% identity and coverage on the DNA sequence level were collapsed into ARG clusters³⁷. Each cluster was represented by the closest hit to known ARGs in the Card³⁵ and ResFinder³⁶ databases (Supplementary Table 6). Donor organisms from which the assembled DNA contig sequences originated were identified by nucleotide sequence similarity search using the DNA contigs as query against the NCBI Reference Prokaryotic database (RefProk, downloaded 21 March 2021) with a threshold e-value of 10⁻¹⁰. The taxonomic hierarchy (kingdom, phylum, class, order, family, genus, species) was acquired using the taxonomizr package in R (v0.8.0).

Mobilization of the isolated ARGs

To create the mobile gene catalogue (that is, a database of recently transferred DNA sequences between bacterial species⁴⁰), we downloaded 1,377 genomes of diverse human-related bacterial species from the Integrated Microbial Genomes and Microbiomes database as done previously⁴⁰ and 1,417 genomes of Gram-negative ESKAPE pathogens from the NCBI RefSeq database (Supplementary Table 8). Using NCBI blastn 2.10.1+⁶², we searched the nucleotide sequences shared between genomes belonging to different species. The parameters for filtering the NCBI blastn 2.10.1+ blast results were the following: minimum percentage of identity, 99%; minimum alignment length, 500; maximum alignment length, 20,000. The blast hits were clustered by cd-hit-est 4.8.1^63,64, with sequence identity threshold of 99%. We predicted the ORFs on the blast hits with prodigal v2.6.3⁶¹, keeping only those longer than 500 nt. Then, to generate the mobile gene catalogue, we compared them with the merged CARD 3.1.0³⁵ and ResFinder (d48a0fe)³⁶ databases using diamond v2.0.4.142⁶⁵. Finally, natural plasmid sequences were identified by downloading 27,939 complete plasmid sequences from the PLSDB database (v2020-11-19)⁴¹. Then, representative sequences of the isolated 114 ARG clusters were BLASTN searched both in the mobile gene catalogue and in natural plasmid sequences, with an identity and coverage threshold of 90%. ARGs that were present in the mobile gene catalogue and/or in natural plasmid sequences were considered as mobile.

Statistical analysis

Statistical analysis was performed using R (v4.1.1). The parametric two-sample t-test was used to assess the differences between the means of the groups of samples. Fisher’s exact test was used to determine significant associations between two variables. Shannon alpha diversity index was used to characterize the diversity of DNA contigs in the libraries using the vegan package (v2.5–7) in R⁶⁶. Data distribution was assumed to be normal, but this was not formally tested.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.