One small-insert library was prepared from genomic DNA obtained from 31 bacterial strains (Table 1) isolated from the sponges A. aerophoba, C. candelabrum and P. ficiformis (Supplementary Materials and Methods). Resistance of these isolates to polymyxin B, erythromycin, ciprofloxacin, cefotaxime, tetracycline, chloramphenicol, rifampicin, ampicillin and imipenem was investigated in a previous study (Supplementary Materials and Methods) (Versluis et al., under review). Additionally, in this study, resistance of the isolates was tested to gentamicin (50 μg/ml), D-cycloserine (50 μg/ml), chlortetracycline (50 μg/ml), amikacin (100 μg/ml) and trimethoprim (50 μg/ml).

An E. coli TOP10 strain harboring the plasmid pZE21 was inoculated into 10 ml LB broth with kanamycin (50 μg/ml). The culture was incubated at 37°C overnight. Plasmid isolation was performed with the GeneJET Plasmid Miniprep Kit (Thermo Fisher Scientific, Waltham, United States). The DNA concentration was measured by Qubit®; (Thermo Fisher Scientific), and the isolated plasmid was electrophoresed on a 1% agarose gel for visual inspection. Next, the plasmid was cut by HincII #R0103C (NEB, Ipswich, Massachusetts) according to manufacturer's instructions. Subsequently, the plasmid ends were dephosphorylated according to the protocol for dephosphorylation of 5′-ends of DNA by rSAP (NEB). Finally, the linearized and dephosphorylated plasmid DNA was purified with the Clean-up Concentrator kit (A&A Biotechnology, Gdynia, Poland).

To prepare a small-insert library with a mixture of genomic DNA from 31 pure culture isolates (Table 1), we individually isolated DNA with the MasterPure^™ DNA Purification Kit (Epicentre, Madison, Wisconsin). The DNA concentration was measured by Qubit® 2.0, and DNA from the different isolates was then pooled at equivalent mass. DNA from sponge tissue samples was isolated with the DNeasy Blood & Tissue Kit (Qiagen, Hilden, Germany). DNA was fragmented by the E210 sonicator (Covaris, Woburn, Massachusetts) where default operating conditions were used to achieve an average fragment size of 2 kbp. The fragmented DNA was electrophoresed on a 1% agarose gel, and DNA ranging from 0.7 to 5 kbp was purified by using the GeneJet Gel Extraction Kit (Thermo Fisher Scientific). In the final gel purification step, DNA was eluted with nuclease-free water instead of elution buffer. Subsequently the DNA ends were repaired with the NEBNext® End Repair Module (NEB) according to manufacturer's instructions. The end-repaired DNA was cleaned with the Clean-up Concentrator kit (A&A Biotechnology), and elution was performed with nuclease-free water. The TaKaRa MightyMix DNA ligation kit (Clontech Laboratories Inc., Mountain View, United States) was used for ligation of the end-repaired fragments into the linearized pZE21 vector. Each ligation mixture consisted of 2 μl Ligation Mix and 2 μl nuclease-free water containing 600 ng insert DNA and 120 ng linearized vector (5:1 ratio). The resulting product was used for two electroporation reactions at 2 μl each, thereby creating two distinct small-insert libraries. At electroporation, 2 μl ligation mixture was added to 50 μl One Shot® TOP10 Electrocomp® E. coli cells (Thermo Fisher Scientific). The sample was transferred to a 1 mm cuvet and subjected to an electric pulse at 1.8 kV, 2.5 μF, and 200 Ω. To recover the E. coli cells, 1 ml SOC medium was added, and the sample was incubated at 37°C for 1 h. After recovery, 2 μl of the cell suspension was ten-fold serially diluted, and plated on LB agar media containing kanamycin (50 μg/ml) in order to estimate library size based on CFU counts. For each library, twelve clones were picked and colony PCRs were performed with primers pZE21_81_104_57C and pZE21_151_174rc_58C (Sommer et al., 2009), which flank the insertion site. The average product size (nt), the CFU count (in no. of colonies per μl) and the total volume (in μl) were multiplied in order to estimate the library size. In order to achieve library amplification, the remainder of the cell suspension after recovery was transferred to 10 ml LB medium with kanamycin (50 μg/ml) and incubated overnight at 37°C. After library amplification, another CFU count was performed. The total number of CFUs pre- and post-amplification was divided in order to estimate the extent of library expansion. Libraries were mixed to contain 20% glycerol and stored in cryotubes at −80°C. We also verified the function of a β-lactamase resistance gene by cloning it into the pZE21 vector as a single gene using the methods described in this paragraph, starting from the ligation step. The insert DNA consisted of amplicon sequences that were obtained by gene-specific PCR.

Libraries were plated at 25X (predicted) library coverage (i.e., every clone is expected to be present 25 times) on LB agar plates containing one of the following antibiotics: ampicillin (20 μg/ml), ciprofloxacin (1 μg/ml), tetracycline (20 μg/ml), chloramphenicol (20 μg/ml), polymyxin B (2 μg/ml), trimethoprim (2 μg/ml), erythromycin (200 μg/ml), cefotaxime (25 μg/ml), rifampicin (25 μg/ml), imipenem (20 μg/ml), gentamicin (20 μg/ml), D-cycloserine (100 μg/ml), chlortetracycline (20 μg/ml) and amikacin (50 μg/ml). E. coli TOP10 cells without a plasmid were used as negative control. Inserts from antibiotic resistant clones were Sanger sequenced from both the 5′ and 3′ flanks with respectively the pZE21_81_104_57C and pZE21_151_174rc_58C primers. Clones were assumed to contain identical inserts if sequences were more than 99% identical over a stretch of >400 bp. One clone per set of identical clones was selected for further analysis. The full-length inserts of these representative clones were sequenced by primer walking. Sanger sequencing was performed by flanking primers custom-designed with Primer3Plus (Untergasser et al., 2007) until GeneStudio version 2.2.0.0 could assemble a contig that spanned the entire insert. Finally, vector sequences were removed by DNA Baser version 3.5.4.2.

ORFs were identified by MetaGeneMark (Zhu et al., 2010), and the corresponding nucleotide and derived amino acid sequences were extracted. Amino acid sequences were used for BLASTp (Altschul et al., 1990) against the NCBI non-redundant protein database and the CARD database (McArthur et al., 2013). BLASTn of the sequences against the NCBI nr/nt and CARD databases were performed as well. Resistance functions were also predicted with HMMER (http://hmmer.org/) using the profile hidden Markov models (pHMMs) of the Resfams database as a reference (Gibson et al., 2014). We defined a resistance gene to be identified with high confidence if it was detected at an E-value of <1E-7 by either BLASTp or BLASTn against the Comprehensive Antibiotic Resistance Database (CARD), or by employing the pHMMs of the Resfams database. In addition, the resistance function needed to match the observed resistance phenotype. InterProScan (Mitchell et al., 2015) was used for classifying proteins into families and for predicting the presence of domains. Clustal Omega (Sievers et al., 2011) was used to globally align resistance genes with their best hit in NCBI's non-redundant protein database and to calculate a global amino acid identity value. ISfinder (Siguier et al., 2006) was used to identify IS elements. The ACLAME server (Leplae et al., 2004) was used to investigate if parts of the insert sequences were previously identified on mobile genetic elements.

Since one small-insert library was prepared with genomic DNA from 31 different sponge-associated bacterial isolates (as opposed to metagenomic DNA), we could use insert sequence-specific PCRs to assign inserts to the bacterium of origin. For this purpose, primers were designed with Primer3Plus (Untergasser et al., 2007). The reaction mixture of a given detection PCR consisted of: 10 μl 5X Green GoTaq Reaction Buffer (Promega, Fitchburg, Wisconsin), 0.2 mM dNTPs (Promega), 1 μM forward primer, 1 μM reverse primer, 0.5 μl GoTaq G2 DNA polymerase (5 U/μl, Promega), 1 μl genomic DNA (10–20 ng/μl), and 32.5 μl nuclease-free water (Promega). The PCR program consisted of: initial denaturation of 30 s at 98°C; 35 cycles of denaturation at 98°C for 10 s, annealing at 57°C for 20 s, and extension at 72°C for 20 s; and final extension at 72°C for 10 min. The PCR product was analyzed on a 1% agarose gel. If neither the initial PCR nor a repetition with a 1°C reduction in annealing temperature yielded a product, then new primers were designed. If a detection PCR yielded a product with genomic DNA from multiple isolates, then the amplification products were sequenced by Sanger sequencing. Inserts were assigned to an isolate if the Sanger sequence was 99–100% identical to the insert sequence.

The metagenomic insert sequences were deposited under accession numbers KU577908 to KU577944. The insert that contained bla_PSV-1 as a single gene was deposited under accession KU926347.

In order to investigate the presence and diversity of functional AR genes in sponge microbiota, small-insert libraries were prepared in E. coli with a mixture of gDNA from 31 bacterial strains isolated from the sponges A. aerophoba, P. ficiformis, and C. candelabrum (library I-31; 0.8 Gb), and with DNA directly isolated from the same three sponges (libraries Aa; 0.2 Gb, Cc; 0.8 Gb, Pf; 1.7 Gb). 16S ribosomal RNA (rRNA) gene amplicon sequence analysis of the bacterial gDNA pool that was used to prepare library I-31 confirmed that DNA from all 31 isolates was present (data not shown). Amongst these 31 sponge-associated bacteria, resistance to 14 out of 16 antibiotics was observed, whereas merely intermediate resistance, signifying a reduced growth rate, to imipenem was observed. None of the isolates were resistant to rifampicin (Table 1). Screening of the four different libraries for AR yielded 37 clones with unique inserts conferring resistance to D-cycloserine (n = 6), gentamicin (n = 1), amikacin (n = 7), trimethoprim (n = 17), chloramphenicol (n = 1), rifampicin (n = 2) and ampicillin (n = 3) (Table 2). No clones were obtained that were resistant to chlortetracycline, polymyxin, erythromycin, ciprofloxacin, cefotaxime, tetracycline or imipenem. The majority of resistant clones (30 of 37) were derived from library I-31, whereas 2 clones, 0 clones and 5 clones were derived from libraries Aa, Cc, and Pf, respectively.

The unique full-length inserts (mean insert size = 2974 bp ± 1627 [s.d.]) were sequenced by Sanger sequencing, which occasionally required multiple iterations of primer walking. We defined a resistance gene to be identified with high confidence if the resistance function was assigned at an E-value of <1E-7, and furthermore the resistance function needed to match the observed resistance phenotype.

In total, 26 resistance genes distributed over 24 different inserts were identified that met our confidence threshold (Table S1). The uniqueness of the resistance genes was evaluated by performing a global alignment at the amino acid level with the closest BLASTp hit in NCBI's non-redundant protein database. The majority of confidently identified resistance genes obtained from library I-31 (17 of 21) were predicted to code for proteins that had >70% identity at the amino acid level with known gene products, whereas proteins encoded by the genes (n = 5) obtained from the metagenomic libraries of environmental sponge DNA had <60% amino acid identity with known gene products (Figures 1A,B). Overall, most confidently identified resistance genes (15/26) were predicted to encode trimethoprim-resistance-conferring dihydrofolate reductases (Figure 1C). The confidently identified AR genes with the lowest similarity to known genes were predicted to code for a glycerol-3-phosphate acyltransferase and a GNAT family acetyltransferase with respectively 32 and 36% amino acid identity with the closest hit in NCBI's non-redundant protein database. Both of these amikacin-resistance-conferring genes were detected in libraries of environmental sponge DNA.

For 13 insert sequences a resistance gene could not be identified meeting our threshold. Those inserts conferred resistance to amikacin (n = 4), D-cycloserine (n = 6), and trimethoprim (n = 3) (Table 3, Table S2). In these cases, besides BLAST searches against the CARD and NCBI nr/nt databases and motif detection by application of Resfams pHMMs, also protein domain and protein family predications by InterProScan were evaluated in order to manually predict AR gene presence. Most (12 of 13) inserts contained multiple open reading frames (ORFs), which prevented the unequivocal assignment of resistance functions to individual genes. For three of four inserts that conferred resistance to amikacin, we predicted the presence of genes encoding aminoglycoside-modifying enzymes. On one of these inserts (clone Env_Ami3) a gene was identified that we predicted to code for an enzyme that putatively modifies amikacin by phosphorylation. On the other two inserts, we predicted genes encoding an aminoglycoside acetyltransferase (clone Iso_Ami2) and an aminoglycoside methyltransferase (clone Iso_Ami3). The fourth insert conferring amikacin resistance (clone Iso_Ami4) contained a single gene that shared 93% amino acid identity with an amino acid transporter, and as such we predicted the gene to confer AR via antibiotic efflux. Six inserts conferred resistance to D-cycloserine. On all of these six inserts we predicted the presence of genes encoding proteins involved in transmembrane transport based on high similarity to known transporters. We predicted that three of the six D-cycloserine resistance-conferring inserts harbored genes encoding transporters that belong to the major facilitator superfamily. Three other inserts on which an AR gene could not be confidently identified conferred resistance to trimethoprim. On one of these inserts (clone Env_Trim4), we predicted that a thymidylate synthethase may be responsible for the resistance phenotype. The two other inserts conferring trimethoprim resistance (clones Iso_trim13 and Iso_trim14) putatively encoded oxidoreductases that could confer trimethoprim resistance.

Resistance genes that were confidently identified in library I-31 were assigned to their strain of origin by gene-specific PCRs, where the PCR amplicon covered at least part of the predicted resistance gene. Gene presence was checked in all 31 strains used for building the library by dedicated PCR reactions, which resulted in confidently identified AR genes being assigned to 15 of 31 bacterial strains in library I-31 (Table 1). Insert sequences for which no resistance gene could be identified with high confidence were assigned to a strain of origin at the loci given in Table S2. We discovered that several inserts may contain hybrid sequences, that is, the inserts consist of two different DNA fragments that were ligated into the vector. For example, the insert of clone Iso_Dcy6 was assigned to Ruegeria spp. in the region between 513 and 1372 bp and assigned to Acinetobacter radioresistens DN138_5C8 in the region between 2240 and 2851 (Table S2). We found that about half of the confidently identified resistance genes (11/21) were assigned to Bacillus spp., and Bacillus was the only taxon to which gentamicin, amikacin, chloramphenicol and rifampicin resistance genes were assigned (Table 1, Figure 1D). Even though gentamicin, amikacin, chloramphenicol, rifampicin and ampicillin resistance genes were assigned to Bacillus spp., the strains themselves were not found to be resistant to these antibiotics. Resistance genes that were assigned to the remaining taxa comprised mostly trimethoprim resistance genes, with an exception being an ampicillin resistance gene that was assigned to Pseudovibrio spp. (Table S1). Strains to which trimethoprim resistance genes were assigned were not trimethoprim resistant themselves in 8 out of 9 cases.

The insert sequence of ampicillin-resistant clone Iso_Amp3 contained a gene predicted to encode a class A β-lactamase. The gene was assigned to Pseudovibrio ascidiaceicola DN64_1D03 and Pseudovibrio ascidiaceicola DN64_8G1, and the predicted protein shared 58% amino acid identity with the closest hit in the CBMAR database (Srivastava et al., 2014). Since this predicted protein displayed high divergence with members of known β-lactamase families it is a candidate to be classified into a new family (Jacoby, 2006). Therefore, we further tested the gene's function by cloning the gene-specific amplicon sequence that was generated with genomic DNA of Pseudovibrio ascidiaceicola DN64_1D03 as template into the pZE21 vector. The new clone, PSV1 (accession KU926347), was resistant to ampicillin, which confirmed the functionality of the β-lactamase gene. Remarkably, the coding sequence of the gene cloned using the genome-derived amplicon sequences was 75 bp longer than the coding sequence of the β-lactamase gene in the original metagenomic insert. Sanger sequencing of the genome-derived amplicon yielded a gene sequence that shared 100% nucleotide identity with the gene sequence in clone PSV1. Hence, we concluded that the longer β-lactamase gene in clone PSV1 is the exact gene that is present in the genome of Pseudovibrio ascidiaceicola DN64_1D03. We predict that even though a frameshift mutation yielded a shorter gene in case of the small-insert library, the encoded enzyme retained its function. We classified the genome-derived (longer) gene into a novel β-lactamase family named PSV after Pseudovibrio. The first member was designated bla_PSV-1. BLASTp of the gene against the CBMAR database showed that LEN and SHV were the closest β-lactamase families. The best hit from the LEN and SHV families in both cases shared 50% amino acid identity with bla_PSV-1. The best hit was bla_LEN−22 (accession AM850912), with which bla_PSV-1shared 41% global amino acid identity. The best hits in the NCBI non-redundant protein sequences database (≥76% amino acid identity) were exclusively genes in sequenced genomes of Pseudovibrio spp. (Figure 2). Additional resistance testing of clone PSV1 showed that bla_PSV-1 conferred resistance to penicillin (50 μg/ml) but not against cefotaxime (20 μg/ml) and imipenem (20 μg/ml).

Horizontal gene transfer is the main mechanism through which resistance genes are disseminated. Therefore, we investigated our sequences for evidence of gene mobilization. Comparison of the insert sequences with the ACLAME plasmid database showed that they were not previously identified on plasmids. Furthermore, ISfinder did not identify insertion sequences. We investigated whether the resistance genes identified in this study had been observed previously in a different genomic context by comparing our sequences with those in the NCBI nr/nt database (Table S3). A β-lactamase and an ADP-ribosyl transferase (those on clones Iso_Amp2 and Iso_Rif2, respectively), both of which were assigned to Bacillus aryabhattai DN67_5C7, shared very high (≥99%) identity with genes previously found in the non-marine Bacillus megaterium QM B1551, which suggests a relatively recent dissemination event. The insert of clone Iso_Amp2 shared 99% identity with the genome of Bacillus megaterium QM B1551 over the full-length insert (1352 bp). This insert produced four other significant alignments with sequences in the NCBI nr/nt database at >90% nucleotide identity, all of which belonged to Bacillus megaterium strains. The insert of clone Iso_Rif2 was 6717 bp in length but only the region from 5198 to 6717 bp aligned with the genome of Bacillus megaterium QM B1551. Besides the rifampicin-resistance-conferring gene, this region encodes a hypothetical protein and a gene that putatively encodes an enzyme involved in site specific recombination of plasmids. The gene putatively encoding a recombination enzyme was the only gene involved in gene mobilization that was predicted on the inserts.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.